Equilibrium and kinetic Si isotope fractionation factors and their implications for Si isotope distributions in the Earth’s surface environments

Several important equilibrium Si isotope fractionation factors among minerals, organic molecules and the H4SiO4 solution are complemented to facilitate the explanation of the distributions of Si isotopes in Earth’s surface environments. The results reveal that, in comparison to aqueous H4SiO4, heavy Si isotopes will be significantly enriched in secondary silicate minerals. On the contrary, quadra-coordinated organosilicon complexes are enriched in light silicon isotope relative to the solution. The extent of 28Si-enrichment in hyper-coordinated organosilicon complexes was found to be the largest. In addition, the large kinetic isotope effect associated with the polymerization of monosilicic acid and dimer was calculated, and the results support the previous statement that highly 28Si-enrichment in the formation of amorphous quartz precursor contributes to the discrepancy between theoretical calculations and field observations. With the equilibrium Si isotope fractionation factors provided here, Si isotope distributions in many of Earth’s surface systems can be explained. For example, the change of bulk soil δ30Si can be predicted as a concave pattern with respect to the weathering degree, with the minimum value where allophane completely dissolves and the total amount of sesqui-oxides and poorly crystalline minerals reaches their maximum. When, under equilibrium conditions, the well-crystallized clays start to precipitate from the pore solutions, the bulk soil δ30Si will increase again and reach a constant value. Similarly, the precipitation of crystalline smectite and the dissolution of poorly crystalline kaolinite may explain the δ30Si variations in the ground water profile. The equilibrium Si isotope fractionations among the quadra-coordinated organosilicon complexes and the H4SiO4 solution may also shed light on the Si isotope distributions in the Si-accumulating plants.


Introduction
The silicon element, only second in abundance to oxygen in Earth's crust, is one of the most dominant constituents in silicates. Since Si is a typical non-redox element and often exists in a stable tetrahedron, the variations of its isotopic composition in Earth's various surface systems are small, no more than 11.8 per mil (Opfergelt and Delmelle 2012). In comparison, there are no or only slight Si isotope fractionations occurring in the formations of the primary silicate minerals (Savage et al. 2011). With the rapid progress of high-resolution analytical instruments including ICP-MS and SIMS, slight variations of Si isotopic compositions can be detected. Therefore, the Si isotope compositions in many of Earth's surface systems, including soils, rivers, oceans, silicon-accumulating algae and high plants, have all been investigated (Basile-Doelsch et al. 2005;Ding et al. 2008Ding et al. , 2009Georg et al. 2007Georg et al. , 2009Hughes et al. 2013;Ziegler et al. 2005). One can refer to the comprehensive reviews of Si isotope geochemistry in Opfergelt and Delmelle (2012). In brief summary, river water and seawater have positive d 30 Si signals, but ground water has negative ones. In addition, d 30 Si values of the ground water decrease with increasing water-rocks interaction and with increasing weathering degree. However, a concave pattern was identified in the Si-accumulating plants including rice, bamboo, and bananas: the d 30 Si values decrease from root to stem and increase from stem to shoots (Ding et al. 2008(Ding et al. , 2009Georg et al. 2009;Ziegler et al. 2005).
There are disagreements on the explanation for the Si isotope signals in Earth's surface systems. For example, it has been argued that the dissolution of primary silicate minerals will preferentially release light Si isotopes to cause the solution to have negative d 30 Si values. This explanation is untenable on a large time scale (He and Liu, 2015). Another point of disagreement is the employment of the Rayleigh model for the interpretation of experimental Si isotope data, especially for complicated and open systems. The deduced Si isotope equilibrium fractionation factors by such models at ambient temperatures between the clay and the solution D 30 Si clay-solution = -1.5 % (Georg et al. 2007) and -2.05 % (Hughes et al. 2013) obviously deviate from the theoretical predictions, which reveals that stiffer chemical bonds will enrich heavier isotopes, i.e., the precipitated minerals will preferentially incorporate heavy isotopes relative to the aqueous H 4 SiO 4 due to their shorter Si-O bonds.
Moreover, Ziegler et al. (2005) suggested that it was the continuous leaching of the 30 Si-entriched pore water during the dissolution-reprecipitation processes that resulted in the gradual decrease of soil bulk d 30 Si with the weathering degree. Considering that this part of the pore water likely discharges into the groundwater, if the leached pore water is always enriched in 30 Si, the sole consequence will be extreme 30 Si-enrichment of the ground water relative to the pore water, which is not in consistent with our current observations. So if in some dissolution-reprecipitation processes the leached pore water converts into 28 Si-enrichment from 30 Si-enrichment, this contradiction could be figured out.
To clarify those issues of contention, the equilibrium Si isotope fractionation factors for the silicate minerals, orthosilicic acid solution and organosilicon complexes were calculated using a newly proposed cluster-modelbased quantum chemistry method (Liu 2013). Albite was chosen as a candidate for its stability in diagenesis. Kaolinite was selected for its ubiquity in weathering environments as a mature weathering product. A few 4-, 5-, or 6-coordinated organosilicon complexes related to the metabolic process of Si in plants were also calculated to explore the Si isotope fractionations in plants. Also, the polymerization of monosilicic acid and dimer when approaching the precipitation of amorphous silica was designed to calculate the kinetic isotope effect.
2 Theory and calculation

