The Absolute Age and Origin of the Giant Gypsum Geode of Pulp í (Almer í a, SE Spain)

: Subaqueous gypsum (CaSO 4 · 2H 2 O) crystals are relatively common in epithermal systems where sulﬁde ore deposits are present. The Giant Geode of Pulp í (Almer í a, SE Spain) hosts some of the largest (up to 2 m in length) subaqueous gypsum crystals discovered to date. Here, we present the ﬁrst U-series ages of its crystals and reconstruct the oxygen and hydrogen isotopic composition ( δ 18 O and δ 2 H) of the Pulp í paleo-aquifer from which the crystals formed by using stable isotopes of gypsum hydration water. We successfully dated the onset of gypsum precipitation in the Geode at 164 ± 15 ka. However, the extremely low U concentration (<11 ppb) and relatively high detrital Th content ( 230 Th/ 232 Th < 3.2) hinder accurate dating other gypsum samples. The δ 18 O and δ D values of the paleo-aquifer during the growth of the crystals aligned with the local meteoric water line, suggesting that the sulfate-enriched mother solution consisted of meteoric water that recharged the aquifer during that period. The mean isotopic composition of the Pulp í paleo-aquifer ( δ 18 O = − 6.5 ± 0.1‰ and δ 2 H = − 42.3 ± 0.5‰) during the formation of the crystals was similar to the current groundwater in this area ( δ 18 O = − 6.1 ± 0.8‰, δ 2 H = − 42 ± 6‰). The isotopic differences observed in samples collected from distinct locations and in individual crystals were probably related to changes in the isotopic composition of the aquifer, as a consequence of varying climate that impacted on the isotopic composition of rainwater during thousands of years in this region. Our results indicated that subaqueous selenite crystals may be useful for paleo-hydrological reconstructions. However, improving the current analytical techniques for dating gypsum with low U concentrations will be essential to obtain accurate and reliable records from Quaternary gypsum cave crystals in the future. Our results suggest that the isotopic composition of subaqueous gypsum crystals can be used for paleo-hydrological reconstructions. However, to achieve this objective, new analytical advances in obtaining accurate and precise ages of hydrothermal gypsum with sub-ppb U levels are required. We anticipate that gypsum crystals formed from unevaporated freshwater will be utilized in the future to reconstruct the isotopic composition of paleo-aquifers or hydrothermal ﬂuids, with important implications for paleoclimate research.


