Unraveling the Sulfate Sources of (Giant) Gypsum Crystals Using Gypsum Isotope Fractionation Factors

We combine newly determined isotope fractionation factors of gypsum precipitated in the laboratory with the isotopic compositions of natural anhydrite and gypsum to unravel the sulfate sources of the giant selenite crystals in the Naica mine (Chihuahua, Mexico). Gypsum was precipitated in the laboratory from CaSO4-NaCl-H2O solutions across a broad temperature range to establish the isotopic fractionation behavior of the sulfate molecule between the solid and dissolved phase. Oxygen isotopes show a significant fractionation dependence on temperature, with the solid phase more depleted in light isotopes with decreasing temperature. Sulfur isotopes display only a weak but similar dependence on temperature. At high salinity (4.5 M NaCl) no temperature dependence was found for the isotope composition. Based on this fractionation behavior, we attempt to elucidate the origin of the sulfate source(s) responsible for the formation of the (giant) gypsum crystals in the Naica mine. Detailed analysis of the isotopic composition of anhydrite, gypsum, and water samples strongly suggests that different types of anhydrite (of hypogenic and sedimentary origin) were dissolved to form these unique gypsum formations. The homogeneous isotopic composition of most gypsum crystals analyzed reveals an effective hydrodynamic mixing and a slow kinetics of precipitation fed by solutions of calcium sulfate from different anhydrite sources.


Introduction
Stable isotopes are valuable markers for tackling a wide variety of geological problems (e.g., Chacko et al. 2001). For example, in ore geology, the isotopic value of gangue or ore minerals (Rye 2005) is frequently used to decipher the genesis of hydrothermal deposits, and in the field of hydrogeology, the variation of isotope ratios is used like a fingerprint to identify different component sources (e.g., Yang et al. 1997;Samborska et al. 2013). On the other hand, isotope variations in speleothems are useful for paleoclimate reconstruction (e.g., McDermott 2004). Key to all these studies is a precise knowledge of the isotopic fractionation factor as a function of the environmental conditions (e.g., temperature [T], salinity, biological activity, etc.).
Despite the prominent role of gypsum (CaSO 4 2H 2 O) in the evolution of the surface of the Earth and that of Mars (e.g., Langevin et al. 2005), no precise values of the fractionation factors are available in literature. To our knowledge, two experimental determinations of the (equilibrium) isotope fractionation for sulfur (Thode and Monster 1965;Raab and Spiro 1991) and one for oxygen (Lloyd 1968), between gypsum crystals and dissolved sulfate, have been reported so far. The values obtained in these works are representative for solutions with high salinities (brines) at room temperature. However, gypsum precipitation can occur in a wide range of solution conditions (e.g., salinities and temperatures; Ossorio et al. 2014), and thus, we take into account that the isotope composition of the sulfate molecules of gypsum crystals might depend significantly on the precipitation conditions due to either kinetic and/or equilibrium fractionation.
One quite spectacular example of precipitation that occurred outside of the conditions studied by the above-mentioned authors comes from the giant gypsum crystals found in the Naica mine (Chihua-hua, Mexico), where crystals have been growing from relatively warm solutions (from 477 5 1.57C to 54.57 5 27C; Garcia-Ruiz et al. 2007;Garofalo et al. 2010;Kruger et al. 2013). The formation of these giant crystals in natural cavities located at different depths (from 2120 to 2590 m) was driven by the gypsumanhydrite solubility disequilibrium. When the temperature of the Naica range dropped below the transition temperature of gypsum-anhydrite (!587C), the latter started to dissolve, creating a solution supersaturated with respect to gypsum, inducing the precipitation of this mineral phase (Garcia-Ruiz et al. 2007). As a consequence, the crystallizing system was maintained close to equilibrium but remained slightly supersaturated with respect to gypsum due to continued dissolution of anhydrite. This resulted in the very slow growth rates of these crystals and their large size (Van Driessche et al. 2011).
Hypogenic anhydrite, related to the later stages of hydrothermal events (Stone 1959;Megaw et al. 1988), was identified as the source for the continuous supply of Ca 21 and SO 4 22 ions (Garcia Ruiz et al. 2007). Nevertheless, inside the Naica mine, a severalhundred-meters-thick lower Cretaceous sedimentary succession, mainly consisting of carbonate material, is found containing intercalations of anhydrite (Aurora Formation; Villasuso 2009). Taking into account that these anhydrite layers are located at the same levels as those of the gypsum caves, it seems reasonable to assume that these materials also acted as sources of SO 4 22 and Ca 21 for the growth of gypsum crystals.
To unravel the contribution of each anhydrite source to the genesis of the giant gypsum deposits, we designed a work strategy divided into two main tasks: (1) experimental determination of the fractionation factor of oxygen and sulfur isotopes between gypsum and dissolved sulfate as a function of the solution temperature, salinity, and saturation index and (2) isotopic and petrographic study of representative samples of gypsum and anhydrite collected from the Naica mine. This included a new set of isotopic analyses but also the reevaluation of previous isotopic data (Garcia-Ruiz et al. 2007). Both new and old isotopic data were used as fingerprints to decipher the geochemical history of the calcium sulfate minerals found at the Naica site.

