Thermodynamics control precipitation in NE Mexico on orbital to millennial timescales

(150 word max) 26

significant discrepancies between climate model projections (Knutti & Sedláček, 2013;Lewis et al., 2010).Improving hydroclimate projections for Northern Mexico is critical given the substantial social, economic, and ecological impacts that shifts in mean precipitation or precipitation extremes can have in the region.For instance, severe droughts in the past have led to agriculture collapse (Liverman, 1990), national food shortages, and surges in international immigration (Johnson, 2011).Although paleoclimate records can contribute to improved understanding of the range and mechanisms of natural precipitation variability in this sensitive region, thus providing a critical test for climate models, few records are currently available.Speleothems are ideally suited for addressing this gap in the paleoclimate record due to their robust U-Th based age models and the multiple hydrologically sensitive proxies they contain.Despite the prevalence of limestone karst landscapes in Northeast (NE) Mexico, they have not yet been studied in this region.The modern climatology of NE Mexico is dominated by the Caribbean Low-Level Jet (CLLJ), which transports moisture from the Atlantic Ocean and Caribbean Sea to most of Mexico and Central America (Mesoamerica) during boreal summer (Méndez & Magaña, 2010a;Wang, 2007).On seasonal timescales, the strength of the CLLJ exhibits a maxima in July driven by increased solar heating and a northward migration of the Intertropical Convergent Zone (ITCZ) (Muñoz et al., 2008).During summer, the CLLJ feeds into the Great Plains Low-Level Jet, which plays a crucial role in modulating summer rainfall in the Great Plains region (Jiang et al., 2007).The CLLJ exhibits another maximum in February, driven by intensified meridional pressure gradient linked to heating over South America.It is important to note, however, only the CLLJ maxima in Boreal summer is associated with increased precipitation over Mexico.On longer timescales, this suggests that insolation variations dominated by orbital precession could impact regional hydroclimate in similar ways, with strengthening of the CLLJ and precipitation increases during Northern Hemisphere insolation maxima, and a strengthening of the CLLJ without a corresponding precipitation increase during Southern Hemisphere insolation maxima.Previous paleoclimate studies have invoked a strong role for insolation driven CLLJ variations in driving moisture variability across Mesoamerica, but most studies have focused on Southern Mexico and Northwest Mexico.The role insolation plays on longer timescales in NE Mexico is less clear.While a sediment record from the El Potosi Basin in NE Mexico demonstrates a strong correlation of runoff to NHSI through a strengthening of the CLLJ (Roy et al., 2016), other paleoclimate studies have found a limited role for NHSI and have suggested other seasons of insolation as more important drivers of hydroclimate (Quiroz-Jiménez et al., 2018;Roy et al., 2015Roy et al., , 2019Roy et al., , 2020)).In contrast, the role of insolation in NW Mexico and Southern Mexico is much better understood, with multiple records demonstrating a strong positive correlation between NHSI and precipitation via alteration of the North American Monsoon (Barron et al., 2012a;Kirby et al., 2006;Metcalfe et al., 2015;Roy et al., 2012Roy et al., , 2013) ) and northward shifts in the ITCZ (Lachniet et al., 2013).Unfortunately, records from NE Mexico do not span multiple orbital cycles, inhibiting our ability to fully constrain the response of precipitation to insolation, and limiting our understanding of the role of external forcing on future regional hydroclimate variability.On millennial timescales, existing paleoclimate records from NE Mexico have linked decreased precipitation during Heinrich Stadials (HS) to a weakening of the CLLJ.For instance, Roy et al. (2016) found decreased Ti concentration in lake sediments from El Potosi Basin, which were interpreted as reflecting reduced precipitation during HS 1.The dry HS1 conditions were assumed to reflect a weakening of the CLLJ as the ITCZ shifted southward in response to a shutdown in the Atlantic Meridional Overturning Circulation (AMOC).However, this interpretation may be inconsistent with modern dynamics of the CLLJ, which actually strengthens during Boreal winter (February) when the ITCZ migrates south and temperatures in northern South America warm (Cook & Vizy, 2010;Mestas-Nuñez et al., 2007).Using the seasonal ITCZ migration as an analogue, this suggests the CLLJ could theoretically strengthen during HS events in response to a southward ITCZ shift.This highlights a potential discrepancy, wherein a strengthened CLLJ may be associated with increased precipitation during interstadial events and decreased precipitation during stadial events, suggesting some other factor may play a more important role in driving regional hydroclimate of NE Mexico.While previous records have not shown a strong SST control on mean precipitation in NE Mexico on orbital or millennial timescales, we hypothesize that decreased CLLJ strength alone cannot explain drying during Heinrich Stadials Unlike most previous work reconstructing climate on millennial and orbital timescales in NE Mexico that have alluded to a strong control of CLLJ on precipitation, previous paleoclimate reconstructions and modeling studies focused on interannual to multi-decadal timescales have linked hydroclimate variations across Mesoamerica to changes in SSTs.For instance, elevated SSTs in the Gulf of Mexico (GOM) and North Atlantic are known to increase the intensity and frequency of tropical storms and hurricanes on seasonal timescales.Tree ring (Stahle et al., 2007(Stahle et al., , 2016;;Villanueva-Diaz et al., 2007) and climate model studies have also shown that both Pacific and Atlantic (Bhattacharya et al., 2017;Bhattacharya & Coats, 2020) SSTs exert a strong precipitation control across Mesoamerica on interannual to multidecadal timescales.Over the Common Era (last 2,000 years), changes in Atlantic SSTs are particularly associated with a strong, out-of-phase, dipole precipitation pattern between northern and southern Mesoamerica (Bhattacharya et al., 2017).The impact of SST variations on Mesoamerican hydroclimate patterns on millennial and orbital timescales is more poorly constrained, largely due to the paucity of records spanning long time periods.For instance, whether the north-south dipole response also dominates on longer timescales is still an open question.While records from southern Mesoamerica consistently show drying during HS events (Lachniet et al., 2013;Medina-Elizalde et al., 2017), lake sediment records from northern Mexico demonstrate a range of responses (Roy et al., 2012(Roy et al., , 2014(Roy et al., , 2020)).Although the inconsistency of the northern Mexico records may simply be driven by uncertainties in the age models (Roy et al., 2020), or the influence of tropical Pacific cyclones (Roy et al., 2014), the sparse paleoclimate record from NE Mexico hinders our ability to assess whether the dipole precipitation pattern dominates on longer timescales.More recent studies have begun to investigate the link between NE Mexico hydroclimate and SSTs on orbital and millennial timescales.For example, increased clay mineral concentration (Al+Si+K+Fe/Ca), interpreted to reflect increased watershed erosion from high intensity rainfall, in lake sediments from the Cieneguilla Basin and the Sandia Basin in NE Mexico have linked periods of increased tropical storm and hurricanes to warm Gulf of Mexico SSTs during the mid-Holocene and Bolling-Allerod (Roy et al., 2019(Roy et al., , 2020)).However, notably these shifts in precipitation extremes were not clearly associated with shifts in mean precipitation.Overall, the correlation between GOM SSTs and NE Mexico precipitation on orbital and millennial timescales has been shown to be inconsistent (Roy et al., 2016).Despite some evidence for SST influence on NE Mexican hydroclimate, most previous paleoclimate studies have concluded that changes in CLLJ strength are more significant than SSTs in driving regional hydroclimate.Resolving this issue and improving our understanding of the dynamic and thermodynamic response of precipitation across Mesoamerica to external forcing and internal climate variability is incredibly important under future climate change (Lim et al., 2018).To evaluate the timing and underlying dynamics of NE Mexico hydroclimate variability more fully, we present a new decadal-resolution, multi-proxy (δ 18 O,  13 C, Mg/Ca) speleothem record from Tamaulipas, Mexico that spans 58.5 to 4.6 ka.Our results show strong hydrologic responses to key millennial-scale events including the Younger Dryas and Heinrich Stadials 1, 3, 4, and 5, and a muted response to NHSI.We utilize results of a freshwater-forcing experiment conducted with a state-of-the-art climate model to investigate the relative importance of dynamic and thermodynamic controls on NE Mexico precipitation.Our results suggest a dominant role for cool SSTs in driving decreased moisture transport to the region during these millennial-scale cold events.This reduction in precipitation occurs despite the strengthening of the Caribbean Low-Level Jet, further highlighting the relative importance of thermodynamics.Finally, our model results combined with a comparison to other high-resolution precipitation records demonstrate that weakened AMOC leads to a spatially uniform drying across Mesoamerica, rather than the dipole response observed on shorter timescales.

Multiproxy reconstruction of hydroclimate variability in NE Mexico
We present a ~55,000 year record of hydroclimate utilizing oxygen isotopes ( 18 O), carbon isotopes ( 13 C) and trace elements (Mg/Ca) from speleothem sample CB2.CB2 is a 78 cm-long, candle-shaped, translucent grey and beige colored stalagmite composed of 100% calcite (Fig 2, see SI).CB2 was retrieved from Cueva Bonita (23N, 99W; 1071 m above sea level), located in the highlands of the Sierra Madre Oriental in the northeastern Mexico state of Tamaulipas (Fig. 1; see SI).The climate of this region is characterized by cool-dry winters and warm-wet summers (Wang & Lee, 2007) (See SI, Fig. S1), with a precipitation maximum in summer (July) driven by warm North Atlantic SSTs and an intensification of the Caribbean Low-Level Jet (Wang & Lee, 2007) (Fig. 1).The stable isotope and trace element speleothem records are tied to a U-series agedepth model, produced by 33 230 Th-234 U ages and the mean of 2000 Monte-Carlo simulations using the age-modeling software COPRA (Breitenbach et al., 2012).The age model indicates the sample formed continuously from 58.5 ± 3.8 ka to 4.6 ± 2.4 ka, with a relatively constant growth rate and an average temporal resolution of 34 years (Fig. 2).This represents the highest resolution continuous paleoclimate proxy record in Mexico over this time period.Previous work in Southern Mexico has consistently interpreted  18 Op to be reflective of precipitation amount (Frappier et al., 2014;Lachniet et al., 2013Lachniet et al., , 2017;;Medina-Elizalde et al., 2016), with greater amounts of rainfall associated with more negative  18 Op values.Risi et al. (2008) argues the amount effect dominates in the tropics due to high rainfall rates, which limits isotopic exchange with near-surface moisture.However in nearby Texas,  18 Op has been interpreted to reflect shifting moisture source, temperature, seasonality and shifts in thunderstorm size and duration (Maupin et al., 2021;Wong et al., 2015).To improve our understanding of modern precipitation isotope systematics, we have used an array of modeling and observational data.We analyzed moisture bearing air trajectories over a 15-year period (see SI) which demonstrates moisture is consistently sourced from the Gulf of Mexico and Caribbean Sea.An observational record of precipitation  18 O from approximately 1 km from the cave, established as part of this study, suggests the  18 O of precipitation is strongly dependent on precipitation amount (p < 0.01; r 2 = 0.88; slope = 2‰) (Fig. S2).This correlation is further supported by isotopeenabled GCM simulations of precipitation over the last 40 years, which suggest lower  18 Op values primarily reflect an increase in precipitation amount (Fig. S2).We therefore interpret  18 Op to be reflective of precipitation amount in NE Mexico and, due to stable cave conditions recorded over a 3-year period (Fig. S3), suggest this signal is likely also reflected, at least in modern, speleothem  18 O.
While  18 O is often a proxy for previously mentioned larger scale atmospheric processes, the controls of  13 C are often more localized and record changes in overlying vegetation (amount and type), soil respiration and CO2 degassing (Fohlmeister et al., 2020).We suggest the most likely driver of  13 C in CB2 is prior calcite precipitation (PCP), which occurs when there is reduced local water balance and is the result of enhanced CO2 degassing in the epikarst (Johnson et al., 2006).As a further constraint on the mechanisms of CB2  13 C variability, we conducted Mg/Ca analyses to test for covariation of  13 C and Mg/Ca, which is a robust indicator of PCP.In addition to the preferential loss of 12 C from the drip water dissolved inorganic carbon pool during CO2 degassing, PCP leads to the preferential uptake of Ca 2+ leaving the remaining drip waters, and the speleothem, enriched in trace elements (Mg 2 + ).CB2  13 C and Mg/Ca ratios strongly moderately correlate on orbital timescales (r = 0.31) and strongly correlate on millennial timescales, especially during Heinrich Stadials (r = 0.69), suggesting they are both dominantly controlled by PCP.The CB2  18 O record is dominated by millennial scale variations during the glacial and deglacial periods, and there is a clear ~1‰ decrease from the last glacial maximum to the Holocene (after correcting for global ice volume).In contrast to speleothem records from Asia and South America (Cheng et al., 2012(Cheng et al., , 2016)), the CB2  18 O record shows no clear precessional signal on orbital timescales, likely reinforcing the notion that this location is not impacted by monsoon dynamics.While increased precipitation in the early-Holocene has been recorded elsewhere in Mesoamerica (Lachniet et al., 2013;Warken et al., 2019), the effects of changing glacial-interglacial cave temperatures could easily explain the 1‰ glacial-interglacial shift observed in our  18 O record (see SI).The lack or orbital variability provides a relatively stable and un-varying background to which large amplitude millennial scale variability is superimposed on.The  18 O time series consistently exhibits large positive excursions during Heinrich Stadials, indicating shifts towards drier conditions.Speleothem  18 O values (n = 1578) range from -1.29‰ to -6.30‰ (VPDB) with the most 18 O-enriched samples occurring during Heinrich Stadial 1 (HS 1) at ~16.8 ka (Fig. 2).However, excursions toward higher  18 O values are also noted during 50-47 ka, 43-42 ka, 31-28 ka, 18-15 ka, and 12-10 ka, corresponding to HS 5, 4, 3, 1 and the Younger Dryas (YD), respectively (Fig. 2).While the timing of HS 4 in our record is slightly older than the expected age (~40-38 ka) (Hemming, 2004), the offset can be attributed to a relatively large uncertainty in our age model around this time.The CB2  13 C data co-varies with  18 O (Pearson's r = 0.52) and is similarly dominated by large amplitude millennial variations during the glacial, and a ~3‰ decrease across the deglaciation (Fig. 2).A lake sediment record from the region suggests the decrease in  13 C could reflect a shift from C4 to C3 dominant vegetation (Roy et al., 2019), but this shift could also potentially reflect some combination of decreased PCP, increased soil respiration or vegetation intensity, and/or decreased water-rock interaction during the Holocene (Fohlmeister, 2020).Millennial-scale shifts in  18 O are reproduced in the  13 C record, consistent with decreased local water balance during the Younger Dryas and Heinrich Stadials.Increases in  13 C were as large as 3.94‰ during HS 5.However, not all changes in  13 C were as dramatic, such as during HS 3 the shift in  13 C was as subtle as 1.04‰.We suggest the varying responses recorded by CB2 proxies may reflect real differences between individual Heinrich Stadials (Lynch-Stieglitz et al., 2014).The CB2 Mg/Ca values exhibit similar variations on millennial timescales as the stable isotopes, but diverge from them during the glacial-interglacial transition, with an increase from mean glacial concentrations of 30 mmol/mol to Holocene concentrations of ~40 mmol/mol (Fig. 2).We suggest this could potentially reflect the influence of temperatures on Mg partitioning into calcite (See SI; Stoll et al., 2012).During the glacial period, Mg/Ca data is consistent with a PCP control, especially during HS 3, 4 and 5, as evidenced by a significant positive correlation between  13 C and Mg/Ca (Pearson's r = 0.8, r = 0.61, r = 0.66, respectively).We therefore interpret higher Mg/Ca ratios and enriched 13 C to reflect enhanced PCP during periods of reduced local water balance (see SI).It is important to note, although local water balance is heavily influenced by precipitation amount, it is also influenced by other factors such as temperature and relative humidity which influence the amount of evapotranspiration.Discussion SST forcing of precipitation on orbital timescales While orbital-driven variations in insolation have been invoked to explain widespread moisture variations in the broader region of Mesoamerica (Lachniet et al., 2013), as well in NE Mexico (Roy et al., 2016), most of these records only span one precession cycle.The CB2 record, which spans ~2.5 precession cycles and extends ~25,000 years beyond the oldest lake record from NE Mexico (Roy et al., 2020), thus offers a unique opportunity to further constrain the impacts of insolation on precipitation and local water balance.While summer insolation has been proposed to drive precipitation in NE Mexico via a northward shift in the ITCZ and a strengthening of the CLLJ (Roy et al., 2016), and in NW Mexico via an intensified NAM (Barron et al., 2012b;Metcalfe et al., 2015), our record does not demonstrate a strong correlation to summer insolation.However, our record is not alone in that a growing number of records across Mesoamerica have also found a weak correlation to NHSI, suggesting the influence of autumn, winter or spring insolation as the dominant driver of hydroclimate variability on orbital timescales.For instance, Roy et al., (2015Roy et al., ( , 2019Roy et al., ( , 2020) ) has attributed a strong correlation between increased watershed erosion and/or runoff in northern Mexico to autumn insolation through increased tropical cyclone and hurricane activity.Furthermore, water scarcity in North Central Mexico recorded by increased authigenic calcite precipitation in a sediment core from ephemeral Lake Santiaguillo, has been linked to peaks in spring insolation (Quiroz-Jiménez et al., 2018).Even winter insolation has been linked to an extended wet season in speleothem  18 O records from Santo Tomás Cave in Cuba and Terciopelo Cave in Costa Rica (Lachniet et al., 2009;Warken et al., 2019).However, comparison of CB2  18 O to Northern Hemisphere insolation from different seasons demonstrates consistently weak correlations over the last ~55 ka (Fig. S4).While CB2  18 O seemingly changes with insolation over the last 20,000 years (Fig. 3, Fig. S4), exhibiting an in-phase response with autumn insolation or a lagged response to summer insolation, this relationship doesn't continue throughout the late-Pleistocene.This lack of a consistent insolation pattern suggests that other boundary conditions such as global ice volume, atmospheric pCO2, or SSTs may play a more direct role in explaining glacial-interglacial hydroclimate variability in NE Mexico than insolation.The CB2  18 O time series shows much more similarity with regional and global temperature records, indicating that precipitation at our study site may be more sensitive to thermodynamic controls on orbital timescales.Comparison with the Greenland ice core temperature record (Andersen et al., 2004; Fig. 3), shows the CB2 record is dominated by relatively cool and/or dry glacial conditions, as evidenced by relatively positive  18 O values from 58.5 to 20 ka, with a shift towards warmer and/or wetter conditions during the deglacial and Holocene.Superimposed on this orbital-scale trend are millennial scale shifts towards more positive  18 O values, indicating even drier conditions during Heinrich Stadials and the Younger Dryas.These trends are reproduced by the CB2  13 C and Mg/Ca records.Compared with insolation, CB2  18 O exhibits a much better correlation (r = -0.59 to -0.65) with regional SSTs, including the Tropical North Atlantic, the Gulf of Mexico, the Caribbean Sea, and the Eastern Equatorial Pacific (Fig. 3c; (Dubois et al., 2011;Lea et al., 2003;Schmidt et al., 2004;Ziegler et al., 2008).Similarly, the CB2 record also exhibits a closer relationship with atmospheric pCO2 than with insolation (r = -0.58,p < 0.05, Fig. 3b).Other records from the region have also shown a strong correlation to changes in local SSTs.For instance, increased rainfall in SE USA indicated by increased Mississippi River discharge into the Gulf of Mexico has been attributed to warmer SSTs (Tripsanas et al., 2013).Also, a stalagmite record from Cuba has attributed variability in  18 O to be driven by Gulf of Mexico and Caribbean Sea SSTs on glacial-interglacial timescales (Warken et al., 2019).Even lake sediment records from NE Mexico have identified SSTs as an important driver of climate on orbital timescales, but this record did not maintain a strong correlation to SSTs into the Holocene, consequently, more emphasis was placed on insolation, shifts in the ITCZ, and a stronger CLLJ (Roy et al., 2016).The CB2 record does not support a direct insolation control on NE Mexico precipitation, but rather suggests that glacial-interglacial variations in sea surface temperatures, indirectly reflecting orbital forcing and associated feedbacks, are more important.This suggests SST variations may be the dominant control on long timescales and future water availability may be more sensitive to SST changes than local atmospheric circulation patterns.Millennial-scale droughts linked to strengthened CLLJ and lower SSTs While previous records have suggested a variable response of Northern Mexico precipitation to millennial-scale AMOC changes (Quiroz-Jimenez et al., 2017;Roy et al., 2020), CB2 proxies consistently record drying in northeast Mexico during Heinrich Stadials, with the exception of HS 2 (Fig. 2).Roy et al. (2016) proposed that a southward shift of the ITCZ drove a weakening of the CLLJ during Heinrich Events, thus leading to dry conditions in NE Mexico.While the seasonal ITCZ migration may not be a perfect analogue for ITCZ shifts during Heinrich Stadials, the CLLJ also strengthens during boreal winter when the ITCZ is positioned farther south, highlighting the need for a better dynamical understanding of hydroclimate variations during Heinrich Stadials.To evaluate the underlying climate dynamics associated with drying during Heinrich Stadials, we analyzed results of an isotope-enabled Earth System Model simulation (iCESM1; Brady et al., 2019) forced with freshwater added to the North Atlantic (see methods).Model results show 10C (or more) of cooling across the Atlantic Basin (Fig. 4a).Cooling is simultaneously driven by a decrease in the warm western Gulf Stream and the positive wind-evaporation-SST feedback via an enhanced Bermuda High (Xie & Philander, 1994).This results in an intensification of the easterlies which funnel into the Caribbean Sea to feed the CLLJ (Fig. 4c).As the jet intensifies, it bifurcates across Mexico, increasing the magnitude of both the southern branch of the CLLJ across the Isthmus of Mexico and the southeasterly branch towards the United States (Fig. 4c).However, unlike the seasonal intensification of the CLLJ leading to increased rainfall, the resulting wind anomalies produce the opposite response.Strong winds combine with considerable SST cooling to reduce the vertically integrated moisture flux and precipitation over much of Mexico (Fig. 4b and 4c).Model results suggest moisture divergence can cause a 1 -2 mm/day reduction in precipitation (Fig. 4b) and a 2‰ increase in precipitation  18 O (Fig. 4d).This very closely matches the 1.5‰ increase we see in speleothem  18 O during HS 1 (Fig. 2).The freshwater added to the model experiment (0.5 Sv for 100 years) is most similar to the freshwater released during HS 1 (Hemming, 2004), upwards of 66% more than during other Heinrich Stadials and explains why we see a stronger drying during HS1 than during other stadials.Although through different mechanisms, the Younger Dryas is estimated to exhibit similar SST cooling as Heinrich Stadials, explaining the similar magnitude of drying indicated by our proxies.
It is important to note, the CB2 record suggests that changes in AMOC work in both directions, with increased precipitation leading to a -2‰ excursion in  18 O, a -3‰ shift in  13 C and, a decrease of 15 mmol/mol of Mg/Ca during the Bølling-Allerød when AMOC was strengthened and SSTs are warmer (Thiagarajan et al., 2014).In juxtaposition to previous work from Northern Mexico which has suggested a weakening of CLLJ in response to Heinrich Stadials and SST cooling, we demonstrate the jet strengthened during shutdowns in AMOC.These new results demonstrate precipitation in Northeast Mexico is potentially inversely correlated to jet strength via enhanced moisture divergence and a positive wind-evaporation-SST feedback, similar to processes controlling decreased rainfall in southern regions of Mesoamerica (Méndez & Magaña, 2010b).However, this interpretation opposes modern climatology which suggests a stronger summer CLLJ leads to increased moisture availability in Northern Mexico and regions of the United States (Wang, 2007).The contrasting response of a strengthening jet in modern climatology and paleoclimate records suggest the precipitation is not necessarily linked to the variability of the CLLJ.Changes in thermodynamics, particular over the ocean, appear to be much more important drivers of hydroclimate variability.The modeling results demonstrated here are consistent with the similarities between CB2  18 O and regional SSTs on orbital timescales.Our work also contradicts previous studies proposing that cooler SSTs drive an out-of-phase North-South precipitation pattern throughout Mesoamerica, at least on millennial timescales.We find CB2 exhibits similar patterns of variability to proxy records from farther south, including speleothem  18 O records in Costa Rica (Lachniet et al., 2009) and Cuba (Warken et al., 2019), Cariaco Basin sediment records (Deplazes et al., 2013), and magnetic susceptibility records from Guatemala (Escobar et al., 2012) (Fig. 5).This comparison to regional records suggests changes in regional SSTs and the CLLJ significantly and uniformly lead to drying across North and Central America and, will likely do so in the future.

Wet conditions in NE Mexico during the Last Glacial Maximum (or HS 2)
In contrast to the other Heinrich Stadials, the CB2  13 C and Mg/Ca records suggest an extended period of wetter conditions at our study site during HS 2 (~24 ka; Fig 2).This is in contrast to other proxy records from Lake Péten Itzá, Guatemala (Escobar et al., 2012;Hodell et al., 2008), Santo Tomás Cave, Cuba (Warken et al., 2019), and Terciopelo Cave, Costa Rica (Lachniet et al., 2009), which all show the expected response of drying during HS 2 (Fig. 5).While the CB2 response could potentially reflect a highly localized signal or be impacted by non-climatic proxy controls, the clear covariation between  13 C and Mg/Ca during this event does suggest it was characterized by increased water balance.There are three potential mechanisms we consider that could explain the wetter conditions in NE Mexico during HS 2, which occurs around the time of the Last Glacial Maximum: 1) Increased winter precipitation derived from the Pacific winter storm track during the LGM, 2) A weaker HS 2 event which NE Mexico hydroclimate is more sensitive to than other regions, and/or 3) increased water balance due to colder temperatures during HS 2 and the LGM.An intensified Pacific winter storm track during the Last Glacial Maximum has been frequently evoked to explain increased precipitation across Southern USA and Northern Mexico, but the spatial extent of the winter storm track has been previously disputed (Bradbury, 1997;Metcalfe et al., 2000).Speleothem  18 O demonstrates no change in precipitation amount during HS 2 while both  13 C and Mg/Ca ratios (Fig. 2) suggest anomalous increases in local water balance, compared to earlier glacial conditions.Furthermore, precipitation sourced from the Pacific would have an isotopic composition heavily depleted in 18 O due to the moisture's distance and terrain traversed (Dansgaard, 1964), which is not observed in CB2  18 O during this time period.Additionally, some of the best performing PMIP3 models (Chevalier et al., 2017) show an increase in the magnitude of the CLLJ (Fig. S5) but no increase in winter precipitation during the LGM (Fig. S6).Therefore, we suggest there was no contribution of rainfall from an enhanced winter storm track in NE Mexico during HS 2 or the LGM.There is still some debate about whether HS 2 resulted in a full AMOC shutdown.Higher Pa/Th ratios in North Atlantic sediment cores are traditionally interpreted to be reflective of a weakened AMOC (Lynch-Stieglitz et al., 2014).However, increases in opal flux as noted during HS 2 and HS 3 could also contribute to higher Pa/Th ratios (Lynch-Stieglitz et al., 2014).Additionally, oxygen isotopes from benthic foraminifera in cores from the Caribbean suggest changes in oceanic circulation during HS 2 were lower in magnitude compared to those during HS 1 and the YD (Lynch-Stieglitz et al., 2014).An incomplete shutdown in the western gulf stream, as potentially indicated from Caribbean benthic foraminifera, would have helped mitigate SST cooling, weakening the precipitation response in NE Mexico.This idea is further supported from SSTs from the Gulf of Mexico, the North Atlantic and the Caribbean, that show elevated temperatures compared to other Heinrich Stadials.We suggest our site is more sensitive to changes in Heinrich Stadial strength than other records of Mesoamerica because Cueva Bonita is located at the divergence of predicted precipitation  18 O change (Fig. 4d), where weaker Heinrich Stadials could result in no change, or even a slight positive excursion in precipitation  18 O.Therefore, the lack of response of CB2  18 O to HS 2 in our proxy record lends further, concrete support to the possibility that HS 2 may have been weaker than other Heinrich Stadials.Finally, we suggest wetter conditions during HS 2 and the LGM are primarily driven by decreased temperature at our study site, with no change in precipitation amount.While this is evident in the response of CB2 proxies, it is also supported by PMIP3 data.For instance, increased local water balance is observed in elevated soil moisture content in the LGM, compared to both the Mid-Holocene and Pre-Industrial period (Fig. S7).This is likely driven by a reduction in evapotranspiration, linked to increased cloudiness (Bush et al., 2009), reduced land temperatures (Lachniet & Vazquez-Selem, 2005), or changes in relative humidity over land (Bush & Philander, 1999).While evaluating the exact mechanism of increased local water balance is beyond the scope of this paper, our record is consistent with PMIP3 data and recent modeling studies (Lowry & Morrill, 2019) demonstrating high lake stands in Mexico and Central America (Escobar et al., 2012) may have been driven by decreased ET.

