Last glacial millennial-scale hydro-climate and temperature changes in Puerto Rico constrained by speleothem fluid inclusion δ¹⁸O and δ²H values

We present speleothem fluid inclusion δ¹⁸O_f and δ²H_f values from Larga Cave, Puerto Rico, that cover the interval between 46.2 and 15.3 ka on the millennial scale, including the Last Glacial Maximum (LGM) and several stadial and interstadial cycles. The data set can be divided in two main clusters of stable isotope compositions of the fluid inclusion water with respect to the global meteoric water line (GMWL), which coincide with strong variations in the water content of the stalagmite. In particular, this clustering is found to be climate related, where one cluster comprises samples from cold and dry periods, such as the Heinrich and Greenland stadials (HSs and GSs), as well as parts of the LGM, which exhibit very high δ¹⁸O_f and δ²H_f values. We interpret this enrichment as being caused by evaporation inside the cave due to enhanced cave ventilation during these colder and drier times. In contrast, in most samples corresponding to warmer and wetter Greenland interstadials (GIs), but also for some from HS 2 and 3, the δ¹⁸O_f and δ²H_f values plot on the meteoric water line and modification of fluid inclusion water due to “in-cave” evaporation are found to be negligible.

Consequently, variations of recent glacial hydro-climate and temperatures in the western tropical Atlantic can be constrained. In general, δ¹⁸O_f values from fluid inclusions are up to 3 ‰ higher than those of modern drip water, which is interpreted as a weaker atmospheric convective activity during the last glacial period. In addition, reconstructed temperatures suggest an average cooling of 2–3 ^∘C during the LGM compared to modern cave temperatures. Reconstructed cave temperatures yield an average cooling of −1.4 ± 2.8 ^∘C for HS 2 and −3.6 ± 2.2 ^∘C for HS 3. Higher δ¹⁸O_f values of these samples further suggest that the drip water was dominated by orographic rainfall and/or cold fronts, along with weak or even absent convective activity. In contrast, during interstadial phases, reconstructed temperatures reached nearly modern values, and convective activity was comparable to or only slightly weaker than today.

1 Introduction

Stable isotopes of oxygen (represented by δ¹⁸O values) are a widely used proxy used to deduce information about past climates from speleothems, since the cave drip water reflects the amount- or recharge-weighted annual mean isotopic signature of the rainfall above the cave (Baker et al., 2019; Lachniet, 2009). This value is in turn modulated by several variables, such as temperature, the type and seasonal distribution of precipitation, and the source moisture trajectory (Lachniet, 2009; Bird et al., 2020; Lases-Hernandez et al., 2020). In the tropics, the isotopic composition of precipitation, δ¹⁸O_p, has been attributed to changes in precipitation amount and convective activity (Dansgaard, 1964; Lachniet and Patterson, 2009; Vieten et al., 2018a). Consequently, tropical speleothem δ¹⁸O values are usually interpreted in terms of this isotopic “amount effect” (Medina-Elizalde et al., 2010; Arienzo et al., 2019; Lases-Hernandez et al., 2019). However, in the soil and karst above caves and at the surface of stalagmites, various (disequilibrium) isotope fractionation effects can obscure the interpretation of stable isotopes in speleothem calcite (Weber et al., 2021; Hansen et al., 2019; Carlson et al., 2020).