Material and Methods
Synthetic Gypsum Crystals Produced in the Lab.
Synthetic isothermal gypsum crystals were grown by mixing 25 mL of equimolar solutions of CaCl 22 NaCl and Na 2 SO 4 NaCl at five different temperatures: 47, 207, 407, 607, and 807C. These crystallization solu-tions (V tot p 50 mL) were stored in sealed glass containers inside a temperature-controlled chamber (50.57C; Memmert). All experiments were run for 1 mo to assure that equilibrium was reached between the solids and the reaction solution, except for experiments at 47C. These experiments were stored up to 7 mo due to the significantly slower growth kinetics of gypsum at low temperatures (Van Driessche et al. 2010), so more time was needed to reach equilibrium of the precipitation reaction. Additionally, at 407C, experiments were stopped at 3 d and 2, 3, and 4 wk, to test the possible influence of experiment duration on the isotopic composition. In table 1 a detailed overview of the experimental conditions of each run is shown. Experiments were divided into three groups with different initial starting conditions of the reaction solution: group A contained 0.035 M CaSO 4 and 0.08 M NaCl, group B contained 0.08 M CaSO 4 and 0.08 M NaCl, and group C contained 0.08 M CaSO 4 and 4.5 M NaCl.
Gypsum crystals were separated from their reaction solution by vacuum filtration. Five milliliters of MilliQ water was added to the collected reaction solution after filtration to avoid any additional precipitation. The crystals were dried at room temperature and stored. At the end of each run, the isotopic compositions of the formed gypsum crystals and the sulfate that remained in solution were measured.

Natural Gypsum and Anhydrite Samples from the Naica Mine
Seven new isotopic analyses of anhydrite from the Aurora Formation and five of gypsum were carried out to complement the isotopic characterization done during our previous work (García-Ruiz et al. 2007). As a result, the total set of data corresponds to nine anhydrite samples from the Aurora Formation (Anhydrite I), six hypogene anhydrite samples (Anhydrite II), and 12 gypsum samples from different parts of the mine. A summary of these natural samples is shown in table 2.

