Equilibrium Sn isotope fractionation between aqueous Sn and Sn-bearing minerals: Constrained by first-principles calculations

Tin is an important metal for science and technology, and it has been identified as a critical metal for the development of decarbonization technology (Moss et al. 2013). Most of the Sn production (99%) throughout the world originates from hydrothermal tin deposits, and, consequently, the ore-forming processes and mineral exploration of tin have attracted close attention (Yuan et al. 2019; Lehmann and Harmanto 1990; Lehmann 2021; Mao et al. 2019, 2021). A great deal of research has shown that Sn (Sn²⁺ and Sn⁴⁺) forms various complexes with different ligands in different geological fluids (e.g., Tobias 1958; Mark 1977; Pettine et al. 1981; Séby et al. 2001; Müller and Seward 2001; Uchida et al. 2002; Duc-Tin et al. 2007; Cigala et al. 2012; Schmidt 2018; Liu et al. 2020; She et al. 2020; Wang et al. 2021). According to limited experimental studies, the Sn isotope fractionations are obvious among aqueous Sn species and Sn-bearing minerals, indicating that Sn isotope geochemistry is a robust tool to trace the Sn-involved geological process (Polyakov et al. 2005; Brugger et al. 2016; Mason et al. 2016; Badullovich et al. 2017; Creech et al. 2017; Mathur et al. 2017; Yao et al. 2018; Wang et al. 2019a; Liu et al. 2021a). The equilibrium Sn isotopic fractionation factors among Sn-bearing substances are a prerequisite to use Sn isotopes to trace geological processes. However, equilibrium Sn isotopic fractionation factors between mineral and fluids and between mineral pairs are mostly lacking because Sn²⁺ is readily oxidized under environmental conditions (Wang et al. 2019a) and dividing Sn²⁺ and Sn⁴⁺ complexes in solutions is challenging (Altunay and Gürkan 2015), which makes it difficult to further decode the evolution process of Sn in geological systems.

Recently, the reduced partition function ratios (β factors) of Sn-Cl species in solution were calculated by using the cluster method (one configuration) (She et al. 2020; Wang et al. 2021; Sun et al. 2022), and the β values of Sn²⁺-Cl species and tin-bearing minerals (cassiterite, megawite, and romarchite) were calculated by two different methods: the cluster method and periodic boundary calculations (Wang et al. 2021), respectively. However, due to the dynamic properties of aqueous solutions, previous studies have demonstrated that the configuration disorder of the aqueous solutions can affect the β values (e.g., Li et al. 2009; Li and Liu 2011; Wang et al. 2019b; Gao and Liu 2021). Therefore, different configurations are required to be sampled when calculating the β factors for the complexes in solutions. First-principles molecular dynamics (FPMD) simulations based on DFT (density functional theory) can generate sufficiently long and stable trajectories for solutions at constant volume, temperature and pressure, which can provide enough configurations to calculate β values of complexes in solution. This approach has successfully predicted the equilibrium isotopic fractionation factors among minerals and aqueous species for many isotope systems (e.g., Blanchard et al. 2009, 2017; Kowalski and Jahn 2011; Dupuis et al. 2015; Qin et al. 2016; Wang et al. 2019b; Gao and Liu 2021).

In the present work, we apply this method to calculate the β values of aqueous Sn species and tin-bearing minerals. All tin-bearing substances are treated at the same theoretical level. We focused on [Sn²⁺$(OH)2](aq)0$⁠, [Sn²⁺$F3](aq)−$⁠, [Sn²⁺Cl₃]⁻_(aq), [Sn²⁺CO₃]⁰_(aq), and [Sn⁴⁺Cl₄(H₂O)₂]_(aq), which have been identified as the predominant stable Sn-bearing species in various geological fluids (e.g., Mark 1977; Heinrich 1990; Edwards et al. 1996; Sherman et al. 2000; Cigala et al. 2012; Schmidt 2018; Liu et al. 2020; She et al. 2020; Wang et al. 2021). As for Sn-bearing minerals, we mainly focused on cassiterite (Sn⁴⁺O₂), malayaite (CaSn⁴⁺SiO₅), and stannite (Cu₂FeSn⁴⁺S₄), which are dominant in various types of tin deposits (Eadington 1983; Eadington and Kinealy 1983; Hosking 1988; Peng and Bromley 1992; Bortnikov 2006; Aleksandrov and Troneva 2007; Yao et al. 2018). The controlling factors of tin isotope fractionation among these minerals and the aqueous Sn species have been further discussed, and the calculated equilibrium Sn isotope fractionation between aqueous Sn species and minerals are key to the utilization of Sn isotopes in geological processes.

