Persistent Link Between Caribbean Precipitation and Atlantic Ocean Circulation During the Last Glacial Revealed by a Speleothem Record From Puerto Rico

The sensitivity of tropical Atlantic precipitation patterns to the mean position of the Intertropical Convergence Zone (ITCZ) at different time scales is well‐known. However, recent research suggests a more complex behavior of the northern hemispheric tropical rain belt related to the ITCZ in the western tropical Atlantic. Here we present a precisely dated speleothem multi‐proxy record from a well‐monitored cave in Puerto Rico, covering the period between 46.2 and 15.3 ka. The stable isotope and trace element records document a pronounced response of regional rainfall to abrupt climatic excursions in the North Atlantic across the Last Glacial such as Heinrich stadials and Dansgaard/Oeschger events. The annual to multidecadal resolution of the proxy time series allows substructural investigations of the recorded events. Spectral analysis suggests that multidecadal to centennial variability persisted in the regional hydroclimate mainly during interstadial conditions but also during the Last Glacial Maximum. In particular, we observe a strong agreement between the speleothem proxy data and the strength of the Atlantic meridional overturning circulation, supporting a persistent link of oceanic forcing to regional precipitation. Comparison to other paleo‐precipitation records enables the reconstruction of past changes in position, strength, and extent of the ITCZ in the western tropical Atlantic in response to millennial‐ and orbital‐scale global climate change.


Introduction
Last Glacial precipitation reconstructions from the tropical Atlantic region have documented pronounced rainfall variability on centennial to orbital timescales (Deplazes et al., 2013;Hodell et al., 2008;Peterson et al., 2000). These observations are interpreted to reflect a variable position of the Intertropical Convergence Zone (ITCZ) in response to changes in cross-equatorial temperature gradients (e.g., Broccoli et al., 2006;Schmidt & Spero, 2011;Stríkis et al., 2015). In particular, pronounced millennial-scale alternations, such as Dansgaard/Oeschger (D/O) oscillations and Heinrich events, were associated with a distinct climate response of warm (interstadial) and cold (stadial) periods in the Northern Hemisphere, including the tropical Atlantic region (Arienzo et al., 2015;Deplazes et al., 2013;Warken et al., 2019). In the case of the prominent Heinrich events, the latitudinal displacements of the ITCZ resulted from a weakening of the Atlantic Meridional Overturning Circulation (AMOC) (Henry et al., 2016;Lynch-Stieglitz et al., 2014;McManus et al., 2004). A weaker AMOC coincides with a reduced northward heat transport, a higher cross-equatorial sea-surface temperature (SST) gradient and, consequently, drier conditions in the tropical Northern Hemisphere (Broccoli et al., 2006;Clark et al., 2001;Waelbroeck et al., 2018). However, the influence of ocean circulation on changes in Last Glacial tropical rainfall, for example, associated with D/O events, is still under debate (Burckel et al., 2015;Henry et al., 2016;Roberts & Hopcroft, 2020;Them Il et al., 2015;Waelbroeck et al., 2018).
In particular, the details of the spatio-temporal structure and evolution of the Atlantic ITCZ are still not well understood . Mounting evidence from both modeling studies and proxy records suggests that the simple explanation of a north-south shifting of the ITCZ is not adequate to explain the observed rainfall patterns, especially because the most northern shifts of the Atlantic ITCZ are linked to a complex set of forcings and feedbacks on different timescales (Asmerom et al., 2020;Oster et al., 2019;Roberts & Hopcroft, 2020;Singarayer et al., 2017). For instance, the magnitude and extent of the ITCZ's latitudinal movements, and the role of external forcing, such as insolation, in modulating this variability especially during pronounced and/or abrupt climate change events during glacial periods, such as Heinrich and D/O events, remain unclear.
Modern tropical Atlantic precipitation variability occurs on a range of time scales, including the interannual to multidecadal periods, where the latter have been associated with the Atlantic Multidecadal Oscillation (AMO) (Bhattacharya et al., 2017;Knight et al., 2006;Moreno-Chamarro et al., 2019). Speleothem stable isotope records have strongly contributed to the understanding of multidecadal to millennial precipitation dynamics in the western tropical Atlantic during the Holocene (e.g., Asmerom et al., 2020;Fensterer et al., 2012Fensterer et al., , 2013Medina-Elizalde et al., 2010). For the Last Glacial, however, the number of studies, whose temporal resolution allows for the investigation of sub-centennial variability, is very limited and provides an inconsistent picture. For instance, Medina-Elizalde et al. (2017) suggest that the Atlantic low-latitude regions respond to internal modes of climate variability on multidecadal to centennial scales similar to today regardless of the global climate state, whereas the analyses by Hertzberg et al. (2012) document that Atlantic multidecadal scale may only operate during the warmer climate periods. This has significant implications for modern and near-future climate variability.
Here we present a stable oxygen isotope record obtained from a speleothem from Larga Cave, Puerto Rico, spanning the time interval from 46.2 to 15.3 ka. Our interpretation of the speleothem δ 18 O values is supported by trace element (Mg/Ca, Sr/Ca, and Ba/Ca) and δ 13 C data, all of which have been shown to constitute a valuable tool to identify potential hydrological processes or even reconstructing past precipitation patterns in the western tropical Atlantic during major climate shifts up to sub-decadal resolution (Arienzo et al., 2017;Cruz et al., 2007;Warken et al., 2018).

