Figure 1. Precipitation data from GPCC (Schneider et al., 2011), showing (a) Mean annual precipitation over Mexico, Central America, and the Circum-Caribbean region, (b) EOF1 of annual precipitation, (c) EOF1 of summer (JJAS) precipitation, and (d) EOF1 of winter (DJFM) precipitation. Panel (e) shows mean monthly precipitation near Cueva Bonita (http://clicom-mex.cicese.mx/).

3 of 14 the more recent past (Common Era) can provide a robust record of past hydroclimate variability to help clarify the role of Atlantic versus Pacific SSTs on regional precipitation, none have been published in this critical region.

To address this gap, we have developed the first continuous inter-annually resolved stalagmite record (CB4) of past hydroclimate spanning the last millennium utilizing four geochemical proxies: stable oxygen isotopes (δ 18 O), carbon isotopes (δ 13 C), trace elements (Mg/Ca), and dead carbon proportion (DCP, based on 14 C). Moreover, novel forced-SST climate modeling experiments are employed to confirm precipitation change measured in speleothem CB4. The CB4 sample was retrieved from Cueva Bonita (23°N, 99°W; 1,071 m above sea level) located in the northern-most tropical cloud forest on the windward side of the Sierra Madre Oriental in the Northeast state of Tamaulipas (Figures 1 and S1 in Supporting Information S1). The climate of NE Mexico is characterized by a warm wet summer and a cool dry winter (JJAS; Figure 1e). The stalagmite age model is constrained by 19 U-Th dates and fluorescent annual lamina counting (Figure S2 in Supporting Information S1), and extends from 833 CE to 2017 CE, when the sample was collected (Supporting Information S1). Previous research has often interpreted δ 18 O as a proxy for weighted mean annual precipitation amount (Baker et al., 2020), which we also demonstrate is the predominant influence on δ 18 O at Cueva Bonita (Figure S3 in Supporting Information S1). However, a growing number of studies have shown that δ 13 C, Mg/Ca, and DCP are also potentially reliable proxies for local water balance (Griffiths et al., 2020), improving our interpretation of hydroclimate when combined with speleothem δ 18 O (see Text S1 in Supporting Information S1).

The CB4 stalagmite was cut, polished and sampled for 15 U-Th dates along its vertical growth axis using a Dremel hand drill with a diamond dental bur. The CB4 sample has uranium concentrations ranging from 37 to 160 ng/g (Table S1 in Supporting Information S1). Calcite powder samples weighing 250-300 mg were prepared at Massachusetts Institute of Technology following methods similar to 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). The corrected ages were calculated using an initial 230 Th/ 232 Th value of 9.8 ± 4.9 ppm to correct for initial 230 Th. The 9.8 ppm initial Th correction value was determined by testing dates corrected with different initial 230 Th corrections for stratigraphic order following methods laid out by Hellstrom (2006) and matching the ages with the radiocarbon bomb peak depth. The uncertainty of 4.9 ppm was scaled proportionally to the normal ±50% correction (4.4 ± 2.2 ppm). U-Th ages range from 78 ± 96 to 2119 ± 162 years before present, however, this study focused on the top 100 mm of the sample with an oldest date of 1189 ± 154 (where present is 1950 CE). All 15 dates fall in stratigraphic order within 2σ uncertainty (Table S1 in Supporting Information S1), but two were identified to be outliers based on low probability of fit for age models (Figure S2 in Supporting Information S1). U-Th ages were combined with fluorescent layer counting to decrease uncertainty. 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).

CB4 was micro-sampled for both stable isotope and trace element analyses using a Sherline micromill at 250 μm increments to a depth of 1 mm, producing 400 samples. The powder for CB4 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. 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.

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 10.1029/2022GL098186 4 of 14 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 CB4 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 (200 total; Table S3 in Supporting Information S1). For plotting/aesthetic purposes, CB2 Mg/Ca, δ 18 O and δ 13 C were smoothed using a moving average. The pandas function DataFrame.rolling().mean() was utilized to smooth the data for plotting only, with the size of the moving window set to 4 years. The full data set reported in Supporting Information S1 is unsmoothed.