Measurement of Oxygen and Sulfur Isotope Ratios
All the sample treatments were performed at the Mineralogia Aplicada i Medi Ambient research group laboratory, and treated samples were analyzed at the Serveis Cientifico-tècnics of Universitat de Barcelona following a standardized protocol for sulfur and oxygen isotopic analysis. Representative samples of 0.1 g of gypsum and anhydrite were dissolved in 50 mL of MilliQ water using stirring to accelerate the dissolution process. When the solid phase was fully dissolved, sulfate was precipitated as BaSO 4 by the addition of excess BaCl 2 2H 2 O and acidifying the sample with HCl (pH ! 2) and boiling it to prevent BaCO 3 precipitation. The sulfur isotope (d 34 S) composition of gypsum and reaction solutions was determined with an elemental analyzer (Carlo Erba 1108) coupled to an isotope-ratio mass spectrometer (IRMS; Delta C Finnigan Mat), and the oxygen isotope (d 18 O) composition was measured with a thermochemical elemental analyzer (TC/EA Thermo-Quest Finnigan) coupled to an IRMS (Delta C Finningan Mat). Notation is expressed in terms of d‰ relative to the Vienna standard mean ocean water and Vienna Canyon Diablo troilite standards. The isotope ratios were calculated using international and internal laboratory standards. For oxygen, two or three analyses were done per sample, and the error was 0.16‰. The average amount of oxygen per sample was 29%, which is the same as the barite standard. For sulfur, one analysis per sample was done, and the average amount of S was 13.2%. Sample weight ranged from 0.180 to 0.190 mg for oxygen and from 0.269 to 0.291 mg for sulfur.
To improve the precision of our measurements, instrumental mass fractionation effects were corrected using two different calibration protocols depending on the origin of the samples, i.e., natural or synthetic. (i) In the case of synthetic crystals, preliminary experiments had already provided us with an isotopic composition of the raw material as well as a notion of the fractionation values. Therefore, when measuring the isotope composition of sulfur,