Equilibrium Sn isotopic fractionation is determined by the vibrational frequency of tin atoms in two distinct phases (Bigeleisen and Mayer 1947; Urey 1947), and more detailed information about the terms and equations of the theory can be found in many reviews (Schauble 2004; Chacko et al. 2001; Liu et al. 2010; Dauphas and Schauble 2016; Blanchard et al. 2017). The isotope exchange reactions of the light tin isotope (¹¹⁶Sn) and the heavy tin isotopes (¹²²Sn or ¹²⁴Sn) in two tin-bearing substances (SnA and SnB) can be expressed as:

where A and B are ligands. The reduced partition function ratios (β_SnA or β_SnB) represent the fractionation factors in the two phases (Sn phase of interest and its ideal gas atom). The isotope fractionation factor (α_SnA-SnB) between the two Sn substances (SnA and SnB) can be described by:

Based on the Bigeleisen-Mayer equation, the β factor is defined as follows (Bigeleisen and Mayer 1947; Urey 1947):

and

where N is the number of vibrational harmonic modes; v_i depicts the vibrational frequency of i^th mode; h represents Planck constant; k_B depicts Boltzmann constant; and T depicts the temperature in Kelvin. The above expressions are applicable to the calculation of isotope fractionation between molecules.

For crystals, the β factors are calculated in a similar way which can be represented as:

where N_q depicts the amount of q-vectors in the Brillouin zone. In this study, the vibrational frequencies of minerals are found by computing supercells. The vibrational frequencies of q = 0 modes can yield adequately precise β values. Three translational modes with vibrational frequencies verging on zero are neglected.

Sn-bearing minerals. We carried out calculations on cassiterite (Sn⁴⁺O₂), malayaite (CaSn⁴⁺SiO₅) and stannite (Cu₂FeSn⁴⁺S₄) (Fig. 1). Based on the density functional theory (DFT), the first-principles periodic method is used to calculate the vibrational frequency which is implemented with the Vienna ab initio simulation package (VASP) (Kresse and Furthmüller 1996). The exchange-correlation functional is defined by the Perdew, Burke, and Ernzerhof (PBE) (Perdew et al. 1996) generalized gradient approximation (GGA) and projector augmented wave (PAW) method. The cut-off energy of the plane wave is set at 600 eV. On the basis of the Monkhorst-Pack scheme (Monkhorst and Pack 1976), the k-point grids in the Brillouin-zone integration are shown in Table 1. Supercells are used for the structural relaxations (Table 1), and both the cell parameters and atomic positions are optimized. The force convergence criterion is set as 10^–3 eV/Å. After relaxation, the vibration frequencies are computed by applying density functional perturbation theory (IBRION = 5).

In addition, to detect the influence of simulation size on Sn isotope fractionation, we compared the calculated β values with different size simulations. For instance, we calculated the β factors of stannite with supercells of 1 × 1 × 1 and 2 × 2 × 2, and the calculated β values are 1.00323 and 1.00326, respectively. The difference is very small, which indicates that 2 × 2 × 2 supercells are already adequately large to remove the size effect.

