# Complexation of Zr and Hf in fluoride-rich hydrothermal aqueous fluids and its significance for high field strength element fractionation

Anselm Loges, Marco Manni, Marion Louvel, Max Wilke, Sandro Jahn, Edmund Welter, Manuela Borchert, Shilei Qiao, Stephan Klemme, Bettina G Keller, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Anselm Loges, Marco Manni, Marion Louvel, Max Wilke, Sandro Jahn, et al.. Complexation of Zr and Hf in fluoride-rich hydrothermal aqueous fluids and its significance for high field strength element fractionation. Geochimica et Cosmochimica Acta, 2024, 366, pp.167-181. 10.1016/j.gca.2023.12.013. insu-04431876

HAL Id: insu-04431876

https://insu.hal.science/insu-04431876

Submitted on 1 Feb 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Complexation of Zr and Hf in fluoride-rich hydrothermal aqueous fluids and its significance for high field strength element fractionation

Anselm Loges, Marco Manni, Marion Louvel, Max Wilke, Sandro Jahn, Edmund Welter, Manuela Borchert, Shilei Qiao, Stephan Klemme, Bettina G Keller, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Anselm Loges, Marco Manni, Marion Louvel, Max Wilke, Sandro Jahn, et al.. Complexation of Zr and Hf in fluoride-rich hydrothermal aqueous fluids and its significance for high field strength element fractionation. Geochimica et Cosmochimica Acta, 2024, 366, pp.167-181. 10.1016/j.gca.2023.12.013. insu-04431876

HAL Id: insu-04431876

https://insu.hal.science/insu-04431876v1

Submitted on 1 Feb 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Complexation of Zr and Hf in fluoride-rich hydrothermal aqueous fluids and its significance for high field strength element fractionation

<figure><img src="image_4.png" /></figure>

Check for
updates

Anselm Loges<sup>a, *,</sup>, Marco Manni<sup>b</sup>, Marion Louvel<sup>c, d</sup>, Max Wilke<sup>e</sup>, Sandro Jahn<sup>f</sup>, Edmund Welter<sup>g</sup>, Manuela Borchert<sup>c, g</sup>, Shilei Qiao<sup>a</sup>, Stephan Klemme<sup>c</sup>, Bettina G. Keller<sup>b</sup>, Timm John<sup>a</sup>

## ARTICLE INFO

Associate editor: Xiandong Liu

**Keywords:**

Zirconium

Hafnium

Aqueous speciation

Hydrothermal fluoride complexation

Geochemical twin fractionation

## ABSTRACT

Zirconium and hafnium behave nearly identically in most geological processes due to their identical nominal ionic charge and similar radius. Some of the most pronounced exceptions from this rule are observed in fluoride-rich aqueous systems, suggesting that aqueous fluoride complexation may be involved in Zr/Hf fractionation. To understand the mechanisms causing this phenomenon, we investigated complexation of Zr<sup>4+</sup> and Hf<sup>4+</sup> in fluoride-rich (1.0 mol/kg HF) aqueous solutions at 40 MPa and 100–400 °C, using synchrotron X-ray absorption spectroscopy (X-ray absorption near edge structure and extended X-ray absorption fine structure) combined with classical and *ab initio* molecular dynamics simulations. The dominant experimentally observed complexes are [Zr(F<sub>2</sub>OH)<sub>4</sub>·2H<sub>2</sub>O]<sup>0</sup> and [Hf(F<sub>2</sub>OH)<sub>4</sub>·2H<sub>2</sub>O]<sup>0</sup>, respectively. The first coordination shell comprises a distorted octahedron, with fluoride and hydroxide ligands at a similar mean radial distance (1.9–2.0 Å) from the central cation, and H<sub>2</sub>O ligands at a slightly greater distance (>2.1 Å). With increasing temperature, the H<sub>2</sub>O ligands move further out, causing first an increasing distortion of the octahedron and subsequently a partial transition to less hydrated complexes as a certain fraction of the H<sub>2</sub>O molecules move to the second shell at > 3 Å. As a consequence, the radial distance of the F<sup>-</sup> and OH<sup>-</sup> anions from the central cation, as well as the overall average radial distance of the first shell decreases due to decreased steric repulsion from the H<sub>2</sub>O ligands. Both experiments and simulations agree in that Hf forms slightly shorter bonds to its nearest neighbors than Zr. The results suggest two hypotheses for the mechanism of Zr/Hf fractionation during precipitation of minerals from fluoride-rich hydrothermal solutions: 1) The heavy twin (Hf) prefers the lower coordination (shorter bonds) and is thus less likely to enter into the higher coordination found in the solids. This mechanism would be analogous to equilibrium isotope fractionation. 2) The change of Hf into a higher coordination environment (e.g., from solution to solid) is slower because it forms stronger ligand-bonds than Zr. This would be analogous to reactive kinetic isotope fractionation. In either case mass dependent fractionation qualitatively matches the observations but mass independent effects on bond strength may also be significant. Quantitative investigations of these effects are needed and may also shed light on the currently still somewhat enigmatic fractionation behavior of Zr isotopes.

## 1. Introduction

Aqueous fluids are among the most important mobile phases in the Earth's crust and are involved in many processes that redistribute and

concentrate certain elements to economically viable deposits. Understanding the transport mechanisms and solubility of metals in aqueous fluids is thus critical for the genetic interpretation of the natural distribution of these metals, and forms the basis for use of element

<footnote><sup>a</sup> Institut für Geologische Wissenschaften, Freie Universität Berlin, Malteserstr. 74-100, 12249 Berlin, Germany</footnote>

<footnote><sup>b</sup> Institut für Chemie und Biochemie, Freie Universität Berlin, Arnimallee 22, 14194 Berlin, Germany</footnote>

<footnote><sup>c</sup> Institut für Mineralogie, Universität Münster, Corrensstrasse 24, 48149 Münster, Germany</footnote>

<footnote><sup>d</sup> Institut des Sciences de la Terre d'Orleans, CNRS UMR7327, 1A rue de la Ferollerie, 45100 Orleans, France</footnote>

<footnote><sup>e</sup> Institut für Geowissenschaften, Universität Potsdam, Karl-Liebknecht-Str. 24, 14476 Potsdam, Germany</footnote>

<footnote><sup>f</sup> Institut für Geologie und Mineralogie, Universität zu Köln, Zülpicher Str. 49b, 50674 Cologne, Germany</footnote>

<footnote><sup>g</sup> Deutsches Elektronen-Synchrotron DESY, A Research Centre of the Helmholtz Association, Notkestrasse 85, 22607 Hamburg, Germany</footnote>

<footnote>* Corresponding author.
E-mail address: anselm.loges@fu-berlin.de (A. Loges).</footnote>
<figure><img src="image_3.png" /></figure>

# Complexation of Zr and Hf in fluoride-rich hydrothermal aqueous fluids and its significance for high field strength element fractionation

<figure><img src="image_4.png" /></figure>

Anselm Loges<sup>a, *,</sup>, Marco Manni<sup>b</sup>, Marion Louvel<sup>c, d</sup>, Max Wilke<sup>e</sup>, Sandro Jahn<sup>f</sup>, Edmund Welter<sup>g</sup>, Manuela Borchert<sup>c, g</sup>, Shilei Qiao<sup>a</sup>, Stephan Klemme<sup>c</sup>, Bettina G. Keller<sup>b</sup>, Timm John<sup>a</sup>

## ARTICLE INFO

Associate editor: Xiandong Liu

**Keywords:**

Zirconium

Hafnium

Aqueous speciation

Hydrothermal fluoride complexation

Geochemical twin fractionation

## ABSTRACT

Zirconium and hafnium behave nearly identically in most geological processes due to their identical nominal ionic charge and similar radius. Some of the most pronounced exceptions from this rule are observed in fluoride-rich aqueous systems, suggesting that aqueous fluoride complexation may be involved in Zr/Hf fractionation. To understand the mechanisms causing this phenomenon, we investigated complexation of Zr<sup>4+</sup> and Hf<sup>4+</sup> in fluoride-rich (1.0 mol/kg HF) aqueous solutions at 40 MPa and 100–400 °C, using synchrotron X-ray absorption spectroscopy (X-ray absorption near edge structure and extended X-ray absorption fine structure) combined with classical and *ab initio* molecular dynamics simulations. The dominant experimentally observed complexes are [Zr(F<sub>2</sub>OH)<sub>4</sub>·2H<sub>2</sub>O]<sup>0</sup> and [Hf(F<sub>2</sub>OH)<sub>4</sub>·2H<sub>2</sub>O]<sup>0</sup>, respectively. The first coordination shell comprises a distorted octahedron, with fluoride and hydroxide ligands at a similar mean radial distance (1.9–2.0 Å) from the central cation, and H<sub>2</sub>O ligands at a slightly greater distance (>2.1 Å). With increasing temperature, the H<sub>2</sub>O ligands move further out, causing first an increasing distortion of the octahedron and subsequently a partial transition to less hydrated complexes as a certain fraction of the H<sub>2</sub>O molecules move to the second shell at > 3 Å. As a consequence, the radial distance of the F<sup>-</sup> and OH<sup>-</sup> anions from the central cation, as well as the overall average radial distance of the first shell decreases due to decreased steric repulsion from the H<sub>2</sub>O ligands. Both experiments and simulations agree in that Hf forms slightly shorter bonds to its nearest neighbors than Zr. The results suggest two hypotheses for the mechanism of Zr/Hf fractionation during precipitation of minerals from fluoride-rich hydrothermal solutions: 1) The heavy twin (Hf) prefers the lower coordination (shorter bonds) and is thus less likely to enter into the higher coordination found in the solids. This mechanism would be analogous to equilibrium isotope fractionation. 2) The change of Hf into a higher coordination environment (e.g., from solution to solid) is slower because it forms stronger ligand-bonds than Zr. This would be analogous to reactive kinetic isotope fractionation. In either case mass dependent fractionation qualitatively matches the observations but mass independent effects on bond strength may also be significant. Quantitative investigations of these effects are needed and may also shed light on the currently still somewhat enigmatic fractionation behavior of Zr isotopes.

## 1. Introduction

Aqueous fluids are among the most important mobile phases in the Earth's crust and are involved in many processes that redistribute and

concentrate certain elements to economically viable deposits. Understanding the transport mechanisms and solubility of metals in aqueous fluids is thus critical for the genetic interpretation of the natural distribution of these metals, and forms the basis for use of element

<footnote><sup>a</sup> Institut für Geologische Wissenschaften, Freie Universität Berlin, Malteserstr. 74-100, 12249 Berlin, Germany</footnote>

<footnote><sup>b</sup> Institut für Chemie und Biochemie, Freie Universität Berlin, Arnimallee 22, 14194 Berlin, Germany</footnote>

<footnote><sup>c</sup> Institut für Mineralogie, Universität Münster, Corrensstrasse 24, 48149 Münster, Germany</footnote>

<footnote><sup>d</sup> Institut des Sciences de la Terre d'Orleans, CNRS UMR7327, 1A rue de la Ferollerie, 45100 Orleans, France</footnote>

<footnote><sup>e</sup> Institut für Geowissenschaften, Universität Potsdam, Karl-Liebknecht-Str. 24, 14476 Potsdam, Germany</footnote>

<footnote><sup>f</sup> Institut für Geologie und Mineralogie, Universität zu Köln, Zülpicher Str. 49b, 50674 Cologne, Germany</footnote>

<footnote><sup>g</sup> Deutsches Elektronen-Synchrotron DESY, A Research Centre of the Helmholtz Association, Notkestrasse 85, 22607 Hamburg, Germany</footnote>

<footnote>* Corresponding author.
E-mail address: anselm.loges@fu-berlin.de (A. Loges).</footnote>
distribution to quantitatively model past geological processes as well as a tool for the assessment of potentially economic deposits. Zirconium and hafnium are typical high field-strength elements (HFSE), meaning they occur naturally in form of small, highly charged cations. In addition to exploration and economic assessment of ore deposits, understanding the behavior of HFSE in fluids is of importance for their use as tracers of processes in large-scale geochemical cycles, particularly the evolution of the mantle-crust differentiation throughout geologic history, and subduction of oceanic crust by plate tectonic processes. This is particularly relevant because poorly understood fluid related HFSE fractionation processes may obscure magmatic fractionation trends, which are used to decipher the formation of Earth's crust and mantle-crust interaction (e. g., Weyer et al., 2003; Pfänder et al., 2007). The characteristic depletion of HFSE (including e.g., Ti, Zr, Nb, Hf, and Ta) in subduction related magmatic systems compared to large-ion lithophile elements (LILE, e.g., Rb, Sr, Cs, Ba, and Pb) is well documented (e.g., Perfit et al., 1980; McCulloch and Gamble, 1991) and helps deciphering the mass transfer from the subducting slab to the overlying mantle wedge and ultimately the magmas of the volcanic arc. Zirconium stable isotopes (δ⁹⁴/⁹⁰Zr) have in recent years been shown to fractionate strongly in magmatic systems. However, the underlying fractionation mechanisms are currently not quantitatively understood, likely due to the prominence of non-equilibrium fractionation (e.g., Inglis et al., 2019; Méheut et al., 2021; Tissot and Ibañez-Mejia, 2021; Bindeman and Melnik, 2022). Remarkably, δ⁹⁴/⁹⁰Zr values can be positively (e.g., Ibañez-Mejia and Tissot, 2019) or negatively (e.g., Zhu et al., 2023) correlated with Zr/Hf ratios, suggesting that the combined investigation of isotopic and elemental ratios may be a particularly valuable diagnostic tool for petrogenesis, once the fractionation mechanisms are fully understood.