Speleothem δ 18 O (δ 18 O speleo ) is primarily controlled by precipitation δ 18 O (δ 18 O precip ) and cave temperature when calcite is deposited under close to isotopic equilibrium conditions (Fairchild et al., 2006;Hendy, 1971;Lachniet, 2009). In low-to-mid latitude regions, cave temperature has a minimal impact and variations in δ 18 O precip dominate the signal. Paleoclimate interpretation of δ 18 O speleo may be complicated, though, by the multiple factors that influence δ 18 O precip. In the tropics, δ 18 O precip has historically reflected the "amount effect" (Dansgaard, 1964), with more precipitation leading to more negative δ 18 O values. However, more recent analyses suggest rainfall δ 18 O is also dependent on processes taking place on more regional scales. For instance, temperature (Lachniet & Patterson, 2009), shifting moisture sources (Aggarwal et al., 2004;Vuille & Werner, 2005), local moisture recycling (Sánchez-Murillo et al., 2016), orographic effects (Sánchez-Murillo et al., 2013), rainout history (Lachniet et al., 2007), microphysical cloud processes (Bony et al., 2008;Konecky et al., 2019;Risi et al., 2008), storm type (Frappier et al., 2007), and relative proportions of stratiform and convective precipitation (Aggarwal et al., 2016) have all been suggested to influence the oxygen isotope composition of rainfall. A comparison between δ 18 O precip to local rainfall amount above Cueva Bonita (Figure S3 in Supporting Information S1), however, suggests precipitation δ 18 O precip is inversely correlated to precipitation amount (r = -0.88, p < 0.05), consistent with the "amount effect" (Dansgaard, 1964). These results are reproduced by IsoGSM data from 1979 to 2015, which demonstrate δ 18 O precip in NE Mexico is closely related to local and regional precipitation amount (Wright et al., 2021;Yoshimura et al., 2008). Lastly, speleothem and precipitation records from Southern Mexico have also demonstrated a strong amount effect (Bernal et al., 2011;Lachniet et al., 2012Lachniet et al., , 2017;;Lases-Hernández et al., 2020;Medina-Elizalde et al., 2010;Pérez Quezadas et al., 2015). While other processes may exert secondary influence on δ 18 O precip , we interpret δ 18 O precip (and consequently δ 18 O speleo ) to dominantly be reflective of precipitation amount. Nevertheless, to minimize the uncertainty of our paleoclimate interpretation, we have applied a multi-proxy approach, using the additional δ 13 C, Mg/Ca, and DCP proxies which are strongly tied to local water balance (See Supporting Information S1).

Calcite samples were analyzed for 14 C at the University of California, Irvine within the Keck Carbon Cycle Accelerator Mass Spectrometry laboratory. Calcite powders were leached with 10% HCl acid, to remove any secondary carbonates, and hydrolyzed with 85% phosphoric acid. Using a modified hydrogen-reduction method (Beverly et al., 2010), samples were then graphitized via Fe catalyzed hydrogen reduction. During data processing calcite powder from a radiocarbon free speleothem was used for blank subtraction. 14 C results were used in combination with the U-Th and lamina counting-based age model to calculate DCP, using IntCal13 data for the atmospheric 14 C activity at the time of speleothem formation (Genty et al., 1999;Reimer et al., 2013). For methods used to convert radiocarbon results to DCP see Supporting Information S1.