Sn-bearing species in aqueous fluids. For aqueous Sn (2+,4+) complexes (Table 1), we carry out the FPMD simulations to provide configurations for calculating the final β factors. This technique has been used to predict Si, Mg, and Br isotope fractionations among various minerals and aqueous species (Kowalski and Jahn 2011; Dupuis et al. 2015; Pinilla et al. 2015; Wang et al. 2019b; Gao and Liu 2021). The simulations of all Sn complexes are carried out in the NVT (canonical ensemble) thermodynamic ensemble at a fixed temperature of 300 K. The cut-off energy of the plane wave is set at 600 eV, and the k-point grids are 1 × 1 × 1. For [Sn²⁺F₃]⁻ and [Sn²⁺Cl₃]⁻ species, the Na atom is added to maintain charge balance. These cell parameters are selected to guarantee that the statistical pressure is verging on zero and the densities are nearly 1 g/cm³. The time step for all simulations is set to 1 fs. The FPMD simulations were performed for 50 ps and the energies converged to a constant value after 30 ps. The Sn-bearing solute molecules and 50 water molecules are embedded in a cubic box (Table 1). Isotope fractionation has previously been shown to be a relatively local effect, and fractionations are closely linked to the local structures of the atoms of interest (e.g., Liu and Tossell 2005; Li and Liu 2011; Fujii and Albarede 2018). It is sufficient to simulate the local structures of Sn (2+, 4+) ions in solutions using a solute molecule surrounded by 50 water molecules (e.g., Gao and Liu 2021; Wang et al. 2019b). The structure at the end of the FPMD simulation contains at least three layers of water molecules surrounding the Sn-bearing solute, and each water molecule is in a location that allows for the formation of two or more hydrogen bonds. They form many four-, five-, and six-member rings (e.g., Xantheas and Dunning 1993; Estrin 1996; Ludwig et al. 1999; Liu and Tossell 2005).

The Sn-bearing species [Sn²⁺$Cl3−$(H₂O)₅₀, Sn²⁺CO₃(H₂O)₅₀, Sn²⁺(OH)₂(H₂O)₅₀, Sn²⁺$F3−$(H₂O)₅₀, and Sn⁴⁺Cl₄(H₂O)₅₀] solution structures are sustained at 300 K for 50 ps and the energies converge to a constant value between 30 and 50 ps. Because of the configuration disorder of the aqueous solution, it is necessary to sample enough configurations to obtain the converged β value. Here, we take 10 configurations for each aqueous solution (Fig. 2). The configurations are extracted from the equilibrated trajectories (>30 ps) with equal intervals (2 ps). Then the atomic positions of the configurations are optimized with a force convergence criterion of 10^–3 eV/Å. Finally, the vibrational frequencies are calculated to predict the β values for each aqueous species.

The relaxed crystal structures of all calculated minerals are shown in Figure 1. The calculated lattice parameters of Sn-bearing minerals, bond lengths, and frequencies are listed in Tables 1 and 2. The β factors of minerals (^124/116Sn and ^122/116Sn) are calculated by Equation 5. In Table 3, we present the polynomial fit parameters for the calculated mineral β factors. The β values for cassiterite, malayaite and stannite are shown in Online Materials¹ Table S2, and their variation with temperature is shown in Figure 2a. The β factors of cassiterite and malayaite are much larger than that of stannite. The 1000ln(β) values of cassiterite (^124/116Sn) decreased from 11.128‰ at 25 °C to 2.666‰ at 350 °C, and the 1000ln(β) values of malayaite (^124/116Sn) decreased from 11.434‰ at 25 °C to 2.766‰ at 350 °C, and the 1000lnβ values of stannite (^124/116Sn) decreased from 3.225‰ at 25 °C to 0.750‰ at 350 °C (Online Materials¹ Table S2). Malayaite has the shortest Sn-O bond length and the largest 1000lnβ value. Though the tin valence states of the three minerals are the same (4+), and the coordination number of Sn in stannite (4) is smaller than those of cassiterite (6) and malayaite (6), the 1000lnβ value of stan-nite is much lower than that of the other two minerals (Online Materials¹ Table S2). This is related to the fact that the bonds formed by Sn and S are much weaker than those formed by Sn and O. Light Sn isotopes are therefore enriched in Sn-S bonds.

To evaluate the structural properties of all Sn-bearing aqueous species, we performed calculations on the partial radial distribution functions (RDFs) of the Sn-O pair of atoms. The coordination numbers (CNs) for O around Sn and the RDFs of various Sn aqueous species are shown in Figure 3. Short-range order and long-range disorder are reflected by the initial peaks and the decrease in magnitude of subsequent fluctuations in liquids (e.g., Karki 2010). For Sn²⁺(OH)_2(aq), Sn²⁺$Cl3(aq)−$⁠, Sn²⁺CO_{3 (aq)}, and Sn⁴⁺Cl_4(aq) at 300 K, the Sn-O RDF peaks occur at ~2.15, ~2.25, ~2.15, and ~1.95 Å, respective (Fig. 3). The calculated coordination numbers (CNs) are ~3.13, ~1.55, ~3.89, and ~1.07, respectively (Fig. 3 inset). For Sn$F3−$⁠, at 300 K, there is no sharp peak in the Sn-O RDF, indicating that F^– dominates the first shell of Sn$F3−$⁠, while water molecules inhabit the second shell.