Conclusion
We have presented the first multi-proxy speleothem record from NE Mexico, that highlights millennial and orbital scale hydroclimate variability from 4.6 to 58.5 ka.In contrast to other speleothem records from tropical and monsoon regions, we find no strong evidence for a direct insolation control on regional hydroclimate.We show instead that glacial-interglacial variations are closely linked to Atlantic SST variations on orbital timescales.We find dry conditions in NE Mexico during the Younger Dryas and HS 1, 3, 4, and 5, and use climate model simulations to demonstrate that this response is driven by cool Atlantic and Caribbean SSTs rather than by a weakening of the CLLJ, as previously thought.Notably, we find evidence for increased local water balance during HS 2 and the LGM, which likely reflects a combination of a relatively weaker Heinrich Stadial 2 and reduced evapotranspiration at our study site.
Our record demonstrates that precipitation in NE Mexico is more sensitive to changes in SSTs than to the strength of the CLLJ, highlighting the thermodynamic response may be more important than the dynamic response in driving precipitation change in response to external forcing and internal climate variability.Furthermore, through comparison with other records and with climate model simulations, we find a uniform drying across Mesoamerica in response to cool SSTs during Heinrich stadials, suggesting that all of Mexico may experience similar hydroclimate shifts in response to future climate forcing.