Variable DCP measurements, such as those measured in sample CB4, have been shown in previous speleothem studies to be indicative of hydrological changes above the cave (Bajo et al., 2017;Griffiths et al., 2012Griffiths et al., , 2020;;Noronha et al., 2015). The exact mechanism of varying DCP is site specific but changes are commonly driven by the response of soil organic matter and open-versus closed-system bedrock dissolution (Fohlmeister et al., 2011;Griffiths et al., 2012). A fully open dissolution system occurs when conditions are dry and cave drip water dissolved inorganic carbon in the epikarst is in complete isotopic equilibrium with modern 14 C values. In this scenario, the DCP would be equal to 0% (Hendy, 1971). Alternatively, closed system dissolution occurs when conditions are wet, and water in the voids and fractures of the epikarst are in isotopic exchange with 10.1029/2022GL098186 5 of 14 both the atmosphere and the bedrock. Therefore, the upper limit of DCP in a completely closed system would be 50% (Bajo et al., 2017;Hendy, 1971;Noronha et al., 2015). Of course, most DCP values in natural cave systems are somewhere between these extremes, with caves demonstrating an average DCP of 15 ± 5% (Genty et al., 1999). The average CB4 DCP values fall within the low side of this range and are 9.5 ± 0.9%. Given the multitude of controls of δ 18 O, δ 13 C and Mg/Ca (Johnson, 2021), DCP provides an additional hydroclimate proxy.

As a complement to the proxy records, we use general circulation model experiments to isolate the atmospheric response to each possible combination of sea surface temperature anomalies associated with AMV and Interdecadal Pacific Variability (IPV). We use the specified chemistry version of the Whole Atmosphere Community Climate Model with Community Atmosphere Model version 4 physics (Marsh et al., 2013;Smith et al., 2015). The model domain extends from the surface up to 145 km over 66 vertical levels with a horizontal resolution of 1.9 by 2.5° latitude and longitude.