In this respect, stable isotopes of oxygen and hydrogen in speleothem fluid inclusions provide a useful proxy since these represent potentially unchanged aliquots of the original drip water, and thus relicts of past precipitation (Schwarcz et al., 1976; van Breukelen et al., 2008; Affolter et al., 2015). The utility of speleothem fluid inclusion analysis improves if the conditions of carbonate precipitation are well understood, i.e. whether or not the minerals precipitated close to “equilibrium conditions”. In the case of equilibrium conditions, palaeo-temperatures can be reconstructed from the combined oxygen isotope ratios of the calcite (δ¹⁸O_c) and the fluid inclusions, δ¹⁸O_f, based on known temperature-dependent isotope fractionation relationships between calcite and water (e.g. Kim and O'Neil, 1997; Tremaine et al., 2011; Hansen et al., 2019; Johnston et al., 2013). In continental mid- and high-latitude locations, estimation of palaeo-temperatures via fluid inclusion δ²H_f values alone is feasible, based on the dominant control of temperature on rainfall isotope composition in such settings (Affolter et al., 2019; Demény et al., 2021). This only works, however, if the contribution of the amount effect is negligible. The analysis of speleothem fluid inclusions has been successfully applied worldwide, including tropical locations (Griffiths et al., 2010; Arienzo et al., 2015; Millo et al., 2017; Meckler et al., 2015). Nevertheless, speleothem fluid inclusion records from the Last Glacial Period (ca. 115–11.7 ka) are still sparse, and in many cases, reliable palaeo-temperature estimates are challenging (Millo et al., 2017; Demény et al., 2021; Nehme et al., 2020).

Here we present a detailed study of fluid inclusions from a Puerto Rican stalagmite (PR-LA-1), spanning the interval between 46.2 and 15.3 ka. The record encompasses a period which was characterized by rapid climatic fluctuations from cold (stadial) to warmer (interstadial) conditions recurring on millennial timescales (Dansgaard et al., 1984; Peterson et al., 2000). In palaeoclimate records from the tropics, the Heinrich stadials (HS) are clearly represented by dramatic cooling and dry conditions in the western tropical Atlantic (Escobar et al., 2012; Arienzo et al., 2015; Grauel et al., 2016). These conditions were likely caused by a southward shift of the Inter-Tropical Convergence Zone (ITCZ) induced by a near shutdown of the Atlantic Meridional Overturning Circulation (AMOC; Broccoli et al., 2006; Lynch-Stieglitz et al., 2014; Deplazes et al., 2013; Escobar et al., 2012). Conversely, during the interstadials, the ITCZ was presumably in a more northerly position, accompanied by warmer and wetter climatic conditions and higher sea surface temperatures (SSTs) in the (sub-)tropical Atlantic (Schmidt et al., 2006a; Peterson et al., 2000; Sachs and Lehman, 1999).

The speleothem PR-LA-1 was selected because its stable isotope values (δ¹⁸O_c and δ¹³C_c) and Mg $/$ Ca ratios demonstrated a pronounced response to the local hydro-climate on the millennial scale (Warken et al., 2020). To further elucidate the millennial-scale climate changes during the last glacial stadials and interstadials, we thus conducted a large number of fluid inclusion measurements following prominent δ¹⁸O_c changes in stalagmite PR-LA-1.

2 Material and Methods

2.1 Site description

Larga Cave (18^∘19^′ N, 66^∘48^′ W; Fig. 1) is located 350 m above sea level in the north central karst region of Puerto Rico. The overburden of Larga Cave is comprised of the Oligocene Lares limestone Fm., which is covered by a thin layer of soil. The host rock is exposed in some places, and the present-day vegetation consists of dense tropical forest. The main passage in the cave is almost horizontal and has a total length of 1400 m with a ceiling height of up to 30 m (Miller, 2010). Stalagmite PR-LA-1 has a total length of 185 cm and was found lying on the cave floor in the main passage, where the rock overburden is 40 to 80 m thick. The stalagmite was deposited between 15.4 and 46.2 ka, with a hiatus between 35.5 and 41.1 ka (Fig. 1; Warken et al., 2020). Microscopic inspection of thin sections shows that PR-LA-1 consists of calcite in columnar fabric and shows no signs of dissolution and/or recrystallization (Fig. S1). Between 2012 and 2017, the mean annual air temperature at the location of the cave was 22.5 ^∘C, and the measured mean local annual rainfall was 2200 mm (Vieten et al., 2018a). Since 2012, cave air parameters, such as temperature (T), pCO₂, and relative humidity (RH), have been recorded at the location of stalagmite PR-LA-1. During this time, the cave temperature has been constant with values of 22.5 ± 0.2 ^∘C. The cave air pCO₂ varies depending on the season, with high values of up to 1800 ppm in summer and low values of around 550 ppm in winter (Vieten et al., 2016). This seasonal pattern results from a temperature-driven, very well-ventilated cave environment in winter and low ventilation in summer. During the winter season, a high buoyancy contrast at night leads to maximal cave ventilation, while it is at a minimum during summer, when cave temperatures are lower than in the outside environment. RH is generally close to 100 %, but has been found to decrease to values to 90 %–95 % for several weeks when ventilation intensifies in winter months (Vieten et al., 2016, 2018a).

Figure 1 Sample locations and previous work on stalagmite PR-LA-1. (a) Location of Larga Cave (1), Puerto Rico, in the western tropical Atlantic region and other sites discussed in this study (map after Warken et al., 2020): (2) Guiana Basin (Rama-Corredor et al., 2015); (3) Cariaco Basin, Venezuela (Deplazes et al., 2013); (4) Lake Petén Itzá, Guatemala (Escobar et al., 2012; Grauel et al., 2016); (5) Florida Strait (Them Il et al., 2015); (6) Abaco, Bahamas (Arienzo et al., 2015). The insert shows a detailed view of the location of Larga Cave (red star), and the position of local meteorological stations: Arecibo Observatory (AO), San Agustín (SA), Luquillo Mountains (LM), and Guánica Dry Forest (GD). Blue arrows in panel (a) indicate the dominant moisture trajectories over the region in boreal summer, when the ITCZ (yellow shaded area) is in its northernmost position. (b) Local rainwater isotope data from Larga Cave (Vieten et al., 2018b), SA (Scholl et al., 2014) and Global Network of Isotopes in Precipitation (GNIP) station San Juan (SJ). In addition, the global meteoric water line (GMWL) and local meteoric water line (LMWLs) are shown from SA (Putman et al., 2019; Scholl et al., 2014); LM (Scholl and Murphy, 2014) and GD (Govender et al., 2013). (c) Image of stalagmite PR-LA-1 and age–depth model from Warken et al. (2020).

As typical for Caribbean locations, the precipitation regime over Larga Cave is characterized by a drier winter season, from December to April and a much wetter season during the rest of the year. During the drier season, occasional light showers occur imbedded in trade winds, which are interrupted by occasional cold fronts emanating from North America (Scholl and Murphy, 2014). In contrast, during the remainder of the year convective rain is formed by easterly moving low pressure systems which occasionally develop into tropical depressions accounting for the majority of total precipitation (Scholl and Murphy, 2014). On inter-annual timescales, the δ¹⁸O values of rainfall in Puerto Rico show a negative correlation with precipitation, suggesting an isotopic amount effect of ca. −0.1 ‰ per 100 mm of rainfall amount (Govender et al., 2013; Scholl et al., 2014; Warken et al., 2020). Over the course of monthly drip water monitoring (2013–2017), the drip water composition at sites with a maximum distance of 10 m to the location where stalagmite PR-LA-1 was found shows only little seasonal variation, with a mean value of δ¹⁸O_w = −2.6 ± 0.3 ‰, and δ²H_w = −9.5 ± 0.2 ‰ Vienna Standard Mean Ocean Water (VSMOW). The seasonal signal in the rainwater is strongly smoothed by the soil and epikarst above the cave with a mixing time of at least several months to years, suggesting that the Larga speleothem time resolution record is at least inter-annual.

A LMWL is available for three stations in Puerto Rico (Fig. 1). At San Agustín, located closest to Larga Cave, a monthly data set covering 1998 to 2012 shows a LMWL of δ²H = 8.4 × δ¹⁸O + 11.6 ‰ (VSMOW; Putman et al., 2019; Scholl et al., 2014). Similar LMWLs (Fig. 1) exist for a site in the Luquillo Mountains in eastern Puerto Rico (Scholl and Murphy, 2014) and for the Guánica Dry Forest, located in south-western Puerto Rico (Govender et al., 2013).

2.2 Fluid inclusion stable isotope composition

Solid samples for fluid inclusion isotope analysis were cut using a band saw along the growth axis of stalagmite PR-LA-1. High growth rates of this specimen of up to 1 mm yr⁻¹ allow for the sampling of distinct millennial-scale climate oscillations evident in the δ¹⁸O_c record, such as D/O or Heinrich events. Replicate samples were taken along the corresponding growth layers to minimize the effect of integrating over too-long time intervals. Individual sample sizes were between 0.21 and 1.34 g. The speleothem fragments were placed in a copper tube and installed in a preparation line connected to a cavity ring down spectrometer (CRDS, L2130-I, Picarro) at the Institute of Environmental Physics, Heidelberg University (Weißbach, 2020). After connecting the sample, the setup requires approximately 90 min until the water vapour background has reached stable conditions. When the standard deviation of the water vapour concentration is below 20 ppmV over 30 min, the sample is crushed with a hydraulic crusher. To avoid the CRDS memory effect, a stable water vapour background was generated to ensure constant saturation of the cavity without condensation (Affolter et al., 2014; De Graaf et al., 2020). This was achieved by mixing dry nitrogen as a carrier gas with water of a known isotopic composition, while the line, including the tube that holds the sample, was held at a constant temperature of 120 ^∘C.

The mixture of background water and sample signal is measured in the CRDS analyser, and evaluated similarly to the method of Affolter et al. (2014) with a Python code (Hemmer, 2022). The isotopic composition and released water volume are calibrated with in-house water standards and glass capillaries with known volumes in the microlitre range respectively. These are inserted into the copper tube and connected to the preparation line, similarly as the speleothem samples. The in-house water standards are analysed regularly against the international primary reference material Vienna Standard Mean Ocean Water 2 (VSMOW2) and Standard Light Antarctic Precipitation 2 (SLAP2). The long-term precision of single stable isotope measurements is 0.5 ‰ for δ¹⁸O and 1.5 ‰ for δ²H, while for released water volumes above 1.0 µL, the precision improves to 0.2 ‰ for δ¹⁸O and 0.4 ‰ for δ²H (Weißbach, 2020). In the case of water volumes < 0.2 µL, the result of a single measurement may become less precise and is therefore reproduced by replicate measurements to verify the result. If not stated differently, the stable isotope composition of fluid inclusions is given relative to VSMOW.

3 Results

3.1 Water content and stable isotope composition of speleothem fluid inclusions

In total, 64 fluid inclusion water isotope measurements were performed at 30 different depth intervals. The results of the single measurements are shown in Table S1 in the Supplement, while the mean values of the replicate measurements are given in Table S2. The mean water content of PR-LA-1 calcite for the different depth intervals shows a variation of close to 0 to 3 µL g⁻¹ (Fig. 2). Sections with a relatively high water content, e.g. around ca. 24 or 28 ka, also show a high (or at least highly variable) growth rate, while comparably low water content, such as at ca. 20, 26, or 32 ka, are mainly found in sections with relatively slow growth rates. The mean fluid inclusion water δ¹⁸O_f values vary between −1.9 and +11.4 ‰, and mean δ²H_f values are between +6.9 and +44.9 ‰ (VSMOW; Fig. 3a). The data exhibit a systematic ¹⁸O and ²H enrichment compared to modern drip water stable isotope values. This general pattern is also reflected in the slope of the linear regression line fitted to the data, which is with 2.42 ± 0.25, significantly lower than expected for a LMWL. This slope agrees with those of globally distributed evaporation lines for lakes and soil waters (Gibson et al., 2008), as well as evaporation experiments (Hu et al., 2009).

Figure 2 Results of fluid inclusion analysis of stalagmite PR-LA-1 (Table S2). From top to bottom: (a) growth rate (logarithmic scale) and (b) δ¹⁸O_c values from Warken et al. (2020). Symbols indicate the averaged δ¹⁸O_c values used for each depth of fluid inclusion isotope analysis respectively. Vertical blue bars (dotted red lines) indicate the timing of cool/dry (warm/humid) periods, such as HS 1 to 3 (Greenland interstadials (GIs) 12c to 6) in the western tropical Atlantic (Warken et al., 2020). The following panels show fluid inclusion isotope values of (c) δ¹⁸O_f, (d) δ²H_f, and (e) measured water content. Symbol colours indicate the samples with high (dark coloured symbols) and low (light coloured stars) water content, with a threshold value of 0.55 µL g⁻¹ respectively. Horizontal lines indicate modern cave conditions. Note that the axes of the δ¹⁸O_c, δ¹⁸O_f, and δ²H_f values are inversed.

Figure 3 (a) Results of stable isotopes of fluid inclusions and drip waters in Larga Cave as indicated in Tables S1 and S2. Lines indicate the GMWL and the LMWL with δ²H = 8.2 × δ¹⁸O+14 ‰ VSMOW by Scholl and Murphy (2014). Single values of the water isotopes in fluid inclusions (small circles) and mean values (large circles) with corresponding standard deviation are shown together with Larga Cave drip water from different drip sites (2013–2019, green triangles) and local rainwater data (compare Fig. 1). The mean values of the replicate measurements are allocated to two groups, with cluster (1) comprising samples on or scattering around the GMWL (indicated by the dark blue circle), and cluster (2) comprising all samples with a clear evaporation signal, plotting off the GMWL (light blue circle). (b) Stable isotope values of fluid inclusions from stalagmite PR-LA-1, categorized according to their allocation to the presumably warm (interstadial) or cold (stadial) periods when they formed. Thereby, a distinction is made between Heinrich stadials, Greenland stadials, and Greenland interstadials. Analysed samples from the LGM are also subdivided into presumably more warm/humid and relatively cool/dry periods. Numbers refer to the sample IDs as shown in Tables S1 and S2.

Furthermore, the data can be allocated to two different clusters, depending on the measured water content of the samples. The first cluster comprises samples which are mainly located on or scatter around the GMWL. Most of these are also characterized by a water content greater than 0.55 µL g⁻¹. In contrast, many samples show more positive δ¹⁸O and δ²H values and clearly plot off the GMWL along the linear regression line. In contrast to the first set of samples, most of these have a water content lower than 0.55 µL g⁻¹. Only one sample clearly does not fit into this pattern and plots away from the GMWL despite a water content of > 0.55 µL g⁻¹ (sample 32b), which is therefore also allocated to cluster (2). Conversely, two samples plot close to the MWLs despite a relatively low water content (samples 5 and 8). The relationship between fluid inclusion stable isotope values and water content is also expressed numerically with a negative correlation of δ¹⁸O_f and the measured water content, with r = −0.72 (r² = 0.51) for logarithmic fitting (Fig. S2).

3.2 Calculation of palaeo-temperatures using the oxygen isotope equilibrium thermometer

In order to assess the applicability of the oxygen isotope equilibrium thermometer based on the mean fluid inclusion δ¹⁸O_f values, the average carbonate δ¹⁸O_c values were calculated from the intervals covered by the depth range of the measured fluid inclusion samples respectively (Table S2). Figure 4a shows that the previously mentioned pattern with two clusters of samples is also reproduced when comparing the fluid inclusion δ¹⁸O_f with the carbonate δ¹⁸O_c values. Samples with high water content, which plot on the GMWL in Fig. 3a, are clearly separated by the samples with lower water content. While cluster (1) samples with high water content show a positive correlation between δ¹⁸O_c and δ¹⁸O_f values, cluster (2) samples do not follow the same trend. Similarly, the fractionation between speleothem calcite and fluid inclusions, ${}^{\mathrm{18}}{\mathit{\epsilon }}_{\text{c/f}}=\mathrm{1000}\times \ln \mathit{\alpha }$ (with $\mathit{\alpha }=(\mathrm{1000}+{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\text{c}})\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}/(\mathrm{1000}+{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\text{f}})$, and δ¹⁸O_c and δ¹⁸O_f in VSMOW), shows a clear dependence on the δ¹⁸O_f value for the evaporated samples of cluster (2) (Fig. 4b), but shows no relationship with carbonate δ¹⁸O_c values (Fig. 4c). In contrast, for cluster (1) samples, no relationship between ¹⁸ε_c/f and δ¹⁸O_f values is observed (Fig. 4b). Figure 4 also suggests that fluid inclusion samples 9b, 32b, 27, 32, and 36 are also enriched in ¹⁸O, and have to be allocated to cluster (2).