Ten configurations are extracted from the equilibrated trajectories at equal intervals of each Sn-bearing species (Fig. 2). The β factors are calculated by averaging the β values of 10 configurations (Online Materials¹ Tables S1 and S2). The β factors of aqueous Sn species are in order of Sn⁴⁺Cl₄(H₂O)₂ > Sn²⁺$F3−$ > Sn²⁺(OH)₂ > Sn²⁺CO₃ > Sn²⁺$Cl3−$ (Fig. 4b; Online Materials¹ Table S2). Sn⁴⁺Cl₄(H₂O)₂ has the largest β factors (^124/116Sn: 7.627 ± 0.04‰ at 25 °C), and the β factors of Sn²⁺$Cl3−$ is the smallest (^124/116Sn: 2.639 ± 0.14‰ at 25 °C) (Online Materials¹ Table S2). In Table 3, we list the polynomial fit parameters for Sn-bearing aqueous species’ β factors.

Previous studies have shown that the nuclear volume effects (NVE), which result from the isotopic change in the nuclear size and shape, is larger for heavy elements and it is recognized as the major origin of the mass-independent isotope fractionation (Bigeleisen 1996; Schauble 2007; Abe et al. 2008, 2010; Fujii et al. 2009, 2011, 2013). For example, the NVE effect is ~1‰/amu (atomic mass unit) for Hg (Schauble 2007), Pb (Fujii et al. 2011; Yang and Liu 2015), TI (Schauble 2007; Yang and Liu 2015), and U (Bigeleisen 1996; Abe et al. 2008, 2010). Here, the method of applying the first-principles calculations obtains the mass-dependent fractionation in Sn-bearing substances resulting from differences in vibrational frequencies, without considering the effect of NVE on Sn isotope fractionation originating from differences in volume and shapes of the atomic nucleus. However, according to previous studies, the influence of NVE on Sn isotope fractionation is usually quite small (0.5‰ in ^124/116Sn for cassiterite and romarchite at 25 °C, 0.14‰ at 727 °C) (Schauble 2013), and therefore the mass-dependent fractionation at low and moderate temperatures is the main factor leading to the fractionation of Sn isotopes.

Our calculated lattice coefficients of Sn-bearing minerals are consistent with previous experimental values (Table 1) (Brockway 1934; Takenouchi 1971; Bolzan et al. 1997). The calculated Sn-O and Sn-S bond lengths of cassiterite and stannite are 2.059 and 2.566 Å, respectively, which are in good agreement with the experimentally measured values, 2.058 and 2.539 Å (Brockway 1934; Takenouchi 1971; Bolzan et al. 1997). The calculated bond lengths of Si-O, Sn-O, and Ca-O of malayaite (1.647, 2.064, and 2.448 Å) are also consistent with the experimentally measured values (1.641, 2.042, and 2.483 Å) (Table 2) (Scott 1970; Higgins and Ribbe 1977; Zhang et al. 1999; Olekseyuk et al. 2002; Fontané et al. 2012; He et al. 2021).

For aqueous Sn species, the variation of 1000ln(β) values exceeded 0.1‰ between different configurations (Fig. 2), indicating that it is essential to sample different configurations when computing β factors in solution, which is similar to the previous calculations of Br and Mg species in solution (Wang et al. 2019b; Gao and Liu 2021). In this study, 10 configurations are sampled for each aqueous Sn (2+, 4+) species to calculate the average values of β, and the values converge to constants (Fig. 2; Online Materials¹ Table S2). The 1000ln^122/116β values of Sn²⁺$Cl3−$·50H₂O using the same theoretical level with the minerals and Sn²⁺$Cl3−$·18H₂O (only one configuration) (Wang et al. 2021) calculated at B3LYP/def2-TZVP level were 2.07 and 2.18‰, respectively (Fig. 4b; Online Materials¹ Table S2). For Sn-bearing minerals, the results show that the1000ln^122/116β values of cassiterite we calculated (PBE) are consistent with the X-ray scattering (NRIXS) values (Polyakov et al. 2005) (Fig. 4a). The consistency of these data indicates that our calculation results are accurate and reliable.