Given that the AMV and IPV can exist in positive (warm tropics), negative (cool tropics), or neutral states, there are nine combinations of SST variability that are used to force the model. The neutral AMV/neutral IPV experiment serves as the control. It is a 200-year continuous integration of the model that is forced with a fixed repeating annual cycle of present-day SST/sea ice concentration variability (1979-2008 average annual cycle from the Hadley Centre Sea Ice and Sea Surface Temperature data set -HadISST (Rayner et al., 2003).

To obtain the anomalous SSTs corresponding to the AMV and IPV, external drivers of climate variability (solar and volcanic) and anthropogenic forcing (greenhouse gasses and aerosols) must be removed from the longterm SST record. The external forcing to the SSTs is obtained by applying signal-to-noise maximizing EOF analysis (Ting et al., 2009) to a global mean SST record derived from a CMIP5 multimodel ensemble of historical and Representative Concentration Pathway 8.5 (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) data. The difference between the spatial pattern of "observed" Extended and Reconstructed SST version 4 (ERSSTv4, Huang et al., 2015) data and the spatial pattern of SSTs corresponding to the aforementioned forced component yields the internal SST variability from which the AMV and IPV patterns are obtained following subsequent filtering for low-frequency (decadal) variability and construction of AMV/IPV indices (see technical note 1 from Boer et al., 2016;Elsbury et al., 2019). The AMV/IPV simulations are identical to the control except that the SST patterns (Figure S4 in Supporting Information S1) corresponding to each combination of AMV and IPV are superimposed on top of the control simulation SSTs and then integrated for 200 years. The precipitation, sea level pressure (SLP) and low-level wind anomalies shown in Figure 4 were produced by differencing the atmospheric fields in these perturbation simulations from that of the control. The AMV index utilized in Figure 3 was created using a detrended, 121-month smoothed, area weighted average of Kaplan SST data set over the N. Atlantic from 0 to 70°N (Enfield et al., 2001).

The most striking centennial-scale feature of this record is the extremely low speleothem δ 18 O values of -6.5‰ at ∼1490 CE (Figure 2), which is also supported by very low δ 13 C values (-11‰) and Mg/Ca ratios (38 mmol/ mol) near the same time. Increased precipitation during the 15th century has been suggested as a dominant driver 10.1029/2022GL098186 6 of 14 of population expansion at major archeological sites, including Tenochtitlan in the Basin of Mexico, and has thus been referred to as the Aztec Pluvial (Sanders et al., 1979). Stahle et al. (2016) provided evidence that the Aztec Pluvial could be a major wet period over the Common Era but lacked older tree ring records to confirm the magnitude and spatial extent of increased precipitation. Our record confirms through multiple geochemical proxies that the Aztec Pluvial represented the wettest conditions in Northeast Mexico over the last millennium. Speleothem δ 18 O from Juxtlahuaca cave in SW Mexico and Ti concentration in lake sediments from Central Mexico also show increased precipitation during this interval (Lachniet et al., 2017;Wogau et al., 2019), suggesting the Aztec Pluvial impacted at least Northeast, Central and Southwest Mexico, if not the entire country. While a strengthening of the North American Monsoon has been invoked to explain increased precipitation in Southwest and Central Mexico (Lachniet et al., 2017), NE Mexico is outside the monsoon's dominant core region, suggesting a strengthened monsoon is not likely the primary mechanism for increased precipitation at this time. A comparison of Cueva Bonita δ 18 O to Eastern and Tropical North Atlantic SSTs instead suggests the 15th century pluvial period was likely forced by anomalously warm Atlantic SSTs (Figures S5 and S6 in Supporting Information S1). Warmer SSTs are known to increase precipitation in NE Mexico by increasing boundary layer moisture convergence, as well as favoring the development of hurricanes (Wang & Lee, 2007). These results are consistent with reconstructions from the Last Millennium Reanalysis Project (Anderson et al., 2019) (Enfield et al., 2001). Both the CB4 and Atlantic SSTs are detrended to account for the impact of anthropogenic warming. The δ 18 O appears to capture both extended periods of positive phases (1920-1960, 2000-2020+) and extended negative phases of the AMV index over the last 150 years. (b) Comparison of AMV phases (positive [2000][2001][2002][2003][2004][2005][2006][2007][2008][2009][2010] minus negative [1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990]) on mean low-level winds and precipitation anomalies. A positive phase leads to increased precipitation in NE Mexico but drying in NW and Southern Mexico.

10.1029/2022GL098186 7 of 14 Notably, the CB4 record also indicates a pattern of decreasing δ 18 O values (from -3.5‰ --6.7‰) beginning near the end of the pre-industrial period, around 1830. This trend is also supported by a shift toward more negative δ 13 C values (-8.9‰ --13‰), decreased Mg/Ca ratios (54-28 mmol/mol), and increased DCP suggesting increased precipitation with anthropogenic warming is robust across multiple geochemical proxies. The decreasing trend of δ 18 O and δ 13 C appears to be in response to Atlantic warming, as Pacific warming is delayed until the early-to mid-20th century (Figure S6 in Supporting Information S1). Interestingly, this trend toward wetter conditions is not obvious in Mexican tree rings (Stahle et al., 2016), but is evident in historical precipitation records, satellite data and re-analysis data from Central Mexico (Martinez-Lopez et al., 2018), suggesting precipitation increases may occur only in summer and early autumn. Although our record is not sub-annually resolved to verify, we suggest the wetting trend may be driven by more extreme pluvial climate events in the late-summer and early-autumn months. This interpretation is supported by evidence for increased tropical cyclone rainfall rates (Knutson et al., 2019) and increased precipitation from more slowly decaying hurricanes on land over the last century (Li & Chakraborty, 2020). Furthermore, extremely wet hurricanes are projected to increase in frequency, making historic flooding (1 in 2000-year) events such as that caused by Hurricane Harvey much more likely (1 c, d, g, andh). Cueva Bonita location is indicated by the star.

in 100-year) by 2100 CE (Emanuel, 2017). Overall, the CB4 record combined with recent analysis of extreme pluvial events suggests NE Mexico could become wetter under anthropogenic climate change.