Chondritic concentrations of Zr and Hf are low, 3.63 µg/g for Zr and 0.11 µg/g for Hf (Palme and O'Neill, 2014), but the upper continental crust is significantly enriched (190 µg/g Zr and 5.80 µg/g Hf; Taylor and McLennan, 1995) due to the lithophile nature of Zr and Hf. Both elements occur exclusively in tetravalent oxidation state in nature and are strongly concentrated in a few refractory minerals that, despite their low overall abundance, are near-ubiquitous in crustal silicate lithologies, except for some mafic and most ultramafic rocks. By far the most important of these minerals is zircon (ZrSiO₄), but baddeleyite (ZrO₂), zirconolite and a wide variety of more or less complex (alkali-)zirconosilicates, including vaslovite, elpidite, and eudialyte, are commonly found in magmatic rock-types that are silica-undersaturated with respect to quartz. Due to the identical charge and very similar ionic radius, 0.72 Å for Zr and 0.71 Å for Hf in octahedral coordination (Shannon, 1976), both elements are usually found together in these minerals in proportions closely resembling those in the host rock and the bulk silicate Earth, and are therefore referred to as a geochemical twin pair. Notable exceptions from this rule are late granitic zircons (e.g., Wang et al., 2010), some pegmatites (e.g., Neves et al., 1974; Yin et al., 2013; Kudryashov et al., 2020), and systems where Zr and Hf are incorporated into or adsorbed onto minerals from an aqueous solution (e.g., Bau, 1996; Schmidt et al., 2014; Inguaggiato et al., 2015; Censi et al., 2017). What mechanisms ultimately control Zr/Hf fractionation in various geological environments in spite of the very similar chemical properties of Zr and Hf remains the subject of ongoing research to which we seek to contribute with the present study. Aqueous complexation is known to be a major factor in fractionating elements (Brugger et al., 2016) and may also be responsible for geochemical twin fractionation. This mechanism is known for the twin pair Y/Ho in fluoride-rich solutions, where Y and Ho are dissolved as different core complexes (YF₂²⁻ difluoride and HoF²⁺ monofluoride, respectively; Loges et al., 2013). However, differences in hydration in the first and second coordination shell of otherwise similar complexes have also been shown to drastically influence the solubility of metal cations (Brugger et al., 2016 and references therein).

Although HFSE in general show low solubility in most natural aqueous fluids, highly mobile aqueous complexes of Zr and Hf can

increase the solubility of these elements by several orders of magnitude compared to a pure H₂O solvent (e.g., Migdisov et al., 2011; Wilke et al., 2012; Aseri et al., 2015). As hard Lewis acids, Zr⁴⁺ and Hf⁴⁺ form complexes most readily with hard bases such as OH⁻ and F⁻ (Pearson, 1963). Hydroxide and fluoride are consequently indeed documented to be the most relevant at 25 °C (Aja et al., 1995). Aqueous complexation of most elements is typically significantly different in geological fluids at elevated temperatures compared to room temperature, due to the strong decrease of the dielectric constant of water with increasing temperature and decreasing density (e.g., Weingärtner and Franck, 2005). Experimental studies on baddeleyite solubility have found that Zr complexation in F-bearing acidic hydrothermal fluids up to 500 °C and 100 MPa is dominated by ZrF(OH)₃ and ZrF₂(OH)₂ complexes and that the solubility of baddeleyite decreases with increasing temperature (Ryzhenko et al., 2008; Migdisov et al., 2011). The solubility of Hf in F-free aqueous solution at room temperature in form of [Hf(OH)₄]⁰ complexes is very limited between pH 4 and 9 (ca. 10⁻¹⁰ mol/kg) but a steep increase in solubility at pH < 4 in form of [Hf(OH)]³⁺ complexes or (presumably hydrated) Hf⁴⁺ ions as well as at pH > 9 as [Hf(OH)₅]⁻ or [Hf(OH)₆]²⁻ complexes has been documented (Rai et al., 2001). Complexation of Zr and Hf with chloride and sodium-silicate has been found to favor 6-fold coordination of the central cation in experimental high-p-T studies (Louvel et al., 2013; Wilke et al., 2012; 2013) as well as with ab initio molecular dynamics (AIMD) modeling (Jahn et al., 2015).

In this study, we investigate the coordination environment of Zr and Hf in HF-rich aqueous solutions under conditions relevant to natural fluorine-rich hydrothermal systems such as greisen or pegmatite fluids (40 MPa, 100–400 °C, 1 mol/kg HF) by combining experimental X-ray absorption near edge structure (XANES) and extended X-ray absorption fine structure (EXAFS) data with ab initio molecular dynamics (AIMD) and classical molecular dynamics (MD). This comprehensive approach is necessary for determining the complex symmetry and ligand ratios in mixed complexes with both F and O in the same coordination shell because the very similar and overall low scattering cross section of these elements do not allow to extract all necessary information from X-ray absorption spectroscopy alone (Farges, 1991; 1996; Louvel et al., 2013). Therefore, EXAFS should not be used to interpret the number of ligands but only the average distance, to which this method is much more sensitive. XANES data can provide insight into the coordination of elements, whose observable X-ray absorption edges show features that can be demonstrated to be indicative of certain coordination symmetries. Most commonly, K edges are used for this purpose. Several experimental studies have demonstrated, using minerals and solutions, that the double peak of the Zr K white line is highly sensitive to coordination symmetry (and thus coordination number) as well as ligand type (e.g., Louvel et al., 2013; Wilke et al., 2012; 2013). For the Hf L₃ edge, no comprehensive data set is available and the data that are published suggest less obvious effects on the XANES region of this edge (e.g., Wilke et al., 2013). Computed XANES spectra of the Zr K and Hf L₃ edge corroborate these experimental findings (Jahn et al., 2015). Classical MD simulations are a very useful complement to XANES and EXAFS spectra analyzed with the help of AIMD. They are computationally less demanding and can therefore sample longer timescales. Here, we use classical MD simulations to simulate the structure and dynamics of the ligands and their exchange with the surrounding solvent between 25 and 400 °C in order to understand the changes in average ligand distance observed with EXAFS. Because classical MD simulations use a parametrized potential energy function, they cannot account for instantaneous changes in the electron cloud caused be changes in the surrounding electrostatic field (dynamic polarization). This drawback becomes acceptable for hard ions, i.e., those with high charge and small radius like the HFSE, whose electron cloud shows little response to changes in the surrounding electrostatic field. We adjust the parameters for the dispersion interactions of the cations to match the experimental data. We propose two possible mechanisms for Zr/Hf fractionation in hydrothermal systems based on our experimental and computational
distribution to quantitatively model past geological processes as well as a tool for the assessment of potentially economic deposits. Zirconium and hafnium are typical high field-strength elements (HFSE), meaning they occur naturally in form of small, highly charged cations. In addition to exploration and economic assessment of ore deposits, understanding the behavior of HFSE in fluids is of importance for their use as tracers of processes in large-scale geochemical cycles, particularly the evolution of the mantle-crust differentiation throughout geologic history, and subduction of oceanic crust by plate tectonic processes. This is particularly relevant because poorly understood fluid related HFSE fractionation processes may obscure magmatic fractionation trends, which are used to decipher the formation of Earth's crust and mantle-crust interaction (e. g., Weyer et al., 2003; Pfänder et al., 2007). The characteristic depletion of HFSE (including e.g., Ti, Zr, Nb, Hf, and Ta) in subduction related magmatic systems compared to large-ion lithophile elements (LILE, e.g., Rb, Sr, Cs, Ba, and Pb) is well documented (e.g., Perfit et al., 1980; McCulloch and Gamble, 1991) and helps deciphering the mass transfer from the subducting slab to the overlying mantle wedge and ultimately the magmas of the volcanic arc. Zirconium stable isotopes (δ⁹⁴/⁹⁰Zr) have in recent years been shown to fractionate strongly in magmatic systems. However, the underlying fractionation mechanisms are currently not quantitatively understood, likely due to the prominence of non-equilibrium fractionation (e.g., Inglis et al., 2019; Méheut et al., 2021; Tissot and Ibañez-Mejia, 2021; Bindeman and Melnik, 2022). Remarkably, δ⁹⁴/⁹⁰Zr values can be positively (e.g., Ibañez-Mejia and Tissot, 2019) or negatively (e.g., Zhu et al., 2023) correlated with Zr/Hf ratios, suggesting that the combined investigation of isotopic and elemental ratios may be a particularly valuable diagnostic tool for petrogenesis, once the fractionation mechanisms are fully understood.

Chondritic concentrations of Zr and Hf are low, 3.63 µg/g for Zr and 0.11 µg/g for Hf (Palme and O'Neill, 2014), but the upper continental crust is significantly enriched (190 µg/g Zr and 5.80 µg/g Hf; Taylor and McLennan, 1995) due to the lithophile nature of Zr and Hf. Both elements occur exclusively in tetravalent oxidation state in nature and are strongly concentrated in a few refractory minerals that, despite their low overall abundance, are near-ubiquitous in crustal silicate lithologies, except for some mafic and most ultramafic rocks. By far the most important of these minerals is zircon (ZrSiO₄), but baddeleyite (ZrO₂), zirconolite and a wide variety of more or less complex (alkali-)zirconosilicates, including vaslovite, elpidite, and eudialyte, are commonly found in magmatic rock-types that are silica-undersaturated with respect to quartz. Due to the identical charge and very similar ionic radius, 0.72 Å for Zr and 0.71 Å for Hf in octahedral coordination (Shannon, 1976), both elements are usually found together in these minerals in proportions closely resembling those in the host rock and the bulk silicate Earth, and are therefore referred to as a geochemical twin pair. Notable exceptions from this rule are late granitic zircons (e.g., Wang et al., 2010), some pegmatites (e.g., Neves et al., 1974; Yin et al., 2013; Kudryashov et al., 2020), and systems where Zr and Hf are incorporated into or adsorbed onto minerals from an aqueous solution (e.g., Bau, 1996; Schmidt et al., 2014; Inguaggiato et al., 2015; Censi et al., 2017). What mechanisms ultimately control Zr/Hf fractionation in various geological environments in spite of the very similar chemical properties of Zr and Hf remains the subject of ongoing research to which we seek to contribute with the present study. Aqueous complexation is known to be a major factor in fractionating elements (Brugger et al., 2016) and may also be responsible for geochemical twin fractionation. This mechanism is known for the twin pair Y/Ho in fluoride-rich solutions, where Y and Ho are dissolved as different core complexes (YF₂²⁻ difluoride and HoF²⁺ monofluoride, respectively; Loges et al., 2013). However, differences in hydration in the first and second coordination shell of otherwise similar complexes have also been shown to drastically influence the solubility of metal cations (Brugger et al., 2016 and references therein).

Although HFSE in general show low solubility in most natural aqueous fluids, highly mobile aqueous complexes of Zr and Hf can

increase the solubility of these elements by several orders of magnitude compared to a pure H₂O solvent (e.g., Migdisov et al., 2011; Wilke et al., 2012; Aseri et al., 2015). As hard Lewis acids, Zr⁴⁺ and Hf⁴⁺ form complexes most readily with hard bases such as OH⁻ and F⁻ (Pearson, 1963). Hydroxide and fluoride are consequently indeed documented to be the most relevant at 25 °C (Aja et al., 1995). Aqueous complexation of most elements is typically significantly different in geological fluids at elevated temperatures compared to room temperature, due to the strong decrease of the dielectric constant of water with increasing temperature and decreasing density (e.g., Weingärtner and Franck, 2005). Experimental studies on baddeleyite solubility have found that Zr complexation in F-bearing acidic hydrothermal fluids up to 500 °C and 100 MPa is dominated by ZrF(OH)₃ and ZrF₂(OH)₂ complexes and that the solubility of baddeleyite decreases with increasing temperature (Ryzhenko et al., 2008; Migdisov et al., 2011). The solubility of Hf in F-free aqueous solution at room temperature in form of [Hf(OH)₄]⁰ complexes is very limited between pH 4 and 9 (ca. 10⁻¹⁰ mol/kg) but a steep increase in solubility at pH < 4 in form of [Hf(OH)]³⁺ complexes or (presumably hydrated) Hf⁴⁺ ions as well as at pH > 9 as [Hf(OH)₅]⁻ or [Hf(OH)₆]²⁻ complexes has been documented (Rai et al., 2001). Complexation of Zr and Hf with chloride and sodium-silicate has been found to favor 6-fold coordination of the central cation in experimental high-p-T studies (Louvel et al., 2013; Wilke et al., 2012; 2013) as well as with *ab initio* molecular dynamics (AIMD) modeling (Jahn et al., 2015).

In this study, we investigate the coordination environment of Zr and Hf in HF-rich aqueous solutions under conditions relevant to natural fluorine-rich hydrothermal systems such as greisen or pegmatite fluids (40 MPa, 100–400 °C, 1 mol/kg HF) by combining experimental X-ray absorption near edge structure (XANES) and extended X-ray absorption fine structure (EXAFS) data with *ab initio* molecular dynamics (AIMD) and classical molecular dynamics (MD). This comprehensive approach is necessary for determining the complex symmetry and ligand ratios in mixed complexes with both F and O in the same coordination shell because the very similar and overall low scattering cross section of these elements do not allow to extract all necessary information from X-ray absorption spectroscopy alone (Farges, 1991; 1996; Louvel et al., 2013). Therefore, EXAFS should not be used to interpret the number of ligands but only the average distance, to which this method is much more sensitive. XANES data can provide insight into the coordination of elements, whose observable X-ray absorption edges show features that can be demonstrated to be indicative of certain coordination symmetries. Most commonly, K edges are used for this purpose. Several experimental studies have demonstrated, using minerals and solutions, that the double peak of the Zr K white line is highly sensitive to coordination symmetry (and thus coordination number) as well as ligand type (e.g., Louvel et al., 2013; Wilke et al., 2012; 2013). For the Hf L₃ edge, no comprehensive data set is available and the data that are published suggest less obvious effects on the XANES region of this edge (e.g., Wilke et al., 2013). Computed XANES spectra of the Zr K and Hf L₃ edge corroborate these experimental findings (Jahn et al., 2015). Classical MD simulations are a very useful complement to XANES and EXAFS spectra analyzed with the help of AIMD. They are computationally less demanding and can therefore sample longer timescales. Here, we use classical MD simulations to simulate the structure and dynamics of the ligands and their exchange with the surrounding solvent between 25 and 400 °C in order to understand the changes in average ligand distance observed with EXAFS. Because classical MD simulations use a parametrized potential energy function, they cannot account for instantaneous changes in the electron cloud caused be changes in the surrounding electrostatic field (dynamic polarization). This drawback becomes acceptable for hard ions, i.e., those with high charge and small radius like the HFSE, whose electron cloud shows little response to changes in the surrounding electrostatic field. We adjust the parameters for the dispersion interactions of the cations to match the experimental data. We propose two possible mechanisms for Zr/Hf fractionation in hydrothermal systems based on our experimental and computational
results. Quantifying their respective influence remains work in progress at this time and has the potential to become an important diagnostic tool for petrologists, especially in combination with Zr stable isotope analysis.

## 2. Methods

### 2.1. Samples

Sintered solids of ZrO₂ and HfO₂, respectively, were equilibrated with 1 mol/kg HF for the two experimental runs. The solutions were acidified to pH = 2 (at room temperature) by adding 0.01 mol/kg HI to suppress hydrolysis and ensure that HF remained associated at low temperatures. Hydrogen iodide was chosen for acidification because Γ is not expected to form complexes with Zr or Hf (Pearson, 1963) and even if it did, the large back-scattering amplitude of iodine would make this immediately obvious in the XAS signal. The ZrO₂ solid was bought as oxide powder (99.7 % metals basis excluding Hf; Hf < 75 µg/g) from Alfa Aesar. The powder was pressed and sintered at 1000 °C in air for 24 h. The HfO₂ was bought as sintered lumps (99,9% metals basis excluding Zr; Zr < 0.5 %) and used directly.

### 2.2. Experimental set-up and XAS acquisition

The speciation of Zr and Hf in HF-bearing solutions was monitored from 25 to 400 °C at 40 MPa conducting X-ray absorption spectroscopy (XAS) measurements at P65 undulator beamline of the PETRA III storage ring, DESY, Hamburg (Welter et al., 2019), using a hydrothermal autoclave. All spectra were collected in transmission geometry. Details of the autoclave are described in Klemme et al. (2021). The aqueous samples resided in a vitreous carbon sample cell with 4 mm inner diameter placed inside a heater assembly and surrounded by a steel pressure vessel. The solids were loaded into the sample cell and sank to the bottom before the experiment. This way, they can interact with the fluids during the experiment but are outside of the beam path. Pressure is transmitted to the aqueous sample by free moving pistons using He as a pressure medium. The sample is heated with a resistive wire-heater (0.1 mm stainless steel wire) that completely surrounds the sample cell, except for the three spectroscopic windows (at 0°, 90°, and 180° relative to the direction of the X-ray beam). Two N-type thermocouples located directly above and below the spectroscopic windows were used to monitor temperature throughout the experiments. The offset between the temperature measured by the thermocouples and the actual temperature of the sample was calibrated by measuring the density of the cell loaded with pure water via total X-ray absorption at 9500 eV under 40 MPa pressure at 25, 100, 200, 300, and 400 °C, and comparing the results to the corresponding theoretical values (1014 to 523 kg/m³) taken from the NIST Chemistry WebBook (Lemmon et al., 2019). The value for the X-ray absorption cross-section of water was taken from Elam et al. (2002). The discrepancy between the measured and calculated temperature was 0 °C at 25 °C and increased linearly to 43 °C at 400 °C. A linear correction was consequently applied. Based on the magnitude of the applied correction and the associated uncertainties, the accuracy of the temperature measurements is estimated to be better than 5 °C at 25 °C and better than 15 °C at 400 °C.

XAS spectra were collected around the Zr K-edge (17998 eV) from 17,850 to 19000 eV and Hf L₃ edge (9561 eV) from 9410 to 10160 eV. The energies were scanned in continuous mode (monochromator and undulator gap) with a Si(311) and Si(111) monochromator crystal for Zr and Hf, respectively. Acquisition time per scan was ~ 5.5 and ~ 4 min, respectively. X-ray intensities in front of and behind the sample were recorded with ionization chambers. Four to 18 scans per temperature were collected and merged to improve signal to noise ratio, with the highest temperatures requiring the largest number of scans due to decreased density and solubility with increasing temperature. Spectra of the empty cell were taken prior to the high p-T measurements to check

for potential absorption edges of impurities in the ceramics of the heater assembly, cell material, or the beryllium windows but none were found in the energy ranges used for Zr or Hf.

### 2.3. EXAFS analysis

For the EXAFS fits, all individual spectra of a given composition at a given temperature were individually inspected and then merged. Analysis of the EXAFS data was performed using the xafsX software package (Winterer, 2020). Extraction of the k-space EXAFS data from the raw spectra was done in the following steps: 1) linear fit to pre-edge background; 2) fit of the sum of an error function and Gaussian to the near-edge region, 100 eV before and after the edge; 3) spline fit to the post-edge region. The EXAFS equation was then fitted to the extracted k-space data directly without filtering in R-space. The k-ranges to be fitted were chosen for each data set individually according to the data quality because the changing concentrations and densities meant that the spectra were much noisier at higher temperatures. A weighting of k³ was used throughout the data processing. Amplitude and phase for single scattering paths for the EXAFS fits were calculated with the FEFF 6.0 l software package's quick first shell algorithm (Rehr et al., 1991).

Two factors limit the EXAFS interpretation in our particular chemical system. Firstly, it is not possible to confidently distinguish F and O as ligands of an element under XAS investigation, due to their nearly identical amplitude and phase-shift functions (Farges, 1991; 1996, Louvel et al., 2013). Secondly, double-electron excitation effects (Filipponi, 1995) cause additional edges in the post-edge background for Zr K at 7.5 Å⁻¹ in k-space (Wang and Chen, 1998). We are unaware of direct descriptions of this effect on Hf L₃ but it has been described for the lanthanides (Chaboy et al., 1994; Ohta et al., 2008) and based on that is expected at around 7 Å⁻¹ in k-space (Roy et al., 2001). This affects the amplitude but not the frequency of the EXAFS oscillations and therefore limits the interpretation of the ligand type and coordination number significantly but interpretation of the radial distance of the ligands from the central atom is less affected. The double-electron excitation effect was not corrected. To account for these two complicating factors, we only fitted the radial distance and mean square displacement of a single shell with a fixed coordination number of 6, based on XANES results. Of the 6 ligands, two were F and four were O, based on published solubility studies (Ryzhenko et al., 2008; Migdisov et al., 2011).

### 2.4. AIMD and XANES calculations

Theoretical XANES spectra were derived from snapshots of AIMD, which were performed using the QUICKSTEP module of the CP2K code (Lippert et al., 1997; VandeVondele et al., 2005; Kühne et al., 2020). Electronic structure calculations were based on Kohn-Sham density functional theory using a mixed Gaussian and planewave approach, double zeta valence plus polarization basis sets (VandeVondele and Hutter, 2007) and Goedecker-Teter-Hutter pseudopotentials (Goedecker et al., 1996). The exchange-correlation functional was treated in the PBE parametrization of the generalized gradient approximation (Perdew et al., 1996) complemented by the D3 method to account for London dispersion interactions (Grimme et al., 2010; Jonchiere et al., 2011). Initial structures were adopted from results of an earlier AIMD study of zirconium chloride speciation (Jahn et al., 2015). The cubic simulation cells contained a total composition of 31 + xH₂O, 6-x HF, and 1 ZrFₓ(OH)₄ₓ or HfFₓ(OH)₄ₓ, where 1 ≤ x ≤ 4, i.e., different initial cells sampled Zr and Hf complexes with between one and four fluoride ligands. Additional simulations with a cell containing 59 H₂O and 1 ZrF₄ or HfF₄ did not show any significant finite size and composition effects concerning the local structures of interest here. AIMD simulations were performed in the canonical ensemble (NVT) at number of particles (N), constant temperature (T) and constant volume (V). At 350 °C, the edge length of the cubic simulation cell was set to 12.22 Å to reproduce the experimental pressure of 40 MPa, approximated by the equation of state
results. Quantifying their respective influence remains work in progress at this time and has the potential to become an important diagnostic tool for petrologists, especially in combination with Zr stable isotope analysis.

## 2. Methods

### 2.1. Samples

Sintered solids of ZrO₂ and HfO₂, respectively, were equilibrated with 1 mol/kg HF for the two experimental runs. The solutions were acidified to pH = 2 (at room temperature) by adding 0.01 mol/kg HI to suppress hydrolysis and ensure that HF remained associated at low temperatures. Hydrogen iodide was chosen for acidification because Γ is not expected to form complexes with Zr or Hf (Pearson, 1963) and even if it did, the large back-scattering amplitude of iodine would make this immediately obvious in the XAS signal. The ZrO₂ solid was bought as oxide powder (99.7 % metals basis excluding Hf; Hf < 75 µg/g) from Alfa Aesar. The powder was pressed and sintered at 1000 °C in air for 24 h. The HfO₂ was bought as sintered lumps (99,9% metals basis excluding Zr; Zr < 0.5 %) and used directly.

### 2.2. Experimental set-up and XAS acquisition

The speciation of Zr and Hf in HF-bearing solutions was monitored from 25 to 400 °C at 40 MPa conducting X-ray absorption spectroscopy (XAS) measurements at P65 undulator beamline of the PETRA III storage ring, DESY, Hamburg (Welter et al., 2019), using a hydrothermal autoclave. All spectra were collected in transmission geometry. Details of the autoclave are described in Klemme et al. (2021). The aqueous samples resided in a vitreous carbon sample cell with 4 mm inner diameter placed inside a heater assembly and surrounded by a steel pressure vessel. The solids were loaded into the sample cell and sank to the bottom before the experiment. This way, they can interact with the fluids during the experiment but are outside of the beam path. Pressure is transmitted to the aqueous sample by free moving pistons using He as a pressure medium. The sample is heated with a resistive wire-heater (0.1 mm stainless steel wire) that completely surrounds the sample cell, except for the three spectroscopic windows (at 0°, 90°, and 180° relative to the direction of the X-ray beam). Two N-type thermocouples located directly above and below the spectroscopic windows were used to monitor temperature throughout the experiments. The offset between the temperature measured by the thermocouples and the actual temperature of the sample was calibrated by measuring the density of the cell loaded with pure water via total X-ray absorption at 9500 eV under 40 MPa pressure at 25, 100, 200, 300, and 400 °C, and comparing the results to the corresponding theoretical values (1014 to 523 kg/m³) taken from the NIST Chemistry WebBook (Lemmon et al., 2019). The value for the X-ray absorption cross-section of water was taken from Elam et al. (2002). The discrepancy between the measured and calculated temperature was 0 °C at 25 °C and increased linearly to 43 °C at 400 °C. A linear correction was consequently applied. Based on the magnitude of the applied correction and the associated uncertainties, the accuracy of the temperature measurements is estimated to be better than 5 °C at 25 °C and better than 15 °C at 400 °C.

XAS spectra were collected around the Zr K-edge (17998 eV) from 17,850 to 19000 eV and Hf L₃ edge (9561 eV) from 9410 to 10160 eV. The energies were scanned in continuous mode (monochromator and undulator gap) with a Si(311) and Si(111) monochromator crystal for Zr and Hf, respectively. Acquisition time per scan was ~ 5.5 and ~ 4 min, respectively. X-ray intensities in front of and behind the sample were recorded with ionization chambers. Four to 18 scans per temperature were collected and merged to improve signal to noise ratio, with the highest temperatures requiring the largest number of scans due to decreased density and solubility with increasing temperature. Spectra of the empty cell were taken prior to the high p-T measurements to check

for potential absorption edges of impurities in the ceramics of the heater assembly, cell material, or the beryllium windows but none were found in the energy ranges used for Zr or Hf.

### 2.3. EXAFS analysis

For the EXAFS fits, all individual spectra of a given composition at a given temperature were individually inspected and then merged. Analysis of the EXAFS data was performed using the xafsX software package (Winterer, 2020). Extraction of the k-space EXAFS data from the raw spectra was done in the following steps: 1) linear fit to pre-edge background; 2) fit of the sum of an error function and Gaussian to the near-edge region, 100 eV before and after the edge; 3) spline fit to the post-edge region. The EXAFS equation was then fitted to the extracted k-space data directly without filtering in R-space. The k-ranges to be fitted were chosen for each data set individually according to the data quality because the changing concentrations and densities meant that the spectra were much noisier at higher temperatures. A weighting of k³ was used throughout the data processing. Amplitude and phase for single scattering paths for the EXAFS fits were calculated with the FEFF 6.0 l software package's quick first shell algorithm (Rehr et al., 1991).

Two factors limit the EXAFS interpretation in our particular chemical system. Firstly, it is not possible to confidently distinguish F and O as ligands of an element under XAS investigation, due to their nearly identical amplitude and phase-shift functions (Farges, 1991; 1996, Louvel et al., 2013). Secondly, double-electron excitation effects (Filipponi, 1995) cause additional edges in the post-edge background for Zr K at 7.5 Å⁻¹ in k-space (Wang and Chen, 1998). We are unaware of direct descriptions of this effect on Hf L₃ but it has been described for the lanthanides (Chaboy et al., 1994; Ohta et al., 2008) and based on that is expected at around 7 Å⁻¹ in k-space (Roy et al., 2001). This affects the amplitude but not the frequency of the EXAFS oscillations and therefore limits the interpretation of the ligand type and coordination number significantly but interpretation of the radial distance of the ligands from the central atom is less affected. The double-electron excitation effect was not corrected. To account for these two complicating factors, we only fitted the radial distance and mean square displacement of a single shell with a fixed coordination number of 6, based on XANES results. Of the 6 ligands, two were F and four were O, based on published solubility studies (Ryzhenko et al., 2008; Migdisov et al., 2011).

### 2.4. AIMD and XANES calculations

Theoretical XANES spectra were derived from snapshots of AIMD, which were performed using the QUICKSTEP module of the CP2K code (Lippert et al., 1997; VandeVondele et al., 2005; Kühne et al., 2020). Electronic structure calculations were based on Kohn-Sham density functional theory using a mixed Gaussian and planewave approach, double zeta valence plus polarization basis sets (VandeVondele and Hutter, 2007) and Goedecker-Teter-Hutter pseudopotentials (Goedecker et al., 1996). The exchange-correlation functional was treated in the PBE parametrization of the generalized gradient approximation (Perdew et al., 1996) complemented by the D3 method to account for London dispersion interactions (Grimme et al., 2010; Jonchiere et al., 2011). Initial structures were adopted from results of an earlier AIMD study of zirconium chloride speciation (Jahn et al., 2015). The cubic simulation cells contained a total composition of 31 + xH₂O, 6-x HF, and 1 ZrFₓ(OH)₄ₓ or HfFₓ(OH)₄ₓ, where 1 ≤ x ≤ 4, i.e., different initial cells sampled Zr and Hf complexes with between one and four fluoride ligands. Additional simulations with a cell containing 59 H₂O and 1 ZrF₄ or HfF₄ did not show any significant finite size and composition effects concerning the local structures of interest here. AIMD simulations were performed in the canonical ensemble (NVT) at number of particles (N), constant temperature (T) and constant volume (V). At 350 °C, the edge length of the cubic simulation cell was set to 12.22 Å to reproduce the experimental pressure of 40 MPa, approximated by the equation of state
<table><caption>**Table 1**

Setup for classical MD simulation boxes. Densities ρ from Wagner and Pruß (2002). L denotes length of the cubic simulation box, each containing 498 H₂O molecules. HF is the number of HF molecules per box.</caption><thead><tr><th>T [°C]</th><th>ρ [kg/m³]</th><th>L [Å]</th><th>HF</th></tr></thead><tbody><tr><td>25</td><td>1011.95</td><td>24.578</td><td>7</td></tr><tr><td>100</td><td>976.09</td><td>24.876</td><td>7</td></tr><tr><td>200</td><td>890.93</td><td>25.671</td><td>8</td></tr><tr><td>300</td><td>764.35</td><td>27.049</td><td>10</td></tr><tr><td>350</td><td>671.84</td><td>28.285</td><td>12</td></tr><tr><td>400</td><td>523.30</td><td>30.821</td><td>16</td></tr></tbody></table>

of pure H₂O (Wagner and Pruß, 2002). A few simulations were performed at 200 °C and 100 °C with a cell length of 11.12 Å to study temperature effects on bond lengths. To control temperature, the canonical sampling through velocity rescaling method was used (Bussi et al., 2007). The equations of motion were integrated with a time step of 0.5 fs. The initial configurations were equilibrated for a few picoseconds, followed by AIMD production runs of 25 ps. For each cell representing different Zr or Hf complexes, 51 snapshots were selected for XANES spectra calculations, i.e., the time between adjacent snapshots was 0.5 ps or 1000 AIMD time steps. As most of the initial Zr and Hf complexes persist over the time scale of the simulations, even if they are not thermodynamically stable, their local structural parameters as well as their theoretical XANES spectra can be readily compared.

XANES spectra simulations were performed using the FDMNES code (Joly, 2001; Bunău and Joly, 2009). This code uses a real space implementation of density functional theory and computes the electronic structure within clusters of a given radius around the absorbing atom. For the present study we used a cluster cutoff distance of 3.35 Å. Convergence of the cutoff was tested by performing a few simulations with cutoffs of 4.5 and 6.0 Å. The 51 individual spectra from one of the AIMD trajectories were averaged to obtain a spectrum representing the complex at the given finite temperature.

## 2.5. Classical MD