Introduction
The Giant Geode of Pulpí (Almería, SE Spain) is a cavity where the surfaces are completely covered with prismatic selenite crystals (CaSO 4 ·2H 2 O) up to 2 m long [1,2] ( Figure 1). The Geode was discovered in 1999 within the Mina Rica of Pilar de Jaravía and is now an important touristic and scientific attraction due to its spectacular characteristics. Indeed, the Giant Geode of Pulpí is currently the only gypsum geode open to the public worldwide [3]. To date, some investigations have addressed the mechanisms that led to precipitation of the massive selenite crystals [1,2,4]. Fluid inclusion microthermometric analyses revealed that subaqueous gypsum precipitation occurred at temperatures from 10 to 35 • C, with the majority of analyses suggesting 20 to 25 • C [4] from a low-salinity solution (from 0.0 to 0.5 wt % eq NaCl; [2]). These results discard gypsum crystallization from evaporated or unevaporated seawater, in addition to precipitation from a high-temperature hydrothermal solution. The possible presence of thermal water at deeper levels did not significantly influence the crystallization of gypsum, but controlled the previous stages of barite, celestine and Fe-oxide precipitation [4]. Sulfur and oxygen isotope analyses suggested that the source of sulfate for the first stages of gypsum precipitation was the dissolution of microcrystalline anhydrite (CaSO 4 ) and interlayered Triassic gypsum hosted in the carbonate bedrock [2,4]. The 87 Sr/ 86 Sr values in gypsum and celestine in the Mina Rica mine (0.7104 ± 0.0001 and 0.7106 ± 0.0001, respectively, [5]) are significantly higher than in seawater during the Phanerozoic and at present ( 87 Sr/ 86 Sr < 0.709; [6]). Dissolution of Triassic gypsum ( 87 Sr/ 86 Sr = 0.707 − 0.708; [7]) and Miocene gypsum ( 87 Sr/ 86 Sr < 0.709; [8]), which are abundant in the surrounding bedrock of the Sierra del Aguilón, cannot explain the strontium isotope values observed in gypsum and celestine. Weathering of silicate minerals in the soils has been proposed as a source of strontium for the relatively high 87 Sr/ 86 Sr values observed in the Giant Gypsum Geode [5], since strontium from the edaphic reservoir commonly has high isotope ratios [9].
U-series analyses of a carbonate coating on a gypsum crystal in the shallower mine passages (~30 m above the Giant Geode level) yielded an age of 60 ka. This indicates that gypsum precipitation, at that level, occurred prior to that time [4]. Despite several attempts being made to date the Geode, an absolute age of its gypsum crystals has not been provided so far. In this study, we performed the first successful U-series analyses of the Giant Gypsum Geode. Furthermore, we studied the oxygen and hydrogen stable isotopes (δ 18 O and δ 2 H) of gypsum hydration water to reconstruct the isotopic composition of the groundwater from which the Geode crystallized, in order to: (1) investigate the characteristics of the paleo-aquifer water (i.e., meteoric water recharge vs. presumed marine intrusion) and (2) to evaluate the potential and drawback of subaqueous selenite crystals as archives to reconstruct the δ 18 O and δ 2 H values of paleo-aquifers.

Sampling
Nineteen gypsum samples (~1-5 g each) from the Mina Rica mine were collected in December 2019 at different altitudes, ranging from 93 to 128 m above sea level (asl) In the particular case of the Giant Gypsum Geode (98 m asl), we collected a sample from the base of a crystal (MI-15) that was in contact with the underlying celestine and, therefore, represented the first stage of gypsum precipitation at the Geode level. In addition, a 22 cm long gypsum crystal (crystal GE-1), which was removed from the Geode entrance during habilitation works, was taken for analysis. Sample MI-15 and the base of the crystal GE-1 were analyzed for U-series dating (Table 1). We attempted to date 3 additional gypsum samples but analyses were unsuccessful because of the low U content (<2 ppb). In addition, 22 subsamples (~50 mg) for stable isotope analyses were extracted along the main growth axis of the crystal GE-1, every centimeter using a dentist tool (Dremel ® ) ( Table 2). Overall, we obtained the age of 2 samples by U-series dating and 42 samples were analyzed for oxygen (δ 18 O) and hydrogen (δ 2 H) isotopes of gypsum hydration water.