Site Description
Larga Cave is located in the north-central karst region of Puerto Rico (18°19′N 66°48′W, Figure 1) at an elevation of 350 m a.s.l. As the region of the Greater Antilles is dominated by sea surface processes, ocean-atmosphere coupling is of great importance for the climate in the Caribbean, which can therefore be classified as predominantly maritime (Granger, 1985;Schellekens et al., 2004). According to the observational record of the closest meteorological station, Arecibo observatory (18°21′N 66°45′W, 323 m a.s.l, data from http://xmacis.rcc-acis.org/), mean annual air temperature from 1980 to 2020 ranges from approximately 20°C to 24°C, and annual rainfall amount varies between 600 and 2,700 mm. In the catchment area of the cave, the mean annual air temperature today is 22.5°C, and annual rainfall amount is high with a mean of 2,200 mm (Vieten, Warken, Winter, Schröder-Ritzrau, et al., 2018). The area above the cave is cov-ered by dense tropical forest with a thin soil cover that is nearly absent on the higher elevated and exposed locations.
The rainfall δ 18 O values at the site show a seasonal pattern with more negative (−3 to −5‰ VSMOW) values occurring during the wet summer season and higher values (−1 to −2‰ VSMOW) in the dry winter season. The interannual anticorrelation of rainfall δ 18 O values with precipitation amount in Puerto Rico suggests an isotopic "amount effect" of about −0.1‰ per 100 mm of rainfall amount (Figure 2), which is in agreement with other regional observations (e.g., Govender et al., 2013;Lases-Hernandez et al., 2019. The observed seasonal pattern of rainfall amount and the associated isotopic signatures reflects the different sources associated with different rainfall types. The dry season is dominated by orographic, low-altitude air masses from the subtropical North Atlantic, while the wet season is sourced from convective, high-altitude air masses from the Caribbean (Govender et al., 2013;Scholl et al., 2009).
Inside Larga Cave, cave air and various drip sites have been monitored since 2012, showing a constant temperature and relative humidity (rH) regime in the main passage with 22.5 ± 0.2°C and close to 100%, respectively. The cave air pCO 2 variability is dictated by the seasonal cave ventilation with a well-ventilated winter (1) Larga Cave (Puerto Rico, this study); (2) Abaco Island, Bahamas (Arienzo et al., 2015(Arienzo et al., , 2017; (3) Bermuda Rise (Böhm et al., 2015;Henry et al., 2016;Lippold et al., 2009;McManus et al., 2004); (4) Cariaco Basin (Deplazes et al., 2013;Peterson et al., 2000); (5) Bonaire Basin (Parker et al., 2015); (6) Terciopelo Cave, Costa Rica ); (7) Juxtlahuaca Cave, Mexico (Lachniet et al., 2013); (8) Santo Tomas Cave, Cuba ; (9) Florida Strait (Them Il et al., 2015). The yellow shaded area indicates the spatial range of the seasonal migration of the mean position of the ITCZ. The red arrow shows the warm surface current in the Caribbean Basin as a part of the AMOC. The blue arrows indicate the mean trajectories of the easterly trade winds and convective systems transporting moisture into the western tropical Atlantic. mode (low pCO 2 values) and a near-stagnant summer mode (high pCO 2 values) due to seasonal temperature differences inside and outside the cave (Vieten et al., 2016).
The seasonal variations in the δ 18 O values of rainwater are smoothed by the soil and karst system acting as a low-pass filter, resulting in a well-mixed seepage water reservoir and a transmission time of atmospheric signals of at least several months up to a few years Vieten, Warken, Winter, Schröder-Ritzrau, et al., 2018). Therefore, the δ 18 O values of drip waters in Larga Cave are relatively constant with an overall mean of −2.6 ± 0.2‰ VSMOW during the monitoring period (2012-2019), which corresponds to the weighted annual mean of the rainfall δ 18 O value of −2.5 ± 0.1‰ VSMOW . The elemental composition of the drip waters also lacks a clear seasonal cycle with rather stable values for the Sr/Ca (0.8-1.0 mmol/mol) and Mg/Ca ratios (20-30 mmol/mol) (Vieten, Warken, Winter, Schröder-Ritzrau, et al., 2018).

Stalagmite PR-LA-1
Speleothem PR-LA-1 was collected lying on the cave floor in the main passage (Figure 3), where the rock overburden is approximately 40-80 m. It has a total length of 1.85 m and was removed in five pieces (L1A to L1E, from top to bottom) with individual lengths between 0.25 and 0.55 m (Figure 3). Stalagmite PR-LA-1 consists of whitish and translucent calcite, which in most parts exhibits convex-shaped lamination. In the bottom part of subsample L1C, between 110 and 118 cm distance from top (dft), the lamination appears more irregular with a concave "dip" in the middle of the stalagmite (Figure 3). In segment L1E, the part closest to the base of the stalagmite, several brownish layers indicate detrital inclusions. Slabs were prepared along the growth axis of each segment, which were subsequently used for 230 Th/U-dating and stable isotope and trace element analyses.

230 Th/U Dating
Samples for 230 Th/U dating were cut along the growth axis using a band saw. Analyses of U and Th isotopes were conducted using an Nu Plasma MC-ICP-MS at the Max Planck Institute for Chemistry (MPIC), Mainz. Chemical preparation of the samples and analytical methods at the MPIC followed Yang et al. (2015) and Obert et al. (2016). All ages and activity ratios used in this study were calculated using the half-lives reported by Cheng et al. (2000) in order to preserve comparability with previous records discussed within this study. Ages are reported relative to CE 2013. Age uncertainties are quoted at the 2σ level and do not include half-life uncertainties.

Proxy Analysis
Dri water samples were taken from different modern drip sites. Samples of speleothem PR-LA-1 were drilled at a spatial resolution of 1 mm. All samples were analyzed at the University of Innsbruck using a ThermoFisher Delta V isotope ratio mass spectrometer equipped with a Gasbench II (Spötl, 2011). Raw data were calibrated against NBS19, and δ values are reported relative to Vienna Pee Dee Belemnite (VPDB) standard. Long-term precision of the δ 13 C and δ 18 O values, estimated as the 1σ-standard deviation of replicate analyses, is 0.06‰ and 0.08‰, respectively (Spötl, 2011;Spötl & Vennemann, 2003).
Element/calcium ratios were measured by laser ablation ICPMS at the MPIC following the procedure of Jochum et al. (2012) and Yang et al. (2015), using the high-resolution sector-field ICPMS Thermo Element2, combined with the New Wave UP-213 Nd:YAG laser ablation system. All analyzed isotopes were measured at low mass resolution (m/Δm~300) and are reported to be interference-free (Jochum et al., 2012). Line scans were performed along the growth axis of PR-LA-1 using a spot size of 110 μm and a scan speed of 10 μm/s, resulting in a spatial resolution of 7 μm per data point (scan time 0.7 s). To avoid potential surface contamination, the scan path was pre-ablated at a scan speed of 80 μm/s. In order to account for matrix effects, data reduction was carried out by using the algorithm of Mischel et al. (2017) by calculating the blank corrected count rates of the analyzed isotopes relative to the intensity of the analyzed isotope 43 Ca of the internal standard element Ca determined in the sample and the reference material. The silicate glass NIST SRM 612 was analyzed for external calibration of the trace element analyses using the reference values determined by Jochum et al. (2011). Element/calcium ratios are given in molar units.

Numerical Methods
Correlation analysis was performed with a test statistic based on Pearson's product moment correlation coefficient r (x, y) following a t distribution with length (x) − 2 degrees of freedom. Reported correlation coefficients are all significant at the 0.05 level. Element/Ca ratios were interpolated to the lower resolution of the stable isotope records for the calculation of the correlation coefficients. The interpretation of the proxy signals is supported by using I-STAL, a model for interpreting Mg/Ca, Sr/Ca, and Ba/Ca ratios in speleothems , and ISOLUTION 1.0, an isotope evolution model describing the stable oxygen (δ 18 O) and carbon (δ 13 C) isotope values of speleothems (Deininger & Scholz, 2019). Spectral analysis was performed using REDFIT (Schulz & Mudelsee, 2002), which fits a first-order autoregressive (AR1) process directly to unevenly spaced time series.

Chronology of PR-LA-1
In total, 74 230 Th/U ages were determined for PR-LA-1 (Figure 3 and supporting information Table S2), showing that this stalagmite covers the period from 46.2 to 15.3 ka with a growth interruption between 41.1 and 35.0 ka. The U content of the samples is about 0.1-0.3 μg/g, and for most samples, ( 230 Th/ 232 Th) activity ratios are >1,000. However, 15 samples are not in stratigraphic order, and most of these showed relatively low ( 230 Th/ 232 Th). The age inversions persist also for the ages corrected for initial Th contribution using the bulk earth ratio. Therefore, the detrital correction is replaced by an isochron method (Beck et al., 2001;Ludwig & Titterington, 1994), in order to generate a correction factor appropriate to the local environment. The calculated ( 230 Th/ 232 Th) ratio of the detritus is 23.7 ± 7.5 (Text S1 and Figure 3). The still remaining nine age inversions, which are limited to the section older than 33 ka, suggest a more complex hydrological system upstream and potentially postdepositional loss of U due to diagenetic effects (Text S1). These inversions are surrounded by overall clean ages, which, in their entirety, constrain a tight and very reliable chronology with average age uncertainties of 0.5-1.5%. Therefore, potentially unreliable ages due to the previously mentioned reasons are not included in the age model (Table S2). The final chronology was constructed using the algorithm StalAge (Scholz & Hoffmann, 2011), which ensures the propagation of the uncertainties of the corrected ages to the final age model. This leads to uncertainties in the age control of the proxy records in the order of 200-400 years before 30 ka and 90-200 years afterwards. In sections with a low density of reliable ages and/or significant detrital Th, mainly segments L1E (45.5 and 44 ka) and L1C (35.5 to 33.5 ka), the uncertainty increases to 500-1,000 years. The average growth rate is approximately 100 μm/a, though the age model suggests very fast rates of more than 1 mm/ a at 43.2-43.0, 41.5-41.1, 33.7-33.5, and 27.8-27.6 ka. Comparably low deposition rates down to 10-30 μm/a occurred at 43.2-41.3, 33.5-30.3, 28.9-27.8, and 17.5-15.3 ka (Figure 4).

Speleothem Climate Proxy Patterns
The mean δ 18 O value of recently precipitated calcite at different drip sites is −3.1 ± 0.1‰ VPDB, which indicates that recent calcite precipitation in the cave occurs near isotopic equilibrium with modern drip water (−2.6 ± 0.2‰ VSMOW) (Hansen et al., 2019;Tremaine et al., 2011). Figure 4 shows the speleothem δ 18 O and δ 13 C values as well as the Mg/Ca, Sr/Ca, and Ba/Ca records. In the earlier part of the record, from

Paleoceanography and Paleoclimatology
46.2 to about 17 ka, the three element ratios show a variability of roughly ±50% around the mean values, with an Mg/Ca baseline of 1.0 ± 0.5 mmol/mol and Sr/Ca (Ba/Ca) values of 0.06 ± 0.03 mmol/mol (0.08 ± 0.04 μmol/mol). After 17 ka, Mg/Ca, Sr/Ca, and Ba/Ca values strongly increase towards their highest values of the whole record. δ 18 O and δ 13 C values show a similar pattern with δ 18 O values around 0 ± 2 ‰ VPDB and reaching the highest values of 1-2‰ VPDB between 17 and 15.4 ka, whereas δ 13 C values vary around −9 ± 3‰ VPDB. On longer timescales, δ 13 C, δ 18 O, and Mg/Ca values follow the general pattern of the speleothem growth rate (Figure 4). Over the whole record between 15.3 and 46.2 ka, δ 18 O values are correlated with δ 13 C values (r δ13C/δ18O = 0.52) as well as Mg/Ca values (r Mg/δ18O = r Mg/δ13C = 0.40). Though Sr/Ca and Ba/Ca are also correlated (r Sr/Ba = 0.49), these elements share a weak relationship with Mg/Ca (r Mg/Sr = 0.32).

Spectral Analysis
For intervals with growth rates of up to 1 mm/a, multidecadal scales are resolved in the stable isotopic and elemental records. The power of sub-millennial variability of the proxies was assessed by multitaper spectral analysis during different boundary conditions across the record. To directly compare stadial versus interstadial conditions, several time slices were selected with sufficiently high and, in particular, relatively constant growth rate to allow for multidecadal resolution. This includes an interval during the Last Glacial Maximum The δ 18 O spectra exhibit maxima in multidecadal to centennial modes during the LGM and GIs ( Figure 5), with concentrations of variance above the 90% confidence levels at periods of about 15-25, 30-50, 60-90, 100-120, and 250-370 a. The spectral power variability of δ 13 C and Mg/Ca is generally weaker but also shows enhanced power at the multidecadal and centennial periods ( Figures S3 and S4). During most analyzed Greenland and Heinrich stadials, the multidecadal modes are less pronounced in all proxies.

Interpretation of the Climate Proxies
The overall positive correlation of δ 18 O and δ 13 C values with Mg/Ca suggests a strong hydrological control of these proxies, as found in other studies Cruz et al., 2007;Sinclair et al., 2012;Wassenburg et al., 2019). The weaker correlation of Sr/Ca and Ba/Ca with Mg/Ca indicates that additional processes modulate the incorporation of Sr and Ba in the speleothem, reducing the suitability of these elements as hydrological proxies.
The main passage of Larga Cave, where PR-LA-1 was collected, is characterized by a pronounced temperature-driven seasonal cave air pCO 2 cycle (Vieten et al., 2016) with a well-ventilated (i.e., low pCO 2 ) regime during the dry winter season. During intervals with weaker and/or shorter summer seasons, the expected higher grade of ventilation would lead to lower pCO 2 in the cave atmosphere and accordingly lower Ca equilibrium concentrations associated with higher values of both stable isotopes and Element/Ca ratios in drip water and speleothem calcite (Baker et al., 2014;Mattey et al., 2010). Since the convective summer rainfall dominates the annual recharge, phases with weaker and/or shorter summer seasons are also likely accompanied by a generally reduced recharge, longer residence times and enhanced water-rock interaction, which would additionally increase both stable isotopic and element/Ca ratios in the drip water (Dreybrodt & Scholz, 2011;Mattey et al., 2010;Sinclair, 2011). The reduced contribution of convective summer rainfall to the reservoir feeding the drip sites in Larga Cave additionally shifts the annual mean δ 18 O value of the infiltrating water towards higher values.
Additional process modulating the elemental and isotopic signatures in Larga Cave is intra-and interannual drip variability in response to recharge, which has been observed at modern drip sites Vieten, Warken, Winter, Schröder-Ritzrau, et al., 2018). Longer drip intervals may coincide with an increase in prior calcite precipitation (PCP) and result in higher Mg/Ca, Sr/Ca, and Ba/Ca ratios (Mattey et al., 2010;Stoll et al., 2012) as well as higher δ 18 O and δ 13 C values (Dreybrodt & Scholz, 2011;Hansen et al., 2019) of the drip water. A higher degree of PCP also decreases the saturation state of the drip water and therefore slows down stalagmite growth, reducing the partitioning of Sr and 10.1029/2020PA003944