Classical molecular dynamics (MD) calculations were performed to investigate the geometric arrangement of the complexes and their development with temperature in detail. Classical MD allows to extract statistical radial distribution functions of ligands (in this case 2000 snapshots per temperature from 2 ns simulation time after equilibration). Systems resembling the starting conditions of the X-ray spectroscopic experiments at 25 °C, were simulated using the GROMACS package, version 2019.4 (Abraham et al., 2015). The parameters of M⁴⁺, F⁻ and OH⁻ were chosen as described below. The water model was TIP4P-2005, which is specifically parametrized for water at high temperatures and pressures (Abascal and Vega, 2005). For HF we used a 3-point model published Orabi and Faraldo-Gómez (2020). The boxes were set up to recreate a concentration of 1.0 mol/L HF, the density of water at 40 MPa and temperatures of 25, 100, 200, 300, and 400 °C. However, this was merely done to keep conditions similar to the experimental conditions, as the model for HF does not allow reactions with the metal complex. The densities were adjusted to the appropriate values taken from Wagner and Pruß (2002) by adjusting the size of the box. An overview over the relevant parameters is shown in Table 1. The systems were generated by placing one manually generated [MF₂(OH)₂·2H₂O] octahedral complex in a cubic box, together with 498 water molecules and the corresponding amount of HF molecules (7 to 16, see Table 1) to adjust to 1.0 mol/L HF. The fluoride anions in the complex had the same Lennard-Jones parameters as F in HF and a charge of -1. The hydroxide ions, were given a partial charge 0.32 e and -1.32 e for the hydrogen and oxygen-atom respectively, following the approach of Chen et al. (2013). The bond distance in the hydroxide was fixed to 0.9572 Å, the O-H bond-length of TIP4P-water. The oxygen had the same Lennard-Jones parameters as in H₂O in the TIP4P/2005 water

<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Raw XAS spectra of the Zr experiments.</figcaption></figure>

model (Abascal and Vega, 2005). The charge of the metal cation was set to +4. The choice of the Lennard-Jones parameters of the metal cation is particularly critical for simulating highly charged transition metal complexes. There are, to our knowledge, currently no published force field parameters for tetravalent cations available that are optimized or tested for use in high-temperature, high-pressure fluids or with fluoride. Therefore, a parameter scan was conducted to derive parameters that reproduce the average ligand distances observed in our experiments, using Steinhardt parameters (Steinhardt et al., 1983) to determine the coordination polyhedron of solvation shells (for details see Supplementary Material). Suitable parameter combinations all lie on a diagonal in semi-logarithmic parameter space (Fig. S1, Supplementary Material). Two sets of parameters, which reproduce the octahedral coordination sphere and the approximate average radial distance observed experimentally with EXAFS, were selected for the metal ion: 1) σ = 1.8 Å, ε = 20.92 kJ/mol; 2) σ = 3.6 Å, ε = 4.184·10⁻³ kJ/mol. The first parameter set was found using the parameter scan (Fig. S1), the second was extrapolated and subsequently tested and shown to reproduce the correct coordination and similar radial distance. These parameters represent a generic tetravalent cation, because within the chosen complexity of our force-field one cannot reliably distinguish between Zr⁴⁺ and Hf⁴⁺.

We used classical simulations to sample and analyze the configurational Boltzmann distribution. We will refer to the central cation in all classical MD results generically as M⁴⁺. Since the atomic mass does not influence the configurational Boltzmann distribution, the mass of the M⁴⁺ cation was set to 178.49 atomic mass units, the mass of hafnium, for all simulations.

In all simulations, the bond length of HF was constrained to the gas phase value of 0.917 Å. The time step was set to 1 fs for all runs. Electrostatics were treated using the Particle Mesh Ewald (PME) method (Darden et al., 1993), all real-space cutoffs were set to 10 Å. The PME-order was set to 6 with a Fourier spacing of 1 Å. Constraints were enforced using the LINCS-algorithm (Hess et al., 1997) with LINCS-order set to 4 and LINCS-iter set to 2. The LINCS-warnangle was set to 30°. Periodic boundary conditions in all directions were employed. After setting up the system, a steepest descent minimization was performed using a displacement step size per energy minimization of 0.01 Å, a maximum force of 0.01 kJ·mol⁻¹·Å⁻¹ and a maximum of 10⁵ steps. The
<table><caption>**Table 1**

Setup for classical MD simulation boxes. Densities ρ from Wagner and Pruß (2002). L denotes length of the cubic simulation box, each containing 498 H₂O molecules. HF is the number of HF molecules per box.</caption><thead><tr><th>T [°C]</th><th>ρ [kg/m³]</th><th>L [Å]</th><th>HF</th></tr></thead><tbody><tr><td>25</td><td>1011.95</td><td>24.578</td><td>7</td></tr><tr><td>100</td><td>976.09</td><td>24.876</td><td>7</td></tr><tr><td>200</td><td>890.93</td><td>25.671</td><td>8</td></tr><tr><td>300</td><td>764.35</td><td>27.049</td><td>10</td></tr><tr><td>350</td><td>671.84</td><td>28.285</td><td>12</td></tr><tr><td>400</td><td>523.30</td><td>30.821</td><td>16</td></tr></tbody></table>

of pure H₂O (Wagner and Pruß, 2002). A few simulations were performed at 200 °C and 100 °C with a cell length of 11.12 Å to study temperature effects on bond lengths. To control temperature, the canonical sampling through velocity rescaling method was used (Bussi et al., 2007). The equations of motion were integrated with a time step of 0.5 fs. The initial configurations were equilibrated for a few picoseconds, followed by AIMD production runs of 25 ps. For each cell representing different Zr or Hf complexes, 51 snapshots were selected for XANES spectra calculations, i.e., the time between adjacent snapshots was 0.5 ps or 1000 AIMD time steps. As most of the initial Zr and Hf complexes persist over the time scale of the simulations, even if they are not thermodynamically stable, their local structural parameters as well as their theoretical XANES spectra can be readily compared.

XANES spectra simulations were performed using the FDMNES code (Joly, 2001; Bunău and Joly, 2009). This code uses a real space implementation of density functional theory and computes the electronic structure within clusters of a given radius around the absorbing atom. For the present study we used a cluster cutoff distance of 3.35 Å. Convergence of the cutoff was tested by performing a few simulations with cutoffs of 4.5 and 6.0 Å. The 51 individual spectra from one of the AIMD trajectories were averaged to obtain a spectrum representing the complex at the given finite temperature.

## 2.5. Classical MD

Classical molecular dynamics (MD) calculations were performed to investigate the geometric arrangement of the complexes and their development with temperature in detail. Classical MD allows to extract statistical radial distribution functions of ligands (in this case 2000 snapshots per temperature from 2 ns simulation time after equilibration). Systems resembling the starting conditions of the X-ray spectroscopic experiments at 25 °C, were simulated using the GROMACS package, version 2019.4 (Abraham et al., 2015). The parameters of M⁴⁺, F⁻ and OH⁻ were chosen as described below. The water model was TIP4P-2005, which is specifically parametrized for water at high temperatures and pressures (Abascal and Vega, 2005). For HF we used a 3-point model published Orabi and Faraldo-Gómez (2020). The boxes were set up to recreate a concentration of 1.0 mol/L HF, the density of water at 40 MPa and temperatures of 25, 100, 200, 300, and 400 °C. However, this was merely done to keep conditions similar to the experimental conditions, as the model for HF does not allow reactions with the metal complex. The densities were adjusted to the appropriate values taken from Wagner and Pruß (2002) by adjusting the size of the box. An overview over the relevant parameters is shown in Table 1. The systems were generated by placing one manually generated [MF₂(OH)₂·2H₂O] octahedral complex in a cubic box, together with 498 water molecules and the corresponding amount of HF molecules (7 to 16, see Table 1) to adjust to 1.0 mol/L HF. The fluoride anions in the complex had the same Lennard-Jones parameters as F in HF and a charge of -1. The hydroxide ions, were given a partial charge 0.32 e and -1.32 e for the hydrogen and oxygen-atom respectively, following the approach of Chen et al. (2013). The bond distance in the hydroxide was fixed to 0.9572 Å, the O-H bond-length of TIP4P-water. The oxygen had the same Lennard-Jones parameters as in H₂O in the TIP4P/2005 water

<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Raw XAS spectra of the Zr experiments.</figcaption></figure>

model (Abascal and Vega, 2005). The charge of the metal cation was set to +4. The choice of the Lennard-Jones parameters of the metal cation is particularly critical for simulating highly charged transition metal complexes. There are, to our knowledge, currently no published force field parameters for tetravalent cations available that are optimized or tested for use in high-temperature, high-pressure fluids or with fluoride. Therefore, a parameter scan was conducted to derive parameters that reproduce the average ligand distances observed in our experiments, using Steinhardt parameters (Steinhardt et al., 1983) to determine the coordination polyhedron of solvation shells (for details see Supplementary Material). Suitable parameter combinations all lie on a diagonal in semi-logarithmic parameter space (Fig. S1, Supplementary Material). Two sets of parameters, which reproduce the octahedral coordination sphere and the approximate average radial distance observed experimentally with EXAFS, were selected for the metal ion: 1) σ = 1.8 Å, ε = 20.92 kJ/mol; 2) σ = 3.6 Å, ε = 4.184·10⁻³ kJ/mol. The first parameter set was found using the parameter scan (Fig. S1), the second was extrapolated and subsequently tested and shown to reproduce the correct coordination and similar radial distance. These parameters represent a generic tetravalent cation, because within the chosen complexity of our force-field one cannot reliably distinguish between Zr⁴⁺ and Hf⁴⁺.

We used classical simulations to sample and analyze the configurational Boltzmann distribution. We will refer to the central cation in all classical MD results generically as M⁴⁺. Since the atomic mass does not influence the configurational Boltzmann distribution, the mass of the M⁴⁺ cation was set to 178.49 atomic mass units, the mass of hafnium, for all simulations.

In all simulations, the bond length of HF was constrained to the gas phase value of 0.917 Å. The time step was set to 1 fs for all runs. Electrostatics were treated using the Particle Mesh Ewald (PME) method (Darden et al., 1993), all real-space cutoffs were set to 10 Å. The PME-order was set to 6 with a Fourier spacing of 1 Å. Constraints were enforced using the LINCS-algorithm (Hess et al., 1997) with LINCS-order set to 4 and LINCS-iter set to 2. The LINCS-warnangle was set to 30°. Periodic boundary conditions in all directions were employed. After setting up the system, a steepest descent minimization was performed using a displacement step size per energy minimization of 0.01 Å, a maximum force of 0.01 kJ·mol⁻¹·Å⁻¹ and a maximum of 10⁵ steps. The
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** XANES spectra of Zr K and Hf L₃, normalized to equal edge step height. Spectra of Zr in solid zircon, baddeleyite, eudialyte, and vlasovite with the respective coordination numbers of Zr in square brackets are shown as dashed gray lines. The data for the solid spectra were taken from Louvel et al. (2013).</figcaption></figure>

minimized system was simulated for 3 ns in the NVT-ensemble using the velocity-rescale thermostat ($\tau_{T} = 0.1$ ps) and the leap-frog integrator. The first nanosecond was discarded as equilibration time. At the beginning of all NVT-runs, the velocities were drawn from the Maxwell

distributions of the target temperatures.

### 3. Results

#### 3.1. XAS raw data and element aqueous concentrations

The K-edge X-ray absorption spectra for Zr in 1 mol/kg HF for sample temperatures from 100 °C to 350 °C are presented in Fig. 1. The continuous decrease in pre-edge absorption intensity with increasing temperature is due to the decrease in density of the solution from 976 kg/m³ (at 100 °C) to 672 kg/m³ (at 350 °C) at a constant pressure of 40 MPa (NIST Chemistry WebBook, Lemmon et al., 2019). For Hf L₃, spectra for sample temperatures from 100 °C to 400 °C were recorded. Due to the high absorption cross section of Hf and significant concentration change with temperature, the Keithley amplifier gain had to be adjusted between temperatures. Therefore, the Hf data are presented in normalized form only (Fig. 2) instead of raw counts. The high intensity of the Hf L₃ white line lead to total absorption of all X-ray photons in the sample at the white line energies around 9563 eV (and therefore zero counts on the I₁ detector) for the isotherm with the highest concentration of Hf, 200 °C (Fig. S2). Therefore, no XANES data exist for this data point (Fig. 2). The EXAFS data are not affected by this because the absorption in the EXAFS region is much weaker than at the white line for Hf L₃. The edge jump was fitted for all recorded spectra individually by fitting a combined error function and Gaussian from the beginning of the spectrum up to 200 eV after the edge. This provides a measure of the concentration of the element in the solutions under investigation and allows to assess whether a steady state was reached at a particular temperature. Using the elemental X-ray mass absorption coefficients (Hubbell and Seltzer, 1996) and the water density at the respective p-T-conditions (Lemmon et al., 2019), element concentrations can be calculated from the height of the edge jumps (for method see Pokrovski et al., 2005). The concentrations and edge jumps are plotted in Fig. S2 For Zr, the 200 °C isotherm converged around 270 µg/g, 250 °C around 220 µg/g and 350 °C around 33 µg/g. The other isotherms did not converge before the temperature was increased. For Hf, only two isotherms converged, 250 °C around 190 µg/g and 400 °C around 24 µg/g. Note also that the density values of the fluid were calculated assuming pure water due to the lack of suitable data for hydrofluoric acid and that the density of the actual fluid is slightly higher due to the dissolved HF and Zr or Hf. However, this error is estimated to be less than 1 % for the compositions in question. The nature of a solid phase was not monitored during the experiments and therefore these data may not necessarily represent solubilities of baddeleyite-type crystalline ZrO₂ and HfO₂ and should not be regarded as such. The purpose of monitoring the apparent concentration was mainly to ensure that no significant precipitation of solid phases has taken place in the beam, which would interfere with the XANES and EXAFS signals from the dissolved species. Experience from past experiments shows that when solid phases precipitate in the beam path, they do so at the walls of the sample chamber because this surface provides nucleation sites, thus lowering the nucleation energy barrier. Spontaneous nucleation in the free fluid instead of at the cell wall would have to overcome this energy barrier unaided and is thus highly unlikely. If precipitation at the cell wall in the beam path were to occur, one would expect a continuous increase of apparent total concentration over time. Because the recorded XANES and EXAFS spectra would then represent a mixture of signal from the dissolved species and the solids (ZrO₂ or HfO₂, in this case), one would also expect a shift in the shape of the spectra due to the continual increase of the proportion of the solid over time as more and more solid precipitates. No such behavior was observed and we therefore conclude that our XANES and EXAFS data represent the dissolved species and not precipitates. However, solid precipitate as well as remnants of the originally added oxide chunks were recovered after the experiments, indicating that the solutions were in continuous contact and chemical exchange with the solids.
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** XANES spectra of Zr K and Hf L₃, normalized to equal edge step height. Spectra of Zr in solid zircon, baddeleyite, eudialyte, and vlasovite with the respective coordination numbers of Zr in square brackets are shown as dashed gray lines. The data for the solid spectra were taken from Louvel et al. (2013).</figcaption></figure>

minimized system was simulated for 3 ns in the NVT-ensemble using the velocity-rescale thermostat ($\tau_{T} = 0.1$ ps) and the leap-frog integrator. The first nanosecond was discarded as equilibration time. At the beginning of all NVT-runs, the velocities were drawn from the Maxwell

distributions of the target temperatures.

### 3. Results

#### 3.1. XAS raw data and element aqueous concentrations

The K-edge X-ray absorption spectra for Zr in 1 mol/kg HF for sample temperatures from 100 °C to 350 °C are presented in Fig. 1. The continuous decrease in pre-edge absorption intensity with increasing temperature is due to the decrease in density of the solution from 976 kg/m³ (at 100 °C) to 672 kg/m³ (at 350 °C) at a constant pressure of 40 MPa (NIST Chemistry WebBook, Lemmon et al., 2019). For Hf L₃, spectra for sample temperatures from 100 °C to 400 °C were recorded. Due to the high absorption cross section of Hf and significant concentration change with temperature, the Keithley amplifier gain had to be adjusted between temperatures. Therefore, the Hf data are presented in normalized form only (Fig. 2) instead of raw counts. The high intensity of the Hf L₃ white line lead to total absorption of all X-ray photons in the sample at the white line energies around 9563 eV (and therefore zero counts on the I₁ detector) for the isotherm with the highest concentration of Hf, 200 °C (Fig. S2). Therefore, no XANES data exist for this data point (Fig. 2). The EXAFS data are not affected by this because the absorption in the EXAFS region is much weaker than at the white line for Hf L₃. The edge jump was fitted for all recorded spectra individually by fitting a combined error function and Gaussian from the beginning of the spectrum up to 200 eV after the edge. This provides a measure of the concentration of the element in the solutions under investigation and allows to assess whether a steady state was reached at a particular temperature. Using the elemental X-ray mass absorption coefficients (Hubbell and Seltzer, 1996) and the water density at the respective p-T-conditions (Lemmon et al., 2019), element concentrations can be calculated from the height of the edge jumps (for method see Pokrovski et al., 2005). The concentrations and edge jumps are plotted in Fig. S2 For Zr, the 200 °C isotherm converged around 270 µg/g, 250 °C around 220 µg/g and 350 °C around 33 µg/g. The other isotherms did not converge before the temperature was increased. For Hf, only two isotherms converged, 250 °C around 190 µg/g and 400 °C around 24 µg/g. Note also that the density values of the fluid were calculated assuming pure water due to the lack of suitable data for hydrofluoric acid and that the density of the actual fluid is slightly higher due to the dissolved HF and Zr or Hf. However, this error is estimated to be less than 1 % for the compositions in question. The nature of a solid phase was not monitored during the experiments and therefore these data may not necessarily represent solubilities of baddeleyite-type crystalline ZrO₂ and HfO₂ and should not be regarded as such. The purpose of monitoring the apparent concentration was mainly to ensure that no significant precipitation of solid phases has taken place in the beam, which would interfere with the XANES and EXAFS signals from the dissolved species. Experience from past experiments shows that when solid phases precipitate in the beam path, they do so at the walls of the sample chamber because this surface provides nucleation sites, thus lowering the nucleation energy barrier. Spontaneous nucleation in the free fluid instead of at the cell wall would have to overcome this energy barrier unaided and is thus highly unlikely. If precipitation at the cell wall in the beam path were to occur, one would expect a continuous increase of apparent total concentration over time. Because the recorded XANES and EXAFS spectra would then represent a mixture of signal from the dissolved species and the solids (ZrO₂ or HfO₂, in this case), one would also expect a shift in the shape of the spectra due to the continual increase of the proportion of the solid over time as more and more solid precipitates. No such behavior was observed and we therefore conclude that our XANES and EXAFS data represent the dissolved species and not precipitates. However, solid precipitate as well as remnants of the originally added oxide chunks were recovered after the experiments, indicating that the solutions were in continuous contact and chemical exchange with the solids.
<caption>**Table 2**</caption>

Cation-ligand distances R from AIMD calculations at different temperatures along with the asymmetry factor α, calculated as the longest bond length of a given complex divided by the shortest.