Chronology
The CB2 stalagmite was cut, polished and sampled for 33 U-Th dates at 2.5 cm intervals along its vertical growth axis using a Dremel hand drill with a diamond dental bur.The CB2 sample has uranium concentrations ranging from 18 to 63 ng/g (Table S1).Calcite powder samples weighing 250-300 mg were prepared at Massachusetts Institute of Technology following methods similar to Edwards et al., (1987) (Lawrence Edwards et al., 1987).Powders were dissolved in nitric acid and spiked with a 229 Th -233 U -236 U tracer, followed by isolation of U and Th by iron co-precipitation and elution in columns with AG1-X8 resin.The isolated U and Th fractions were analyzed using a Nu Plasma II-ES multi-collector inductively coupled plasma mass spectrometer (MC-ICP-MS) equipped with an Aridus 2 desolvating nebulizer, following methods described in Burns et al. (2016) (Burns et al., 2016).The corrected ages were calculated using an initial 230 Th/ 232 Th value of 10.5 ± 2 ppm to correct for detrital 230 Th.This value and its uncertainty were determined by testing dates corrected with different initial 230 Th corrections for stratigraphic order.U-Th ages range from 5550 ± 1400 to 58200 ± 1600 years before present (where present is 1950 CE), and all 33 dates fall in stratigraphic order within 2σ uncertainty (Table S1).The 95% confidence interval for the age-depth model was constructed using 2000 Monte-Carlo simulations through the age-depth modeling software COPRA (Breitenbach et al., 2012).