The strong positive correlation of NE Mexico rainfall to Atlantic SSTs is surprising, considering tree ring records have provided robust evidence of spatially widespread drying in response to positive phases of the AMV (Stahle et al., 2016). However, our record demonstrates a strong positive correlation to Atlantic SSTs not only on centennial timescales but on multidecadal timescales as well. This is evident in a direct comparison of proxies over the last ∼800 years (Figures S5 andS6 in Supporting Information S1), wavelet power spectrum analysis demonstrating a periodicity of 66 years for δ 18 O and 55 years for δ 13 C (Figure S8 in Supporting Information S1), which is close to the previously suggested periodicity of 65 years for the AMV (Schlesinger & Ramankutty, 1994), and a direct comparison of the AMV index to CB4 δ 18 O over the last century (Figure 3a). From lake sediment reconstructions of runoff in NE Mexico, Roy et al. (2016) speculated that the AMV likely altered regional hydroclimate during the early Holocene and Late-Pleistocene but the records did not retain the temporal resolution required to verify the impact of the AMV on regional precipitation. Our record provides the first multiproxy evidence of AMV influence on NE Mexico precipitation over the last millennium.

Interestingly, CB4 proxies not only stand in contrast to tree ring interpretations of the role of Atlantic SSTs, but speleothem δ 18 O, δ 13 C, and Mg/Ca ratios also record wetter conditions during periods of cool Eastern Pacific SSTs, a response also reflected in additional speleothem records from Southern Mexico and Central America (Lachniet et al., 2004(Lachniet et al., , 2017)). The similarity in speleothem records across both Northern and Southern Mexico suggests precipitation may not be out-of-phase in these two regions as previously thought (Bhattacharya & Coats, 2020;Méndez & Magaña, 2010;Stahle et al., 2016). Modern instrumental data also suggests precipitation is mostly in-phase during seasonal (summer) and annual timescales (Figure 1) and only shows a strong dipole precipitation pattern for winter rainfall (Figures 1 and S1 in Supporting Information S1). While changes in winter-spring soil moisture, as typically recorded by tree ring chronologies, are closely linked to changes in early summer soil moisture, they can be poorly correlated with late summer and autumn rainfall (Stahle et al., 2016;St. George et al., 2010). Furthermore, most low-elevation, sub-tropical and tropical tree ring reconstructions from North America are biased toward recording dry extremes while completely missing wet extremes, especially during the summer (Wise & Dannenberg, 2019). Given that summer precipitation in NE Mexico accounts for ∼70% of total annual rainfall (Figure 1e), we suggest the discrepancies between the speleothem data presented here and tree ring-based interpretations are driven by the winter-spring moisture bias of tree rings. While analysis of instrumental data supports this notion, with a positive phase of the AMV leading to increased precipitation in NE Mexico (Figure 3b), historical rainfall is complicated as it is also impacted by Pacific variability and complex forcing of Atlantic variability on the Pacific, and vice versa (Bhattacharya & Coats, 2020).

To test the seasonality and spatial pattern of rainfall in response to SST variability, we utilized a state-of-the-art general circulation model with prescribed patterns of Atlantic and Pacific SST variability. This experimental design allows us to disentangle the influence of the Pacific variability on the Atlantic, and vice-versa. Control runs of this model reliably capture global patterns of observational precipitation and low-level winds (Smith et al., 2015), including Mexico and Central America (Figures S9 and S10 in Supporting Information S1). While our analysis includes a full range of forced-SST conditions (Figures S11-S14 in Supporting Information S1), the natural environment on interannual to decadal timescales is most likely to exhibit an Atlantic-Pacific out-ofphase warming or cooling via changes in the strength of the Walker Circulation (Fosu et al., 2020), and are therefore the focus of this discussion.