<table><thead><tr><th>T [°C]</th><th>complex</th><th>R(Zr-F) [Å]</th><th>R(Zr-OH) [Å]</th><th>R(Zr-H<sub>2</sub>O) [Å]</th><th>α</th><th>R(Hf-F) [Å]</th><th>R(Hf-OH) [Å]</th><th>R(Hf-H<sub>2</sub>O) [Å]</th><th>α</th></tr></thead><tbody><tr><td>100</td><td>[MF<sub>4</sub>(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>1.97(1)</td><td>–</td><td>2.25(3)</td><td>1.14</td><td>1.95(1)</td><td>–</td><td>2.18(3)</td><td>1.11</td></tr><tr><td>200</td><td>[MF<sub>4</sub>(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>1.96(1)</td><td>–</td><td>2.26(1)</td><td>1.15</td><td>1.94(1)</td><td>–</td><td>2.22(2)</td><td>1.14</td></tr><tr><td rowspan="3">350</td><td>[MF<sub>4</sub>(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>1.96(1)</td><td>–</td><td>2.30(2)</td><td>1.17</td><td>1.93(1)</td><td>–</td><td>2.23(1)</td><td>1.16</td></tr><tr><td>[MF<sub>3</sub>(OH)(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>1.97(1)</td><td>1.97(1)</td><td>2.33(2)</td><td>1.18</td><td>1.94(1)</td><td>1.93(1)</td><td>2.27(3)</td><td>1.18</td></tr><tr><td>[MF<sub>2</sub>(OH)<sub>2</sub>(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>1.99(1)</td><td>1.96(1)</td><td>2.32(2)</td><td>1.18</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td></td><td>[MF(OH)<sub>3</sub>(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>–</td><td>–</td><td>–</td><td>–</td><td>1.96(1)</td><td>1.94(2)</td><td>2.24(4)</td><td>1.15</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Calculated XANES spectra for Zr and Hf hydroxy-fluoride complexes and the corresponding experimental spectra. Complexes with different F/OH ratios at 350 °C show little difference in line shape. The Hf spectra show a rather more pronounced difference between 100 °C and 350 °C, compared to the Zr spectra. The calculated spectra are calibrated to the same edge positions as the experimental spectra. All intensities are normalized for easier comparison.</figcaption></figure>

## 3.2. XANES and AIMD calculations

The intensity normalized XANES regions of Zr K and Hf L3 edges are shown in Fig. 2. The white line of the Zr K edge in our data consists of two distinct peaks at ca. 18015 eV and 18030 eV, with the latter having the higher intensity. This shape is typical for 6-fold coordination of Zr with O and clearly distinct from 7- or 8-fold coordination, in which the lower energy peak is more intense (Louvel et al., 2013; Wilke et al., 2012; 2013). Due to the similarity of F and O, we would not expect the exchange of one or several O with F in the coordination sphere to have a pronounced influence on the XANES spectrum. Any potential change in coordination is not obvious from the Zr K edge data but a broadening of the features becomes apparent at higher temperatures. The Hf L3 white line consists of two distinct peaks at 100 °C. The 200 °C data for Hf are not included because the white line absorbed all photons in this instance, rendering those XANES data unusable. At 250 to 400 °C, the splitting of the white line is not observed.

Structural analysis of the AIMD simulations showed that all complexes up to 350 °C prefer 6-fold coordination. Transient deviations to lower coordination were only observed in few cases where ligand

exchange took place. Note that the main purpose of the AIMD simulations here was to generate reasonable structures for calculation of theoretical XANES spectra rather than obtaining information on the thermodynamic stability of the different complexes. The average radial distances of the ligands to the central Zr⁴⁺ and Hf⁴⁺ ions were obtained from the maxima of the respective radial distribution functions and are summarized in Table 2. The Zr/Hf-F and Zr/Hf-OH distances are between 1.93 Å and 1.99 Å. While Zr-F/OH distances are systematically longer than Hf-F/OH distances by 0.02–0.05 Å, no systematic difference in bond length is observed between the metal cations and fluoride or hydroxyl ligands. There seems to be a small decrease of these distances with increasing temperature but the changes are of similar size as the errors. The cation distances to H₂O ligands were consistently larger, i.e. between 2.18 and 2.33 Å, with a significant positive temperature dependence. This leads to an overall increase of asymmetry with increasing temperature (Table 2).

Derived structures for the equilibrated Zr and Hf complexes at 100 °C and 350 °C were then used to calculate theoretical XANES spectra. These spectra are shown in Fig. 3 together with experimental spectra collected in 1 mol/kg HF solution at similar p-T conditions. The calculated spectra
<caption>**Table 2**</caption>

Cation-ligand distances R from AIMD calculations at different temperatures along with the asymmetry factor α, calculated as the longest bond length of a given complex divided by the shortest.

<table><thead><tr><th>T [°C]</th><th>complex</th><th>R(Zr-F) [Å]</th><th>R(Zr-OH) [Å]</th><th>R(Zr-H<sub>2</sub>O) [Å]</th><th>α</th><th>R(Hf-F) [Å]</th><th>R(Hf-OH) [Å]</th><th>R(Hf-H<sub>2</sub>O) [Å]</th><th>α</th></tr></thead><tbody><tr><td>100</td><td>[MF<sub>4</sub>(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>1.97(1)</td><td>–</td><td>2.25(3)</td><td>1.14</td><td>1.95(1)</td><td>–</td><td>2.18(3)</td><td>1.11</td></tr><tr><td>200</td><td>[MF<sub>4</sub>(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>1.96(1)</td><td>–</td><td>2.26(1)</td><td>1.15</td><td>1.94(1)</td><td>–</td><td>2.22(2)</td><td>1.14</td></tr><tr><td rowspan="3">350</td><td>[MF<sub>4</sub>(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>1.96(1)</td><td>–</td><td>2.30(2)</td><td>1.17</td><td>1.93(1)</td><td>–</td><td>2.23(1)</td><td>1.16</td></tr><tr><td>[MF<sub>3</sub>(OH)(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>1.97(1)</td><td>1.97(1)</td><td>2.33(2)</td><td>1.18</td><td>1.94(1)</td><td>1.93(1)</td><td>2.27(3)</td><td>1.18</td></tr><tr><td>[MF<sub>2</sub>(OH)<sub>2</sub>(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>1.99(1)</td><td>1.96(1)</td><td>2.32(2)</td><td>1.18</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td></td><td>[MF(OH)<sub>3</sub>(H<sub>2</sub>O)<sub>2</sub>]<sup>0</sup></td><td>–</td><td>–</td><td>–</td><td>–</td><td>1.96(1)</td><td>1.94(2)</td><td>2.24(4)</td><td>1.15</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Calculated XANES spectra for Zr and Hf hydroxy-fluoride complexes and the corresponding experimental spectra. Complexes with different F/OH ratios at 350 °C show little difference in line shape. The Hf spectra show a rather more pronounced difference between 100 °C and 350 °C, compared to the Zr spectra. The calculated spectra are calibrated to the same edge positions as the experimental spectra. All intensities are normalized for easier comparison.</figcaption></figure>

## 3.2. XANES and AIMD calculations

The intensity normalized XANES regions of Zr K and Hf L3 edges are shown in Fig. 2. The white line of the Zr K edge in our data consists of two distinct peaks at ca. 18015 eV and 18030 eV, with the latter having the higher intensity. This shape is typical for 6-fold coordination of Zr with O and clearly distinct from 7- or 8-fold coordination, in which the lower energy peak is more intense (Louvel et al., 2013; Wilke et al., 2012; 2013). Due to the similarity of F and O, we would not expect the exchange of one or several O with F in the coordination sphere to have a pronounced influence on the XANES spectrum. Any potential change in coordination is not obvious from the Zr K edge data but a broadening of the features becomes apparent at higher temperatures. The Hf L3 white line consists of two distinct peaks at 100 °C. The 200 °C data for Hf are not included because the white line absorbed all photons in this instance, rendering those XANES data unusable. At 250 to 400 °C, the splitting of the white line is not observed.

Structural analysis of the AIMD simulations showed that all complexes up to 350 °C prefer 6-fold coordination. Transient deviations to lower coordination were only observed in few cases where ligand

exchange took place. Note that the main purpose of the AIMD simulations here was to generate reasonable structures for calculation of theoretical XANES spectra rather than obtaining information on the thermodynamic stability of the different complexes. The average radial distances of the ligands to the central Zr⁴⁺ and Hf⁴⁺ ions were obtained from the maxima of the respective radial distribution functions and are summarized in Table 2. The Zr/Hf-F and Zr/Hf-OH distances are between 1.93 Å and 1.99 Å. While Zr-F/OH distances are systematically longer than Hf-F/OH distances by 0.02–0.05 Å, no systematic difference in bond length is observed between the metal cations and fluoride or hydroxyl ligands. There seems to be a small decrease of these distances with increasing temperature but the changes are of similar size as the errors. The cation distances to H₂O ligands were consistently larger, i.e. between 2.18 and 2.33 Å, with a significant positive temperature dependence. This leads to an overall increase of asymmetry with increasing temperature (Table 2).

Derived structures for the equilibrated Zr and Hf complexes at 100 °C and 350 °C were then used to calculate theoretical XANES spectra. These spectra are shown in Fig. 3 together with experimental spectra collected in 1 mol/kg HF solution at similar p-T conditions. The calculated spectra
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** EXAFS spectra of the experiments with fits in k-space. For data treatment and fitting procedure see text.</figcaption></figure>

reproduce the general shape of the white line of the experimental data well, although the lines are somewhat broadened in the calculated spectra. The height ratio of the first and second peak of the white line in the Zr spectra changes slightly with the F/OH ratio of the complex, with higher ratios reproducing the experimental data more closely. However, the correlation is too weak to confidently use it to determine the true F/OH ratio in the experimental data. The splitting of the Hf L₃ white line at 100 °C in the experiments is reproduced by the simulated spectrum at the same temperature. At 350 °C, this split is only visible as a shoulder in the simulated spectra and no longer clearly identifiable in the experimental spectrum. This is likely an effect of thermal broadening of the features.

<figure><img src="image_9.png" /><figcaption>**Fig. 5.** EXAFS spectra of the experiments in R-space. For data treatment see text.</figcaption></figure>

### 3.3. EXAFS

The EXAFS data of both Zr and Hf experiments show a single oscillation frequency in k-space at all temperatures, indicating a single dominant coordination shell (Fig. 4). The Fourier-transformed spectra in R-space also show no evidence of a well-defined second shell (Fig. 5). A deviation between fit and data of the Zr data can be seen at 7.5 Å⁻¹ in k-space, particularly in the 100, 200, and 250 °C data sets (Fig. 4). The position and appearance of this feature correspond to published observations of the dual-electron excitation effect (e.g., Filipponi, 1995; Wang and Chen, 1998). This affects the apparent amplitude of the oscillation before and after this secondary edge. A similar but weaker edge is observed in the Hf data around 7 Å⁻¹ (Fig. 4). However, the edge caused by multi-electron excitation does not affect the interpretation of the
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** EXAFS spectra of the experiments with fits in k-space. For data treatment and fitting procedure see text.</figcaption></figure>

reproduce the general shape of the white line of the experimental data well, although the lines are somewhat broadened in the calculated spectra. The height ratio of the first and second peak of the white line in the Zr spectra changes slightly with the F/OH ratio of the complex, with higher ratios reproducing the experimental data more closely. However, the correlation is too weak to confidently use it to determine the true F/OH ratio in the experimental data. The splitting of the Hf L₃ white line at 100 °C in the experiments is reproduced by the simulated spectrum at the same temperature. At 350 °C, this split is only visible as a shoulder in the simulated spectra and no longer clearly identifiable in the experimental spectrum. This is likely an effect of thermal broadening of the features.

<figure><img src="image_9.png" /><figcaption>**Fig. 5.** EXAFS spectra of the experiments in R-space. For data treatment see text.</figcaption></figure>

### 3.3. EXAFS

The EXAFS data of both Zr and Hf experiments show a single oscillation frequency in k-space at all temperatures, indicating a single dominant coordination shell (Fig. 4). The Fourier-transformed spectra in R-space also show no evidence of a well-defined second shell (Fig. 5). A deviation between fit and data of the Zr data can be seen at 7.5 Å⁻¹ in k-space, particularly in the 100, 200, and 250 °C data sets (Fig. 4). The position and appearance of this feature correspond to published observations of the dual-electron excitation effect (e.g., Filipponi, 1995; Wang and Chen, 1998). This affects the apparent amplitude of the oscillation before and after this secondary edge. A similar but weaker edge is observed in the Hf data around 7 Å⁻¹ (Fig. 4). However, the edge caused by multi-electron excitation does not affect the interpretation of the
<figure><img src="image_10.png" /><figcaption>Fig. 6. a) Mean radial distances (R) of ligands around Zr and Hf, calculated from EXAFS data, plotted against inverse absolute temperature. The dashed trend lines are linear fits with their associated statistical R² value (not to be confused with the radial distance R). The error bars represent random relative experimental uncertainty between individual points and are estimated from the maximum residuals between the data points and the fits. b) Total displacement range around the mean values (R ± σ) as calculated from EXAFS data (not related to experimental error) for all experimental temperatures. Zr data are offset in the graph by -5 °C and Hf by +5 °C for visual clarity. The total displacement includes statistical displacement (difference between individual ligands, i.e. asymmetry of the complex) as well as thermal displacement (change in R over time due to thermal vibration). Mind the difference in y-axis scaling between plots a) and b).</figcaption></figure>

radial distances, when a fixed coordination number is chosen (Fig. S3). The distortion of the amplitude at k-values below ca. 5.5 Å⁻¹ in the Hf data most likely stems from the high-energy end of the XANES effects and the background correction but did not significantly affect the radial ligand distances extracted from the processed EXAFS data at k > 3 Å⁻¹, which was chosen as the lower bound of the range interpreted by the EXAFS fitting. Adjusting the data range over which the near edge region is fitted influences the effect on the amplitude, making interpretation of this parameter more reflective of the fitting procedure than of the actual data. It was therefore decided to use a fixed range of ± 100 eV around the edge for near edge fitting for better comparability of the data.

Due to the double electron excitation effect and the uncertainties in background subtraction, the significance of the EXAFS amplitude for interpretation of complexes is limited. Because the amplitude reduction term S₀² is unknown and cannot be derived from fits independently of the coordination number N, with which it is correlated, we will assume S₀² = 1 for all fits.

Previous studies have found that ZrF₂(OH)₂ and HfF₂(OH)₂ complexes are expected to be dominant under our experimental conditions (Rhyzhenko et al., 2008; Migdisov et al., 2011). However, these studies both used methods that cannot observe potential additional molecular water in the first hydration shell. Based on the comparison of the XANES spectra with AIMD modeling, the complexes were assumed to have close to six-fold coordination but not perfect octahedral symmetry. However, the systematic differences in R between the three different ligand sub-shells were too small to be resolved with EXAFS fitting, because the scattering cross sections of F and O are too similar to each other. Attempts to fit two shells with separately optimized R parameters resulted in inconsistent results between the various experiments and strongly correlated parameters. Such results indicate that the data were over-fitted in these attempts, meaning too many parameters were fitted to the available amount of experimental information. It is therefore important to note that the radial distances R for all fits presented in Fig. 6a are mean values for all inner ligands, with the σ parameter (Fig. 6b) accounting for the overall mean displacement between all ligands, including statistical and thermal displacement. The ligands are, by nature of the EXAFS equation, weighted in the analysis with the inverse square of their radial distance (R⁻²). Based on the considerations above, the complexes ZrF₂(OH)₂·2H₂O and HfF₂(OH)₂·2H₂O were fitted to the Zr and Hf data, respectively, by defining a single ligand shell containing 2F and 4O at an identical distance. The radial distance R and displacement σ were refined, the coordination number was kept fixed.

<caption>**Table 3**</caption>

EXAFS fitting results as reported by xafsX. No error can be calculated for the radial distance R because it is totally correlated with the mean displacement σ.

<table><thead><tr><th></th><th>T [°C]</th><th>R [Å]</th><th>σ² [Å²]</th><th>±</th><th>E₀ [eV]</th><th>±</th></tr></thead><tbody><tr><td rowspan="5">Zr</td><td>100</td><td>2.008</td><td>0.004</td><td>0.001</td><td>-2.952</td><td>0.0001</td></tr><tr><td>200</td><td>1.990</td><td>0.005</td><td>0.001</td><td>-3.331</td><td>0.0001</td></tr><tr><td>250</td><td>1.990</td><td>0.005</td><td>0.002</td><td>-4.810</td><td>0.0002</td></tr><tr><td>300</td><td>1.994</td><td>0.006</td><td>0.005</td><td>-4.495</td><td>0.0003</td></tr><tr><td>350</td><td>1.974</td><td>0.005</td><td>0.008</td><td>-2.852</td><td>0.0006</td></tr><tr><td rowspan="5">Hf</td><td>100</td><td>2.013</td><td>0.003</td><td>0.002</td><td>0.190</td><td>0.0001</td></tr><tr><td>200</td><td>1.975</td><td>0.004</td><td>0.002</td><td>3.850</td><td>0.0002</td></tr><tr><td>250</td><td>1.970</td><td>0.005</td><td>0.002</td><td>0.791</td><td>0.0002</td></tr><tr><td>300</td><td>1.969</td><td>0.005</td><td>0.004</td><td>0.573</td><td>0.0003</td></tr><tr><td>350</td><td>1.965</td><td>0.006</td><td>0.007</td><td>0.960</td><td>0.0005</td></tr><tr><td></td><td>400</td><td>1.969</td><td>0.01</td><td>0.01</td><td>1.473</td><td>0.001</td></tr></tbody></table>

The results are summarized in Table 3. The nature of the fitting procedure precludes the independent determination of an experimental error or uncertainty on R in addition to σ for individual data points. It is, however, possible to estimate if changes in R with temperature are significant by taking into account the R values for all temperatures for either the Zr or Hf experiments. Temperature dependent properties tend to follow simple trends when plotted against reciprocal absolute temperature. The simplest such dependency being linear, thus providing the most conservative error estimate. We can fit a linear function in reciprocal absolute temperature space to the data (Fig. 6a) and use the largest observed residuals between the data points and this fit as an estimate of the relative experimental uncertainty between the individual data points. These estimated uncertainties of 0.01 Å are shown as error bars in Fig. 6a. It is worth noting that this is an estimate for a purely random relative error between the data points at different temperatures. There is no systematic temperature dependent experimental error for XANES and EXAFS data because no temperature dependent parameters are used in the interpretation of the data and all of the spectroscopic equipment is far away from the autoclave and at a stable temperature. We can, therefore, be confident that any differences beyond the error bars in Fig. 6a are significant and the observed apparent contraction of the complexes with increasing temperature is real. Because no variables that could introduce a systematic error between Zr and Hf are used in the EXAFS analysis and considering that all individual R values of Zr complexes at ≥ 200 °C are larger than the corresponding values of Hf complexes, we would argue that there is also a significant difference
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** a) Mean radial distances (R) of ligands around Zr and Hf, calculated from EXAFS data, plotted against inverse absolute temperature. The dashed trend lines are linear fits with their associated statistical R² value (not to be confused with the radial distance R). The error bars represent random relative experimental uncertainty between individual points and are estimated from the maximum residuals between the data points and the fits. b) Total displacement range around the mean values (R ± σ) as calculated from EXAFS data (not related to experimental error) for all experimental temperatures. Zr data are offset in the graph by -5 °C and Hf by +5 °C for visual clarity. The total displacement includes statistical displacement (difference between individual ligands, i.e. asymmetry of the complex) as well as thermal displacement (change in R over time due to thermal vibration). Mind the difference in y-axis scaling between plots a) and b).</figcaption></figure>

radial distances, when a fixed coordination number is chosen (Fig. S3). The distortion of the amplitude at k-values below ca. 5.5 Å⁻¹ in the Hf data most likely stems from the high-energy end of the XANES effects and the background correction but did not significantly affect the radial ligand distances extracted from the processed EXAFS data at k > 3 Å⁻¹, which was chosen as the lower bound of the range interpreted by the EXAFS fitting. Adjusting the data range over which the near edge region is fitted influences the effect on the amplitude, making interpretation of this parameter more reflective of the fitting procedure than of the actual data. It was therefore decided to use a fixed range of ± 100 eV around the edge for near edge fitting for better comparability of the data.

Due to the double electron excitation effect and the uncertainties in background subtraction, the significance of the EXAFS amplitude for interpretation of complexes is limited. Because the amplitude reduction term S₀² is unknown and cannot be derived from fits independently of the coordination number N, with which it is correlated, we will assume S₀² = 1 for all fits.

Previous studies have found that ZrF₂(OH)₂ and HfF₂(OH)₂ complexes are expected to be dominant under our experimental conditions (Rhyzhenko et al., 2008; Migdisov et al., 2011). However, these studies both used methods that cannot observe potential additional molecular water in the first hydration shell. Based on the comparison of the XANES spectra with AIMD modeling, the complexes were assumed to have close to six-fold coordination but not perfect octahedral symmetry. However, the systematic differences in R between the three different ligand sub-shells were too small to be resolved with EXAFS fitting, because the scattering cross sections of F and O are too similar to each other. Attempts to fit two shells with separately optimized R parameters resulted in inconsistent results between the various experiments and strongly correlated parameters. Such results indicate that the data were over-fitted in these attempts, meaning too many parameters were fitted to the available amount of experimental information. It is therefore important to note that the radial distances R for all fits presented in Fig. 6a are mean values for all inner ligands, with the σ parameter (Fig. 6b) accounting for the overall mean displacement between all ligands, including statistical and thermal displacement. The ligands are, by nature of the EXAFS equation, weighted in the analysis with the inverse square of their radial distance (R⁻²). Based on the considerations above, the complexes ZrF₂(OH)₂·2H₂O and HfF₂(OH)₂·2H₂O were fitted to the Zr and Hf data, respectively, by defining a single ligand shell containing 2F and 4O at an identical distance. The radial distance R and displacement σ were refined, the coordination number was kept fixed.

<caption>**Table 3**</caption>

EXAFS fitting results as reported by xafsX. No error can be calculated for the radial distance R because it is totally correlated with the mean displacement σ.

<table><thead><tr><th></th><th>T [°C]</th><th>R [Å]</th><th>σ² [Å²]</th><th>±</th><th>E₀ [eV]</th><th>±</th></tr></thead><tbody><tr><td rowspan="5">Zr</td><td>100</td><td>2.008</td><td>0.004</td><td>0.001</td><td>-2.952</td><td>0.0001</td></tr><tr><td>200</td><td>1.990</td><td>0.005</td><td>0.001</td><td>-3.331</td><td>0.0001</td></tr><tr><td>250</td><td>1.990</td><td>0.005</td><td>0.002</td><td>-4.810</td><td>0.0002</td></tr><tr><td>300</td><td>1.994</td><td>0.006</td><td>0.005</td><td>-4.495</td><td>0.0003</td></tr><tr><td>350</td><td>1.974</td><td>0.005</td><td>0.008</td><td>-2.852</td><td>0.0006</td></tr><tr><td rowspan="5">Hf</td><td>100</td><td>2.013</td><td>0.003</td><td>0.002</td><td>0.190</td><td>0.0001</td></tr><tr><td>200</td><td>1.975</td><td>0.004</td><td>0.002</td><td>3.850</td><td>0.0002</td></tr><tr><td>250</td><td>1.970</td><td>0.005</td><td>0.002</td><td>0.791</td><td>0.0002</td></tr><tr><td>300</td><td>1.969</td><td>0.005</td><td>0.004</td><td>0.573</td><td>0.0003</td></tr><tr><td>350</td><td>1.965</td><td>0.006</td><td>0.007</td><td>0.960</td><td>0.0005</td></tr><tr><td></td><td>400</td><td>1.969</td><td>0.01</td><td>0.01</td><td>1.473</td><td>0.001</td></tr></tbody></table>

The results are summarized in Table 3. The nature of the fitting procedure precludes the independent determination of an experimental error or uncertainty on R in addition to σ for individual data points. It is, however, possible to estimate if changes in R with temperature are significant by taking into account the R values for all temperatures for either the Zr or Hf experiments. Temperature dependent properties tend to follow simple trends when plotted against reciprocal absolute temperature. The simplest such dependency being linear, thus providing the most conservative error estimate. We can fit a linear function in reciprocal absolute temperature space to the data (Fig. 6a) and use the largest observed residuals between the data points and this fit as an estimate of the relative experimental uncertainty between the individual data points. These estimated uncertainties of 0.01 Å are shown as error bars in Fig. 6a. It is worth noting that this is an estimate for a purely random relative error between the data points at different temperatures. There is no systematic temperature dependent experimental error for XANES and EXAFS data because no temperature dependent parameters are used in the interpretation of the data and all of the spectroscopic equipment is far away from the autoclave and at a stable temperature. We can, therefore, be confident that any differences beyond the error bars in Fig. 6a are significant and the observed apparent contraction of the complexes with increasing temperature is real. Because no variables that could introduce a systematic error between Zr and Hf are used in the EXAFS analysis and considering that all individual R values of Zr complexes at ≥ 200 °C are larger than the corresponding values of Hf complexes, we would argue that there is also a significant difference
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** a) Mean distances of the 6F or O atoms closest to the central cation, weighted with R⁻² for two sets of Lennard-Jones parameters: σ = 1.8 Å, ε = 20.92 kJ/mol (blue); σ = 3.6 Å, ε = 4.184·10⁻³ kJ/mol (red). Experimental EXAFS results for Zr and Hf (symbols) shown for comparison. b) Same for the 4 closest F and O atoms. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

between the complexes of Zr and Hf at these temperatures.

The bars in Fig. 6b depict the σ value of the total radial displacement of the ligands. Despite the name σ, this is not an error or uncertainty but a distribution of ligand distances in real space! Each ligand has its own thermal distribution around its rest position due to vibration, which broadens with increasing temperature. The rest position is also different for each ligand (complex asymmetry). The bars in Fig. 6b therefore denote the range in which any particular ligand is expected at any instant in time with a 1σ probability. This range is large compared to the differences in R (Fig. 6a) but describes actual physical properties of the complexes and has no bearing on the significance of differences in R. To understand the meaning of the EXAFS results, we employ classical MD in the following.

### 3.4. Classical MD

In an octahedral complex with three types of ligands, each of which occurs twice, as in the [MF₂(OH)₂·2H₂O] complexes simulated here, five symmetrically distinct configurations of ligands are possible: all-trans (all identical ligands lie opposite each other), all-cis (all ligands lie at 90° to their identical counterpart), and three permutations with one of the ligand types in trans and the other two in cis configuration. Calculations were performed for all five possible configurations. However, the all-trans and the three mixed configurations yield very similar results, whereas the all-cis configuration is unstable. Therefore, only the results of the all-trans configuration will be presented hereinafter. A comparison of the main results for all different configurations can be found in the

### Supplementary Material (Fig. S4 and S5).

The calculated mean distances of the six (and four) F or O atoms closest to the central cation are shown in Fig. 7a (and b, respectively), weighted with R⁻² and calculated for two sets of Lennard-Jones parameters for the central cation: σ = 1.8 Å, ε = 20.92 kJ/mol (blue in Fig. 7); and σ = 3.6 Å, ε = 4.184·10⁻³ kJ/mol (red). The R⁻² weighting is applied for better comparability with the experimental EXAFS results (Fig. 6 and also shown in Fig. 7) because the amplitude of the EXAFS effect also includes the same weighting factor. The full density functions for all three ligand types and both Lennard-Jones parameter sets, as calculated with classical MD, show the expected thermal broadening of the distributions with increasing temperature (Fig. 8). The weighted average cation-ligand distances for the closest six F and O ligand atoms show thermal expansion (Fig. 7a), whereas the closest four alone show thermal contraction (Fig. 7b). This is due to the two H₂O ligands moving out from < 2.2 Å on average at ≤ 200 °C to > 2.3 Å on average at ≥ 300 °C, and finally an increasing loss of H₂O ligands from the octahedron to an outer and much less well-defined ligand position > 3 Å between 300 and 400 °C (Fig. 8 bottom row). This effect is observable for both simulated sets of Lennard-Jones parameters but much more pronounced for the set shown in red in Fig. 7 and in the right column in Fig. 8, which also displays a more pronounced asymmetry of the octahedral complex even at low temperatures. The other ligands, F⁻ and OH⁻, show an opposing trend compared to H₂O, with a slight decrease of the distance to the central cation (Fig. 8) with increasing temperature above 200 °C, likely due to the decreased repulsion by the H₂O ligands.

The H₂O ligand distance distribution becomes bimodal with peaks at ~ 2.3 and ~ 3.7 Å above 300 °C (Fig. 8 bottom right). The relocation of the H₂O ligands from ~ 2.3 to ~ 3.7 Å is in effect a loss of these ligands from the first coordination shell, which is thus reduced from 6- to 5- and ultimately 4-fold coordination. The remaining inner ligands reorganize and assume first an irregular and finally a tetrahedral configuration.

Computation of the classical interaction potential (Lennard-Jones potential plus Coulomb potential) as a function of cation-fluoride distance also reproduces the negative correlation between the cation distance of H₂O and fluoride (Fig. 9). If the cation-ligand distance for H₂O is increased from 2.0 Å to 4.0 Å, the minimum of the cation-ligand potential for the fluoride ions shifts to smaller values. The absolute values for Rₓ are offset compared to the full MD calculations because the interaction potential calculations were carried out *in vacuo*, thus not accounting for solvent interactions.

## 4. Discussion

### 4.1. Complex stoichiometry and first-shell coordination

Previous experimental studies on the solubility of ZrO₂ in HF-bearing solutions have found that neutral complexes of the stoichiometry [ZrF₂(OH)₂]¹⁰ should be stable under conditions investigated here (Ryzhenko et al., 2008; Migdisov et al., 2011). Because mixed-ligand coordination with F and O as nearest neighbor atoms of a central cation cannot be unequivocally determined with X-ray spectroscopy due to the similar behavior of F and O as X-ray scatterers (Farges, 1991; 1996, Louvel et al., 2013), we used the published [ZrF₂(OH)₂]¹⁰ complex to provide the fluoride and hydroxide stoichiometry for our XAS data analysis. However, because Ryzhenko et al. (2008) and Migdisov et al. (2011) utilized the dissolution reactions of a solid to derive the complex stoichiometry, they could not observe electrically neutral molecules of the solvent, such as H₂O, in the complexes. Comparison of our experimental XANES data (Fig. 2) with spectra calculated from AIMD calculations (Fig. 3) as well as published experimental and theoretical data (Louvel et al., 2013; Wilke et al., 2012; 2013) suggest that the complex coordination is predominantly 6-fold. Wilke et al., (2012, 2013), using X-ray absorption/fluorescence spectroscopy in hydrothermal diamond anvil cells (HDAC), found that the solubility of zircon (ZrSiO₄) and hafnon (HfSiO₄) in solutions with up to 30 % sodium silicate at
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** a) Mean distances of the 6F or O atoms closest to the central cation, weighted with R⁻² for two sets of Lennard-Jones parameters: σ = 1.8 Å, ε = 20.92 kJ/mol (blue); σ = 3.6 Å, ε = 4.184·10⁻³ kJ/mol (red). Experimental EXAFS results for Zr and Hf (symbols) shown for comparison. b) Same for the 4 closest F and O atoms. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

between the complexes of Zr and Hf at these temperatures.

The bars in Fig. 6b depict the σ value of the total radial displacement of the ligands. Despite the name σ, this is not an error or uncertainty but a distribution of ligand distances in real space! Each ligand has its own thermal distribution around its rest position due to vibration, which broadens with increasing temperature. The rest position is also different for each ligand (complex asymmetry). The bars in Fig. 6b therefore denote the range in which any particular ligand is expected at any instant in time with a 1σ probability. This range is large compared to the differences in R (Fig. 6a) but describes actual physical properties of the complexes and has no bearing on the significance of differences in R. To understand the meaning of the EXAFS results, we employ classical MD in the following.

### 3.4. Classical MD

In an octahedral complex with three types of ligands, each of which occurs twice, as in the [MF₂(OH)₂·2H₂O] complexes simulated here, five symmetrically distinct configurations of ligands are possible: all-trans (all identical ligands lie opposite each other), all-cis (all ligands lie at 90° to their identical counterpart), and three permutations with one of the ligand types in trans and the other two in cis configuration. Calculations were performed for all five possible configurations. However, the all-trans and the three mixed configurations yield very similar results, whereas the all-cis configuration is unstable. Therefore, only the results of the all-trans configuration will be presented hereinafter. A comparison of the main results for all different configurations can be found in the

### Supplementary Material (Fig. S4 and S5).

The calculated mean distances of the six (and four) F or O atoms closest to the central cation are shown in Fig. 7a (and b, respectively), weighted with R⁻² and calculated for two sets of Lennard-Jones parameters for the central cation: σ = 1.8 Å, ε = 20.92 kJ/mol (blue in Fig. 7); and σ = 3.6 Å, ε = 4.184·10⁻³ kJ/mol (red). The R⁻² weighting is applied for better comparability with the experimental EXAFS results (Fig. 6 and also shown in Fig. 7) because the amplitude of the EXAFS effect also includes the same weighting factor. The full density functions for all three ligand types and both Lennard-Jones parameter sets, as calculated with classical MD, show the expected thermal broadening of the distributions with increasing temperature (Fig. 8). The weighted average cation-ligand distances for the closest six F and O ligand atoms show thermal expansion (Fig. 7a), whereas the closest four alone show thermal contraction (Fig. 7b). This is due to the two H₂O ligands moving out from < 2.2 Å on average at ≤ 200 °C to > 2.3 Å on average at ≥ 300 °C, and finally an increasing loss of H₂O ligands from the octahedron to an outer and much less well-defined ligand position > 3 Å between 300 and 400 °C (Fig. 8 bottom row). This effect is observable for both simulated sets of Lennard-Jones parameters but much more pronounced for the set shown in red in Fig. 7 and in the right column in Fig. 8, which also displays a more pronounced asymmetry of the octahedral complex even at low temperatures. The other ligands, F⁻ and OH⁻, show an opposing trend compared to H₂O, with a slight decrease of the distance to the central cation (Fig. 8) with increasing temperature above 200 °C, likely due to the decreased repulsion by the H₂O ligands.

The H₂O ligand distance distribution becomes bimodal with peaks at ~ 2.3 and ~ 3.7 Å above 300 °C (Fig. 8 bottom right). The relocation of the H₂O ligands from ~ 2.3 to ~ 3.7 Å is in effect a loss of these ligands from the first coordination shell, which is thus reduced from 6- to 5- and ultimately 4-fold coordination. The remaining inner ligands reorganize and assume first an irregular and finally a tetrahedral configuration.

Computation of the classical interaction potential (Lennard-Jones potential plus Coulomb potential) as a function of cation-fluoride distance also reproduces the negative correlation between the cation distance of H₂O and fluoride (Fig. 9). If the cation-ligand distance for H₂O is increased from 2.0 Å to 4.0 Å, the minimum of the cation-ligand potential for the fluoride ions shifts to smaller values. The absolute values for Rₓ are offset compared to the full MD calculations because the interaction potential calculations were carried out *in vacuo*, thus not accounting for solvent interactions.

## 4. Discussion

### 4.1. Complex stoichiometry and first-shell coordination

Previous experimental studies on the solubility of ZrO₂ in HF-bearing solutions have found that neutral complexes of the stoichiometry [ZrF₂(OH)₂]¹⁰ should be stable under conditions investigated here (Ryzhenko et al., 2008; Migdisov et al., 2011). Because mixed-ligand coordination with F and O as nearest neighbor atoms of a central cation cannot be unequivocally determined with X-ray spectroscopy due to the similar behavior of F and O as X-ray scatterers (Farges, 1991; 1996, Louvel et al., 2013), we used the published [ZrF₂(OH)₂]¹⁰ complex to provide the fluoride and hydroxide stoichiometry for our XAS data analysis. However, because Ryzhenko et al. (2008) and Migdisov et al. (2011) utilized the dissolution reactions of a solid to derive the complex stoichiometry, they could not observe electrically neutral molecules of the solvent, such as H₂O, in the complexes. Comparison of our experimental XANES data (Fig. 2) with spectra calculated from AIMD calculations (Fig. 3) as well as published experimental and theoretical data (Louvel et al., 2013; Wilke et al., 2012; 2013) suggest that the complex coordination is predominantly 6-fold. Wilke et al., (2012, 2013), using X-ray absorption/fluorescence spectroscopy in hydrothermal diamond anvil cells (HDAC), found that the solubility of zircon (ZrSiO₄) and hafnon (HfSiO₄) in solutions with up to 30 % sodium silicate at
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Calculated radial distributions functions of F⁻, OH⁻, and H₂O ligands at various temperatures for two sets of Lennard-Jones parameters: σ = 1.8 Å, ε = 20.92 kJ/mol (left column); σ = 3.6 Å, ε = 4.184·10⁻³ kJ/mol (right column). Each line shows the envelope of a histogram (bin width = 0.01 Å) of the distances of both ligands of each type from the central cation in 2000 snapshots of the classical MD simulation. Note that F⁻ and OH⁻ ligands are shown on common x and y scales but the H₂O ligands have their own scales due to their much wider distribution. The sum of the distributions of all three ligands at each temperature corresponds to the total displacement determined with EXAFS (Fig. 6b).</figcaption></figure>

temperatures up to 700 °C and pressures up to 1.3 GPa, is much higher
compared to NaOH solutions, and interpreted this to be due to the for-
mation of 6-fold coordinated zirconosilicate-like complexes. Vlasovite-
like alkali-zirconosilicate clusters with 6-fold coordination were also
found by Louvel et al. (2013) in solutions containing 35 % Na₂Si₂O₅
(NS2). The characteristic shape of our experimental XANES spectra is
well in the variation range of the various investigated and published 6-
fold coordinated complexes and allows us to rule out higher coordina-
tion. Based on these published data and our own results, we therefore
conclude that the dominant complex stoichiometry including neutral
solvent molecules under the physicochemical conditions investigated
here, is [Zr(F,OH)₄·2H₂O]⁰. Since crystalline solids of Zr that could be
used as reference material all have 6-fold or higher Zr coordination, we
computed the Zr XANES spectra of a transient 5-fold coordinated Zr
species existing for about 5 ps in an AIMD simulation at 350 °C. As can
be seen in Fig. S6, the white line is only slightly affected by the coor-
dination change but the intensity of the pre-edge is significantly

enhanced. Such a change is not observed in the experimental spectra
suggesting that a 5-fold coordinated environment is not abundant at the
experimental conditions.

The potential formation of polynuclear complexes needs to be
considered because (1) it may have significant implications for complex
diffusion and mineral precipitation from fluids and (2) Zr and Hf are
known to have a propensity for complex polymerization. Hagfeldt et al.
(2004) demonstrated spectroscopically that hydrated Zr⁴⁺ and Hf⁴⁺ ions
polymerize to edge-sharing tetramers by hydrolysis at any pH greater
than zero in pure water. A tendency of Zr to form polymers is also known
from more complex chemical systems. Zirconium oxychloride forms
tetramers whereas zirconosulfate trimers occur in sulfate-rich solutions
(Kanazhevskii et al., 2006). Solubility based studies (Ryzhenko et al.,
2008; Migdisov et al., 2011) cannot distinguish between mono- and
polynuclear complexes of the same dissolution reaction stoichiometry.
However, heavy elements such as Zr or Hf at a distance that would be
expected for a polynuclear complex, should be visible in EXAFS spectra,
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Calculated radial distributions functions of F⁻, OH⁻, and H₂O ligands at various temperatures for two sets of Lennard-Jones parameters: σ = 1.8 Å, ε = 20.92 kJ/mol (left column); σ = 3.6 Å, ε = 4.184·10⁻³ kJ/mol (right column). Each line shows the envelope of a histogram (bin width = 0.01 Å) of the distances of both ligands of each type from the central cation in 2000 snapshots of the classical MD simulation. Note that F⁻ and OH⁻ ligands are shown on common x and y scales but the H₂O ligands have their own scales due to their much wider distribution. The sum of the distributions of all three ligands at each temperature corresponds to the total displacement determined with EXAFS (Fig. 6b).</figcaption></figure>

temperatures up to 700 °C and pressures up to 1.3 GPa, is much higher
compared to NaOH solutions, and interpreted this to be due to the for-
mation of 6-fold coordinated zirconosilicate-like complexes. Vlasovite-
like alkali-zirconosilicate clusters with 6-fold coordination were also
found by Louvel et al. (2013) in solutions containing 35 % Na₂Si₂O₅
(NS2). The characteristic shape of our experimental XANES spectra is
well in the variation range of the various investigated and published 6-
fold coordinated complexes and allows us to rule out higher coordina-
tion. Based on these published data and our own results, we therefore
conclude that the dominant complex stoichiometry including neutral
solvent molecules under the physicochemical conditions investigated
here, is [Zr(F,OH)₄·2H₂O]⁰. Since crystalline solids of Zr that could be
used as reference material all have 6-fold or higher Zr coordination, we
computed the Zr XANES spectra of a transient 5-fold coordinated Zr
species existing for about 5 ps in an AIMD simulation at 350 °C. As can
be seen in Fig. S6, the white line is only slightly affected by the coor-
dination change but the intensity of the pre-edge is significantly

enhanced. Such a change is not observed in the experimental spectra
suggesting that a 5-fold coordinated environment is not abundant at the
experimental conditions.

The potential formation of polynuclear complexes needs to be
considered because (1) it may have significant implications for complex
diffusion and mineral precipitation from fluids and (2) Zr and Hf are
known to have a propensity for complex polymerization. Hagfeldt et al.
(2004) demonstrated spectroscopically that hydrated Zr⁴⁺ and Hf⁴⁺ ions
polymerize to edge-sharing tetramers by hydrolysis at any pH greater
than zero in pure water. A tendency of Zr to form polymers is also known
from more complex chemical systems. Zirconium oxychloride forms
tetramers whereas zirconosulfate trimers occur in sulfate-rich solutions
(Kanazhevskii et al., 2006). Solubility based studies (Ryzhenko et al.,
2008; Migdisov et al., 2011) cannot distinguish between mono- and
polynuclear complexes of the same dissolution reaction stoichiometry.
However, heavy elements such as Zr or Hf at a distance that would be
expected for a polynuclear complex, should be visible in EXAFS spectra,
<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Lennard-Jones potential $V_{LJ}$ plus Coulomb potential $V_C$ as a function of cation-fluoride distance $R_F$ for several cation $H_2O$ distances ($R_{H2O}$). The vertical lines indicate the minima of the functions for $R_{H2O} = 2.0$ Å and $4.0$ Å, respectively.</figcaption></figure>

as has been observed in solids as well as aqueous complexes (e.g., Hagfeld 2004; Louvel et al., 2013). However, this was not observed in our EXAFS data. The radial distribution functions do not show a clear second shell in any of our experiments (Fig. 5). Attempts to fit a second ligand shell containing one or more Zr at distances that would represent face-, edge-, or corner-sharing octahedra also did not fit the results better than without the additional Zr and were consequently abandoned. We therefore conclude that no polynuclear complexes occurred in our experiments.

Our AIMD simulations indicate that a range of different F/OH ratios in the complex are at least metastable over the course of the simulations (25 ps). It is therefore likely that the average complex formula represents a mixture of various different individual complex configurations. Neither the experimental methods employed in the present study nor those in the literature can determine the individual complexes at a single point in time but represent averages over many individual complexes over a time span (minutes to days) that is many orders of magnitude greater than the time scale of ligand exchange (pico- to milliseconds). It is worth noting that the AIMD simulations that were started with tetrahedral Zr or Hf complexes, equilibrated quickly to octahedral coordination by incorporating nearby $H_2O$ molecules into the first coordination shell, suggesting that purely tetrahedral coordination is not stable under the investigated conditions. Consequently, the classical MD calculations were set up to start from octahedrally coordinated $[MF_2(OH)_2 \cdot 2H_2O]^0$ complexes. Classical MD allows to track the development of systems over much longer time periods than AIMD (3 ns vs. 25 ps) and therefore provides more sound statistical information about the radial distribution of each ligand (Fig. 8). The radial distances of both F⁻ and OH⁻ ligands decrease slightly with increasing temperature, whereas those of the $H_2O$ ligands increase. Up to 200 °C, this causes increasing distortion in the octahedral coordination in the MD

calculations. Above 300 °C, however, the partial loss of $H_2O$ ligands from the first shell (<3 Å) to the second (>3 Å) leads to an apparent equilibrium between 6-, 5-, and 4-fold first shell coordination. The formation and persistence (at least for a few picoseconds) of a 5-fold coordinated complex at 350 °C in the AIMD simulations during ligand exchange (see above and Fig. S6) also suggests that lower coordinated species may become stable at high temperatures. The two different parameter sets used for the classical MD calculations both qualitatively show the increasing asymmetry and tendency to transition to lower than 6-fold coordination with increasing temperature but yield substantially different proportions of $H_2O$ beyond 3 Å at high temperature (bottom two panels in Fig. 8) and thus different proportions of 6-, 5-, and 4-fold coordination. This indicates that the increasing asymmetry of the octahedron likely does lead to a transition to lower coordination at higher temperature, but the coordination proportions should not be used to extrapolate free energies of transition from our classical MD calculations. For this purpose, dedicated AIMD calculations with longer simulation times and free energy thermodynamic integration need to be performed in the future.

Comparison of the R values obtained with AIMD and classical MD with those from the EXAFS experiments (Figs. 6 and 7) strongly suggest that the apparent contraction of complexes observed in the experiments in fact constitutes an increase in asymmetry of the first shell or even a partial transition from octahedral to lower coordination. Because the amplitude of $\chi(k)$, which is defined by the EXAFS function used in the fitting procedure (Newville, 2014), scales with the inverse square of the radial distribution, ligands that are further away from the central cation have a much lower influence on the overall EXAFS function than ligands that are closer. This means that light ligands with a large R and $\sigma$, such as $H_2O$ in this case, become much less visible to the EXAFS analysis. Beyond ca. 3 Å, light ligands become virtually invisible to EXAFS and cannot be properly fitted as a second shell.

Unlike in the MD results, the transition to lower average radial distance is already visible at 200 °C in the EXAFS data (Fig. 6a). It is particularly interesting that this effect seems to be more pronounced for the Hf complexes than for those of Zr between 200 and 300 °C, suggesting that Hf complexes develop a somewhat more distorted coordination environment and potentially even a slightly lower average coordination than Zr complexes under these conditions. It would therefore be useful to have direct evidence on the Hf coordination. The Hf $L_3$-edge does indeed show a change in shape of the white line from 100 to 250 °C, which may correspond to a change in coordination (no XANES data are available for 200 °C, see Results). At 100 °C, the white line consists of two clearly identifiable peaks. At 250 °C and above, this is not observed and the white line seems to consist of one single peak. The theoretical XANES spectra obtained from the AIMD calculations reproduce this feature (Fig. 3). The disappearance of the line splitting with increasing temperature may be the result of increasing asymmetry or simply due to thermal line broadening. Apart from the line splitting in the Hf spectra, the experimental XANES spectra (Fig. 2) show no obvious evidence for a coordination change in the investigated temperature range. The AIMD simulations and theoretical XANES spectra (Table 2, Fig. 3) suggest that the increase in asymmetry only causes minor changes to the XANES spectra of Zr.

## 4.2. Implications for Zr/Hf fractionation

Zirconium and hafnium behave almost identically in the vast majority of geochemical environments and are therefore considered geochemical twin elements. However, extreme enrichment of Hf and formation of high-Hf zircon to near-end-member hafnon ($HfSiO_4$) is known from a few rare-metal pegmatite locations across the globe (e.g., Neves et al., 1974; Yin et al., 2013; Kudryashov et al., 2020). The Hf enrichment in these examples is tied to the late magmatic to hydrothermal phase transition in the pegmatite evolution when fluxing, i.e., depolymerizing, elements such as F, Li, B, Na, K, and water fractionate
<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Lennard-Jones potential $V_{LJ}$ plus Coulomb potential $V_C$ as a function of cation-fluoride distance $R_F$ for several cation $H_2O$ distances ($R_{H2O}$). The vertical lines indicate the minima of the functions for $R_{H2O} = 2.0$ Å and $4.0$ Å, respectively.</figcaption></figure>

as has been observed in solids as well as aqueous complexes (e.g., Hagfeld 2004; Louvel et al., 2013). However, this was not observed in our EXAFS data. The radial distribution functions do not show a clear second shell in any of our experiments (Fig. 5). Attempts to fit a second ligand shell containing one or more Zr at distances that would represent face-, edge-, or corner-sharing octahedra also did not fit the results better than without the additional Zr and were consequently abandoned. We therefore conclude that no polynuclear complexes occurred in our experiments.

Our AIMD simulations indicate that a range of different F/OH ratios in the complex are at least metastable over the course of the simulations (25 ps). It is therefore likely that the average complex formula represents a mixture of various different individual complex configurations. Neither the experimental methods employed in the present study nor those in the literature can determine the individual complexes at a single point in time but represent averages over many individual complexes over a time span (minutes to days) that is many orders of magnitude greater than the time scale of ligand exchange (pico- to milliseconds). It is worth noting that the AIMD simulations that were started with tetrahedral Zr or Hf complexes, equilibrated quickly to octahedral coordination by incorporating nearby $H_2O$ molecules into the first coordination shell, suggesting that purely tetrahedral coordination is not stable under the investigated conditions. Consequently, the classical MD calculations were set up to start from octahedrally coordinated $[MF_2(OH)_2 \cdot 2H_2O]^0$ complexes. Classical MD allows to track the development of systems over much longer time periods than AIMD (3 ns vs. 25 ps) and therefore provides more sound statistical information about the radial distribution of each ligand (Fig. 8). The radial distances of both F⁻ and OH⁻ ligands decrease slightly with increasing temperature, whereas those of the $H_2O$ ligands increase. Up to 200 °C, this causes increasing distortion in the octahedral coordination in the MD

calculations. Above 300 °C, however, the partial loss of $H_2O$ ligands from the first shell (<3 Å) to the second (>3 Å) leads to an apparent equilibrium between 6-, 5-, and 4-fold first shell coordination. The formation and persistence (at least for a few picoseconds) of a 5-fold coordinated complex at 350 °C in the AIMD simulations during ligand exchange (see above and Fig. S6) also suggests that lower coordinated species may become stable at high temperatures. The two different parameter sets used for the classical MD calculations both qualitatively show the increasing asymmetry and tendency to transition to lower than 6-fold coordination with increasing temperature but yield substantially different proportions of $H_2O$ beyond 3 Å at high temperature (bottom two panels in Fig. 8) and thus different proportions of 6-, 5-, and 4-fold coordination. This indicates that the increasing asymmetry of the octahedron likely does lead to a transition to lower coordination at higher temperature, but the coordination proportions should not be used to extrapolate free energies of transition from our classical MD calculations. For this purpose, dedicated AIMD calculations with longer simulation times and free energy thermodynamic integration need to be performed in the future.

Comparison of the R values obtained with AIMD and classical MD with those from the EXAFS experiments (Figs. 6 and 7) strongly suggest that the apparent contraction of complexes observed in the experiments in fact constitutes an increase in asymmetry of the first shell or even a partial transition from octahedral to lower coordination. Because the amplitude of $\chi(k)$, which is defined by the EXAFS function used in the fitting procedure (Newville, 2014), scales with the inverse square of the radial distribution, ligands that are further away from the central cation have a much lower influence on the overall EXAFS function than ligands that are closer. This means that light ligands with a large R and $\sigma$, such as $H_2O$ in this case, become much less visible to the EXAFS analysis. Beyond ca. 3 Å, light ligands become virtually invisible to EXAFS and cannot be properly fitted as a second shell.

Unlike in the MD results, the transition to lower average radial distance is already visible at 200 °C in the EXAFS data (Fig. 6a). It is particularly interesting that this effect seems to be more pronounced for the Hf complexes than for those of Zr between 200 and 300 °C, suggesting that Hf complexes develop a somewhat more distorted coordination environment and potentially even a slightly lower average coordination than Zr complexes under these conditions. It would therefore be useful to have direct evidence on the Hf coordination. The Hf $L_3$-edge does indeed show a change in shape of the white line from 100 to 250 °C, which may correspond to a change in coordination (no XANES data are available for 200 °C, see Results). At 100 °C, the white line consists of two clearly identifiable peaks. At 250 °C and above, this is not observed and the white line seems to consist of one single peak. The theoretical XANES spectra obtained from the AIMD calculations reproduce this feature (Fig. 3). The disappearance of the line splitting with increasing temperature may be the result of increasing asymmetry or simply due to thermal line broadening. Apart from the line splitting in the Hf spectra, the experimental XANES spectra (Fig. 2) show no obvious evidence for a coordination change in the investigated temperature range. The AIMD simulations and theoretical XANES spectra (Table 2, Fig. 3) suggest that the increase in asymmetry only causes minor changes to the XANES spectra of Zr.

## 4.2. Implications for Zr/Hf fractionation

Zirconium and hafnium behave almost identically in the vast majority of geochemical environments and are therefore considered geochemical twin elements. However, extreme enrichment of Hf and formation of high-Hf zircon to near-end-member hafnon ($HfSiO_4$) is known from a few rare-metal pegmatite locations across the globe (e.g., Neves et al., 1974; Yin et al., 2013; Kudryashov et al., 2020). The Hf enrichment in these examples is tied to the late magmatic to hydrothermal phase transition in the pegmatite evolution when fluxing, i.e., depolymerizing, elements such as F, Li, B, Na, K, and water fractionate
from the melt into an exsolved aqueous phase, thus drastically altering the properties of the remaining melt, which becomes highly polymerized and also increasingly aluminous. This indicates that the relative solubility of Zr and Hf in the exsolving aqueous fluids, or in the melt (which has undergone drastic changes to its chemical composition), or both, may be responsible for Zr/Hf fractionation in these systems. A similar but less extreme decline in Zr/Hf ratio has been observed in late magmatic and especially hydrothermal zircons in granites (Wang et al., 2010). Fractionation of Zr and Hf has also been experimentally demonstrated in peraluminous melts (Linnen and Keppler, 2002) as well as aqueous fluids (Shikina et al., 2015) and is known from a range of different natural aqueous systems, where aqueous and surface complexation are likely involved (e.g., Bau, 1996; Schmidt et al., 2014; Inguaggiato et al., 2015; Censi et al., 2017). This complicates the interpretation of Zr/Hf fractionation especially in settings where more than two phases (fluid, melt, solids) are present.

Despite the abundance of studies describing Zr/Hf fractionation (e.g., Bau, 1996; Linnen and Keppler, 2002; Wang et al., 2010; Schmidt et al., 2014; Inguaggiato et al., 2015; Shikina et al., 2015; Censi et al., 2017), the fundamental cause of the effect is currently unclear. That means it is unknown, which physical properties of Zr and Hf affect element partitioning on a molecular level, and via which mechanism they affect it.

While nominal ionic charge and ionic radius of Zr⁴⁺ and Hf⁴⁺ are very similar (~1.4 % difference; Shannon, 1976), the mass of Hf is almost twice that of Zr. Hafnium also has additional electron shells (68 electrons in total for Hf⁴⁺ compared to 36 for Zr⁴⁺), which do not alter the valence electron configuration directly but may nevertheless play a role in modifying bond strength with anions.

The difference in atomic mass (Wieser and Coplen, 2011) between Zr (91.224) and Hf (178.49) can be expected to play at least some role in Zr/Hf fractionation, if we assume that the same mechanisms that are responsible for isotope fractionation, are involved. Mass dependent fractionation in aqueous systems is known to occur via adsorption of stable isotopes of, among others, B (Lemarchand et al., 2007), Zn (Juillot et al., 2008), and Mo (Wasylenki et al., 2011) onto minerals, in cases where the adsorption is accompanied by a coordination change between dissolved complex and surface complex. If we were, for the sake of argument to assume that Zr and Hf do not differ from each other in complexation behavior, we could treat them like “pseudo isotopes” of the same element with a mass difference of almost a factor of 2. The fractionation effect produced by such a mechanism can be expected to be much more drastic for Zr/Hf than the fractionation observed for ⁹⁷Mo/⁹⁵Mo (Wasylenki et al., 2011), at only 2 % relative mass difference. In the most common Zr- and Hf-bearing minerals, Zr and Hf are 6-fold (zirconosilicates), 7-fold (baddeleyite), or 8-fold (zircon, hafnon) coordinated. If Zr and Hf are 6-fold or less coordinated on average in complexes occurring in an aqueous fluid, a significant increase of the coordination number must take place during either adsorption of a dissolved complex onto a mineral surface or during incorporation of the surface complex into the mineral structure. The change from a highly asymmetric complex (as we found in solution) to a more symmetric configuration with the same coordination number can also increase partitioning because the shorter bonds in the asymmetric complex contribute more to the overall energy budget (e.g., Wunder et al., 2011). The increase in effective bond length during precipitation as a result of the coordination or symmetry change would be expected to cause preferential precipitation of the lighter component (Zr) and retain the heavier component (Hf) in solution in analogy to Mo isotope fractionation (Wasylenki et al., 2011). This is due to the fact that the heavier cation forms bonds with a lower vibrational energy. As a result, heavy isotopes or elements are favored in low coordination environments, which have shorter bond lengths and stronger bonds.

Fig. 10 schematically demonstrates the development of Zr/Hf in minerals precipitating from a fluid in a Rayleigh fractionation process (as opposed to equilibrium fractionation) plotted against the proportion

of Zr that has been precipitated for several different equilibrium fractionation factor ratios. It shows that even low differences in fractionation factors for Zr and Hf can produce very low Zr/Hf ratios, if fractionation is allowed to progress sufficiently. This means that minerals precipitating from highly evolved F-rich aqueous solutions (e.g., hydrothermal vein fluorite or hydrothermal zircon), which have undergone prolonged Rayleigh fractionation should have subchondritic Zr/Hf ratios, which is precisely what studies on natural samples have observed (e.g., Bau, 1996; Wang et al., 2010). If the same fractionation mechanism also affects zircon precipitation from melts, this would also explain the extreme Hf enrichment in pegmatitic zircons or, indeed, the formation of near end-member hafnon (e.g., Neves et al., 1974; Yin et al., 2013; Kudryashov et al., 2020). This mechanism has also been proposed for Zr stable isotope fractionation (e.g., Inglis et al., 2019). However, several newer studies found that equilibrium Rayleigh fractionation is not sufficient to explain the large positive as well as negative δ⁹⁴/⁹⁰Zr values found in some magmatic systems (e.g., Méheut et al., 2021; Tissot and Ibañez-Mejia, 2021; Bindeman and Melnik, 2022; Tompkins et al., 2023). These studies report strong mass dependent kinetic isotope fractionation effects. The fact that δ⁹⁴/⁹⁰Zr values can be positively (e.g., Ibañez-Mejia and Tissot, 2019) or negatively (e.g., Zhu et al., 2023) correlated with Zr/Hf ratios, may suggest that Zr/Hf fractionation is either controlled by strong mass independent effects (unlike Zr stable isotopes) or that mass dependent equilibrium fractionation plays a greater role than for Zr stable isotopes due to the larger mass difference between Zr and Hf. It is worth noting that mass dependent diffusion kinetic effects can be expected to be less relevant in aqueous compared to magmatic systems due to the lower viscosity of aqueous fluids, whereas (mass dependent or mass independent) reaction kinetic effects may be of similar magnitude in both kinds of systems.

To the best of our knowledge, Zr stable isotope studies are at the moment limited to magmatic systems and do not cover fluid-mineral interaction. This is unfortunate because the same reasoning applied to magmatic systems may also explain the Zr/Hf fractionation trend observed in deep oceanic hydrogen genetic ferromanganese crusts (Schmidt et al., 2014). However, in this case the precipitating phase has a lower Zr/Hf ratio than the aqueous phase, which should be the case if coordination decreases as Zr or Hf from solution get adsorbed onto the mineral surface. Currently, neither the complex chemistry of Zr and Hf in deep sea water nor their surface complexation on ferromanganese crusts are comprehensively understood. In water without ligands that form strong complexes with these cations, 8-fold coordination of Zr has been predicted with ab initio modeling (Messner et al., 2011). If this is valid for Zr and Hf in deep sea water, and if we postulate 6-fold surface coordination of Zr and Hf on ferromanganese crusts, then the observed trend of lower Zr/Hf ratios in the precipitate compared to the solution could be explained by mass dependent “pseudo isotope” fractionation. However, experimental investigation of the involved coordination environments is needed to test the plausibility of this hypothesis.

Our data also allow a second interpretation that is based on mass independent differences in bond strength rather than mass. According to our interpretation of the EXAFS data, the Hf complexes have slightly but significantly stronger and consequently shorter bonds of Hf to the fluoride and hydroxide ligands, which would be expected to push the H₂O ligands out more easily due to steric hindrance. This may be caused by the higher reduced mass of the ligands' bonds with Hf compared to Zr, but may also be the result of electronic effects. Our observations of different bond lengths to individual ligands between Zr and Hf as central atoms are supported by recent *ab initio* density functional theory (DFT) studies that have calculated higher cohesive energies for oxides and fluorides of Hf compared to the respective Zr counterparts (Anders et al., 2022). The formation of HfO₂ from the elements was calculated to be 0.261 eV more favorable than that of ZrO₂. The formation of HfF₄ was 0.274 eV more favorable than that of ZrF₄, meaning that Hf has a higher affinity for O and F compared to Zr and that this effect is even more pronounced for F than for O. Applied to our complex data, this means
from the melt into an exsolved aqueous phase, thus drastically altering the properties of the remaining melt, which becomes highly polymerized and also increasingly aluminous. This indicates that the relative solubility of Zr and Hf in the exsolving aqueous fluids, or in the melt (which has undergone drastic changes to its chemical composition), or both, may be responsible for Zr/Hf fractionation in these systems. A similar but less extreme decline in Zr/Hf ratio has been observed in late magmatic and especially hydrothermal zircons in granites (Wang et al., 2010). Fractionation of Zr and Hf has also been experimentally demonstrated in peraluminous melts (Linnen and Keppler, 2002) as well as aqueous fluids (Shikina et al., 2015) and is known from a range of different natural aqueous systems, where aqueous and surface complexation are likely involved (e.g., Bau, 1996; Schmidt et al., 2014; Inguaggiato et al., 2015; Censi et al., 2017). This complicates the interpretation of Zr/Hf fractionation especially in settings where more than two phases (fluid, melt, solids) are present.

Despite the abundance of studies describing Zr/Hf fractionation (e.g., Bau, 1996; Linnen and Keppler, 2002; Wang et al., 2010; Schmidt et al., 2014; Inguaggiato et al., 2015; Shikina et al., 2015; Censi et al., 2017), the fundamental cause of the effect is currently unclear. That means it is unknown, which physical properties of Zr and Hf affect element partitioning on a molecular level, and via which mechanism they affect it.

While nominal ionic charge and ionic radius of Zr⁴⁺ and Hf⁴⁺ are very similar (~1.4 % difference; Shannon, 1976), the mass of Hf is almost twice that of Zr. Hafnium also has additional electron shells (68 electrons in total for Hf⁴⁺ compared to 36 for Zr⁴⁺), which do not alter the valence electron configuration directly but may nevertheless play a role in modifying bond strength with anions.

The difference in atomic mass (Wieser and Coplen, 2011) between Zr (91.224) and Hf (178.49) can be expected to play at least some role in Zr/Hf fractionation, if we assume that the same mechanisms that are responsible for isotope fractionation, are involved. Mass dependent fractionation in aqueous systems is known to occur via adsorption of stable isotopes of, among others, B (Lemarchand et al., 2007), Zn (Juillot et al., 2008), and Mo (Wasylenki et al., 2011) onto minerals, in cases where the adsorption is accompanied by a coordination change between dissolved complex and surface complex. If we were, for the sake of argument to assume that Zr and Hf do not differ from each other in complexation behavior, we could treat them like “pseudo isotopes” of the same element with a mass difference of almost a factor of 2. The fractionation effect produced by such a mechanism can be expected to be much more drastic for Zr/Hf than the fractionation observed for ⁹⁷Mo/⁹⁵Mo (Wasylenki et al., 2011), at only 2 % relative mass difference. In the most common Zr- and Hf-bearing minerals, Zr and Hf are 6-fold (zirconosilicates), 7-fold (baddeleyite), or 8-fold (zircon, hafnon) coordinated. If Zr and Hf are 6-fold or less coordinated on average in complexes occurring in an aqueous fluid, a significant increase of the coordination number must take place during either adsorption of a dissolved complex onto a mineral surface or during incorporation of the surface complex into the mineral structure. The change from a highly asymmetric complex (as we found in solution) to a more symmetric configuration with the same coordination number can also increase partitioning because the shorter bonds in the asymmetric complex contribute more to the overall energy budget (e.g., Wunder et al., 2011). The increase in effective bond length during precipitation as a result of the coordination or symmetry change would be expected to cause preferential precipitation of the lighter component (Zr) and retain the heavier component (Hf) in solution in analogy to Mo isotope fractionation (Wasylenki et al., 2011). This is due to the fact that the heavier cation forms bonds with a lower vibrational energy. As a result, heavy isotopes or elements are favored in low coordination environments, which have shorter bond lengths and stronger bonds.

Fig. 10 schematically demonstrates the development of Zr/Hf in minerals precipitating from a fluid in a Rayleigh fractionation process (as opposed to equilibrium fractionation) plotted against the proportion

of Zr that has been precipitated for several different equilibrium fractionation factor ratios. It shows that even low differences in fractionation factors for Zr and Hf can produce very low Zr/Hf ratios, if fractionation is allowed to progress sufficiently. This means that minerals precipitating from highly evolved F-rich aqueous solutions (e.g., hydrothermal vein fluorite or hydrothermal zircon), which have undergone prolonged Rayleigh fractionation should have subchondritic Zr/Hf ratios, which is precisely what studies on natural samples have observed (e.g., Bau, 1996; Wang et al., 2010). If the same fractionation mechanism also affects zircon precipitation from melts, this would also explain the extreme Hf enrichment in pegmatitic zircons or, indeed, the formation of near end-member hafnon (e.g., Neves et al., 1974; Yin et al., 2013; Kudryashov et al., 2020). This mechanism has also been proposed for Zr stable isotope fractionation (e.g., Inglis et al., 2019). However, several newer studies found that equilibrium Rayleigh fractionation is not sufficient to explain the large positive as well as negative δ⁹⁴/⁹⁰Zr values found in some magmatic systems (e.g., Méheut et al., 2021; Tissot and Ibañez-Mejia, 2021; Bindeman and Melnik, 2022; Tompkins et al., 2023). These studies report strong mass dependent kinetic isotope fractionation effects. The fact that δ⁹⁴/⁹⁰Zr values can be positively (e.g., Ibañez-Mejia and Tissot, 2019) or negatively (e.g., Zhu et al., 2023) correlated with Zr/Hf ratios, may suggest that Zr/Hf fractionation is either controlled by strong mass independent effects (unlike Zr stable isotopes) or that mass dependent equilibrium fractionation plays a greater role than for Zr stable isotopes due to the larger mass difference between Zr and Hf. It is worth noting that mass dependent diffusion kinetic effects can be expected to be less relevant in aqueous compared to magmatic systems due to the lower viscosity of aqueous fluids, whereas (mass dependent or mass independent) reaction kinetic effects may be of similar magnitude in both kinds of systems.

To the best of our knowledge, Zr stable isotope studies are at the moment limited to magmatic systems and do not cover fluid-mineral interaction. This is unfortunate because the same reasoning applied to magmatic systems may also explain the Zr/Hf fractionation trend observed in deep oceanic hydrogen genetic ferromanganese crusts (Schmidt et al., 2014). However, in this case the precipitating phase has a lower Zr/Hf ratio than the aqueous phase, which should be the case if coordination decreases as Zr or Hf from solution get adsorbed onto the mineral surface. Currently, neither the complex chemistry of Zr and Hf in deep sea water nor their surface complexation on ferromanganese crusts are comprehensively understood. In water without ligands that form strong complexes with these cations, 8-fold coordination of Zr has been predicted with ab initio modeling (Messner et al., 2011). If this is valid for Zr and Hf in deep sea water, and if we postulate 6-fold surface coordination of Zr and Hf on ferromanganese crusts, then the observed trend of lower Zr/Hf ratios in the precipitate compared to the solution could be explained by mass dependent “pseudo isotope” fractionation. However, experimental investigation of the involved coordination environments is needed to test the plausibility of this hypothesis.

Our data also allow a second interpretation that is based on mass independent differences in bond strength rather than mass. According to our interpretation of the EXAFS data, the Hf complexes have slightly but significantly stronger and consequently shorter bonds of Hf to the fluoride and hydroxide ligands, which would be expected to push the H₂O ligands out more easily due to steric hindrance. This may be caused by the higher reduced mass of the ligands' bonds with Hf compared to Zr, but may also be the result of electronic effects. Our observations of different bond lengths to individual ligands between Zr and Hf as central atoms are supported by recent *ab initio* density functional theory (DFT) studies that have calculated higher cohesive energies for oxides and fluorides of Hf compared to the respective Zr counterparts (Anders et al., 2022). The formation of HfO₂ from the elements was calculated to be 0.261 eV more favorable than that of ZrO₂. The formation of HfF₄ was 0.274 eV more favorable than that of ZrF₄, meaning that Hf has a higher affinity for O and F compared to Zr and that this effect is even more pronounced for F than for O. Applied to our complex data, this means
<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Top: Schematic depiction of the incorporation pathways of complexes from solution into minerals. Zirconium is shown in a regular octahedral coordination, Hf in asymmetric octahedral coordination. The more significant rearrangement of ligands necessary for the more asymmetrical complexes (Hf) may hinder their incorporation compared to regular octahedral complexes and lead to mass-dependent or mass-independent fractionation of Zr and Hf. Bottom: Simulation of the development of Zr/Hf ratios in the precipitating minerals during a Rayleigh fractionation process as a function of the proportion of Zr that has already been precipitated (only the final 5 % of the precipitation process are shown). The five curves depict the Zr/Hf ratios in minerals precipitating from an initial solution with a chondritic Zr/Hf ratio of 33 by mass for different ratios of fractionation coefficients between mineral and solution ($D_{Hf}/D_{Zr}$) from 0.5 to 1. Note that all curves except $D_{Hf}/D_{Zr}$ = 1 eventually trend to zero, meaning extreme relative Hf enrichment, as Zr precipitation approaches 100 %.</figcaption></figure>

that the bond strength controls the distortion of the octahedral complexes and eventually the proportion of lower coordinated complexes in solution. Surface complexes on and metal coordination in minerals on the other hand are partially or fully controlled by the crystal structure and are typically 6-fold or higher (see above). Hafnium would require more significant changes to its coordination environment than Zr in order to be incorporated into a mineral from solution, meaning it is more likely to be retained in solution. Therefore, a mass independent bond-strength related fractionation mechanism would cause the same trend as the mass dependent fractionation discussed above (Fig. 10). The different effect of Zr and Hf on the bonds of their coordinating anions

can also be seen in the different vibrational frequencies observable with Raman spectroscopy (Grüneberger et al., 2016).

Which of the above discussed effects is more relevant for Zr/Hf fractionation in natural systems is not obvious from our data and both may indeed occur simultaneously. The observed effect of bond-strength is likely to affect also mineral surface complexation in the temperature range of this study. Mass dependent isotope fractionation is known to diminish with increasing temperature, following an Arrhenius-type 1/$T_{absolut}$ [1/K] law, typical of equilibrium fractionation processes, but absolute fractionation factors for Zr/Hf are currently unknown. One may speculate then that mass independent bond-strength plays a more significant role compared to mass difference for Zr/Hf fractionation at higher temperatures. This is likely even more relevant for Zr/Hf and Nb/$T_{Ta}$ fractionation in magmatic or high-grade metamorphic systems (e.g., Weyer et al., 2003; Pfänder et al., 2007; Lucassen et al., 2010; John et al., 2011), which may also be explained by the mechanisms proposed here if these elements change coordination during exchange between melt and mineral (e.g., zircon) or between two minerals (e.g., rutile and titanite).

## 5. Conclusions

Our data show that Zr and Hf occur in fluoride-rich hydrothermal solutions predominantly as neutral complexes, at 40 MPa and up to 400 °C. We conclude that these complexes take the form $[Zr(F, OH)_{4} \cdot 2H_{2}O]^{0}$ and $[Hf(F,OH)_{4} \cdot 2H_{2}O]^{0}$, respectively, in the experimentally investigated temperature range. As temperature increases, the octahedral complexes become more asymmetrical, with the fluoride and hydroxide ligands moving closer to the central cation and the $H_{2}O$ ligands moving further away. This may ultimately lead to partial transition to less than octahedral coordination, at high temperature. These findings suggest two possible mechanistic hypotheses for Zr/Hf fractionation in hydrothermal systems during the transfer of $Zr^{4+}$ and $Hf^{4+}$ in complexes in aqueous solution into the structure of minerals, both of which would produce the same general fractionation trend and are equally compatible with published observation from natural system: 1) Mass dependent equilibrium or kinetic fractionation due to coordination increase, as is known for isotope fractionation; or 2) mass independent differences in average coordination number due to different bond strength, which lead to different likelihood of incorporation into minerals due to different degrees of necessary bond reorganization (reaction kinetic effect). Due to the inverse temperature dependence of mass dependent fractionation, the bond-strength effect may dominate Zr/Hf fractionation at higher temperatures but further research is required to investigate the relative importance of the two effects under various conditions.

## Declaration of competing interest

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

## Acknowledgments

We acknowledge DESY (Hamburg, Germany), a member of the Helmholtz Association HGF, for the provision of experimental facilities and support within the Center for Molecular Water Science (CMWS). Parts of this research were carried out at DESY beamline P65, and we would like to thank W. Caliebe (P64), M. Hermann, and R. Nemausat for their excellent support and A. Ciobanu for her generous help and sterling assistance in the chemical laboratories at DESY. This project was funded by the German Research Foundation (Gefördert durch die Deutsche Forschungsgemeinschaft (DFG) – Projektnummer 387284271 – SFB 1349), projects A05 and A07, and by Bundesministerium für Bildung und Forschung (BMBF) (Grant No. FKZ 05K16PMA), DESY Center for
<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Top: Schematic depiction of the incorporation pathways of complexes from solution into minerals. Zirconium is shown in a regular octahedral coordination, Hf in asymmetric octahedral coordination. The more significant rearrangement of ligands necessary for the more asymmetrical complexes (Hf) may hinder their incorporation compared to regular octahedral complexes and lead to mass-dependent or mass-independent fractionation of Zr and Hf. Bottom: Simulation of the development of Zr/Hf ratios in the precipitating minerals during a Rayleigh fractionation process as a function of the proportion of Zr that has already been precipitated (only the final 5 % of the precipitation process are shown). The five curves depict the Zr/Hf ratios in minerals precipitating from an initial solution with a chondritic Zr/Hf ratio of 33 by mass for different ratios of fractionation coefficients between mineral and solution ($D_{Hf}/D_{Zr}$) from 0.5 to 1. Note that all curves except $D_{Hf}/D_{Zr}$ = 1 eventually trend to zero, meaning extreme relative Hf enrichment, as Zr precipitation approaches 100 %.</figcaption></figure>

that the bond strength controls the distortion of the octahedral complexes and eventually the proportion of lower coordinated complexes in solution. Surface complexes on and metal coordination in minerals on the other hand are partially or fully controlled by the crystal structure and are typically 6-fold or higher (see above). Hafnium would require more significant changes to its coordination environment than Zr in order to be incorporated into a mineral from solution, meaning it is more likely to be retained in solution. Therefore, a mass independent bond-strength related fractionation mechanism would cause the same trend as the mass dependent fractionation discussed above (Fig. 10). The different effect of Zr and Hf on the bonds of their coordinating anions

can also be seen in the different vibrational frequencies observable with Raman spectroscopy (Grüneberger et al., 2016).

Which of the above discussed effects is more relevant for Zr/Hf fractionation in natural systems is not obvious from our data and both may indeed occur simultaneously. The observed effect of bond-strength is likely to affect also mineral surface complexation in the temperature range of this study. Mass dependent isotope fractionation is known to diminish with increasing temperature, following an Arrhenius-type 1/T<sub>absolut</sub> [1/K] law, typical of equilibrium fractionation processes, but absolute fractionation factors for Zr/Hf are currently unknown. One may speculate then that mass independent bond-strength plays a more significant role compared to mass difference for Zr/Hf fractionation at higher temperatures. This is likely even more relevant for Zr/Hf and Nb/Ta fractionation in magmatic or high-grade metamorphic systems (e.g., Weyer et al., 2003; Pfänder et al., 2007; Lucassen et al., 2010; John et al., 2011), which may also be explained by the mechanisms proposed here if these elements change coordination during exchange between melt and mineral (e.g., zircon) or between two minerals (e.g., rutile and titanite).

## 5. Conclusions

Our data show that Zr and Hf occur in fluoride-rich hydrothermal solutions predominantly as neutral complexes, at 40 MPa and up to 400 °C. We conclude that these complexes take the form [Zr(F, OH)₄·2H₂O]⁰ and [Hf(F,OH)₄·2H₂O]⁰, respectively, in the experimentally investigated temperature range. As temperature increases, the octahedral complexes become more asymmetrical, with the fluoride and hydroxide ligands moving closer to the central cation and the H₂O ligands moving further away. This may ultimately lead to partial transition to less than octahedral coordination, at high temperature. These findings suggest two possible mechanistic hypotheses for Zr/Hf fractionation in hydrothermal systems during the transfer of Zr⁴⁺ and Hf⁴⁺ in complexes in aqueous solution into the structure of minerals, both of which would produce the same general fractionation trend and are equally compatible with published observation from natural system: 1) Mass dependent equilibrium or kinetic fractionation due to coordination increase, as is known for isotope fractionation; or 2) mass independent differences in average coordination number due to different bond strength, which lead to different likelihood of incorporation into minerals due to different degrees of necessary bond reorganization (reaction kinetic effect). Due to the inverse temperature dependence of mass dependent fractionation, the bond-strength effect may dominate Zr/Hf fractionation at higher temperatures but further research is required to investigate the relative importance of the two effects under various conditions.

## Declaration of competing interest

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

## Acknowledgments

We acknowledge DESY (Hamburg, Germany), a member of the Helmholtz Association HGF, for the provision of experimental facilities and support within the Center for Molecular Water Science (CMWS). Parts of this research were carried out at DESY beamline P65, and we would like to thank W. Caliebe (P64), M. Hermann, and R. Nemausat for their excellent support and A. Ciobanu for her generous help and sterling assistance in the chemical laboratories at DESY. This project was funded by the German Research Foundation (Gefördert durch die Deutsche Forschungsgemeinschaft (DFG) – Projektnummer 387284271 – SFB 1349), projects A05 and A07, and by Bundesministerium für Bildung und Forschung (BMBF) (Grant No. FKZ 05K16PMA), DESY Center for
Molecular Water Science—Early Science Project. Computing time for AIMD simulations was provided by Gauss Centre for Supercomputing e. V. through the John von Neumann Institute for Computing (NIC) on the JUWELS supercomputer at Jülich Supercomputing Centre (JSC). Computing for the force-field based MD simulations was provided by the Zentraleinrichtung für Datenverarbeitung (ZEDAT) of Freie Universität Berlin.

The authors thank Artas Migdisov, Mathieu Chassé, and one anonymous reviewer for their very constructive comments and suggestions, as well as Xiandong Liu for editorial handling of the manuscript.

## Appendix A. Supplementary material

Supplementary material contains method details for the classical MD calculations, as well as figures showing a Lennard-Jones parameter scan to demonstrate the dependence of calculated mean ligand distances on the choice of these parameters (Fig. S1), the concentration of Zr and Hf in solution during the XAS measurements (Fig. S2), example EXAFS fits demonstrating the influence of ligand number (Fig. S3), the effect of ligand configuration on the MD calculation results (Fig. S4 and S5), and a comparison of simulated spectra of Zr in 5- and 6-fold coordination.

Supplementary material to this article can be found online at [https://doi.org/10.1016/j.gca.2023.12.013](https://doi.org/10.1016/j.gca.2023.12.013).

## References

Abascal, J.L.F., Vega, C., 2005. A general purpose model for the condensed phases of water: TIP4P/2005. J. Chem. Phys. 123, 234505.

Abraham, M.J., Murtola, T., Schulz, R., Päll, S., Smith, J.C., Hess, B., Lindahl, E., 2015. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. Software X 1–2, 19–25.

Aja, S.U., Wood, S.A., Williams-Jones, A.E., 1995. The aqueous geochemistry of Zr and the solubility of some Zr-bearing minerals. Appl. Geochem. 10, 603–620.

Anders, J., Göritz, F., Loges, A., John, T., Paulus, B., 2022. Stability of hydroxo/oxo/ fluoro zirconates vs. hafniates—a DFT study. Inorganics 10, 259.

Aseri, A.A., Linnen, R.L., Che, X.D., Thibault, Y., Holtz, F., 2015. Effects of fluorine on the solubilities of Nb, Ta, Zr and Hf minerals in highly fluxed water-saturated haplogranitic melts. Ore Geol. Rev. 64, 736–746.

Bau, M., 1996. Controls on the fractionation of isovalent trace elements in magmatic and aqueous systems: evidence from Y/Ho, Zr/Hf, and lanthanide tetrad effect. Contr. Mineral. and Petrol. 123, 323–333.

Bindeman, I.N., Melnik, O.E., 2022. The rises and falls of zirconium isotopes during zircon crystallisation. Geochem. Persp. Lett. 24, 17–21.

Brugger, J., Liu, W., Etschmann, B., Mei, Y., Sherman, D.M., Testemale, D., 2016. A review of the coordination chemistry of hydrothermal systems, or do coordination changes make ore deposits? Chem. Geol. 447, 219–253.

Bunäu, O., Joly, Y., 2009. Self-consistent aspects of x-ray absorption calculations. J. Phys. Condens. Matter 21, 345501.

Bussi, G., Donadio, D., Parrinello, M., 2007. Canonical sampling through velocity rescaling. J. Chem. Phys. 126, 014101.

Censi, P., Raso, M., Yechieli, Y., Ginat, H., Saiano, F., Zuddas, P., Brusca, L., D'Alessandro, W., Inguaggiato, C., 2017. Geochemistry of Zr, Hf, and REE in a wide spectrum of Eh and water composition: The case of Dead Sea Fault system (Israel). Geochem. Geophys. Geosyst. 18, 844–857.

Chaboy, J., Marcelli, A., Tyson, T.A., 1994. Influence of double-electron transitions on the EXAFS L edges of rare-earth systems. Phys. Rev. B 49, 11652–11661.

Chen, W., Wallace, J.A., Yue, Z., Shen, J.K., 2013. Introducing titratable water to all-atom molecular dynamics at constant pH. Biophys. J. 105, L15–L17.

Darden, T., York, D., Pedersen, L., 1993. Particle mesh Ewald: An N-log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092.

Elam, W.T., Ravel, B.D., Sieber, J.R., 2002. A new atomic database for X-ray spectroscopic calculations. Rad. Phys. Chem. 63, 121–128.

Farges, F., 1991. Structural environment around Th4+ in silicate glasses: Implications for the geochemistry of incompatible Me4+ elements. Geochem. Cosmochim. Acta 55, 3303–3319.

Farges, F., 1996. Does Zr-F “complexation” occur in magmas? Chem. Geol. 127, 253–268.

Filipponi, A., 1995. Double-electron excitation effects above inner shell X-ray absorption edges. Physica B Condens. Matter 208–209, 29–32.

Goedecker, S., Teter, M., Hutter, J., 1996. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B 54, 1703–1710.

Grimme, S., Antony, J., Ehrlich, S., Krieg, H., 2010. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 132, 154104.

Grüneberger, A.M., Schmidt, C., Jahn, S., Rhede, D., Loges, A., Wilke, M., 2016. Interpretation of Raman spectra of the zircon-hafnion solid solution. Eur. J. Mineral. 28, 721–733.

Hagfeldt, C., Kessler, V., Persson, I., 2004. Structure of the hydrated, hydrolysed and solvated zirconium(IV) and hafnium(IV) ions in water and aprotic oxygen donor solvents. A crystallographic, EXAFS spectroscopic and large angle X-ray scattering study. Dalton Trans. 2004, 2142–2151.

Hess, B., Bekker, H., Berendsen, H.J.C., Fraaije, J.G.E.M., 1997. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 18, 1463–1472.

Hubbell, J.H., Seltzer, S.M., 1996. NIST standard reference database 126. National Institute of Standards and Technology, Gaithersburg, MD.

Ibanez-Mejia, M., Tissot, F.L.H., 2019. Extreme Zr stable isotope fractionation during magmatic fractional crystallization. Sci. Adv. 5, eaaax684.

Inglis, E.C., Moynier, F., Creech, J., Deng, Z., Day, J.M.D., Teng, F.-Z., Bizzarro, M., Jackson, M., Savage, P., 2019. Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth. Geochem. Cosmochim. Acta 250, 311–323.

Inguaggiato, C., Censi, P., Zuddas, P., Londoño, J.M., Chacón, Z., Alzate, D., Brusca, L., D'Alessandro, W., 2015. Geochemistry of REE, Zr and Hf in a wide range of pH and water composition: The Nevado del Ruiz volcano-hydrothermal system (Colombia). Chem. Geol. 417, 125–133.

Jahn, S., Dubrail, J., Wilke, M., 2015. Complexation of Zr and Hf monomers in supercritical aqueous solutions: Insights from ab initio molecular dynamics simulations. Chem. Geol. 418, 30–39.

John, T., Klemd, R., Klemme, S., Pfänder, J.A., Elis, H.J., Gao, J., 2011. Nb-Ta fractionation by partial melting at the titanite-rutile transition. Contrib. Mineral. Petrol. 161, 35–45.

Joly, Y., 2001. X-ray absorption near-edge structure calculations beyond the muffin-tin approximation. Phys. Rev. B 63, 125120.

Jonchiere, R., Seitsonen, A.P., Ferlat, G., Saitta, A.M., Vuilleumier, R., 2011. Van der Waals effects in ab initio water at ambient and supercritical conditions. J. Chem. Phys. 135, 154503.

Juillot, F., Maréchal, C., Ponthieu, M., Calacay, S., Morin, G., Benedetti, M., Hazemann, J. L., Proux, O., Guyot, F., 2008. Zn isotopic fractionation caused by sorption on goethite and 2-Lines ferrihydrite. Geochem. Cosmochim. Acta 72, 4886–4900.

Kanazhevskii, V.V., Shmachkova, V.P., Kotsarenko, N.S., Kolomiichuk, V.N., Kochubei, D.I., 2006. Changes in the zirconium local surrounding on ligand substitution in solutions. J. Struct. Chem. 47, 860–868.

Klemme, S., Feldhaus, M., Potapkin, V., Wilke, M., Borchert, M., Louvel, M., Loges, A., Rohrbach, A., Weitkamp, P., Welter, E., Kokh, M., Schmidt, C., Testemale, D., 2021. A hydrothermal apparatus for x-ray absorption spectroscopy of hydrothermal fluids at DESY. Rev. Sci. Instrum. 92, 063903.

Kudryashov, N.M., Skublov, S.G., Galankina, O.L., Udaratina, O.V., Voloshin, A.V., 2020. Abnormally high-hafnium zircon from rare-metal pegmatites of the Vasin-Mylk deposit (the northeastern part of the Kola Peninsula). Geochem. 80, 125489.

Kühne, T.D., Iannuzzi, M., Del Ben, M., Rybkin, V.V., Seewald, P., Stein, F., Laino, T., Khaliullin, R.Z., Schütt, O., Schiffmann, F., Golze, D., Wilhelm, J., Chulkov, S., Bani-Hashemian, M.H., Weber, V., Borstnik, U., Taillefumier, M., Jakobovits, A.S., Lazzaro, A., Pabst, H., Müller, T., Schade, R., Guidon, M., Andermatt, S., Holberg, N., Schenter, G.K., Hehn, A., Bussy, A., Belcelflamme, F., Tabacchi, G., Gloß, A., Lass, M., Bethune, I., Mundy, C.J., Plessl, C., Watkins, M., VandeVondele, J., Krack, M., Hutter, J., 2020. CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations. J. Chem. Phys. 152, 194103.

Lemarchand, E., Schott, J., Gaillardet, J., 2007. How surface complexes impact boron isotope fractionation: Evidence from Fe and Mn oxides sorption experiments. Earth Planet. Sci. Lett. 260, 277–296.

Lemmon E. W., McLinden M. O. and Friend D. G. 2019. Thermophysical Properties of Fluid Systems. In: Linstrom, P.J. and Mallard, W.G. (Eds.), NIST Chemistry Workbook. NIST Standard Reference Database. National Institute of Standards and Technology, Gaithersburg MD.

Linnen, R.L., Keppler, H., 2002. Melt composition control of Zr/Hf fractionation in magmatic processes. Geochem. Cosmochim. Acta 66, 3293–3301.

Lippert, G., Hutter, J., Parrinello, M., 1997. A hybrid Gaussian and plane wave density functional scheme. Mol. Phys. 92, 477–487.

Louvel, M., Sanchez-Valle, C., Malfait, W.J., Testemale, D., Hazemann, J.-L., 2013. Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones. Geochem. Cosmochim. Acta 104, 281–299.

Lucassen, F., Dulski, P., Abart, R., Franz, G., Rhede, D., Romer, R.L., 2010. Redistribution of HFSE elements during rutile replacement by titanite. Contrib. Mineral. Petrol. 160, 279–295.

McCulloch, M.T., Gamble, J.A., 1991. Geochemical and geodynamical constraints on subduction zone magmatism. Earth Planet. Sci. Lett. 102, 358–374.

Méheut, M., Ibañez-Mejia, M., Tissot, F.L.H., 2021. Drivers of zirconium isotope fractionation in Zr-bearing phases and melts: The roles of vibrational, nuclear field shift and diffusive effects. Geochem. Cosmochim. Acta 292, 217–234.

Messner, C.B., Hofer, T.S., Randolf, B.R., Rode, B.M., 2011. Structure and dynamics of the Zr4+ ion in water. Phys. Chem. Chem. Phys. 13, 224–229.

Migdisov, Art. A., Williams-Jones, A. E., van Hinsberg, V. and Salvi, S. 2011. An experimental study of the solubility of baddeleyite (ZrO₂) in fluoride-bearing solutions at elevated temperature. Geochem. Cosmochim. Acta 75, 7426–7434.

Neves, J.M.C., Nunes, J.E.L., Sahama, T.G., 1974. High hafnium members of the zircon-hafnion series from the granite pegmatites of Zambézia, Mozambique. Contrib. Mineral. Petrol. 48, 73–80.

Newville, M., 2014. Fundamentals of XAFS. Rev. Min. Geochem. 78, 33–74.

Ohta, A., Kagi, H., Tsuno, H., Nomura, M., Kawabe, I., 2008. Influence of multi-electron excitation on EXAFS spectroscopy of trivalent rare-earth ions and elucidation of change in hydration number through the series. Am. Mineral. 93, 1384–1392.
Molecular Water Science—Early Science Project. Computing time for AIMD simulations was provided by Gauss Centre for Supercomputing e. V. through the John von Neumann Institute for Computing (NIC) on the JUWELS supercomputer at Jülich Supercomputing Centre (JSC). Computing for the force-field based MD simulations was provided by the Zentraleinrichtung für Datenverarbeitung (ZEDAT) of Freie Universität Berlin.

The authors thank Artas Migdisov, Mathieu Chassé, and one anonymous reviewer for their very constructive comments and suggestions, as well as Xiandong Liu for editorial handling of the manuscript.

## Appendix A. Supplementary material

Supplementary material contains method details for the classical MD calculations, as well as figures showing a Lennard-Jones parameter scan to demonstrate the dependence of calculated mean ligand distances on the choice of these parameters (Fig. S1), the concentration of Zr and Hf in solution during the XAS measurements (Fig. S2), example EXAFS fits demonstrating the influence of ligand number (Fig. S3), the effect of ligand configuration on the MD calculation results (Fig. S4 and S5), and a comparison of simulated spectra of Zr in 5- and 6-fold coordination.

Supplementary material to this article can be found online at https://doi.org/10.1016/j.gca.2023.12.013.

## References

Abascal, J.L.F., Vega, C., 2005. A general purpose model for the condensed phases of water: TIP4P/2005. J. Chem. Phys. 123, 234505.

Abraham, M.J., Murtola, T., Schulz, R., Päll, S., Smith, J.C., Hess, B., Lindahl, E., 2015. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1–2, 19–25.

Aja, S.U., Wood, S.A., Williams-Jones, A.E., 1995. The aqueous geochemistry of Zr and the solubility of some Zr-bearing minerals. Appl. Geochem. 10, 603–620.

Anders, J., Göritz, F., Loges, A., John, T., Paulus, B., 2022. Stability of hydroxo/oxo/fluoro zirconates vs. hafniates—a DFT study. Inorganics 10, 259.

Aseri, A.A., Linnen, R.L., Che, X.D., Thibault, Y., Holtz, F., 2015. Effects of fluorine on the solubilities of Nb, Ta, Zr and Hf minerals in highly fluxed water-saturated haplogranitic melts. Ore Geol. Rev. 64, 736–746.

Bau, M., 1996. Controls on the fractionation of isovalent trace elements in magmatic and aqueous systems: evidence from Y/Ho, Zr/Hf, and lanthanide tetrad effect. Contr. Mineral. and Petrol. 123, 323–333.

Bindeman, I.N., Melnik, O.E., 2022. The rises and falls of zirconium isotopes during zircon crystallisation. Geochem. Persp. Lett. 24, 17–21.

Brugger, J., Liu, W., Etschmann, B., Mei, Y., Sherman, D.M., Testemale, D., 2016. A review of the coordination chemistry of hydrothermal systems, or do coordination changes make ore deposits? Chem. Geol. 447, 219–253.

Bunäu, O., Joly, Y., 2009. Self-consistent aspects of x-ray absorption calculations. J. Phys. Condens. Matter 21, 345501.

Bussi, G., Donadio, D., Parrinello, M., 2007. Canonical sampling through velocity rescaling. J. Chem. Phys. 126, 014101.

Censi, P., Raso, M., Yechieli, Y., Ginat, H., Saiano, F., Zuddas, P., Brusca, L., D'Alessandro, W., Inguaggiato, C., 2017. Geochemistry of Zr, Hf, and REE in a wide spectrum of Eh and water composition: The case of Dead Sea Fault system (Israel). Geochem. Geophys. Geosyst. 18, 844–857.

Chaboy, J., Marcelli, A., Tyson, T.A., 1994. Influence of double-electron transitions on the EXAFS L edges of rare-earth systems. Phys. Rev. B 49, 11652–11661.

Chen, W., Wallace, J.A., Yue, Z., Shen, J.K., 2013. Introducing titratable water to all-atom molecular dynamics at constant pH. Biophys. J. 105, L15–L17.

Darden, T., York, D., Pedersen, L., 1993. Particle mesh Ewald: An N-log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092.

Elam, W.T., Ravel, B.D., Sieber, J.R., 2002. A new atomic database for X-ray spectroscopic calculations. Rad. Phys. Chem. 63, 121–128.

Farges, F., 1991. Structural environment around Th4+ in silicate glasses: Implications for the geochemistry of incompatible Me4+ elements. Geochem. Cosmochim. Acta 55, 3303–3319.

Farges, F., 1996. Does Zr-F “complexation” occur in magmas? Chem. Geol. 127, 253–268.

Filipponi, A., 1995. Double-electron excitation effects above inner shell X-ray absorption edges. Physica B Condens. Matter 208–209, 29–32.

Goedecker, S., Teter, M., Hutter, J., 1996. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B 54, 1703–1710.

Grimme, S., Antony, J., Ehrlich, S., Krieg, H., 2010. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 132, 154104.

Grüneberger, A.M., Schmidt, C., Jahn, S., Rhede, D., Loges, A., Wilke, M., 2016. Interpretation of Raman spectra of the zircon-hafnion solid solution. Eur. J. Mineral. 28, 721–733.

Hagfeldt, C., Kessler, V., Persson, I., 2004. Structure of the hydrated, hydrolysed and solvated zirconium(IV) and hafnium(IV) ions in water and aprotic oxygen donor solvents. A crystallographic, EXAFS spectroscopic and large angle X-ray scattering study. Dalton Trans. 2004, 2142–2151.

Hess, B., Bekker, H., Berendsen, H.J.C., Fraaije, J.G.E.M., 1997. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 18, 1463–1472.

Hubbell, J.H., Seltzer, S.M., 1996. NIST standard reference database 126. National Institute of Standards and Technology, Gaithersburg, MD.

Ibanez-Mejia, M., Tissot, F.L.H., 2019. Extreme Zr stable isotope fractionation during magmatic fractional crystallization. Sci. Adv. 5, eaxa8648.

Inglis, E.C., Moynier, F., Creech, J., Deng, Z., Day, J.M.D., Teng, F.-Z., Bizzarro, M., Jackson, M., Savage, P., 2019. Isotopic fractionation of zirconium during magmatic differentiation and the stable isotope composition of the silicate Earth. Geochem. Cosmochim. Acta 250, 311–323.

Inguaggiato, C., Censi, P., Zuddas, P., Londoño, J.M., Chacón, Z., Alzate, D., Brusca, L., D'Alessandro, W., 2015. Geochemistry of REE, Zr and Hf in a wide range of pH and water composition: The Nevado del Ruiz volcano-hydrothermal system (Colombia). Chem. Geol. 417, 125–133.

Jahn, S., Dubrail, J., Wilke, M., 2015. Complexation of Zr and Hf monomers in supercritical aqueous solutions: Insights from ab initio molecular dynamics simulations. Chem. Geol. 418, 30–39.

John, T., Klemd, R., Klemme, S., Pfänder, J.A., Elis, H.J., Gao, J., 2011. Nb-Ta fractionation by partial melting at the titanite-rutile transition. Contrib. Mineral. Petrol. 161, 35–45.

Joly, Y., 2001. X-ray absorption near-edge structure calculations beyond the muffin-tin approximation. Phys. Rev. B 63, 125120.

Jonchiere, R., Seitsonen, A.P., Ferlat, G., Saitta, A.M., Vuilleumier, R., 2011. Van der Waals effects in ab initio water at ambient and supercritical conditions. J. Chem. Phys. 135, 154503.

Juillot, F., Maréchal, C., Ponthieu, M., Calacay, S., Morin, G., Benedetti, M., Hazemann, J. L., Proux, O., Guyot, F., 2008. Zn isotopic fractionation caused by sorption on goethite and 2-Lines ferrihydrite. Geochem. Cosmochim. Acta 72, 4886–4900.

Kanazhevskii, V.V., Shmachkova, V.P., Kotsarenko, N.S., Kolomiichuk, V.N., Kochubei, D.I., 2006. Changes in the zirconium local surrounding on ligand substitution in solutions. J. Struct. Chem. 47, 860–868.

Klemme, S., Feldhaus, M., Potapkin, V., Wilke, M., Borchert, M., Louvel, M., Loges, A., Rohrbach, A., Weitkamp, P., Welter, E., Kokh, M., Schmidt, C., Testemale, D., 2021. A hydrothermal apparatus for x-ray absorption spectroscopy of hydrothermal fluids at DESY. Rev. Sci. Instrum. 92, 063903.

Kudryashov, N.M., Skublov, S.G., Galankina, O.L., Udaratina, O.V., Voloshin, A.V., 2020. Abnormally high-hafnium zircon from rare-metal pegmatites of the Vasin-Mylk deposit (the northeastern part of the Kola Peninsula). Geochem. 80, 125489.

Kühne, T.D., Iannuzzi, M., Del Ben, M., Rybkin, V.V., Seewald, P., Stein, F., Laino, T., Khaliullin, R.Z., Schütt, O., Schiffmann, F., Golze, D., Wilhelm, J., Chulkov, S., Bani-Hashemian, M.H., Weber, V., Borštnik, U., Taillefumier, M., Jakobovits, A.S., Lazzaro, A., Pabst, H., Müller, T., Schade, R., Guidon, M., Andermatt, S., Holberg, N., Schenter, G.K., Hehn, A., Bussy, A., Belflamme, F., Tabacchi, G., Gloß, A., Lass, M., Bethune, I., Mundy, C.J., Plessl, C., Watkins, M., VandeVondele, J., Krack, M., Hutter, J., 2020. CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations. J. Chem. Phys. 152, 194103.

Lemarchand, E., Schott, J., Gaillardet, J., 2007. How surface complexes impact boron isotope fractionation: Evidence from Fe and Mn oxides sorption experiments. Earth Planet. Sci. Lett. 260, 277–296.

Lemmon E. W., McLinden M. O. and Friend D. G. 2019. Thermophysical Properties of Fluid Systems. In: Linstrom, P.J. and Mallard, W.G. (Eds.), NIST Chemistry Workbook. NIST Standard Reference Database. National Institute of Standards and Technology, Gaithersburg MD.

Linnen, R.L., Keppler, H., 2002. Melt composition control of Zr/Hf fractionation in magmatic processes. Geochem. Cosmochim. Acta 66, 3293–3301.

Lippert, G., Hutter, J., Parrinello, M., 1997. A hybrid Gaussian and plane wave density functional scheme. Mol. Phys. 92, 477–487.

Louvel, M., Sanchez-Valle, C., Malfait, W.J., Testemale, D., Hazemann, J.-L., 2013. Zr complexation in high pressure fluids and silicate melts and implications for the mobilization of HFSE in subduction zones. Geochem. Cosmochim. Acta 104, 281–299.

Lucassen, F., Dulski, P., Abart, R., Franz, G., Rhede, D., Romer, R.L., 2010. Redistribution of HFSE elements during rutile replacement by titanite. Contrib. Mineral. Petrol. 160, 279–295.

McCulloch, M.T., Gamble, J.A., 1991. Geochemical and geodynamical constraints on subduction zone magmatism. Earth Planet. Sci. Lett. 102, 358–374.

Méheut, M., Ibanez-Mejia, M., Tissot, F.L.H., 2021. Drivers of zirconium isotope fractionation in Zr-bearing phases and melts: The roles of vibrational, nuclear field shift and diffusive effects. Geochem. Cosmochim. Acta 292, 217–234.

Messner, C.B., Hofer, T.S., Randolf, B.R., Rode, B.M., 2011. Structure and dynamics of the Zr4+ ion in water. Phys. Chem. Chem. Phys. 13, 224–229.

Migdisov, Art. A., Williams-Jones, A. E., van Hinsberg, V. and Salvi, S. 2011. An experimental study of the solubility of baddeleyite (ZrO₂) in fluoride-bearing solutions at elevated temperature. Geochem. Cosmochim. Acta 75, 7426–7434.

Neves, J.M.C., Nunes, J.E.L., Sahama, T.G., 1974. High hafnium members of the zircon-hafnion series from the granite pegmatites of Zambézia, Mozambique. Contrib. Mineral. Petrol. 48, 73–80.

Newville, M., 2014. Fundamentals of XAFS. Rev. Min. Geochem. 78, 33–74.

Ohta, A., Kagi, H., Tsuno, H., Nomura, M., Kawabe, I., 2008. Influence of multi-electron excitation on EXAFS spectroscopy of trivalent rare-earth ions and elucidation of change in hydration number through the series. Am. Mineral. 93, 1384–1392.
Orabi, E.A., Faraldo-Gómez, J.D., 2020. New molecular-mechanics model for simulations of hydrogen fluoride in chemistry and biology. J. Chem. Theory Comput. 16, 5105–5126.

Palme, H. and O'Neill, H. St. C. 2014. Cosmochemical estimates of mantle composition. In: Treatise on Geochemistry Elsevier. pp. 1–39.

Pearson, R.G., 1963. Hard and soft acids and bases. J. Am. Chem. Soc. 85, 3533–3539.

Perdew, J.P., Burke, K., Ernzerhof, M., 1996. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868.

Perfit, M.R., Gust, D.A., Bence, A.E., Arculus, R.J., Taylor, S.R., 1980. Chemical characteristics of island-arc basalts: Implications for mantle sources. Chem. Geol. 30, 227–256.

Pfänder, J.A., Münker, C., Stracke, A., Mezger, K., 2007. Nb/Ta and Zr/Hf in ocean island basalts — Implications for crust–mantle differentiation and the fate of Niobium. Earth Planet. Sci. Lett. 254, 158–172.

Pokrovski, G.S., Roux, J., Hazemann, J.-L., Testemale, D., 2005. An X-ray absorption spectroscopy study of argutite solubility and aqueous Ge(IV) speciation in hydrothermal fluids to 500 °C and 400 bar. Chem. Geol. 217, 127–145.

Rai, D., Xia, Y., Hess, N.J., Strachan, D.M., McGrail, B.P., 2001. Hydroxo and Chloro Complexes/Ion Interactions of Hf4+ and the Solubility Product of HfO2(am). J. Solution Chem. 19.

Rehr, J.J., Mustre de Leon, J., Zabinsky, S.I., Albers, R.C., 1991. Theoretical x-ray absorption fine structure standards. J. Am. Chem. Soc. 113, 5135–5140.

Roy, M., Lindsay, J.D., Louch, S., Gurman, S.J., 2001. Multiple-electron excitation in X-ray absorption: a simple generic model. J. Synchrotron Rad. 8, 1103–1108.

Ryzhenko, B.N., Kovalenko, N.I., Prisyagina, N.I., Starshinova, N.P., Krupskaya, V.V., 2008. Experimental determination of zirconium speciation in hydrothermal solutions. Geochem. Int. 46, 328–339.

Schmidt, K., Bau, M., Hein, J.R., Koschinsky, A., 2014. Fractionation of the geochemical twins Zr–Hf and Nb–Ta during scavenging from seawater by hydrogenetic ferromanganese crusts. Geochim. Cosmochim. Acta 140, 468–487.

Shannon, R.D., 1976. Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides. Acta Crystallogr. A 32, 751–767.

Shikina, N.D., Tagirov, B.R., Volchenkova, V.A., Bychkova Ya, V., 2015. Zr/Hf ratio in supercritical chloride fluids: Experimental study of zirconium and hafnium complexation at 450°C and 0.6–1 kbar. Petrol. 23, 93–101.

Steinhardt, P.J., Nelson, D.R., Ronchetti, M., 1983. Bond-orientational order in liquids and glasses. Phys. Rev. B 28, 784–805.

Taylor, S.R., McLennan, S.M., 1995. The geochemical evolution of the continental crust. Rev. Geophys. 33, 241.

Tissot, F.L.H., Ibañez-Mejia, M., 2021. Unlocking the Single-Crystal Record of Heavy Stable Isotopes. Elements 17, 389–394.

Tomkins, H.G.D., Ibañez-Mejia, M., Tissot, F.L.H., Bloch, E., Wang, Y., Trail, D., 2023. Zircon growth experiments reveal limited equilibrium Zr isotope fractionation in magmas. Geochem. Persp. Lett. 25, 25–29.

VandeVondele, J., Hutter, J., 2007. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 127, 114105.

VandeVondele, J., Krack, M., Mohamed, F., Parrinello, M., Chassaing, T., Hutter, J., 2005. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 167, 103–128.

Wagner, W., Pruß, A., 2002. The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. J. Phys. Chem. Ref. Data 31, 387–535.

Wang, W.-C., Chen, Y., 1998. An improved method for the extraction of the X-ray absorption atomic background at the Zr K-edge. Phys. Stat. Sol. (a) 168, 351–357.

Wang, X., Griffin, W.L., Chen, J., 2010. Hf contents and Zr/Hf ratios in granitic zircons. Geochem. J. 44, 65–72.

Wasylenki, L.E., Weeks, C.L., Bargar, J.R., Spiro, T.G., Hein, J.R., Anbar, A.D., 2011. The molecular mechanism of Mo isotope fractionation during adsorption to birnessite. Geochim. Cosmochim. Acta 75, 5019–5031.

Weingärtner, H., Franck, E.U., 2005. Supercritical water as a solvent. Angew. Chem. Int. Ed. 44, 2672–2692.

Welter, E., Chernikov, R., Herrmann, M., Nemausat, R., 2019. A beamline for bulk sample x-ray absorption spectroscopy at the high brilliance storage ring PETRA III. AIP Conf. Proc. 2054, 040002.

Weyer, S., Münker, C., Mezger, K., 2003. Nb/Ta, Zr/Hf and REE in the depleted mantle: implications for the differentiation history of the crust–mantle system. Earth Planet. Sci. Lett. 205, 309–324.

Wieser, M.E., Coplen, T.B., 2011. Atomic weights of the elements 2009 (IUPAC Technical Report). Pure Appl. Chem. 83, 359–396.

Wilke, M., Schmidt, C., Dubrail, J., Appel, K., Borchert, M., Kvashnina, K., Manning, C.E., 2012. Zircon solubility and zirconium complexation in H2O+Na2O+SiO2+Al2O3 fluids at high pressure and temperature. Earth Planet. Sci. Lett. 349–350, 15–25.

Wilke, M., Jahn, S., Schmidt, C., Dubrail, J., Appel, K., Borchert, M., Kvash-nina, K., Pascarelli, S., Manning, C.E., 2013. Insights from X-ray absorption/fluorescence spectroscopy and ab-initio molecular dynamics on concentration and complexation of Zr and Hf in aqueous fluids at high pressure and temperature. J. Phys. Conf. Ser. 430, 012122.

Winterer M. 2020. xafsX: a program to process, analyse and reduce X-ray absorption fine structure spectra (XAFS). In International Tables for Crystallography.

Wunder, B., Meixner, A., Romer, R.L., Jahn, S., 2011. Li-isotope fractionation between silicates and fluids: Pressure dependence and influence of the bonding environment. Eur. J. Mineral. 23, 333–342.

Yin, R., Wang, R.C., Zhang, A.-C., Hu, H., Zhu, J.C., Rao, C., Zhang, H., 2013. Extreme fractionation from zircon to hafnon in the Koktokay No. 1 granitic pegmatite, Altai, northwestern China. Am. Mineral. 98, 1714–1724.

Zhu, Z., Zhang, W., Wang, J., Wang, Z., Guo, J.-L., Hoffmann, E.J., Feng, L., Luo, T., Hu, Z., Liu, Y., Moynier, F., 2023. Magmatic crystallization drives zircon Zr isotopic variations in a large granite batholith. Geochim. Cosmochim. Acta 342, 15–30.
Orabi, E.A., Faraldo-Gómez, J.D., 2020. New molecular-mechanics model for simulations of hydrogen fluoride in chemistry and biology. J. Chem. Theory Comput. 16, 5105–5126.

Palme, H. and O'Neill, H. St. C. 2014. Cosmochemical estimates of mantle composition. In: Treatise on Geochemistry Elsevier. pp. 1–39.

Pearson, R.G., 1963. Hard and soft acids and bases. J. Am. Chem. Soc. 85, 3533–3539.

Perdew, J.P., Burke, K., Ernzerhof, M., 1996. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868.

Perfit, M.R., Gust, D.A., Bence, A.E., Arculus, R.J., Taylor, S.R., 1980. Chemical characteristics of island-arc basalts: Implications for mantle sources. Chem. Geol. 30, 227–256.

Pfänder, J.A., Münker, C., Stracke, A., Mezger, K., 2007. Nb/Ta and Zr/Hf in ocean island basalts — Implications for crust–mantle differentiation and the fate of Niobium. Earth Planet. Sci. Lett. 254, 158–172.

Pokrovski, G.S., Roux, J., Hazemann, J.-L., Testemale, D., 2005. An X-ray absorption spectroscopy study of argutite solubility and aqueous Ge(IV) speciation in hydrothermal fluids to 500 °C and 400 bar. Chem. Geol. 217, 127–145.

Rai, D., Xia, Y., Hess, N.J., Strachan, D.M., McGrail, B.P., 2001. Hydroxo and Chloro Complexes/Ion Interactions of Hf4+ and the Solubility Product of HfO2(am). J. Solution Chem. 19.

Rehr, J.J., Mustre de Leon, J., Zabinsky, S.I., Albers, R.C., 1991. Theoretical x-ray absorption fine structure standards. J. Am. Chem. Soc. 113, 5135–5140.

Roy, M., Lindsay, J.D., Louch, S., Gurman, S.J., 2001. Multiple-electron excitation in X-ray absorption: a simple generic model. J. Synchrotron Rad. 8, 1103–1108.

Ryzhenko, B.N., Kovalenko, N.I., Prisyagina, N.I., Starshinova, N.P., Krupskaya, V.V., 2008. Experimental determination of zirconium speciation in hydrothermal solutions. Geochem. Int. 46, 328–339.

Schmidt, K., Bau, M., Hein, J.R., Koschinsky, A., 2014. Fractionation of the geochemical twins Zr–Hf and Nb–Ta during scavenging from seawater by hydrogenetic ferromanganese crusts. Geochim. Cosmochim. Acta 140, 468–487.

Shannon, R.D., 1976. Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides. Acta Crystallogr. A 32, 751–767.

Shikina, N.D., Tagirov, B.R., Volchenkova, V.A., Bychkova Ya, V., 2015. Zr/Hf ratio in supercritical chloride fluids: Experimental study of zirconium and hafnium complexation at 450°C and 0.6–1 kbar. Petrol. 23, 93–101.

Steinhardt, P.J., Nelson, D.R., Ronchetti, M., 1983. Bond-orientational order in liquids and glasses. Phys. Rev. B 28, 784–805.

Taylor, S.R., McLennan, S.M., 1995. The geochemical evolution of the continental crust. Rev. Geophys. 33, 241.

Tissot, F.L.H., Ibañez-Mejia, M., 2021. Unlocking the Single-Crystal Record of Heavy Stable Isotopes. Elements 17, 389–394.

Tomkins, H.G.D., Ibañez-Mejia, M., Tissot, F.L.H., Bloch, E., Wang, Y., Trail, D., 2023. Zircon growth experiments reveal limited equilibrium Zr isotope fractionation in magmas. Geochem. Persp. Lett. 25, 25–29.

VandeVondele, J., Hutter, J., 2007. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 127, 114105.

VandeVondele, J., Krack, M., Mohamed, F., Parrinello, M., Chassaing, T., Hutter, J., 2005. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 167, 103–128.

Wagner, W., Pruß, A., 2002. The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. J. Phys. Chem. Ref. Data 31, 387–535.

Wang, W.-C., Chen, Y., 1998. An improved method for the extraction of the X-ray absorption atomic background at the Zr K-edge. Phys. Stat. Sol. (a) 168, 351–357.

Wang, X., Griffin, W.L., Chen, J., 2010. Hf contents and Zr/Hf ratios in granitic zircons. Geochem. J. 44, 65–72.

Wasylenki, L.E., Weeks, C.L., Bargar, J.R., Spiro, T.G., Hein, J.R., Anbar, A.D., 2011. The molecular mechanism of Mo isotope fractionation during adsorption to birnessite. Geochim. Cosmochim. Acta 75, 5019–5031.

Weingärtner, H., Franck, E.U., 2005. Supercritical water as a solvent. Angew. Chem. Int. Ed. 44, 2672–2692.

Welter, E., Chernikov, R., Herrmann, M., Nemausat, R., 2019. A beamline for bulk sample x-ray absorption spectroscopy at the high brilliance storage ring PETRA III. AIP Conf. Proc. 2054, 040002.

Weyer, S., Münker, C., Mezger, K., 2003. Nb/Ta, Zr/Hf and REE in the depleted mantle: implications for the differentiation history of the crust–mantle system. Earth Planet. Sci. Lett. 205, 309–324.

Wieser, M.E., Coplen, T.B., 2011. Atomic weights of the elements 2009 (IUPAC Technical Report). Pure Appl. Chem. 83, 359–396.

Wilke, M., Schmidt, C., Dubrail, J., Appel, K., Borchert, M., Kvashnina, K., Manning, C.E., 2012. Zircon solubility and zirconium complexation in H2O+Na2O+SiO2+Al2O3 fluids at high pressure and temperature. Earth Planet. Sci. Lett. 349–350, 15–25.

Wilke, M., Jahn, S., Schmidt, C., Dubrail, J., Appel, K., Borchert, M., Kvash-nina, K., Pascarelli, S., Manning, C.E., 2013. Insights from X-ray absorption/fluorescence spectroscopy and ab-initio molecular dynamics on concentration and complexation of Zr and Hf in aqueous fluids at high pressure and temperature. J. Phys. Conf. Ser. 430, 012122.

Winterer M. 2020. xafsX: a program to process, analyse and reduce X-ray absorption fine structure spectra (XAFS). In International Tables for Crystallography.

Wunder, B., Meixner, A., Romer, R.L., Jahn, S., 2011. Li-isotope fractionation between silicates and fluids: Pressure dependence and influence of the bonding environment. Eur. J. Mineral. 23, 333–342.

Yin, R., Wang, R.C., Zhang, A.-C., Hu, H., Zhu, J.C., Rao, C., Zhang, H., 2013. Extreme fractionation from zircon to hafnon in the Koktokay No. 1 granitic pegmatite, Altai, northwestern China. Am. Mineral. 98, 1714–1724.

Zhu, Z., Zhang, W., Wang, J., Wang, Z., Guo, J.-L., Hoffmann, E.J., Feng, L., Luo, T., Hu, Z., Liu, Y., Moynier, F., 2023. Magmatic crystallization drives zircon Zr isotopic variations in a large granite batholith. Geochim. Cosmochim. Acta 342, 15–30.