Precipitation  18 O and Cave Monitoring
The closest precipitation stations of the International Atomic Energy Agency (Veracruz and San Salvador) are over 600 km from our field site and likely do not represent local patterns at the cave.Therefore, we established a local precipitation collection station directly at our field site (Alta Cima).In order to reduce kinetic effects due to evaporation, evaporation-limiting precipitation collectors were built and deployed during fieldwork (Fig. S2).Samples were collected after rainfall events in 1.5 ml glass vials with conical inserts, sealed with parafilm and refrigerated to limit evaporation.In total, 48 samples were collected from June 2018 to May 2019 and analyzed for D and  18 O using a Picarro L2130i cavity ringdown spectrometer at Chapman University (Table S2).The long-term standard deviation of an independent quality control sample is 0.51 ‰  2 H and 0.11 ‰  18 O (VSMOW).
Moisture source analysis at Cueva Bonita was conducted using air parcel back-trajectory simulations using the NOAA HYSPIT model (Stein et al., 2015) (Fig. S2).Air parcel trajectories were launched every 6 hours at an elevation of 1500m between 2005 and 2018 using GDAS weekly data.Each trajectory evaluated the air parcel's location over the previous 72 hours from launch.In total, 3600 rain-bearing trajectories were produced for the combined summer (JJAS) and winter (DJFM) months, only rain-bearing trajectories were included in analysis (Fig. S2).