During summer, in response to a warm Pacific and cold Atlantic, precipitation decreases across almost all of Mexico and Central America (Figure 4a). Anomalous convection in the Pacific in response to warmer conditions has previously been attributed to drying via an enhanced Walker Circulation, which is also simulated in this study (Figures 4a and S11c in Supporting Information S1). This results in a southward migration of the Atlantic ITCZ and stronger easterly trade winds (Bhattacharya & Coats, 2020;Bhattacharya et al., 2017;Chiang & Sobel, 2002;Giannini et al., 2000Giannini et al., , 2001)). While the contraction of the ITCZ is known to decrease precipitation in Southern Mexico and Central America (Asmerom et al., 2020), stronger easterly trade winds are thought to increase precipitation in Northern Mexico via an intensification of easterlies and the Caribbean Low-Level Jet (CLLJ) (Wang & Lee, 2007). Although low-level wind anomalies in model simulations correctly replicate the intensification of the CLLJ (Figure S10 in Supporting Information S1), models demonstrate a stronger CLLJ instead leads to decreased precipitation over much of Mexico (Figures 4b and4f). This response is also replicated when SST conditions are reversed, which drives a weakening of the CLLJ and increased precipitation (Figured 4b and 4c). On longer orbital to interannual timescales, a stronger CLLJ has been linked to drier conditions in NE Mexico through Atlantic SST cooling and an enhanced wind-evaporation-SST feedback loop (Wright et al., 2021). However, this experiment utilizes prescribed SSTs and we therefore cannot attribute observed precipitation changes to this mechanism. We instead suggest that warmer Atlantic SSTs lead to a reduction in the strength of the CLLJ and, consequently, a reduction in vertical wind shear. Decreased vertical wind shear appears to be further amplified by cooler Pacific SSTs (Figure S14 in Supporting Information S1), which fosters the formation of deep convective storms and increases precipitation throughout most of Mexico (Figure 4b). This mechanism is further supported by observational records and previous modeling results (Wang, 2007;Wang & Lee, 2007), which have linked decreased vertical wind shear to more frequent and larger magnitude hurricanes in the Tropical North Atlantic.

Another notable result of the fixed SST simulations is the contrasting response of precipitation throughout most of Mexico during summer and winter (Figures 4a and4b vs. 4c and 4d, Figure S13 in Supporting Information S1). The only region that appears to respond consistently in both seasons is Northwest Mexico, which is strongly influenced by the North American Monsoon and is likely to be more sensitive to variability in Pacific SSTs. In response to a warm Pacific/cold Atlantic, winter precipitation throughout much of Mexico slightly increases (Figure 4c). Increased winter precipitation in Northwest and Central Mexico has been attributed to a strengthening of the North Pacific storm-track extending further south and west in response to anomalous waves (Seager & Hoerling, 2014). While this is not likely to drive increased precipitation as far east as NE Mexico (Wright et al., 2021), warmer Pacific SSTs could drive more frequent and intensified cold fronts from the North, increasing light, low-intensity rainfall to the region (Magana et al., 2003). This interpretation is consistent with results presented here, where anomalous northerly low-level winds are produced by forced-SST experiments (Figure 4g).

When SST conditions are reversed during winter (warm Atlantic/cold Pacific), model results suggest there is a reversal in the low-level winds pattern, leading to drier conditions (Figure 4h). This is consistent with the response recorded in tree rings which have suggested drier conditions during a positive phase of the AMV (Méndez & Magaña, 2010;Stahle et al., 2016). Previous work has suggested warmer Atlantic SSTs reduce the interbasin SST and SLP gradient, which weakens the CLLJ leading to decreased moisture transport to Northern Mesoamerica and increased precipitation in Southern Mesoamerica (Bhattacharya & Coats, 2020;Bhattacharya et al., 2017;Mestas-Nuñez et al., 2007;Stahle et al., 2016). Although a decreased and southward diverging CLLJ is reproduced in this study (Figure 4h), producing the north-south dipole precipitation pattern (Figure 4d), we believe the decrease in northern frontal storms is a more important mechanism for reducing precipitation in NE Mexico since a weaker CLLJ tends to increase precipitation during summer months (Figures 4b and4f). Therefore, winter simulations of rainfall support tree ring interpretations of: (a) an out-of-phase dipole spatial pattern, and (b) a dominant role of Pacific SSTs in controlling regional precipitation. However, winter precipitation (DJFM) only accounts for a small fraction of total annual rainfall, with winter contributing less than 8% of annual rainfall in NE Mexico (Figure 1e). We therefore suggest relative changes in Atlantic SSTs are much more important in controlling total annual precipitation amount in the region.