Bigeleisen and Mayer equation
The works by Bigeleisen and Mayer (1947) and Urey (1947) about the equilibrium constant calculation for isotopic exchange reactions have established the base of theoretical and computational geochemistry. Specifically, for an ideal single-atom-substitution system, where the compound is exchanging isotopes with the atom, the reaction equation is shown as: Element X has two isotopes, with the primed one denoting the light isotope. The equilibrium constant K can be depicted as the quotient of the partition functions of the products and reactants as follows: The partition function comprises of the transitional partition function Q trans , rotational partition function Q rot , vibrational partition function Q vib and electronic partition function Q elec . For light elements, the electronic energy of a compound after isotope substitution does not change. So the electronic partition function can be canceled out. The theorem of Teller and Redlich (1935) is often employed to simplify the calculation of the equilibrium constant, K. K can be expressed as, where u = hcx/kT, with h denoting the Plank constant; r, the symmetry factor; k, the Boltzmann constant; T, absolute temperature in K; c, velocity of light; and x, the harmonic frequencies in cm -1 . The Reduced Partition Function Ratio (RPFR) is defined as: Herein, K is equal to the isotopic fractionation factor a. Given the RPFRs of a pair of minerals in equilibrium, a can be derived from the ratio of their RPFRs. Bigeleisen and Wolfsberg (1958) introduced the transition state theory to tackle the kinetic isotope effect in chemical reactions. It is assumed that the isotope exchange equilibrium between the reactants and the transition state complex is reached and that the transition state complexes with different isotopes have the identical probability of forming products. Therefore, the isotope fractionation between the products and reactants is determined by the isotope exchange equilibrium between the reactants and the transition state complex. The Kinetic Isotope effect (KIE) is defined as the ratio of rate constants for light and heavy isotopes in formation of the transition state complex:

Kinetic isotope effect
where k L and k H are the rate constants for light and heavy isotopes, respectively, and RPFR react and RPFR ts are the reduced partition function ratios for the reactant and transition state complex. In the calculation for RPFR ts , one of the vibrational degrees of freedom along the path of decomposition is missing from the vibrational partition function. Finally, RPFR ts takes the form: where m L = is the imaginary frequency along the path of decomposition.