Results and Discussion
Isotope Fractionation of Synthetic Gypsum Crystals as a Function of T and Salinity. Figure 1 summarizes the sulfur and oxygen isotopic compositions of gypsum precipitated from solution as a function of temperature and for different saturation indexes. As expected, due to the larger mass difference, oxygen isotopes show the highest fractionation factor, and a significant dependence on temperature is found for two experimental conditions ( fig. 1a, 1b). On the other hand, the sulfur isotope composition of precipitated gypsum shows only a weak dependence on temperature for one experimental condition ( fig. 1b).
In order to verify how long it takes to reach a stable isotope composition, we determined the fractionation factor as a function of reaction time (i.e., the time the precipitated product was in contact with the mother solution) at 47 and 407C. For the latter, we observed that after 3 d already a stable isotopic fractionation was reached, although condition A did not yet attain chemical equilibrium. Since the reaction time has no significant effect on the isotopic composition of the precipitated phase, these experiments provide us with an estimation of the "natural" scatter of the oxygen and sulfur isotopic composition at 407C for most of the starting conditions ( fig. 1d-1f), which is approximately 0.5‰. The scatter is higher, up to 1‰, for oxygen in the experiments conducted at high ionic strength ( fig. 1f). At 47C, for oxygen fractionation, a steep dependence on time was found ( fig. 1g-1i). After 1 mo, no isotopic equilibrium was reached, and at least 2 mo of reaction time was necessary to reach a stable isotope composition. After 7 mo of reaction, no significant difference was observed. Although the growth of gypsum at low temperature (!207C) is rather slow (Van Driessche et al. 2010), both conditions B and C were nearly finished precipitating out after 1 mo, but condition A was not (a significant amount of excess SO 4 22 could still be detected in solution). Noteworthy also is the fact that for sulfur, no time dependence was found. Currently, more experiments are being conducted to further study the isotope fractionation behavior at low temperature.
No significant difference was observed between the fractionation behavior for different supersaturations (saturation index p 0.2-0.8; fig. 1a, 1b), the opposite of that found for calcium carbonate, where a strong dependence on supersaturation (i.e., growth rate) is observed (e.g., Lemarchand et al. 2004;Dietzel et al. 2009;Watkins et al. 2013). On the other hand, the lack of temperature dependence for high-salinity solutions ( fig. 1c) indicates that solution composition can influence the fractionation behavior. In this case, the salinity effect (e.g., Horita 2005; Oi and Morimoto 2013) is most likely responsible for the absence of a clear temperature influence at high salinity. In any case, our results demonstrate a significant dependence on temperature, especially for oxygen isotopes, but more specific experiments are needed to fully elucidate the influence of other factors, such as salinity, on the fractionation behavior of the sulfate molecule during gypsum precipitation. A detailed understanding could be helpful for exploring new proxies for paleoclimatology reconstruction of gypsum growth environments.
For our high-salinity experiments, we found an average value of 3.3‰ 5 0.3‰ for the oxygen fractionation factor. This is in fair agreement with the value of 3.6‰ reported by Lloyd (1968), determined from gypsum and associated brines of evaporating pans. Hence, the value of 3.3‰ 5 0.3‰ is valid for evaporitic environments with high salinity. But, our results also indicate that, depending on the precipitation conditions (temperature, ionic strength) of gypsum, different values of the oxygen fractionation factor should be used. In the case of sulfur, even though Feely and Kulp (1957) reported the absence of fractionation during gypsum precipitation, Thode and Monster (1965) found a fractionation factor of 1.6‰. This value was obtained from gypsum precipitating from a slowly evaporating saturated solution of CaSO 4 under reduced pressure (16-20 mm Hg) at room temperature (T not specified). Raab and Spiro (1991) found a similar value, 1.59‰, during the stepwise evaporation of seawater at 23.57C, up to the middle of the halite precipitation field. In this work we found a constant, slightly higher, value of ∼2.0‰, except for lower-salinity experiments conducted at higher temperatures (1607C), with a value of ∼1.7‰.
Anhydrite in the Naica Mine. The Naica mine is an Ag-Pb-Zn-(Cu) carbonate-hosted manto massive sulfide deposit located in the Sierra Madre Occidental, on the northern flank of a NW-SE 12-kmlong and 7-km-wide dome structure. The structure is mainly composed of Cretaceous carbonate rocks intruded by a set of sills and dikes of Oligocene age probably related to an igneous body located at a depth of 11100 m (Alva- Valdivia et al. 2003). For a more extended review of the geology of the Naica area, the interested reader is referred to the works of Stone (1959), Megaw et al. (1988), and Villasuso (2009).
Inside the mine, the host carbonates belong to the Aurora Formation (low Cretaceous). They are located at depths between 2290 and 2520 m and contain some anhydrite beds and pods (hereafter referred to as anhydrite I; fig. 2a). Using optical microscopy, two types of anhydrite can be distinguished within the Aurora Formation: (1) anhydrite made by euhedral prismatic crystals (anhydrite Ia; fig. 2d) and (2) diagenetic anhydrite that occurs as lenticular gypsum pseudomorphs (anhydrite Ib; fig. 2c). Crystal sizes of both types of anhydrites are between tens to hundreds of micrometers. In addition, a third type of anhydrite corresponding to latestage hypogenic anhydrite is also abundantly present in Naica (anhydrite II; fig. 2b). These hypogenic crystals are characterized by their euhedral prismatic shape and light blue color, sometimes exposed as beautiful fans reaching up to 20 cm in size. Considering that some of the Aurora Formation carbonates underwent dolomitization, it is likely that anhydrite Ia, enclosed within the Aurora Formation, could be of the same origin as anhydrite II, both related to the ore-mineralization episode (Jones and Xiao 2005 and references therein).
In figure 2e, 2f, petrographic evidence of isovolumetric replacement of anhydrite by gypsum is shown in samples from the Aurora Formation. This replacement process should create an excess amount of Ca 21 and SO 4 22 ions in the aqueous solution because the molar volume of gypsum exceeds that of anhydrite by ∼40% (77,440 and 46,103 cm 3 /mol, respectively), and thus, locally, an increase of supersaturation with respect to gypsum is created. If this mechanism occurred at a large scale, a significant increase in Ca 21 and SO 4 22 concentration should be expected, fostering crystal nucleation. Gypsum to anhydrite transformation, and vice versa, does not produce significant isotope fractionation of the sulfur and oxygen atoms of the sulfate molecule (Worden et al. 1997). Therefore, diagenetic anhydrites should preserve the original isotope composition of the primary sedimentary gypsum (e.g., Pierre 1988;Utrilla et al. 1992;Clark and Fritz 1997). Within the d 18 O values of anhydrite I, two groups (A and B, fig. 3) can be distinguished. According to the estimated uncertainty of the seawater isotopic variation Claypool curve (Paytan et al. 2004;Bottrell and Newton 2006), the   fig. 3). Hence, isotope compositions indicate, as was already inferred from the petrography, that anhydrite I contains sulfates from different sources, most likely one being sedimentary and the other hydrothermal.
If we consider the isotopic composition of hypogenic anhydrite (type II) and those of group A from the Aurora Formation, it is interesting to note the narrow d 18 O variation (3.0‰) and the wide d 34 S range (5‰). This pattern is also observed for anhydrites from other hydrothermal deposits. This is interpreted in terms of the disproportionation of SO 2 to SO 4 22 and H 2 S and the subsequent precipitation of sulfate as anhydrite (Seal et al. 2000) in a magmatic-hydrothermal environment. In El Teniente, a porphyry copper deposit in Chile, the variation of d 18 O of anhydrite is 2‰, whereas the d 34 S range is 8‰ (Kusakabe et al. 1984).
A total of 13 gypsum samples (see table 2 . 3). The only apparent difference between these samples and the other two samples (N10-2 and N10-21) from the CS is their size, the latter being much larger (see fig. 4a).
The isotopic composition of dissolved SO 4 22 of present-day groundwater has also been plotted in figure 3 (d 34 S p 14.6‰ and d 18 O p 14.4‰; Garcia-Ruiz et al. 2007). In previous works, the composition of d 18 O of H 2 O was determined and falls in the range of 27.5‰ and 27.8‰ (Dames and More 1977;Garcia-Ruiz et al. 2007).
In the present-day cave system, the microbial biomass in the thermal water is very low, and attempts to amplify diagnostic functional genes for sulfate reduction were unsuccessful, suggesting that this activity, if present, is not important in the aquifer (Ragon et al. 2013). There are no obvious indications to think this was different in past times, and thus, we assume that no d 18 O exchange between the water and the sulfate occurred due to biological activity.
Unraveling the Sulfate Sources of Gypsum at Naica.
The isotopic composition of dissolved sulfate in present-day Naica waters ( fig. 3) falls somewhere in between the isotopic composition of anhydrite I and II. Thus, taking into account that (1) anhydrite dissolution does not produce isotopic fractionation (Seal et al. 2000) and (2) this mineral is the main sulfate phase present in Naica (Stone et al. 1959), we can safely assume that the sulfate isotopic compo- sition of the actual Naica waters is the result of the dissolution of anhydrite IB (although small contributions of IA cannot be ruled out) and II. Taking into account that the isotopic composition of the Aurora Formation anhydrites cannot be correlated with their location and that both types of anhydrites act as a source for the gypsum crystals, it seems reasonable to assume that a homogeneously mixed solution prevailed inside the cave system.
If we apply the fractionation factor determined from our low-salinity experiments at 507C to the composition of dissolved sulfates of actual Naica waters, we obtain a theoretical value for gypsum (dashed arrow, fig. 3) that matches very well with the isotopic values obtained from present-day gypsum crystals (table 2). The oxygen composition of this sample, within the margin of error, overlaps with the values obtained for the giant gypsum crystals. Therefore, we propose that the oxygen isotopic values of the dissolved sulfate in Naica during past gypsum precipitation were similar to what they are today. This, in turn, strongly suggests that the overall water composition (i.e., salt content and pH) and temperature of Naica have been very stable during gypsum crystallization. This was also indirectly inferred from our previous results on their formation mechanism (Garcia-Ruiz et al. 2007) and growth history (Van Driessche et al. 2011). Now that we have a detailed understanding of the isotopic composition of the Naica system with respect to calcium sulfate, and disposing of a precise value of the fractionation factor as a function of temperature, we should be able to isotopically differentiate between crystals grown at different depths inside the Naica mine. Previous to any mining activities, a significant difference in average growth temperatures existed between caves located closer to the surface (e.g., CS p 477 5 1.57C) and those found at greater depth (e.g., CC p 54.57 5 27C; Kruger et al. 2013). The temperature difference between both cavities is small but in principle significant enough to differentiate isotopically (for oxygen), because it falls within the range where the fractionation factor is most dependent on temperature. Thus, a shift toward higher isotopic values for lower temperatures is expected; i.e., for the crystals grown in the CS, we should see higher d 18 O. But, if we look at figure 3, no significant and consistent difference in isotopic composition of crystals formed in different caves can be detected.
As  . 4b), confirm the T h values obtained (477 5 1.57C) from previously analyzed crystals from the CS (Kruger et al. 2013). Thus, the difference in isotopic composition is not induced by a difference in growth temperature compared to the other samples from CS. So the only plausible option left to explain the difference in isotopic composition of the gypsum crystals is to contemplate that they were precipitated from solutions containing sulfate from a different source(s) (i.e., the anhydrites from which SO 4 22 ions were released had a different isotopic composition).
When applying the inverse of the fractionation factors determined in this work (i.e., from crystal to solution composition) to the values of both outliers from the CS, we find that the calculated isotopic composition (dash-dotted arrow, fig. 3) of the sulfate source fits very well with that of group B of anhydrite type I (Aurora Formation). Thus, taking into account the variability of the isotopic composition of the anhydrites and the rather high homogeneity in isotopic composition of most gypsum crystals formed at Naica, it seems reasonable to assume that the CS experienced different hydrological regimes; i.e., waters slightly supersaturated in CaSO 4 accessed the cavity through different pathways (i.e., waters got saturated in CaSO 4 from dissolving different anhydrite sources). Thus, the small crystals from the CS seemed to have formed when the whole cave complex was experiencing a common hydrological system, and the large crystals from CS would have formed in a more localized regime, with, most likely, a larger influence of the anhydrites from the Aurora Formation. Precipitation in a closed system containing only sulfate from anhydrite IA could also explain the isotope composition of the large crystals of the CS. However, this scenario seems unlikely given the large total volume of gypsum crystals. A more plausible mechanism would be a change in the hydrological regime due, perhaps, to seismic activity, which abruptly altered the flow path, shifting the principal anhydrite source to IB.
Also noteworthy is the fact that if we apply our experimentally determined fractionation factor to the isotopic values of the anhydrite sample from the same borehole as the gypsum sample (open square and open triangle, fig. 3), we obtain a theoretical value for gypsum far from the experimental value obtained for the gypsum from that location. This indicates that the sulfate source, or at least part, was located away from the crystallization location. This observation also indicates that hydrological regimes inside the mine play an important role in the formation of these giant gypsum crystals.