The new CB4 speleothem record from NE Mexico combined with forced-SST climate model results highlights the precipitation dipole pattern is far more spatially complex in Mexico than previously thought (Figures S15 and S16 in Supporting Information S1). We suggest previous reconstructions of the dipole pattern may have utilized tree rings that are potentially biased toward winter, spring, or early-summer rainfall (Stahle et al., 2016(Stahle et al., , 2020)), minimizing the role of the Atlantic in modulating late-summer and early-autumn precipitation throughout most of the region. While variability in Pacific SSTs can still play an important secondary role in regulating precipitation, we suggest warmer Pacific SSTs predominantly drive decreased precipitation over the region. We instead propose that warmer Atlantic SSTs drive increased precipitation in the region by increasing summer precipitation, which accounts for 70% of total annual rainfall in the region (Figure 1e).

Anthropogenic carbon emissions will continue to warm Tropical Atlantic SSTs in the future (Chen et al., 2018). Given the observed trend of increased precipitation in NE Mexico over the industrial period indicated by decreasing CB4 δ 18 O, this suggests NE Mexico could become wetter in the future, in contrast to current model projections. Importantly, though, forced-SST experiments show that precipitation in Mexico is most sensitive to the warming of the Tropical Atlantic relative to the Tropical Pacific. Unfortunately, models currently exhibit significant discrepancies in determining this interbasin SST gradient, greatly decreasing confidence in regional rainfall projections (Bhattacharya & Coats, 2020). In part, this is further underscored by the fact that models produce different SST responses to external forcing (volcanic vs. orbital) even over the last millennium (Bhattacharya & Coats, 2020). Furthermore, SST biases in the eastern Tropical Pacific and the Tropical Atlantic, driven by poor model representation of topography in Mesoamerica (Baldwin et al., 2021) and oceanic and atmospheric circulation (Imbol Nkwinkwa et al., 2021), further confound the ability to accurately reproduce the SST interbasin gradient. Lastly, the large inter-model spread in precipitation is also partially driven by the paucity of paleoclimate records and uncertainties in the response of rainfall and SSTs to internal variability (Deser et al., 2020). By using a model with prescribed SSTs based on observational data, this study was able to circumnavigate some of these issues and highlight conditions and mechanisms that drive precipitation change in Mexico. However, accurate projections of future rainfall will necessitate improvements in model SST sensitivity to various forcings, more accurate topography, better representation of atmospheric circulation, and a robust understanding of the response of precipitation to internal variability. In addition, improved constraints on precipitation isotope controls, through increased modern sampling and detailed analyses of isotope-enabled climate model simulations of the region, are an important area for future study that could further strengthen paleoclimatic interpretations.

We suggest robust proxy records of past hydroclimate such as CB4 be utilized in future studies to isolate models that best reproduce both the amount and temporal pacing of precipitation variability. Given the demonstrated dependence of precipitation on regional SSTs, we believe these models will also accurately simulate large scale interactions between the Pacific and Atlantic (i.e., the interbasin gradient), which has important implications for predicting hydroclimate in other regions as well. Particularly, annual to sub-annual speleothem records would be ideal to confirm precipitation during the dominantly wet summer season and for the direct evaluation of ENSO, which unfortunately was not possible given the temporal resolution of CB4. Overall, however, CB4 is an important contribution to the paleoclimate record in this critically understudied region. We have provided four independent hydrological proxies of past precipitation change over the last millennium and have validated the proposed mechanisms of precipitation change utilizing forced-SST experiments. We hope this study will serve as the foundation for future work in Northern Mexico. Zhang, Jessica Wang, Chris Wood, and Elizabeth Patterson for assistance with lab work. We thank Crystal Tulley-Cordova for sharing the precipitation collector design.