Stable Isotope and Trace Element Analysis
CB2 was micro-sampled for both stable isotope and trace element analyses using a Sherline micromill at 500 m increments to a depth of 1 mm, producing 1578 samples.The powder for CB2 was collected, weighed out to 40 -80 μg and analyzed on a Kiel IV Carbonate Preparation Device coupled to a Thermo Scientific Delta V-IRMS at the UC Irvine Center for Isotope Tracers in Earth Sciences (CITIES) following methods described by McCabe-Glynn et al. (2013) to determine  18 O and  13 C (McCabe-Glynn et al., 2013).Every 32 samples of unknown composition were analyzed with 14 standards which included a mix of NBS-18, IAEA-CO-1, and an in-house standard.The analytical precision for  18 O and  13 C is 0.08‰ and 0.05‰, respectively.Speleothem  18 O values were ice-volume corrected using mean ocean  18 O values (Fig. S8).1578 samples were analyzed (Table S3).
For trace element analysis, 20 -60 μg calcite powder samples were dissolved in 500 μL of a double distilled 2% nitric acid solution.The samples were analyzed using a Nu Instruments Attom High Resolution Inductively Coupled Plasma Mass Spectrometer (HR-ICP-MS) at the CITIES laboratory.Mg/Ca ratios were calculated from the intensity ratios using a bracketing technique with five standards of known concentration and an internal standard (Ge) added to all samples to correct for instrumental drift.Trace element analysis of CB2 serves to complement the interpretation of speleothem  18 O and  13 C; therefore, a lower-resolution (multi-decadal to centennial) analysis was conducted over the complete record by analyzing every other sample (789 total; Table S3).For plotting/aesthetic purposes, CB2 Mg/Ca,  18 O and  13 C were smoothed using a moving average.The pandas function DataFrame.rolling().mean()wasutilized to smooth the data, with the size of the moving window set to 2. The full data set reported in the supplementary materials is unsmoothed.