Conclusions
From our laboratory experiments we can infer that fractionation of isotopes during gypsum precipitation is influenced by both the solution temperature and the initial solution speciation: at low salinity, a clear dependence on temperature is detected, while at high salinity, similar to evaporitic environments, no dependence on temperature is found, and the obtained value is in fair agreement with the one used in literature. Gypsum formed at low temperature (47C) needs more time than other experiments in this study to reach a stable isotopic composition, but, even so, the time (!2 mo) is very fast compared to the geological timescale. These observations imply that the use of the previously reported value for oxygen fractionation is justified only for highly saline solutions; when gypsum precipitation is studied in low-salinity conditions, significantly different values should be used, depending on the reaction temperature. Also for sulfur fractionation, we found a significantly different value compared to those previously reported. In addition, specific experiments are needed to fully comprehend the influence of salinity on the fractionation behavior of the sulfate molecule during gypsum precipitation.
The isotopic study of gypsum and anhydrite samples collected from the Naica mine revealed that (i) several sulfate sources from different types of anhydrite contributed to the growth of gypsum crystals in the Naica mine; (ii) these sources, or at least part of them, were located away from the crystallization site; (iii) hydrological regimes inside the mine played an important role in the formation of these giant gypsum crystals; and (iv) Naica's groundwater had a stable composition during gypsum growth, confirming previously postulated hypotheses (Garcia-Ruiz et al. 2007;Van Driessche et al. 2011).

A C K N O W L E D G M E N T S
This work has been carried out within the framework of project CGL2010-16882 of the Spanish Ministerio de Economía y Competitividad and cofunded with Fondo Europeo de Desarrollo Regional. A. E. S. Van Driessche is grateful for FWO grant 1523115N. We acknowledge C. Ayora (Consejo Superior de Investigaciones Científicas) and C. Johnson (US Geological Survey) for useful insights on the manuscript and S. Bottrell and an anonymous reviewer for their constructive review. M. Ossorio acknowledges the Junta para la Ampliación de Estudios-Predoc fellowship.