Based on the calculations, the β values of Sn⁴⁺Cl₄(H₂O)₂ are larger than those of Sn²⁺ species, indicating that heavy Sn isotopes are preferentially enriched in higher Sn oxidation states, and the enrichment trend is similar to that of Br, Se, and V species in solution (Li and Liu 2011; Wu et al. 2015; Gao and Liu 2021). For the Sn²⁺ aqueous species, the bond lengths of Sn²⁺-F_(aq), Sn²⁺-O_(aq), and Sn²⁺-Cl_(aq) are 2.13, 2.14, and 2.77 Å, respectively, and the β factors decrease with the increase of bond length (Fig. 2b; Online Materials¹ Table S2), indicating that bond length is the major controlling factor for the Sn isotope fractionation of aqueous Sn²⁺ species. It conforms to the rule that shorter bonds have a higher vibrational frequency, which favors enrichment in heavier isotopes over longer and weaker bonds (Urey 1947; Schauble 2004; Young et al. 2009; Huang et al. 2013, 2014). When the valence states and bond lengths of Sn-bearing species are different, the valence state has a stronger effect on Sn isotopic fractionation than that of the bond length among these Sn complexes. For instance, at a constant temperature, the bond length of Sn⁴⁺Cl₄(H₂O)₂ (2.44 Å) is longer than those of Sn²⁺$F3−$ (2.13 Å) and Sn²⁺(OH)₂(2.14 Å), but the β factor of Sn⁴⁺Cl₄(H₂O)₂ is larger than those of Sn²⁺$F3−$ and Sn²⁺(OH)₂ (Fig. 2b; Online Materials¹ Table S2).

In previous theoretical studies of other elements, a strong correlation has been demonstrated between β values and cell parameters of minerals during equilibrium isotope fractionation, such as the valence state of an element (e.g., Schauble et al. 2001), the coordination number (e.g., Huang et al. 2013; Liu et al. 2018), and the average bond length (e.g., Huang et al. 2019; Liu et al. 2018; Méheut and Schauble 2014). For cassiterite (malayaite) and stannite with different coordination numbers and coordinating elements, heavy tin isotopes are preferentially enriched in cassiterites with high CNs (6), and light tin isotopes are enriched in stannite with low CNs (4), deviating from the general rule that minerals with low CNs are preferentially enriched in heavy isotopes. This is because the coordination element of Sn is S in stannite (Sn-S), while the coordination element of Sn is O in cassiterite and malayaite (Sn-O). The bond strength of Sn-O is stronger than that of Sn-S, which makes the Sn-O bond length of cassiterite (or malayaite) shorter than the Sn-S bond length of stannite, suggesting that the properties of the coordination element have a stronger effect on isotopic fractionation than the coordination number in these three minerals. The general principle of coordination in qualitative fractionation factors is also unsuitable for Cu- and Fe-bearing minerals. For instance, the CNs of Cu in cuprite is two, and that of Cu in carbonates and sulfates is four, but cuprite tends to be enriched in lighter Cu isotopes than carbonates and sulfates (Liu et al. 2021b); the CN of Fe in chalcopyrite is four, but it is more enriched in lighter Fe isotopes than Fe in hematite and pyrite, in which the CN is eight (Blanchard et al. 2017). Although stiff bonds tend to enrich heavy isotopes in various isotopic systems, it is unclear which of these crystal-chemical parameters is the most important factor in determining the interatomic bonds (Schauble 2004). This demonstrates the complexity of factors affecting equilibrium stable isotope fractionation (Liu et al. 2021b).