Paleoceanography and Paleoclimatology
WARKEN ET AL.
Ba, which may explain weak or even absent correlations of these elements as also observed in stalagmite PR-LA-1 Warken et al., 2018). The weak correlation of Mg/Ca with Sr/Ca and Ba/Ca observed both in the speleothem and modern drip water (Vieten, Warken, Winter, Schröder-Ritzrau, et al., 2018) indicates only minor PCP in the epikarst. However, it cannot be excluded that PCP exerts an additional control on different timescales, which would exacerbate the impact of dry/wet conditions on Mg/Ca as well as δ 18 O and δ 13 C in the speleothem. Numerical simulations with I-STAL  using input parameters closely oriented on the observed variations during the cave monitoring and in the speleothem support the hydrological interpretation of the speleothem proxies (Text S2). In particular, the (1) false-alarm levels of 80% (purple), 90% (violet), 95% (dark blue), and 99% (light blue), respectively. For all shown spectra, the AR(1) passed the REDFIT runs test, which checks the equality of theoretical AR(1) and data spectrum.

Paleoceanography and Paleoclimatology
analyses show that the observed elemental variability may be simulated by drip rate changes and a variable effective cave air pCO 2 (i.e., the pCO 2 recorded by the stalagmite during the season of main calcite precipitation), which is likely linked to seasonally variable ventilation ( Figure S1). Therefore, we regard the combination of changes of drip rate and cave ventilation in response to weaker and/or shorter wet summer seasons as possible controls of speleothem Mg/Ca, δ 18 O, and δ 13 C values in Larga Cave, which may have acted on annual up to multi-centennial timescales.
The relative contribution of kinetic isotope fractionation effects to the observed speleothem δ 18 O and δ 13 C variations was numerically explored using the ISOLUTION model (Deininger & Scholz, 2019). The highest speleothem δ 18 O and δ 13 C values of about + 2‰ and −4‰ VPDB, respectively, are observed at 17-15.3 ka. The simulations reveal that only about half of the observed magnitude of the positive excursions in the stable isotope record can be explained by cave-internal isotope fractionation effects due to long drip intervals, cooler temperatures, and stronger evaporative ventilation effects ( Figure S2). Hence, the additional increase can only be explained by drip water δ 18 O values significantly higher than the modern values, which, in turn, supports that drip water, and therefore, speleothem δ 18 O values are directly related to changes of the isotopic composition of the infiltrating water.
Modern drip water δ 18 O values in Larga Cave predominantly represent the mean annual rainfall δ 18 O value Vieten, Warken, Winter, Schröder-Ritzrau, et al., 2018). In the tropical Atlantic region and also at the site of Larga Cave, δ 18 O values of precipitation are generally linked to summer rainfall amount (Figure 2), a relationship which is based on the type and source of wet season (convective) versus drier season (orographic) rainfall (Govender et al., 2013;Lases-Hernandez et al., 2020;Scholl et al., 2009). At the (former) drip site of PR-LA-1, this signature is then amplified by kinetic isotope fractionation effects in response to drip interval and ventilation effects, an effect which is also observed in Mg/Ca and δ 13 C values. Since both cave-internal and external effects point towards the same interpretation, the combination of local (Mg/Ca, δ 13 C) and regional (δ 18 O) hydroclimate proxies allows to interpret the PR-LA-1 record as a robust reconstruction of past (summer) rainfall amount and convective activity in the western tropical Atlantic.