Geometry optimization
The Volume Variable Cluster Model (VVCM) method (Liu 2013) was employed to optimize the mineral structures. One can refer to the detailed procedure in He and Liu (2015). Herein, to express it concisely, the cluster is built from three dimensional crystal structures established by X-ray diffraction or neutron diffraction data. Virtual point charges were added to the outermost darling bonds to keep the neutrality of the whole cluster. By adjusting the distance between point charges and the outermost atoms, different stable mineral structures were sought out, from which the calculated structures with the least deviation from the experimental structures and the lowest electronic energy were selected for the calculation of the RPFR values. All geometry optimization was implemented with the GAUSSIAN09 software package (Frisch et al. 2010). The hybrid DFT method at B3LYP/6-311G(2df) level was employed to optimize the 3-dimensional structures of the silicate minerals (Becke 1993). As for the hydrogen-containing minerals and orthosilicic acid solution, an extra set of p functions was added to this level.

Optimized mineral structures
The optimized mineral structures tabulated in Table 1 and the Supplementary data are consistent with experimental data. The mismatches are in the order of per mil.

D 30 Si mineral-mineral
We used weighted-RPFR values for kaolinite and albite because they have more than one Si site in their structures. Apparently, quartz is enriched in the heavy Si isotope relative to kaolinite and albite ( Fig. 1), which is qualitatively consistent with results of Méheut et al. (2009) and Méheut and Schauble (2014). The enrichment order quartz [ albite [ kaolinite is consistent with the order of the gradually prolonged average Si-O bond length in these minerals. Although in good agreement with the previous results between quartz and albite, there are large discrepancies in the cases of ''quartz vs. kaolinite'' and ''albite vs. kaolinite''; much smaller Si isotope fractionations are predicted here (Fig. 1).

D 30 Si solid-solution
Based on the RPFR value of the H 4 SiO 4 solution calculated by the identical level (i.e., 1.0718, from He and Liu (2015)), kaolinite tends to be enriched in 30 Si if it is in equilibrium with the solution (Fig. 2). However, the field observation results diverge. The smectite result (Georg et al. 2009) is smaller than our kaolinite result, which is reasonable as the average Si-O bond length in smectite is longer. However, the field observation results of kaolinite  Rastsvetaeva et al. 2009). Initial structures of the kaolinite are extracted from the neutron powder diffraction data (Bish 1993 (Georg et al. 2007) and 2:1 clay (Hughes et al. 2013) are in completely opposite directions. We will address these disagreements in Discussions. The temperature dependencies of the equilibrium Si isotope fractionations between the minerals and the H 4 SiO 4 solution can be found in Table 2. The results of Méheut et al. (2007) and Méheut and Schauble (2014) are also listed for comparison.

D 30 Si organosilicon complexes-H4SiO4
We arbitrarily selected propanetriol to build the four-coordinated organosilicon complexes. In addition, the structures of the hyper-coordinated organosilicon complexes from Kubicki and Heaney (2003) were rebuilt. Figure 3 shows the optimized geometries of the organosilicon complexes we studied. Extremely large equilibrium Si isotope fractionations were found between the hyper-coordinated organosilicon complexes and the aqueous H 4 SiO 4 at 298.15 K, such as -9.1 % for the [5] Si-disorbitol complex and -19.2 % for the [6] Si-trisorbitol complex, whilst there are only moderate fractionations for the four-coordinated complexes, such as -0.36 % for the [4] Si-propanetriol complex and -0.92 % for the [4] Sidipropanetriol complex.

KIE in the polymerization of monosilicic acid and dimer
He and Liu (2015) has calculated the KIE in the dimerization of monosilicic acid, 1.005. Recently, McIntosh (2013) indicated that there were two routes for monosilicic acid and dimer in basic solution to form trimer. Given the same mechanism, there must also be two routes for monosilicic acid and dimer to form trimer in neutral solution (Fig. 4). We have successfully searched out the transition states by synchronization transition quasi-Newton searching methods (i.e. the QST3 method in Gaus-sian09). The calculated activation energies and KIEs are shown in Table 3. It should be noted that route 1 discriminates 28 Si in the formation of the transition state, with a KIE value, 0.999548, whilst route 2 prefers 28 Si in the formation of transition state, with a KIE value, 1.010618. Given the same pre-exponential factor, the slightly low activation energy of route 2 indicates that, in the incorporation of monosilicic acid into the crystal lattice, the monosilicic acid molecule is selected to adjust geometric structure and  to lose a water molecule. At ambient temperature, the reaction rate of route 2 is 2.24 times of that of route 1. The weighted average KIE value for the formation of amorphous silica is 1.0072. Although the experimental data deviate from the numerical calculations, they both show the positive temperature dependence (Fig. 5).

Plane-wave DFT method vs. VVCM method
There are significant discrepancies between our Volume Variable Cluster Model (VVCM) calculations and some of the pseudo-potential based plane-wave DFT results (Méheut et al. 2007;Méheut and Schauble 2014) (Fig. 1).
Here we compare the results of the geometry optimization and frequency of these two methods.
In the case of geometry optimization, pseudopotential plane-wave DFT methods can effectively reproduce cell parameters of silicates despite the large elongation of the average Si-O bond; for example, 1.2 % error for quartz and 1.4 % error for phyllosilicates (i.e., clay minerals) (Balan et al. 2001), due to the shortcoming of the PBE functional employed in the calculations (Méheut et al. 2007;Méheut and Schauble, 2014). In comparison, the largest average Si-O bond length error in our VVCM calculations is only 0.43 % in all studied minerals (Table 1), showing the power of a large full-electron wavefunction (i.e., B3LYP/6-311G(2df,p)) used in the geometry optimization. In fact, the isotope effect highly depends on the local kinetic energy differences surrounding the atom of interest before/after isotope substitution; so to effectively reproduce the local experimental structure surrounding the atom of interest it is very important to calculate the isotope fractionation. In these cases, large deviations of the calculated structures of a pair of minerals from the experimental ones, especially disproportionate deviations, may exert a considerable influence on the Si equilibrium isotope fractionation calculations (cf. quartzkaolinite, albite-kaolinite, Méheut et al. 2007, Méheut andSchauble 2014). However, if the deviations from the experimental structures are by the identical extent, the influence of the structural uncertainty on the Si isotope equilibrium fractionation may be canceled as a systematic error (cf. quartz-albite, Méheut and Schauble 2014).
The essence of structural uncertainty is the uncertainty of harmonic frequency, which is an important parameter in calculating RPFR values. Méheut et al. (2009) pointed out that about 5 % systematic errors in frequencies would bring about 10 % relative errors in the calculation of 1000lna at low temperatures. Typically, the calculated harmonic frequencies are overestimated and a scaling factor smaller than unity is needed to calibrate the calculated frequencies. For a pair of minerals, whose structures deviate disproportionately from the experimental ones, it is difficult to determine a sole scaling factor. Therefore, in the calculations of the Si isotope fractionation factors, it seems safe to keep the calculated structures within the uncertainty of experimental data. Based on the discussion above, our results are reliable.

Si isotope fractionation between clay and solution
We provide equilibrium Si isotope fractionations between kaolinite and the solution at different temperatures (Fig. 2). Our calculation results are significantly larger than the field observations of Georg et al. (2007), (2009) and Hughes et al. (2013). Because there was a much shorter average Si-O bond length (1.615-1.620Å ) (e.g., Bish and Von Dreele 1989) compared to that of smectite (1.635-1.645 Å ) (e.g., Gournis et al. 2008) for the smectite case (Georg et al.  Georg et al. (2007) and Hughes et al. (2013) suggested that aqueous H 4 SiO 4 enriched heavier Si isotopes than the clay minerals (Fig. 2), which were in obvious disagreement with the theoretical predictions. According to either the experimental data (Rastsvetaeva et al. 2009) or our quantum chemistry calculations, the average Si-O bond length in aqueous H 4 SiO 4 is longer than those in clay minerals (e.g., 1.639 vs. 1.62-1.63Å ), meaning aqueous H 4 SiO 4 cannot enrich heavier Si isotopes than clay minerals if under an equilibrium condition. Those field observations actually reflected the ''net'' or ''apparent'' Si isotope fractionation between the river water and the clay minerals. However, it is not an equilibrium Si isotope fractionation factor for a single isotope exchange reaction. The clay mineral precipitating processes are very complicated, including the non-equilibrium processes and the formation of other Si-bearing compounds (e.g., organosilicon complexes). These affecting factors will lead to a large negative ''apparent'' isotope fractionation factor, as observed in field. With the equilibrium Si isotope fractionation factors provided here, we could further explore the molecular-level clay mineral precipitation mechanisms by distinguishing the equilibrium isotope fractionations from the other processes.

Si isotope distributions in soil profile
With the precise equilibrium Si isotope fractionation factors, distributions of the Si isotopes in the Earth's surface reservoirs can be constrained. For example, Ziegler et al. (2005) found that in the chronosequence of Hawaii soil, soil bulk d 30 Si decreased to a more negative value from young to old soil in the lower soil layer (30-100 cm). Based on our theoretical calculations, we can provide reasonable explanations for their observations and predict a secular variation of the Si isotope in the soil profile. As a response to the precipitation of amorphous allophane, which is formed following the dissolution of primary silicate minerals and is enriched with light Si isotopes, the bulk soil d 30 Si will decrease. Since allophane is metastable, it will dissolve and transform into poorly crystalline minerals-including clays-in which the light Si isotope will continuously be enriched in products and the bulk d 30 Si in the soil profile will gradually decrease. Meanwhile, Fe and Mn form sesqui-oxides and accumulate in old soils. Enhanced 28 Si adsorption on Fe-oxides drives the silicon isotope of the bulk soil to a more negative value (Delstanche et al. 2009). We can predict that bulk soil d 30 Si in relation to the weathering degree has a concave pattern, with a minimum where allophane completely dissolves and the total amounts of sesqui-oxides and poorly crystalline minerals reaches their maximum.
When well crystalline clays precipitate from the pore solution under the isotope exchange equilibrium, bulk soil d 30 Si increases and begins to level off. To substantiate our prediction we need more studies relevant to the crystallinity and Si isotope composition of clays in old soils. In conclusion, the decrease of d 30 Si during kinetic precipitation will reasonably explain the stepwise decrease of d 30 Si in the soil profile.
Similarly, the variation of ground water d 30 Si with gradually enhanced water-rocks interaction can also be explained. Ground waters can effectively escape the perturbation of plants and will show distinct d 30 Si variations during diagenesis. For example, Georg et al. (2009) determined the d 30 Si of the ground water in the Navajo Sandstone aquifer and attributed the gradual decrease of the d 30 Si along the flow path to the dissolution of kaolinite and amorphous silica, evidenced by the saturation index calculation. It should be noted that smectite also had precipitated from the ground water, which was shown by the saturation index calculation and electron microscopic characterization (Zhu et al. 2006). Smectite was found to be enriched in heavier silicon isotopes relative to the ground water, which matched our theoretical predictions. Therefore, the dissolution of poorly crystalline kaolinite and the precipitation of well crystalline smectite may cocontrol the observed d 30 Si trend in ground water.

Si isotope fractionations in plants
In the metabolic process of plants, silicon (in the form of H 4 SiO 4 ) is taken up by the root and accumulated in the xylem sap, where, even if the concentration of silicon is well above the saturation of the amorphous silica, polymerization does not occur. This point was confirmed by NMR experiments; the 29 Si NMR peak at approximately -71 ppm of xylem sap extraction was assigned to the 4-coordinated Si either as a monomer or as an organosilicon complex (Mitani et al. 2005;Sahai 2004). Polyol carbohydrates were suggested to be able to stabilize the Sibearing monomer species to a considerably elevated concentration in the solution (Kinrade et al. 1999(Kinrade et al. , 2001. Moreover, the coordination number of Si in the organosilicon complexes increases with the pH value. The five-membered ring [5] Si-disorbitol complex and the [6] Sitrisorbitol complex had successfully reproduced the peaks near -102 and -144 ppm on the 29 Si NMR spectra and the heptet in the 1 H-29 Si coupled spectra in solutions with high pH values (Kinrade et al. 1999(Kinrade et al. , 2001Kubicki and Heaney 2003;Sahai 2004). Comparison between experimental data and calculated Si kinetic isotope fractionation during amorphous silica precipitation. The red line is the result of dimerization (He and Liu 2015); the blue line is the result of the polymerization of monosilicic acid and dimer; the blue squares are the experimental data of Geilert et al. (2014) and the black diamonds are the experimental data of Roerdink et al. (2015) A unique d 30 Si variation was observed in different Siaccumulating plant organs including rice and bamboo: a concave d 30 Si pattern for the sequence root ? stem ? branch ? leaf with relative positive d 30 Si in the shoots and the minimum of d 30 Si in the stem (Ding et al. 2008(Ding et al. , 2009). Based on the Rayleigh fractionation model, these unique d 30 Si distributions observed by Ding and his co-authors were attributed to the different ratios of the biogenic 28 Sienriched silica gel to remnant the isotopically heavy H 4 SiO 4 solution in different organs. However, rice and bamboo are not strictly closed systems for Si. In other words, Si can be excreted out of the rice and the bamboo through various mechanisms including guttation (Yamaji et al. 2008;Tian 2008).
Recently, a new transporter Lsi6, which can unload Si from the xylem sap and facilitates the subsequent intracellular transportation to silicon ''pools'', has been identified in the shoots (Yamaji et al. 2008). The molecular-level details of how Lsi6 transports the Si is still not clear. However, evidence shows that Lsi6 only transports H 4 SiO 4, instead of the organosilicon complexes in the pool of xylem sap. With the help of the equilibrium Si isotope fractionation between H 4 SiO 4(aq) and the organosilicon complexes provided here, we can discuss the Si isotope fractionations associated with such a process. Because large Si isotope fractionations existed between H 4 SiO 4(aq) and the 5-or 6-coordinated organosilicon complexes, the pH value of the xylem-sap pool will be important for the Si isotope distributions in the plants. A higher pH value of the xylem-sap pool will let the transported Si become isotopically lighter.

Conclusion
Several important equilibrium Si isotope fractionation factors among minerals, organic molecules and the H 4 SiO 4 solution are complemented to facilitate the explanation of the distributions of Si isotope in Earth's surface environments. Similar to quartz, precipitated minerals are enriched with heavier Si isotopes than other Si-bearing species or compounds if they are at equilibrium. The order of Si heavy isotope enrichment is found to be albite [ kaolinite [ aqueous H 4 SiO 4 [ organosilicon complexes. Discrepancies related to field observations or previous theoretical calculations are carefully discussed. That the large kinetic isotope effect associated with the formation of amorphous quartz precursor results in the light Si isotope enrichment of field observed quartz is further confirmed.
Based on equilibrium Si isotope fractionation factors provided here, the distribution of Si isotopes in a few Earth's surface reservoirs can be constrained. For example, d 30 Si variations in the soil can be predicted as a concave pattern in relation to weathering degree. Similarly, the precipitation of crystalline clay minerals and the dissolution of poorly crystalline kaolinite may explain the d 30 Si variation in the ground water profile. Finally, equilibrium Si isotope fractionations between H 4 SiO 4(aq) and the 4-, 5or 6-coordinated organosilicon complexes can facilitate further discussion of the Si isotope distribution in plants.