Temperature is an important factor influencing the formation of ore deposits. As a function of temperature, equilibrium stable isotope fractionation among minerals is a robust tool to reconstruct the thermal evolution process of hydrothermal deposits. The first-principles calculation provides a reliable method to predict equilibrium isotope fractionation factors (1000lnα). Recently, a few mineral-pair isotope geothermometer have been obtained by first-principles calculation, such as Ag (argentite-stephanite), Mg (garnet-pyroxene, garnet-olivine) and Ti (ilmenite-clinopyroxene) (Huang et al. 2013; Wang et al. 2020, 2022). Here, we calculated the equilibrium Sn isotope fractionation factors among the three minerals (Fig. 5) and found that the Sn isotope fractionation of cassiterite-malayaite varies little (1000lnα_{cassiterite-malayaite} ≈ 0) with temperature (Fig. 5), which indicates that it is not suitable for an isotopic geothermometer. The 1000lnα values of cassiterite-stannite and malayaite-stannite show a strong dependence on temperature and vary widely (1000lnα^124/116malayaite-stannite = 8.209‰ and 1000lnα^124/116cassiterite-stannite = 7.903‰ at 25 °C) (Fig. 5). Therefore, cassiterite-stannite and malayaite-stannite have potential as Sn isotope geothermometers. The tin isotope fractionation of cassiterite-stannite and malayaite-stannite as functions of temperature are as follows.

malayaite-stannite:

cassiterite-stannite:

The precipitation of various Sn-bearing complexes from hydrothermal solution into tin-bearing minerals is important for tin mineralization. According to the calculated results of this study (Fig. 6), the oxidation of Sn²⁺ complexes in solution into Sn⁴⁺-bearing minerals produces significant tin isotopic fractionation, (Online Materials¹ Tables S2 and S3), which indicates that the oxidation state of Sn (Sn²⁺ aqueous species → Sn⁴⁺ minerals) (except stannite) is the major controlling factor on Sn isotopic fractionation during precipitation. Therefore, the precipitation of cassiterite or malayaite from aqueous Sn²⁺ complexes (Figs. 6a and 6b) can lead to the enrichment of the lighter Sn isotope in residual hydrothermal fluids and in the later-precipitated Sn-bearing minerals. This is consistent with the Sn isotope fractionation in the typical tin deposits investigated by Yao et al. (2018), which show the precipitation of heavy-Sn-enriched isotope cassiterite results in residual dissolved Sn with a lighter isotopic composition in different types of tin ores.

There are two cases of Sn isotope fractionation in the process of precipitation of these Sn²⁺ complexes into stannite (Fig. 6c). When tin transports dominantly as [Sn²⁺(OH)₂]⁰, [Sn²⁺F₃]⁻, and [Sn²⁺CO₃]⁰ the precipitation of stannite will lead to the enrichment of heavy tin isotopes in the residual hydrothermal and later-deposited minerals. On the contrary, the precipitation of stannite from hydrothermal solution dominated by [Sn²⁺Cl₃]⁻ leads to the enrichment of light tin isotopes in the residual hydrothermal solution.

In addition, recent studies indicate that Sn can also be transported by Sn⁴⁺ aqueous species (Wang et al. 2021; Schmidt 2018). For the precipitation process of Sn⁴⁺ species, our results show that when Sn⁴⁺ precipitates as cassiterite (or malayaite), the values of 1000ln(α)[Sn(IV)Cl₄(H₂O)₂]⁰-SnO₂ (^124/116Sn) range from –3.479 to –0.834‰ (from 25 to 350 °C) (Online Materials¹ Tables S2 and S3), which indicate that direct precipitation of Sn⁴⁺ (Sn⁴⁺ aqueous species → Sn⁴⁺ mineral) from solution can also cause significant Sn isotope fractionation. When Sn⁴⁺ precipitates into stannite, it will lead to the enrichment of heavy tin isotopes in residual hydrothermal fluids and late minerals. The direct precipitation of Sn⁴⁺ species may also be a fractionation mechanism that cannot be ignored in the process of tin mineralization. The 1000ln(α) values between Sn-bearing aqueous species and minerals at different temperatures (Fig. 5; Online Materials¹ Table S3) provide a guideline to understand Sn isotope fractionation in different geological processes.

Our data on equilibrium Sn isotope fractionation factors among Sn-bearing aqueous species and minerals are important to better understand the transport, enrichment, and mineralization processes of tin when applying tin isotopic geochemistry. Cassiterite-stannite and malayaite-stannite geothermometers appear to be robust tools to reconstruct the thermal evolution process of hydrothermal Sn deposits. This study is crucial for the further application of Sn isotope in tracing various geological processes involving tin.

We thank Editor Don Baker and Associate Editor Jianwei Wang for handing the manuscript, providing insightful comments, and helping to polish the writing. Constructive comments form three anonymous reviews are appreciated. This study was assisted financially by the National Natural Science Foundation of China (92062218, 41822304).