Figure 4 Analysis of the fractionation between speleothem calcite and fluid inclusions. Panel (a) shows the measured δ¹⁸O_c values in speleothem carbonate (Warken et al., 2020) vs. measured δ¹⁸O_f in speleothem fluid inclusion water (this study). (b) Fractionation δ_c/f=1000ln α between carbonate and fluid inclusion water vs. the δ¹⁸O_f values. (c) Fractionation ¹⁸ε_c/f=1000ln α vs. δ¹⁸O_c values. Numbers indicate the sample IDs of the fluid inclusion measurements (Table S2). The uncertainty of the calculated fractionation is estimated from propagating the analytical error of the fluid inclusion measurements, and from the range of the δ¹⁸O_c calcite value due to dating and sampling uncertainty. Blue symbols indicate samples with water content > 0.55 µL g⁻¹, and light blue symbols water content < 0.55 µL g⁻¹. The grey lines indicate the water–calcite oxygen isotope fractionation relationship after Tremaine et al. (2011) for cave temperatures of 15, 20 and 25 ^∘C respectively. The ellipses indicate the samples which are allocated to cluster (2) due to large fractionation effects in δ¹⁸O_f. All of these samples also plot off the water–calcite oxygen isotope fractionation curve for 25 ^∘C.

For a first-order temperature calculation from ¹⁸ε_c/f, three different isotope fractionation relationships have been tested (Kim and O'Neil, 1997; Tremaine et al., 2011; Johnston et al., 2013). All results are listed in Table S2, where the given uncertainties of the palaeo-temperatures result from propagating the analytical error of the fluid inclusion measurements, as well as the range of the δ¹⁸O_c calcite value due to dating and sampling uncertainty. All palaeo-temperatures reported and discussed in the main text are calculated using the Tremaine et al. (2011) equation, which has been developed from cave monitoring and modern speleothem calcite and yields realistic temperatures compared to the modern cave air temperature (Table S2). Note that we discard samples from the temperature discussion where replicate measurements exhibit a poor reproducibility in terms of largely variable water content and large isotopic differences, in order to focus on the most robust palaeoclimatic inferences (Table S2). Note also that we do not report temperatures for samples of cluster (2), which clearly show significant influence from fractionation through evaporation (Table S2). We regard the temperature results obtained from the 12 replicate samples from cluster (1) as reliable within the large uncertainty of 1.5–2.5 ^∘C. The calculated temperatures scatter around the modern cave temperature of 22.8 ^∘C, with palaeo-temperatures varying between 16 ± 1.8 ^∘C (sample 8, 21.39 ± 0.25 ka) and 25.7 ± 1.4 ^∘C (sample 24, 33.52 ± 0.52 ka), and a mean value of 21.1 ± 2.7 ^∘C.

4 Discussion

4.1 Influence of evaporation on PR-LA-1 fluid inclusion isotope data

4.1.1 Categorization of the fluid inclusion isotope data

For the fluid inclusion isotope data from cluster (1), δ¹⁸O_f and δ²H_f values scatter around the LMWL (Govender et al., 2013; Putman et al., 2019; Scholl and Murphy, 2014; Scholl et al., 2014). Thus, we suggest that these samples represent pristine fossil drip water. The second cluster of fluid inclusion isotope values mainly contains samples with water content below 0.55 µL g⁻¹, which all show increased δ¹⁸O_f values that are noticeably more fractionated than hydrogen isotopes. In addition to categorization with respect to water content, the fluid inclusion samples can be further associated with stadial and interstadial phases during the Last Glacial according to their stable isotope and elemental signatures in speleothem calcite (Fig. 3b). The fluid inclusion samples representing warm GIs (red triangle) plot close to the MWL and can be all associated with cluster (1). In contrast, the samples formed during the colder intervals (HS, GS, LGM) are mainly attributed to cluster (2) and trend away from the MWL. In addition, for Last Glacial Maximum (LGM) samples, the differentiation between the two clusters appears to be valid, with relatively warm/humid phases plotting closer to the MWL, and relatively dry/cool phases rather clustering with the GS and HS samples (Fig. 3b). This demonstrates that our data set allows us to distinguish cooler and/or drier climatic background states during the Heinrich and Greenland stadials, represented by cluster (2) samples, from the presumably warmer and more humid phases during the Greenland interstadials, which comprise cluster (1). Note that evaporation-driven fractionation must indicate a general annual decrease of humidity in the cave spanning decades that is required to trap the water in the calcite. Therefore, we suggest that samples of cluster (2) reflect periods with stronger cave ventilation as recorded presently in the thermal contrast of cave and ambient air in winter. This categorization suggests that the underlying process, which causes this clustering and affects the fluid inclusion isotopic signatures, is purely climate related.

Isotope effects due to the partial evaporation of fluid inclusion water can in principle also occur during the warm-up phase in the crusher when fluid inclusions are insufficiently sealed to withstand the pressure that builds up during warm-up. Based on our experience, water loss in the warm-up phase causes noticeable disruptions in the water background of the analyses, but this has been observed mostly for diagenetically compromised samples so far. We note that the water vapour background was continuously monitored before crushing, and no anomalies were observed in the samples from PR-LA-1 (Weißbach, 2020). Figure S1 also shows no clear relationship between shape and distribution of fluid inclusions in PR-LA-1, nor do we see signs of diagenetic alteration in stalagmite PR-LA-1, which could have favoured inclusion leaking during heating. Previous studies have also invoked recrystallization-induced oxygen isotope changes in inclusion-hosted water of speleothems (Demény et al., 2016). This process is unlikely to be systematically related with climate and would shift the fluid inclusion δ¹⁸O values towards more negative values, which is opposite to the observed pattern in this study.

Other studies suggest that in situ oxygen isotope exchange of fluid inclusion water with the host carbonate is unlikely to cause significant isotope fractionation (Uemura et al., 2020; Demény et al., 2021). The systematic enrichment of cluster (2) samples indicates that the stable isotope composition of the low water content samples has been affected by fractionation processes in the gas phase as ¹⁸O is enriched over H isotopes, just as in non-equilibrium evaporation of thin water films (Gat et al., 1994). Thus, all information combined suggests evaporation inside the cave on the surface of the stalagmite as the scenario most likely to cause the lower slope of the water isotope relationship of 2.42 ± 0.25 (Fig. 3) via temperature- and humidity-related climate change (Gibson et al., 2008; Hu et al., 2009).

4.1.2 Climate-related evaporation effects on fluid inclusion isotope values in Larga Cave

Recent studies have suggested that in-cave evaporation may affect the isotopic composition of speleothem fluid inclusions (Nehme et al., 2020; Weißbach, 2020). A potential climate-related driver of this evaporative fractionation that has not been described previously is temperature-driven ventilation. This effect has been observed in Larga Cave, especially during the dry winter season (Vieten et al., 2016). Cave air ventilation may provoke evaporation of remaining water at the stalagmite's surface, which in turn influences the δ¹⁸O_f and δ²H_f values of the inclusion water seconds (Dreybrodt and Scholz, 2011; Deininger et al., 2012; Day and Henderson, 2011). At the same time, evaporation would also reduce the water content of the inclusions, explaining the systematically lower water content of cluster (2) samples (Sect. 3.1 and Fig. S2). Since speleothems are covered only by a thin film of water, evaporation can be very effective, especially given the high average cave temperature of Larga Cave and assuming a persistently lowered humidity. Calculation of the evolution of both water content and fluid inclusion δ¹⁸O values with a Rayleigh fractionation approach (Fig. S2) supports this interpretation and shows that a moderate evaporation of 50 %–75 % of the initial water content (between 0.5 and 3 µL g⁻¹) is sufficient to explain the observed enrichment of more than 12 ‰ of the fluid inclusion δ¹⁸O_f values.

Figure 4 shows that the fractionation ¹⁸ε_c/f between speleothem carbonate and fluid inclusion water is dominated by the enrichment of the fluid inclusion water, which means that the process resulting in the large fractionation of the δ¹⁸O_f values indeed only has a minor effect on the δ¹⁸O_c values of CaCO₃. When regarding both cluster (1) and cluster (2), ¹⁸ε_c/f is largely independent from δ¹⁸O_c values (Fig. 4c), while a very strong correlation with δ¹⁸O_f values is evident (Fig. 4b). This argues for the significant influence of evaporation, which would have a strong effect on the fluid inclusion water, but only a minor effect on the precipitated carbonate. The reason for the different effect of evaporation on the δ¹⁸O_f values of the fluid inclusion water and the precipitated calcite is related to the very different timescales of precipitation and oxygen isotope exchange between water and dissolved bicarbonate. Precipitation of calcite occurs on a timescale of a few hundred seconds (Dreybrodt and Scholz, 2011), depending – amongst other parameters – on cave temperature and film thickness. The timescale of oxygen isotope exchange between water and dissolved bicarbonate, in contrast, is in the range of hours to days (Beck et al., 2005; Dreybrodt and Scholz, 2011), and thus orders of magnitude longer.

When the drip water enters the cave, it is reasonable to assume that the δ¹⁸O value of the dissolved bicarbonate is in oxygen isotope equilibrium with the δ¹⁸O value of the drip water. When the drip water then falls onto the surface of a stalagmite, degassing of excess CO₂ occurs within a few seconds (Dreybrodt and Scholz, 2011), with a negligible effect on the δ¹⁸O value of the dissolved bicarbonate. Then, precipitation of CaCO₃ starts immediately, and the δ¹⁸O value of the bicarbonate that was previously equilibrated with the drip water is imprinted on the δ¹⁸O value of the precipitated calcite. We note that oxygen isotope fractionation during precipitation of speleothem calcite is a complex process (e.g. Dreybrodt and Scholz, 2011; Guo and Zhou, 2019; Hansen et al., 2019; Scholz et al., 2009), and oxygen isotope fractionation does not, in most cases, occur under conditions of oxygen isotope equilibrium. After rapid precipitation of speleothem calcite, the drip water remaining on the surface of the speleothem, which forms the fluid inclusion water when it is eventually trapped, is progressively affected by evaporation, which leads to increasing δ¹⁸O and δ²H values of the water, and finally to fluid inclusions that are not in isotope equilibrium with the host-calcite. However, due to the very long time constant of oxygen isotope exchange between water and bicarbonate, it takes hours to imprint the oxygen isotope signature of the (evaporated) water on the dissolved HCO${}_{\mathrm{3}}^{-}$, and hence, the precipitated calcite. In addition, since the Ca concentration of the drip water drops exponentially, only a minor component of calcite could still be precipitated from the remaining solution (Scholz et al., 2009; Mühlinghaus et al., 2007; Deininger et al., 2012; Dreybrodt and Scholz, 2011). Given the smaller amount of Ca typically left in the solution, even under evaporative conditions, only an insignificant amount of calcite may be formed from the solution affected by evaporation. Thus, in the case of evaporation inside the cave, the δ¹⁸O values of the precipitated calcite would still reflect the δ¹⁸O value of the original drip water, whereas the trapped fluid inclusion water may be strongly affected by evaporation and exhibit elevated δ¹⁸O values (Fig. 3). These processes have been modelled by Deininger et al. (2012), who indeed predicted an enrichment of the δ¹⁸O values of the water film on the surface of a stalagmite of several per mil due to evaporation, whereas the effect on the δ¹⁸O values of the precipitated calcite was below 1 ‰.

Ventilation in Larga Cave is mainly driven by temperature differences between the atmosphere and the cave (Vieten et al., 2016). During the dry winter season, over the course of the day, the part of the cave hosting the studied speleothem is warmer than the atmosphere, inducing a buoyancy-driven air flow. According to our results, evaporation effects are strongest during the comparably cooler and drier periods of the record. This in turn suggests that cave ventilation was particularly pronounced in these times, which would argue for a stronger seasonal temperature contrast during stadial intervals. This observation is consistent with the hypothesis of enhanced Northern Hemisphere seasonality during abrupt cooling events (Denton et al., 2005). Similarly, Ziegler et al. (2008) reported persistently high summer temperatures of the Atlantic warm pool also during North Atlantic cold events, such as the Heinrich stadials, suggesting that these seem to influence mainly the winter climate conditions in the Caribbean. It is therefore plausible that a pronounced and persistently strong seasonal temperature difference with particularly cool temperatures during the winter season enhanced winter type ventilation conditions in Larga Cave. Hence, ventilation, and thus, the occurrence of evaporation on the stalagmite's surface would have been significant.

4.1.3 Applicability of the oxygen isotope equilibrium thermometer

The previously discussed processes have implications for the applicability of the oxygen isotope equilibrium thermometer using fluid inclusion and carbonate δ¹⁸O values from PR-LA-1. Post-depositional (but pre-entrapment) enrichment of ¹⁸O in fluid inclusion water due to evaporation apparently results in high ¹⁸ε_c/f fractionation values between fluid inclusion water and calcite, because such influenced samples do not preserve the original δ¹⁸O_f signal. Samples plotting off the MWL are thus discarded from further inference. Note that evaporation effects would increase reconstructed temperatures by a factor of 2 to even 3; hence, temperature reconstruction is very sensitive to evaporation, which must be considered when interpreting meaningful temperature values as likely being maximum temperature estimates. Figure 4 shows that samples plotted on the GMWL in Fig. 3 might also show signs of evaporation. This is the case for the sample taken from HS 1 (15.9 ka). Here, other palaeoclimate records, as well as PR-LA-1 calcite δ¹⁸O values and elemental proxies, suggest that this phase was the coolest and driest in Puerto Rico and the Caribbean realm (Warken et al., 2020; Escobar et al., 2012; Grauel et al., 2016). The temperature calculated from the PR-LA-1 fluid inclusion value, however, is considerably warmer than at present day. As outlined above, we have to assume that evaporative effects also influenced speleothem and fluid inclusion formation during HS 1, shifting the δ¹⁸O_f to heavier (more positive) values. This demonstrates that samples may be influenced from evaporation despite a position on or close to the MWL. We therefore consider all calculated temperatures as a maximum estimate, and results less than +3 ^∘C greater than modern temperatures are considered, at first, plausible.

4.2 Reconstruction of hydro-climate and palaeo-temperatures in the western tropical Atlantic region

4.2.1 Fluid inclusion δ¹⁸O values suggest reduced convective activity during the last glacial

Despite their position close to the MWL, δ¹⁸O_f values in cluster (1) are significantly higher, by +2.6 ± 1.3 ‰, compared to modern drip water δ¹⁸O_w. At the time of deposition of speleothem PR-LA-1 (46.2 to 15.3 ka), the δ¹⁸O signature of the (sub)tropical Atlantic surface water was higher by about 1.5 ‰ due to larger global continental ice volume and increased salinity of the western tropical Atlantic (Hagen and Keigwin, 2002; Lea et al., 2003; Schmidt and Spero, 2011). Consequently, glacial drip water represented by fluid inclusion δ¹⁸O_f values must be corrected for this 1.5 ‰ offset from the present day. Millo et al. (2017) reported a similar enrichment of about 4 ‰ compared to modern drip water in last glacial fluid inclusion water in a subtropical stalagmite from Brazil, and a similar glacial–interglacial offset has been observed in δ¹⁸O values of stalagmites from Guatemala (Winter et al., 2020) and Florida (van Beynen et al., 2017).

We interpret the offset in cluster (1) as reflecting generally weaker convective activity during the last glacial period compared to modern conditions and associated reduced summer rainfall (e.g. Arienzo et al., 2019; Winter et al., 2020; Lases-Hernandez et al., 2019, 2020; Bird et al., 2020). Today, rainfall in Puerto Rico is isotopically enriched when it is formed by occasional northerly cold fronts and the orographic processes of trade winds (Scholl and Murphy, 2014). In contrast, isotopically lighter precipitation is associated with convective rain, such as hurricanes and tropical storm events. The type of precipitation, i.e. convective vs. orographic, depends to some extent on the meridional location of the ITCZ and related Caribbean SSTs (Winter et al., 2020). In this study, the mean of the reconstructed fluid inclusion temperatures for the LGM (ca. 25–18 ka) is 20.6 ± 3.2 ^∘C. Since our reconstructed temperatures can be regarded as upper limits, due to an underlying imperceptible amount of evaporation, this supports generally cooler temperatures than today over the course of the LGM. The observation of a glacial–interglacial temperature difference of 2–3 ^∘C highlights previous results of Caribbean SST records (Rama-Corredor et al., 2015; Lea et al., 2003; Schmidt et al., 2006b). Such low temperatures were presumably too cold to sustain strong convective activity during the LGM. Accordingly, Winter et al. (2020) suggested that deep atmospheric convection during the Holocene was triggered by exceeding an SST threshold of 27.5 ^∘C in the tropical ocean, leading to a shift of the hydro-climatic regime from a drier state and weak convection during the last glacial to a humid state and strong convection after 9 ka over Central America.

The mean value of Puerto Rican fluid inclusion temperatures obtained during MIS 3 (45.5–29 ka) is 21.7 ± 2.6 ^∘C. Despite the relatively large uncertainties, this may support a warmer MIS 3 than LGM, which has been found in the Guiana Basin (Rama-Corredor et al., 2015) and Brazil (Millo et al., 2017). For central Mexico, Caballero et al. (2019) report a mean annual temperature of MIS 3, comparable to the present day. From a global perspective, MIS 3 was ca. 2 ^∘C warmer than the LGM, due to a northern hemispheric summer insolation maximum and the resulting smaller ice sheets (Van Meerbeeck et al., 2009; Svendsen et al., 2004). According to Rama-Corredor et al. (2015), western tropical Atlantic SSTs were close to the threshold of 27.5 ^∘C during certain intervals of MIS 3 (Fig. 5), allowing sustainment of convective activity at least temporarily and/or over a longer time interval over the course of the year. Consequently, relatively wet environmental conditions – even though still drier than those of the Holocene – are supported by records from Mexico, Guadeloupe and the Bahamas (Lozano-García et al., 2015; Royer et al., 2017; Arienzo et al., 2017).

Figure 5 Palaeo-temperatures and the hydro-climate in the western tropical Atlantic region during the LGM and MIS2. From top to bottom: (a) Cariaco Basin reflectance record representing the mean meridional position of the ITCZ over northern South America (Deplazes et al., 2013); (b) speleothem PR-LA-1 δ¹⁸O_c values from Warken et al. (2020), representing the local high-resolution hydro-climate; (c) Guiana Basin Uk37 SST reconstruction (Rama-Corredor et al., 2015); (d) Ostracod δ¹⁸O-derived temperature record from Lake Petén Itzá (Guatemala; Escobar et al., 2012); (e) speleothem fluid inclusion-based temperature reconstruction from the Bahamas (Arienzo et al., 2015); (f) fluid inclusion-derived temperatures from PR-LA-1 (this study). The horizontal dashed line indicates modern cave air temperatures at Larga Cave. (g) Green triangles indicate the timing of Larga Cave ventilation events as suggested by fluid inclusion stable isotopic composition influenced by evaporation. In the tropical Atlantic, HS 1 to 4 and the Greenland stadials are associated with cold and dry conditions (light blue bars), and Greenland interstadials are associated with warmer and more humid climatic conditions (light red lines).

4.2.2 Millennial-scale climate variability in the western tropical Atlantic

On the millennial scale, fluid inclusion δ¹⁸O_f (and δ²H_f) values exert pronounced swings between very high values of up to 2 ‰ (30 ‰) during stadial phases and values down to 0 ‰ (15 ‰) during interstadial phases between 45 and 27 ka. The pattern largely follows the evolution of the high-resolution calcite stable isotope record (Fig. 2). Vertical blue bars in Figs. 2 and 5 illustrate the intervals of elevated δ¹⁸O_c values during the Heinrich stadials, which are associated with relatively cold and dry conditions in the circum-Caribbean area. The upper limits calculated for palaeo-temperatures are between 17.4 ± 1.5 and 22.8 ± 1.6 ^∘C for HS 2 and between 18.6 ± 1.9 and 19.8 ± 1.5 ^∘C for HS 3 (Fig. 5, Table S2). Accordingly, we estimate an minimum average cooling compared to present-day cave air temperature of −1.4 ± 2.8 ^∘C for HS 2 and −3.6 ± 2.2 ^∘C for HS 3. These values are consistent with the results of Arienzo et al. (2015), who found an average temperature decrease of 4 ^∘C across HS 1 to 3 for the Bahamas. In addition, many studies document that HS 1 exhibited the most pronounced cooling of all Heinrich stadials (Arienzo et al., 2015; Grauel et al., 2016). HS 1 was also more pronounced in the isotopic and elemental calcite signals of our sample compared to HS 2 or 3 (Warken et al., 2020). Escobar et al. (2012) and Grauel et al. (2016) found an even stronger cooling, down to 10 ^∘C for HS 1 in Lake Petén Itzá, Guatemala (Fig. 5). We note, however, that these studies used bottom dwelling ostracods, which may have recorded winter-biased lake temperatures. In addition, the strong maritime climate of Puerto Rico, located in the eastern Caribbean, makes less pronounced temperature variations appear reasonable. In addition to absolute temperature estimates, the samples which show signs of evaporation from cluster (2) suggest that during these times an increased seasonal temperature difference between persistently warm summer temperatures but cooler winter temperatures lead to enhanced cave ventilation. The timing of these samples, which includes not only the Heinrich stadials but also Greenland stadials and certain intervals of the LGM, is highlighted in Fig. 5. This further illustrates that our data set suggests that enhanced seasonality may have occurred during these intervals, which would be indicative of a persistent summer expansion of the Atlantic warm pool, but relatively cool winter temperatures.

Besides the temperature drops associated with Heinrich stadials, the expression of Greenland interstadials is still not well constrained in the American tropics (Warken et al., 2019). The mean temperature calculated from samples associated with interstadial events in speleothem PR-LA-1 is 22.8 ± 1.9 ^∘C (samples 15, 24, 28, 35), which is in the range of the present-day cave temperature. In comparison to the average stadial temperatures of 19.7 ± 0.9 ^∘C (samples 16, 18, 19), this indicates a difference of up to 3 ^∘C between stadial and interstadial phases in Puerto Rico. In the tropical Guiana Basin, at 7^∘ N latitude, abrupt glacial SST variations on the millennial scale are in the range of about ± 1 ^∘C (Rama-Corredor et al., 2015). Further north, Them Il et al. (2015) report a shift of ca. 2 ^∘C associated with stadial/interstadial transitions in the Florida Strait (24^∘ N), whereas shifts by 2 to 5 ^∘C have been observed in subtropical SSTs, such as at the onsets of GI8 and GI12 at the Bermuda Rise at 33^∘ N (Sachs and Lehman, 1999). Despite the large uncertainties associated with the oxygen isotope equilibrium thermometer, our results are thus consistent with these previous reconstructions and suggest that the interstadial phases in the western tropical Atlantic are also associated with relatively warm temperatures for Puerto Rico, and which may have been even comparable to the present day. Such warm temperatures favoured convective activity, and thus, higher precipitation amounts in the region. Marine reconstructions showed that interstadials were associated with reduced upper water column salinity, while stadial events were associated with increased salinity (Them Il et al., 2015). These observations support our results that meridional shifts in the position of the ITCZ in response to changes in the strength of the AMOC affected the evaporation–precipitation ratio in the western tropical Atlantic (Schmidt and Spero, 2011; Arienzo et al., 2017; Peterson et al., 2000; Schmidt et al., 2006a; Warken et al., 2020).

5 Conclusions

We demonstrate that the stable isotopic composition of fluid inclusions can be strongly influenced by evaporative effects. In Larga Cave, the temperature-driven ventilation potentially favoured such conditions, and our data suggest that this effect may have been especially pronounced during cold and dry conditions, such as during Heinrich and Greenland stadial phases and parts of the LGM. During warmer Greenland interstadials, most samples suggest that isotopic modification due to evaporation is negligible. In Larga Cave, this process is apparently climate related, making it possible to omit the potentially affected sub-samples for palaeo-temperature calculation.

Our data set thus provides valuable information about last glacial hydro-climate and temperature variability in the western tropical Atlantic and allows for constraining climatic changes on the millennial to orbital scale. For the presumably coldest phases of the last glacial, the Heinrich stadials, we observe a minimum cooling of −1.4 ± 2.8 ^∘C during HS 2 and −3.6 ± 2.2 ^∘C during HS 3. In contrast, during interstadial phases, temperatures were on average 23.4 ± 1.7 ^∘C, which is comparable to modern values. The analysis of stable water isotope values in speleothem fluid inclusions further suggests that, compared to the present day, the precipitation regime in the western tropical Atlantic was characterized by generally weaker convective activity during the glacial period as compared to the present day. Especially during the LGM, average temperatures were 2–3 ^∘C cooler, and thus too cold to sustain deep atmospheric convection. The intermediate background climate state during MIS 3 was potentially slightly warmer than the LGM, and favoured an increased sensitivity of western tropical Atlantic hydro-climate to perturbations of the AMOC and northern hemispheric temperature changes on millennial scales. In this period, the region experienced significant shifts between a slightly weaker – but still active – convective regime during the warmer interstadial phases and a dry/non-convective regime during the colder stadial phases.