Inferences on Western Tropical Atlantic Rainfall Patterns During the Last Glacial
Speleothem PR-LA-1 grew from 46.2 to 15.3 ka, during the second half of the Last Glacial period. Compared to modern calcite (δ 18 O = −3.1‰ VPDB, this study), the δ 18 O values are about 3-4‰ higher with an average level of about 0 ± 1‰ during the LGM. This suggests that Puerto Rico-in agreement with reconstructions from the Caribbean area-experienced a generally drier climate during glacial times and possibly a weaker convective summer regime (Correa-Metrio et al., 2012;Hodell et al., 2008;Peterson et al., 2000;Schmidt & Spero, 2011). Increased continental ice volume, lower tropical SSTs by 2-4°C, and salinity values account for a change in the δ 18 O value of seawater of up to 2‰ during glacial periods for the (sub)tropical Atlantic (Hagen & Keigwin, 2002;Lea et al., 2003;Schmidt & Spero, 2011). These effectively reduce the difference between modern calcite and speleothem PR-LA-1 δ 18 O values during the LGM to about 2‰. Accompanied by increased Mg/Ca and δ 13 C values, peak speleothem δ 18 O values reach 2-3‰ VPDB during the presumably coldest and driest intervals of the Last Glacial in the region, that is, the Heinrich stadials (Arienzo et al., 2015(Arienzo et al., , 2017. Arienzo et al. (2015) reconstructed a temperature decrease of about 4°C during these periods compared to recent times, which would only explain about 1‰ change in cave temperature-related isotope fractionation during carbonate precipitation. The ISOLUTION simulations show that disequilibrium isotope fractionation effects would lead to a maximum increase of 1-2‰ (section 4.1), still leaving a residual shift in speleothem δ 18 O values of about 1-2‰, which indicates both a drop in temperatures and decreasing convective activity accompanied with lower rainfall amounts (Arienzo et al., 2017;Escobar et al., 2012).
Compared to the modern Puerto Rican "amount effect," an increase of 1‰ in δ 18 O values roughly corresponds to a rainfall decrease of about 50% (Figure 2). However, a number of factors can affect the δ 18 O values of drip water and speleothems (Baker & Fritz, 2015;. For example, partially increased d-excess values of rainwater at the location of Larga Cave suggest enhanced moisture recycling, especially during the dry winter season . In addition, the modern isotopic amount effect also reflects the seasonally varying moisture sources and rainfall types, that is, convective activity (Govender et al., 2013;Scholl et al., 2009). During drier glacial periods, convective activity may have been weaker (Winter et al., 2020) and moisture recycling more dominant, which would both shift the weighted annual mean of the δ 18 O values of the drip water towards higher values. This would imply that the amount effect was less evident as it is today as suggested by the modern relationship. The difficulty of quantifying these effects therefore prohibits a direct, robust calibration of glacial speleothem δ 18 O values with rainfall amount.  (Böhm et al., 2015;Henry et al., 2016;Lippold et al., 2009Lippold et al., , 2019McManus et al., 2004). Vertical blue bars indicate intervals of high PR-LA-1 δ

O values indicative of cool/dry periods in Puerto
Rico coinciding with a weak AMOC, such as Heinrich stadials (HS) 1-5.

Millennial-Scale Spatio-Temporal Patterns
The PR-LA-1 δ 18 O record exhibits a pattern of millennial-scale variability reminiscent of D/O oscillations and Heinrich stadials recorded in Greenland ice cores, marine sediments, and other archives ( Figure 6). This structure is also visible in the δ 13 C and Mg/Ca values (Figure 4), indicating significant changes in regional precipitation amount and convective activity. In agreement with previous hydroclimate reconstructions from the tropical Atlantic realm (Arienzo et al., 2017;Winter et al., 2020) and SST reconstructions from the Caribbean and the Florida Strait (Parker et al., 2015;Them Il et al., 2015), lower (higher) speleothem δ 18 O values during interstadial (stadial) phases are associated with wetter and warmer (drier and cooler) conditions and greater (less) convective rainfall intensity ( Figure 6).
Among these events, Heinrich stadial 1 (HS1, 17 to 15.3 ka) stands out, where the speleothem proxies δ 18 O, δ 13 C, and Mg/Ca indicate a significant reduction of rainfall amounts along with strongly enhanced disequilibrium fractionation at the drip site (section 4.1). The onset of HS1 in PR-LA-1 occurred at about 17.5 ka with a gradual increase in δ 18 O values indicating a decrease of rainfall amount and temperature, interrupted by a short warm and wet period at about ∼17 ka (Figures 6 and 7). Comparing PR-LA-1 δ 18 O values with the speleothem records from Juxtlahuaca Cave, Mexico (Lachniet et al., 2013) and Abaco Island, Bahamas (Arienzo et al., 2017) shows that the timing of the increase in δ 18 O during the substages of HS1 was nearly synchronous (Figure 7). A three-stage structure of HS1 is also recorded by lacustrine sediments in Guatemala (Escobar et al., 2012), supporting its region-wide expression. The observation that HS1 was the severest stadial in the region is consistent with other reconstructions (Arienzo et al., 2015;Escobar et al., 2012) and is supported by reconstructions showing that HS1 was accompanied by a near shutdown of the AMOC and an abrupt and synchronous response of the ITCZ to changes in ocean circulation (Lachniet et al., 2013;Stríkis et al., 2015). In contrast, there were only muted changes of the AMOC during HS2 and HS3 (Böhm et al., 2015;Lynch-Stieglitz et al., 2014).
HS2 shows a different picture in Puerto Rico characterized by an abrupt, one-stage increase in δ 18 O between 24.3 and 23.8 ka (Figure 8), which is also clearly imprinted in the records from the Bahamas and Cuba (Arienzo et al., 2017;Warken et al., 2019). The records from the north-eastern Caribbean appear to follow the pattern of 231 Pa/ 230 Th on the Bermuda Rise more closely, indicating a temporary weakening of the AMOC and a synchronous reduction of regional rainfall. In contrast, the speleothem record from the Yucatan Peninsula (Medina-Elizalde et al., 2017) further west is characterized by a continuous precipitation reduction after 23.8 ka (Figure 8), synchronous and in-phase with precipitation records from South America. The difference in the climate response in the north-eastern Caribbean compared to reconstructions from further south-west indicates a diminishing influence of the Atlantic from the east towards the south-west, in contrast to the basin-wide impact of HS1.

10.1029/2020PA003944
Paleoceanography and Paleoclimatology may indeed reflect an expression of D/O events in the Larga record. These appear to be characterized by a hydrological regime similar as (or only slightly drier and/or cooler than) today when taking into account glacial background conditions in the tropical Atlantic source regions (Hagen & Keigwin, 2002;Schmidt & Spero, 2011).
The amplitude of the response of speleothem δ 18 O values to some Greenland stadial and interstadials appears to be attenuated in the northernmost records from Cuba and the Bahamas (Arienzo et al., 2017;Warken et al., 2019) but also towards Costa Rica ) compared to precipitation records in the Cariaco Basin (Deplazes et al., 2013). Similarly, the δ 18 O values of PR-LA-1 show decreasing amplitudes of centennial to millennial oscillations after 30 ka compared to NGRIP and Cariaco with weak or even absent excursions across GIs 2-5 ( Figure 6). This suggests a progressive weakening of the connection between Caribbean precipitation and North Atlantic forcing as previously suggested by other studies (Oster et al., 2019;Warken et al., 2019), highlighting the complex climatological patterns in this region and arguing for regionally diverse expansions/contractions of the zonal mean of the ITCZ combined with meridional shifts.

Persistent Link of Multidecadal to Millennial-Scale Rainfall to Ocean Circulation
A remarkable feature of the PR-LA-1 record is that speleothem δ 18 O values closely track changes in the reconstructed strength of the AMOC ( Figure 6). This similarity becomes even more obvious when taking into account the uncertainty of the chronology of the marine sediment records of about 0.5 ka, which adds to the uncertainty of the LA-1 age model of 0.1-0.4 ka. Within the uncertainty of the age models, major events of freshwater input into the North Atlantic and subsequent slowdown of ocean circulation, such as during Heinrich stadials 1-3, are recorded by elevated speleothem δ 18 O values (Figures 6 and 7), reflecting pronounced transitions to drier and cooler conditions in the western tropical Atlantic during HS1 (17.2-15.5 ka), HS2 (24.3-23.8 ka), and HS3 (30-29 ka) (Arienzo et al., 2015(Arienzo et al., , 2017Correa-Metrio et al., 2012;Escobar et al., 2012;Lachniet et al., 2013).
The Larga δ 18 O record strikingly agrees with sedimentary 231 Pa/ 230 Th from the Bermuda Rise as well as regional Mg/Ca-SST reconstructions (Parker et al., 2015;Them Il et al., 2015), not only during Heinrich stadials but also in the course of millennial to centennial variations reminiscent of D/O main and substages ( Figure 6). The role of the AMOC as the driver of D/O variability has recently also been corroborated by marine sediment cores from the Caribbean and the Florida Strait (Parker et al., 2015;Them Il et al., 2015). The agreement of the PR-LA-1 record with these marine sediment cores provides additional support for a strong and persistent link between AMOC variability and precipitation changes in the western tropical Atlantic. Hence, we regard AMOC strength as a dominant regional control of atmospheric circulation changes (Burckel et al., 2015;Henry et al., 2016;Waelbroeck et al., 2018). In agreement with recent modeling studies suggesting that asymmetric extratropical forcing, such as ice sheets and freshwater input into the North Atlantic, produce meridional shifts in the zonal mean rain belt, these observations indicate a southward shift of the ITCZ associated with abrupt decreases of AMOC strength (Broccoli et al., 2006;Clark et al., 2001;Parker et al., 2015;Singarayer et al., 2017;Them Il et al., 2015).
To explore the nature of decadal-to centennial-scale variability previously reported in the tropical Atlantic for certain intervals during a time when boundary conditions were very different from today (

10.1029/2020PA003944
Paleoceanography and Paleoclimatology proxy data to identify significant modes of variability. Spectral analysis of the speleothem PR-LA-1 δ 18 O data suggests relatively high power on decadal (20-30 a), multidecadal (40-90 a), and centennial periods (100-120 and 200-300 a), whereby the latter may be poorly defined due to the finite lengths of the individual sections. The decadal and interdecadal modes of variability have been also observed in the Cariaco Basin and represent periodicities that have no clear modern analogs. Due to the different boundary conditions from today, it is likely that there were modes of climate variability during glacial times that are not known from the modern climate system (Hertzberg et al., 2012).
The multidecadal variability in the speleothem proxies may be related to a process similar to the modern AMO, which has a period of ∼65-80 years (Enfield et al., 2001;Knight et al., 2005). The present-day AMO involves surface warming and cooling of the inter-American Seas and is strongly correlated with Atlantic hurricane variability and migrations of the ITCZ (del Monte-Luna et al., 2015;Enfield et al., 2001;Knight et al., 2006;Moreno-Chamarro et al., 2019). Because variations in the strength of the AMOC are linked to meridional redistribution of ocean heat within the North Atlantic, whose spatial SST pattern is superimposed on that of the AMO, it has been suggested that multidecadal variability in the North Atlantic is related to fluctuations of the AMOC (Knight et al., 2005;Latif et al., 2004;Oelsmann et al., 2020;Zhang & Zhang, 2015).
Instrumental, paleoclimate, and modeling data support a positive correlation between the AMO and hydroclimate over the western tropical Atlantic including Puerto Rico during the Late Holocene (Bhattacharya et al., 2017;Fensterer et al., 2012;Winter et al., 2011). Previous studies have demonstrated the presence of interannual and multidecadal variability also during the Last Glacial, such as the lacustrine record from Lake Petén Itzá (Correa-Metrio et al., 2012;Escobar et al., 2012;Hodell et al., 2008), a high-resolution speleothem δ 18 O record from the Yucatan Peninsula covering HS 2 (Medina-Elizalde et al., 2017), and foraminifera data from Cariaco Basin sediments across GI 12 (Hertzberg et al., 2012). In our record, the multidecadal variability is less pronounced or even absent during the colder and drier intervals, the Heinrich and most Greenland stadials. Stadial conditions represent a state, where the AMOC is relatively weak or even in an "off" mode (Böhm et al., 2015;McManus et al., 2004). This would provide support for the hypothesis that instabilities in the coupled ocean-atmosphere system are introduced when AMOC intensifies and that multidecadal climate fluctuations similar to the AMO operate primarily during warmer climate intervals, which would have significant implications for modern and near-future climate variability (Hertzberg et al., 2012). However, in order to test if this observation is indeed a persistent feature of the tropical climate system or just an artifact in our record, further reconstructions of a similar or even higher temporal resolution across different climate states would be required.
The presence of multidecadal, centennial, and millennial scale precipitation variability in the same proxy record further suggests a common mechanism of heat redistribution in an AMO-analogue mode on these timescales (Bhattacharya et al., 2017;Ting et al., 2009;Zhang & Zhang, 2015). The existence of multidecadal variability similar to the AMO may be a function of relatively warm climate states with a strong AMOC, while being dampened or even absent during colder, stadial, and near-stadial conditions and weak overturning circulation, although this idea requires further testing against higher-resolution proxies of the AMOC (Hertzberg et al., 2012). This argues for a dominant role of the strength of ocean circulation in tropical Atlantic precipitation and adds support to the hypothesis that the AMOC is a key element for the link between northern high latitudes and tropical climate.

Conclusions
The speleothem multi-proxy record from Larga Cave provides valuable insights into regional paleoclimate variability on multidecadal to orbital timescales, whereby δ 18 O values are the most sensitive proxy of past precipitation amount. This study highlights the close connection of western tropical Atlantic precipitation variability and Northern Hemisphere climate variability. In particular, our record underlines the central and persistent role of the strength of the AMOC in abrupt glacial climate change, adding valuable information to recent debates and underscoring that a more complete understanding of the underlying mechanisms can only be achieved by integrating observations and models (Roberts & Hopcroft, 2020;Them Il et al., 2015). The compilation of regional records from speleothems and other archives suggests the existence of spatially and temporally varying precipitation patterns during the Last Glacial period. These observations suggest that meridional shifts in ITCZ variability occurred along with expansions and contractions of the associated rain belt on different timescales depending on the nature of the forcing. Thus, the concept of north-south shifts of the ITCZ causing wet-dry conditions in the region is insufficient to explain the observed tropical precipitation variability and requires further in-depth research on the variability and drivers of the western tropical Atlantic rain belt.

Data Availability Statement
The data presented in this paper were uploaded to the open PANGAEA data library (https://doi.pangaea.de/ 10.1594/PANGAEA.915000) and will be available after acceptance of this paper. In order to ensure comparability with future studies, we also provide the ages and activity ratios calculated with the half-lives of Cheng et al. (2013) in Table S2.