Earth System Model Simulations
We use a water isotope tracer enabled version of the Community Earth System Model, iCESM1 (Brady et al., 2019).Model physics are consistent with CESM1, which simulates present-day and historical climate change quite well (Hurrell et al., 2013).Here, we run a fullycoupled configuration of iCESM1 with 1.9 x 2.5° atmosphere (CAM5) and land (CLM4), and nominal 1° ocean (POP2) and sea ice (CICE4).The model tracks stable water isotopes of oxygen and hydrogen through all model components.Previous studies demonstrate that the water isotope tracer components of iCESM1 can accurately simulate the δ 18 O distribution of both present (Nusbaumer et al., 2017) and past climates (He et al., 2021).
We configured our 21 ka simulation with period-appropriate orbital parameters and greenhouse gas concentrations, as in the PMIP4 protocol (Kageyama et al., 2017), and ICE-6G ice sheets (Peltier et al., 2015).Initial isotopic distribution in the ocean comes from the GISS interpolated ocean δ 18 O dataset (LeGrande & Schmidt, 2006) with globally uniform enrichment of +1 ‰ for δ 18 O (Duplessy et al., 2002).All other ocean initial conditions come from a previously performed LGM simulation (Di Nezio et al., 2016).We run this model configuration for 500 years to reach quasi-equilibrium, with our analyses coming from the final 50 years of simulation.We then branch this simulation to explore the effects of melt water flux in the North Atlantic at the LGM.Starting from year-500 of the 21 ka simulation, we add 0.50 Sv of freshwater with a δ 18 O of -30‰ into the North Atlantic (50°N -70°N), sufficient to rapidly and substantially weaken AMOC (Zhu et al., 2017).After 100 years, we shutoff the freshwater flux and extend the simulation for another 50 years.Analyses come from the final 50 years of the 150year simulation.5) Costa Rica (Lachniet, 2009), additional ocean sediment cores from the (6) Caribbean Sea (Schmidt et al., 2004), and (7) Cariaco Basin (Deplazes et al., 2013), and a lake sediment core (8) from Guatemala (Escobar et al., 2012).et al., 2013) (Escobar et al., 2012), Juxtlahuaca (light green, Lachniet et al., 2013), Costa Rica (blue; Lachniet, 2009), Cuba (dark green; (Warken et al., 2019) and Lake Petén-Itzá (black; Escobar et al., 2012).

Figures
FiguresThis document includes the following figures: Fig. 1.Summer (JJAS) climatology and nearby paleoclimate records Fig. 2. Stalagmite CB2, Age-Depth Model and Mg/Ca,  18 O and  13 C results Fig. 3. Comparison of CB2  18 O to various forcings Fig. 4. Earth System Model simulations of Heinrich Stadials Fig. 5. Comparison of CB2 with other paleoclimate records

Fig
Fig. 2. Stalagmite CB2, Age-Depth Model and Mg/Ca,  18 O and  13 C results.A) Results of 1578 stable isotope and 789 trace element measurements. 18 O is top and blue,  13 C is central and in orange, Mg/Ca is in bottom and green.Dates with associated uncertainties are below Mg/Ca.Heinrich Stadials highlighted in light red.B) CB2 Age-Depth Model constructed using 2000 Monte-Carlo simulations via the age-depth modeling software COPRA and 33 U-Th ages.Uncertainty in age-depth model indicated by grey shading, uncertainty in U-Th ages indicated by red error bars.C) Sample CB2 after being cut and polished.

Fig. 4 .
Fig. 4. Earth System Model simulations of Heinrich Stadials.Summary dynamical changes associated with high latitude freshwater input under glacial boundary conditions.a) Temperature and Sea Level Pressure changes; b) Changes in precipitation; c) Changes in total column precipitable water and 850 mb winds; d) Changes in the stable oxygen isotope ratio of precipitation.