U-Th Dating
U-series analyses were performed at the Institute of Geology and Mineralogy (University of Cologne, Germany). The gypsum was ground to a fine powder that was divided into aliquots (3 for sample GE-1 and 4 for sample MI-15). U-Th separation was conducted by column chemistry. For that, gypsum samples (~150-200 mg) were dissolved in 10 mL 7 M HNO 3 and placed on a hotplate at 120 • C overnight. A mixed 229 Th-233 U-236 U tracer was added to the solution, which was again put on the hotplate at 120 • C overnight for sample spike equilibration and subsequently dried down. The 229 Th-233 U-236 U mixed tracer spike at Cologne is a combination of the IRMM-3636 U double spike material with a nominal 233 U/ 236 U ratio of 1.01906 [10] and 229 Th spike material resulting in a 229 Th/ 236 U ratio of 1.1909 ± 0.0044 × 10 −2 dissolved in HCl-HF. Accuracy of the tracer calibrations was further ensured by repeated measurements of silicate rock standards that were close to or near secular equilibrium In the next step, the samples were taken up again in 10 mL 7 M HNO 3 , and 2 mL concentrated HCl and 0.5 mL H 2 O 2 were added to remove any organic material. The samples were dried down again and re-dissolved in 10 mL 7 M HNO 3 , ready for chemical separation. The samples were passed through a two column chemical separation procedure. We used BioRad AG1-X8 resin (100-200 mesh) in custom-made columns with a column volume (CV) of 1 mL. The resin was cleaned using 6 M HCl, H 2 O and 7 M HNO 3 in several steps and conditioned, before the sample solution was passed through the column in 7 M HNO 3 media in individual 1 mL steps, until the whole sample was loaded onto the column. The matrix was rinsed from the resin with two CV of 7 M HNO 3 . Subsequently, the Th fraction was collected in 5 individual CV 6 M HCl, then U was eluted using 3 CV 1 M HBr. Both fractions were dried down, Th was taken up again in 1 mL 7 M HNO 3 and passed through the same AG1-X8 column to purify the Th cut from any residual Ca. Both the purified Th and the U fraction were dried down, treated with H 2 O 2 to remove any organic residue from the resin, re-dissolved in 0.2 M HCl + 0.01 M HF (Th) and 0.14 M HNO 3 (U) for measurement. The U-Th separation of each aliquot was conducted in independent columns and the concentrations and isotope ratios of each aliquot were measured individually (Table 1). By replicating the column chemistry separation, the external reproducibility of the method can be assessed. This is especially important for samples with low U and Th content, which are more prone to contaminations or background effects.
Concentrations and isotope ratios were measured using a Thermo Scientific Neptune MC-ICPMS with a central SEM. We used an Aridus II desolvator system following a standard sample bracketing procedure. For U, CRM112A reference material was used, and for Th, IRMM35 and IRMM36 standards doped with CRM112A were used. Measured U isotope ratios are corrected for mass bias using the known ratio of the 233 U/ 236 U double spike. Thorium samples were doped with the CRM112A U standard, using the known 238 U/ 235 U ratio (137.88) for mass bias correction. The 'Faraday Cup-SEM' yield was corrected using the known 234 U/ 238 U or 230 Th/ 232 Th ratio of the respective bracketing standard. We used a mathematical correction in which the U-Th isotopic composition of the detrital contamination was arbitrary estimated, in a manner resembling the approach by Ludwig and Paces [11]. 232 Th was used as the index for the correction and assumed a typical crystal Th/U ratio of 3.5 (atomic ratio), with 234 U/ 238 U and 230 Th/ 238 U activity ratios near secular equilibrium. In our model, we specifically used the initial activities ratios 232 Th/ 238 U = 1.21 ± 50%, 234 U/ 238 U = 1 ± 10% and 230 Th/ 238 U = 1 ± 10%. Ages and initial ( 234 U/ 238 U) ratios were calculated using an iterative approach using the decay constants of Cheng et al. [12] and Jaffey et al. [13]. The respective age uncertainties of each aliquot were derived from Monte Carlo simulations taking the uncertainties of the activity ratios into account. The ages of all the aliquots of each sample (n = 4 for MI-15 and n = 3 for GE-1) were averaged and 1 standard deviation (1σ) was calculated.

Stable Isotopes of Gypsum Hydration Water
The powdered gypsum samples were dried at 45 • C for 24 h, and then placed under vacuum to a pressure of~5 × 10 −3 mbar for at least 3 h at room temperature to remove adsorbed water. This low vacuum pumping is effective at removing adsorbed water with no detectable loss of hydration water [14]. The isotopic analyses of gypsum hydration water were conducted at University of Almería (Almería, Spain) using a Heat Induction Module (IM-CRDS, Picarro©) coupled to a Cavity Ring-Down Spectrometer (CRDS, Pi-carro© L2140i). The water vapor released after heating the powdered samples (~10 mg) by the IM to a final temperature of 250 • C was transferred on-line to the CRDS analyzer that measured the isotopic values (δ 18 O and δ 2 H) of the vapor with a frequency of 1 s. The Gaussian peak generated from the samples was integrated by the Picarro software that used an approach similar to that described in Bauska et al. [15]. The results were standardized to the Vienna Standard Mean Ocean Water (V-SMOW), Standard Light Antarctic Precipitation (SLAP) and Greenland Ice Sheet Precipitation (GISP) by analysing four gypsum standards before each set of 30-35 samples. The gypsum standards were previously calibrated against liquid water using the cryogenic extraction method by Gázquez et al. [14]. Each sample was analyzed 3-4 consecutive times. The mean analytical precision (1SD) of the 3-4 repetitions was 0.1‰ for δ 18 O and 0.9‰ for δ 2 H. The drift of the CRDS instrument was monitored and corrected (if needed) by measuring one of the gypsum standards every 10-12 samples.
The isotopic composition of the gypsum original solution (called paleo-aquifer in this study) was calculated by using isotope fractionation factor at 25 • C, since gypsum in the Geode crystallized at 20-25 • C [4]; 1.0034 for α 18 O gypsum-water and 0.981 for α 2 H gypsum-water [16], where: and, These fractionation factors, especially α 18 O gypsum-water , are very insensitive to temperature and have been previously used to reconstruct the isotopic composition of ancient water bodies [8,[17][18][19].

The Age of the Giant Geode
The U concentration in samples MI-15 and GE-1 were~11 ppb and~3 ppb, respectively, ( Table 1). These concentrations are lower than in subaerial gypsum stalactites analyzed in previous studies (i.e., 94-150 ppb U; [19]), but are in agreement with the U content reported in the subaqueous selenite crystals from the caves of Naica (2 to 10 ppb U, [20]). We attempted to analyze other gypsum samples from the Mina Rica mine, but the U contents were even lower (<2 ppb), impeding the obtainment of additional ages. As 230 Th/ 232 Th activity ratios in the analyzed samples were <20 (2.0 to 3.2), a correction for 230 Th contamination was needed to obtain accurate ages. The corrected ages were 21-56 ka younger than the uncorrected ages ( Table 1). The age of sample MI-15, that represents the early stage of gypsum precipitation at the Geode level, was 164 ± 15 ka (n = 4). This indicated that gypsum precipitation in the Giant Geode started during the Marine Isotope Stage 6 (MIS 6) (ca. 191 to ca. 131 ka; [21]). The age of sample GE-1 was 108 ± 87 ka (n = 3), more imprecise because of the lower U content, but within analytical errors with the age of sample MI-15. These results indicated that gypsum samples with U content > 10 ppb can be dated by the method used here, however lower U concentrations result in more imprecise and unreliable ages.
Previous studies reported U-series ages of a calcite coating that covers selenite crystals in the upper mine levels (30 m above the Giant Geode level) [4]. Although these authors did not provide information about the origin of this coating (i.e., subaerial/subaqueous) they interpreted that conditions for calcite precipitation prevailed at 60 ka; this may indicate that gypsum precipitation in the Sierra del Aguilón mountain range, and probably in the Giant Geode, ceased before that time. Thus, we concluded that gypsum precipitation in the Giant Geode extended during a maximum period of 100 ka, from 164 ± 15 ka to sometime before 60 ka.

Paleoclimatic Significance of Subaqueous Gypsum Crystals
The selenite crystals of the Giant Geode and the Mina Rica mine formed in subaqueous conditions [2]. During crystallization, the two molecules of structurally-bound hydration water of gypsum (CaSO 4 ·2H 2 O) are incorporated from the mother solution (i.e., paleoaquifer). In consequence, the δ 18 O and δ 2 H values of gypsum hydration water (δ 18 O gypsum and δ 2 H gypsum hereafter) reflect those of the mother solution (δ 18 O aquifer and δ 2 H aquifer hereafter). The oxygen and hydrogen isotope fractionation factors for the gypsum-solution pair (α 18 O gypsum-water and α 2 H gypsum-water ) are well-known and very insensitive to temperature [16]. Thus, by using the δ 18 O gypsum and δ 2 H gypsum the isotopic composition of the paleo-aquifer at the time of gypsum precipitation (δ 18 O aquifer and δ 2 H aquifer ) can be calculated ( Table 2). Here we use this property of gypsum to reconstruct the isotopic composition of the paleo-aquifer of Pulpí during gypsum crystallization.
The δ 18 O aquifer during the formation of the investigated samples ranged from −5.3 to −7.1‰ (mean ± 1 standard deviation = −6.2 ± 0.4‰) and the δ 2 H aquifer ranged from −36.6 to −45.0‰ (−40.5 ± 2.2‰). The d-excess of the paleo-aquifer (d-excess = δ 2 H aquifer − 8 × δ 18 O aquifer ) ranged from 3.4 to 14.4‰ (9.1 ± 2.3‰). These values are close to those of modern rainwater and groundwater in this region (δ 18 O = −6.1 ± 0.8‰, δ 2 H = −42 ± 6‰, d-excess~6.3 ± 4.8‰; [22,23]) ( Figure 2). In particular, the samples collected from the Giant Geode (sample MI-15 and 22 subsamples from the crystal GE-1) crystallized from a solution with a δ 18 O aquifer value of −6.0 ± 0.2‰, δ 2 H aquifer value of −39.1 ± 1.0‰ and d-excess of 8.5 ± 2.1‰. The δ 18 O aquifer and δ 2 H aquifer pairs are aligned with modern rainwater values in this region ( Figure 2). This indicates that seawater (with δ 18 O~0‰ and δ 2 H~0‰) or fluids from deeper magmatic/metamorphic processes (usually with higher δ 18 O values; [24]) did not contribute significantly to the solution from which gypsum precipitated. In contrast, the gypsum crystals recorded the isotopic composition of meteoric water seepage that recharged the paleo-aquifer of Sierra del Aguilón at ca. 164 ka and during MIS 6. Indirectly, these results indicate that during the formation of the Geode the Mediterranean Sea level was at a similar or lower position than at present (Figure 3). Otherwise, seawater would have influenced the solution by increasing salinity and isotopic values, which was not seen in our results of gypsum hydration water, or in previous results of microthermometry of fluid inclusions [2]. A relatively low sea level position during the formation of the Geode is supported by Mediterranean Sea level reconstructions from subaqueous carbonate speleothems that recorded a lowstand of~−60 m relative to present during MIS 6 [25].
Given the uncertainties in the ages of the analyzed samples and the lack of additional dates as a result of the difficulties we found in the U-Th analyses, robust interpretations about the evolution of the Pulpí paleo-aquifer are not possible. However, the range of isotopic values observed between the various groups of samples from Mina Rica (Giant Geode vs. the rest of the mine levels) can be due to the fact that gypsum formed under different climate conditions. The isotopic variability of the aquifer water was likely controlled by changes in the δ 18 O and δ 2 H values of rainwater, as a result of long-term variations in atmospheric temperature and changes in the source of moisture. Indeed, both parameters govern the modern seasonal variations in δ 18 O, δ 2 H and d-excess values of rainwater in this region [26,27], and also in the past [19,28]. Lower δ 18 O and δ 2 H values are expected for meteoric water seepage during glacial (colder) periods [28][29][30]. In contrast, higher isotopic values of the recharge indicate rainwater infiltration to the aquifer during an interglacial (warmer) period (e.g., modern conditions) [28][29][30]. Figure 3. Conceptual sketch of gypsum precipitation in the Giant Geode of Pulpí from at least 164 ± 15 ka to sometime before 60 ka. Our results of stable isotopes in gypsum hydration water suggest a prevailing source of sulfate-enriched water from the paleo-aquifer, with the sea level located below the Geode level (currently at 97 m asl). Meteoric water seepage that recharged the aquifer and the solution from which the gypsum crystals were formed was freshwater, with insignificant contributions of seawater. Dissolution of Triassic gypsum and anhydrite associated to phyllites was the main source of sulfate to the solution. Note that the scale is exaggerated and a precise topography omitted for a better display.
Studies of stable isotopes in fluid inclusions in subaqueous carbonate speleothems from the central Mediterranean reported a change in δ 18 O by 1.7‰ in paleo-groundwater waters from the glacial MIS8 to the interglacial MIS5c [29]. Similarly, a δ 18 O difference by 2‰ in paleo-groundwater during the transition from MIS6 to MIS5 was found in a subaqueous flowstone from Devil's Hole (southwest USA) [30]. Considering that the range of reconstructed δ 18 O values of the Pulpí paleo-aquifer ranged by 1.7‰ during crystallization of gypsum, it is possible that the gypsum precipitation extended over a long period of time and that the crystals recorded major climate changes.
In particular, the δ 18 O values of the Pulpí paleo-aquifer at 164 ± 15 ka during the onset of gypsum precipitation in the Giant Geode (−6.5‰, sample MI-15) are more positive than during the crystallization of crystals stratigraphically younger inside the geode (maximum δ 18 O aquifer of −5.3‰), but more negative than in the fluids that generated gypsum in other mine levels (minimum δ 18 O aquifer of −7.1‰). A paleoclimate reconstruction from carbonate speleothems (i.e., Gitana Cave [21],~50 km from Pulpí) suggests that dry/cold conditions during MIS6 prevailed between 175-154 ka and peaked at~157 ka, while there is evidence for a wetter/warmer interstadial around 151 ka. Sea surface temperatures reconstructed from Alboran Sea sediments (core ODP 077; [31]) show that temperature declined gradually by~4 • C from~180 to~140 ka. Altogether, these paleoclimatic reconstructions demonstrate that despite climatic conditions during MIS6 in southern Iberia being generally cold and dry, there were significant variations in rainfall amount and temperature that probably resulted in changes in the δ 18 O of rainwater in this region. It is possible that the gypsum crystals in the Aguilón Mountain Range formed during different climatic phases of MIS6, or even during MIS5 in the case of the younger samples with higher δ 18 O and δ 2 H values of gypsum hydration water.

Conclusions
U-series dating of gypsum has revealed that the Giant Geode of Pulpí grew during the Upper Pleistocene, starting at least 164 ± 15 ka and probably until 60 ka. The Pulpí paleo-aquifer was occupied by freshwater at that time, rather than by brackish solution or seawater. The δ 18 O and δ 2 H values of the paleo-aquifer during the formation the gypsum crystals varied by 1.7‰ for δ 18 O and by~9‰ for δ 2 H. This indicates that gypsum in the Aguilón Mountain range crystallized probably over a long period of time, in which a succession of wetter/drier and colder/warmer climatic conditions took place. Climate variability led to changes in the isotopic composition of rainfall and resulted in secular variations of the isotopic composition of the Pulpí paleo-aquifer that were ultimately recorded by the isotopic composition of the gypsum crystals. Table 1. U-Th results and calculated ages of samples MI-15 (4 repetitions) and GE-1 (3 repetitions) from the Giant Gypsum Geode. We used a mathematical correction in which the U-Th isotopic composition of the detrital contamination was arbitrarily estimated, in a manner resembling the approach by Ludwig and Paces [11]. 232 Th was used as index for the correction and assumed a typical crustal Th/U ratio of 3.5 (atomic ratio), with 234 U/ 238 U and 230 Th/ 238 U activity ratios near secular equilibrium. In our model, we specifically used the initial activity ratios 232 Th/ 238 U = 1.21 ± 50%, 234 U/ 238 U = 1 ± 10% and 230 Th/ 238 U = 1 ± 10%.   Our results suggest that the isotopic composition of subaqueous gypsum crystals can be used for paleo-hydrological reconstructions. However, to achieve this objective, new analytical advances in obtaining accurate and precise ages of hydrothermal gypsum with sub-ppb U levels are required. We anticipate that gypsum crystals formed from unevaporated freshwater will be utilized in the future to reconstruct the isotopic composition of paleo-aquifers or hydrothermal fluids, with important implications for paleoclimate research.