Atlantic‐Pacific Gradients Drive Last Millennium Hydroclimate Variability in Mesoamerica

The paleoclimatic record from Mexico and Central America, or Mesoamerica, documents dramatic swings in hydroclimate over the past few millennia. However, the dynamics underlying these past changes remain obscure. We use proxy indicators of hydroclimate to show that last millennium hydroclimate variability was dominated by opposite‐signed moisture anomalies in northern and southern Mesoamerica. This pattern results from changes in moisture convergence driven by Atlantic‐Pacific interbasin temperature gradients. While this pattern is reproduced by several models and multiple experiments with a single model, models appear to disagree about the underlying dynamics of this interbasin gradient. Moreover, disagreement about the interbasin gradient, and associated hydroclimate pattern, dominates spread in 21st century regional hydroclimate projections. These results emphasize the role of interbasin asymmetries in governing past and future regional climate change. They also demonstrate that paleoclimate studies can elucidate mechanisms directly relevant to projecting future hydroclimate in climate change hot spots like Mesoamerica.


Introduction
In 2018, a drought stretching from southern Mexico to Panama created food shortages for over two million people (Moloney, 2018). Events like this highlight the profound impact of drought on livelihoods in Mexico and Central America, hereafter called Mesoamerica. Even small shifts in rainfall can place stress on regional water resources, since the mean hydroclimatology of the region is highly seasonal. For example, the majority of regional rainfall occurs in summer, and much of Mesoamerica experiences a midsummer drought, or canìcula, with limited rainfall during July and August (Magaña et al., 1999).
Instrumental data and historical model simulations provide robust evidence of a region-wide twentieth century drying trend (Anderson et al., 2019;Marvel et al., 2019;Neelin et al., 2006). Observed drying has been attributed to both higher temperatures and reductions in precipitation (Fuentes-Franco et al., 2015;Herrera et al., 2018;Rauscher et al., 2011). Models project additional drying by the end of the 21st century (Rauscher et al., 2011). However, the limited instrumental record hinders our ability to understand how the magnitude of historical or projected drying compares to natural climate variability, which in many areas is larger than observed drying trends (Anderson et al., 2019). Longer-term records of hydroclimate are therefore essential to adequately characterize regional climate variability and contextualize future projections.
Paleoclimatic data from Mesoamerica provides ample evidence of past hydroclimate variability across a range of timescales. For instance, lake and marine sediments, as well as speleothems, record century-long dry intervals over the last few millennia (Akers et al., 2016;Bhattacharya et al., 2015;Curtis et al., 1996; ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.  Haug et al., 2003;Hodell et al., 1995Hodell et al., , 2005aHodell et al., , 2005bJones et al., 2015a;Kennett et al., 2012;Lachniet et al., 2012;Medina-Elizalde & Rohling, 2012;Medina-Elizalde et al., 2010;Metcalfe et al., 2010;Park et al., 2010;Stansell et al., 2013;Park et al., 2019). Researchers disagree, however, about the spatial extent and magnitude of these droughts (Bhattacharya et al., 2017;Evans et al., 2018;Medina-Elizalde & Rohling, 2012). On shorter timescales, tree-ring chronologies provide evidence of historic droughts with antiphased hydroclimate anomalies between northern Mexico and southern Mexico and Guatemala, although data from the southern region is limited (Stahle et al., 2016). Moreover, comparing patterns of drought across proxy archives has been difficult given the distinct timescales recorded by different archives (Metcalfe et al., 2010). As a result, the spatiotemporal characteristics of past Mesoamerican hydroclimate variability are poorly understood, and the relationship between regional hydroclimate, external forcing, and internal variability remains ambiguous (Bhattacharya et al., 2017;Cook et al., 2012;Hodell et al., 2001;Lachniet et al., 2017;Oglesby et al., 2010).

10.1029/2020GL088061
Here, we identify the dominant patterns of hydroclimate variability over the last millennium at decadal and longer timescales. Our analysis incorporates multiple types of proxy archives, including lakes and speleothems ( Figure S1 and Table S1 in the supporting information) and the Mexican Drought Atlas (MXDA), a tree-ring-based reconstruction of the summer (June-August, JJA) Palmer Drought Severity Index, or PDSI (Stahle et al., 2016). We leverage the Last Millennium Ensemble (LME) of the Community Earth System Model Version 1.2 (CESM1.2), a suite of last millennium simulations run with different forcing combinations, as well as last millennium simulations from other models, to investigate the dynamics underlying past hydroclimate variability and assess models' ability to simulate the patterns found in proxy archives.

Proxy and Instrumental Data
We assembled a network of hydroclimate-sensitive proxy records, including carbonate records from closed-basin lakes that are sensitive to precipitation minus evaporation (P-E) (Bhattacharya et al., 2015;Curtis et al., 1996;Hodell et al., 2005a;Park et al., 2019) and records of elemental concentration interpreted as a signal of total annual rainfall (Haug et al., 2003;Jones et al., 2015;Metcalfe et al., 2010). We also include speleothem oxygen isotope records, interpreted as a signal of rainfall amount (Akers et al., 2016;Kennett et al., 2012;Lachniet et al., 2012;Medina-Elizalde et al., 2010). We chose multidecadally sampled records constrained by at least five age markers over the past 2,000 years and at least four over the last 1,000 years. We focus on the interval between 850 and 1850 CE rather than earlier periods as this is contemporaneous with model simulations of the last millennium (see Text S1).
To account for the role of age uncertainty in the timing of events in proxy records, we propagated age uncertainty through each record using Bayesian age models (Blaauw & Christen, 2011). Highly skewed data (e.g., Juanacatlàn and Chichancanab) were log-transformed to decrease skewness (Wilks, 2011). For precisely dated records (e.g., speleothem or varved records that are dated via layer-counting or high precision dates), we add an error term of ±20 years at each depth. Given the large age uncertainties associated with the records exclusively dated with radiocarbon, our analysis is insensitive to the magnitude of age error assigned to these higher resolution proxy archives. Following the methods of Anchukaitis and Tierney (2013), for each iteration of our proxy empirical orthogonal function (EOF) analysis, we sampled one iteration of each sites' age model, interpolated all the records to a decadal time step, and performed EOF analysis on the correlation matrix of standardized data (z scores). This prevents sites with higher amplitudes of variability from having an outsized influence on the resultant EOFs (Wilks, 2011). We performed 1,000 iterations of this analysis to create an ensemble of loadings, time series, and variance estimates ( Figures S2 and S4 and Table S2). We use scree plots to assess significance, where eigenvalues are plotted by decreasing value. Because the variance associated with each eigenvector is directly proportional to the eigenvalue, we plot the variance explained by each eigenvector for ease of interpretation ( Figures S3 and S4). The "break" in the slope on a scree plot divides significant modes of variability from those indistinguishable from EOFs of random Gaussian noise (Wilks, 2011).
We also analyze an independent, tree-ring-based hydroclimate reconstruction. The MXDA is a reconstruction of JJA PDSI between 1400 and 2012 CE and features over 250 climate sensitive tree-ring chronologies (Stahle et al., 2016). While this data set may be more sensitive to spring and early summer hydroclimate variability than the lake and speleothem data, potentially complicating interpretations, it is a valuable comparison since the MXDA features annually resolved, spatially explicit data. We applied a 10-year low-pass Butterworth filter to emphasize decadal and centennial variability (Wilks, 2011).
We contrast the patterns found in proxy data with the spatiotemporal variability found in the Global Precipitation Climatology Centre (GPCC) data product (Schneider et al., 2017) in the domain from 5-40°N and 45-135°W, limiting our analysis to the period from 1930 to 2017. We also explore the relationship between GPCC rainfall and sea surface temperatures (SSTs) (Ishii et al., 2005) and large-scale atmospheric fields from the National Centers for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) reanalysis (Kalnay et al., 1996). Given limitations in data coverage over Mesoamerica and the short duration of the instrumental record, these analyses are primarily exploratory, but they provide an important complement to proxy data and models.

Model Simulations
We leverage the CESM1.2 LME to investigate the large-scale dynamics underlying hydroclimate variability in Mesoamerica. The LME features transient simulations spanning 850 to 1850 CE, including 13 ensemble members forced with all forcings, and several smaller ensembles with single forcings (e.g., volcanic forcing only) (Otto-Bliesner et al., 2016). Simulations are run with a 2°atmosphere and land model and a roughly 1°o cean. We analyzed 12 ensemble members with the full suite of forcing factors, 5 volcanic-only ensemble members, and an unforced control simulation. This large ensemble of model simulations allows us to analyze the interactions of climate variability and climate forcing (Otto-Bliesner et al., 2016). Although the model is run at a relatively course resolution and overestimates precipitation in northern Mexico, it captures the broad features of the region's climatology, including the seasonal northward progression of rainfall ( Figure S3).
Proxy data likely reflect lower frequency variability, though on very different timescales: Lacustrine proxy records capture centennial to millennial-scale trends, while speleothem and tree-ring data reflect decadal to centennial-scale variations. We analyze decadally filtered (10-year low-pass Butterworth filter) fields of modeled precipitation as these data will contain variability across all of the timescales reflected in the various proxy data. Proxies likely also reflect different seasonal signals: In this region, speleothems and lakes may predominantly capture summer rainfall (Bhattacharya et al., 2017), while the MXDA may reflect spring and early summer rainfall. To allow for the varying seasonal signals archived in the proxy records, we analyzed patterns associated with annual average rainfall as well as different seasonal averages.
Model EOFs of precipitation were assessed over land for the domain from 5-40°N and 45-135°W, using standardized data to prevent regions with higher climatological rainfall values from skewing the analysis. We investigate large-scale dynamics using model fields of surface pressure, temperature, zonal and meridional winds, and moisture convergence calculated from monthly averaged zonal and meridional moisture transports. We also analyzed individual moisture budget terms that decompose moisture convergence into contributions from dynamical changes, changes in moisture, and transient eddies (Text S2).
To assess whether similar modes of variability are produced by other coupled models, we analyze last millennium simulations from CCSM4, IPSL-LR, MPI-LR, and HadCM3 (Landrum et al., 2013;Schmidt et al., 2011). We complement this analysis by investigating sources of intermodel variation in future model projections for Mesoamerica. For this latter analysis, we compiled 17 simulations of the high-emissions RCP8.5 scenario from the fifth Coupled Model Intercomparison Project (CMIP5) archive (Table S3), and analyzed fields of precipitation and SST, regridded to a common 1°by 1°resolution via bilinear interpolation. We quantify precipitation change in each model as the difference between rainfall at 2070-2100 and rainfall from a preindustrial control simulation. We then subtract the multimodel mean precipitation change from that of each individual model before conducting an EOF analysis over the region spanning 5-40°N and 45-135°W. This EOF analysis reflects the dominant patterns of disagreement in rainfall projections across models-in short, an analysis of the spatial patterns that best characterize intermodel spread in future regional hydroclimate projections.

Leading Modes of Hydroclimate Variability
Our results reveal that last millennium tree-ring and lake and speleothem proxies show consistent spatial patterns of hydroclimate variation but resolve different timescales of variability. The leading mode in the lake and speleothem proxy network is a statistically significant dipole that explains 25% of the variance in the data set, with wet conditions in southern Mesoamerica accompanied by drier conditions in highland and desert regions of Mexico (Figures 1 and S4). For sites in southern Mesoamerica, this mode is characterized by wetter conditions during the Medieval Climate Anomaly (MCA), between approximately 900-1300 CE, and drier conditions starting at 1450 CE. In contrast, three sites near 20°N show wetter conditions after 1450 CE. Due to high age uncertainty for radiocarbon-dated records, the dipole mode primarily exhibits centennial variability ( Figure S2 and Table S2). A similar, statistically significant north-south dipole mode, with a change in sign near 19°N, explains 47% of the variance in the 10-year low-pass-filtered MXDA ( Figure S5). In contrast to the centennial-scale fluctuations in the lake and speleothem network, the dipole mode in the MXDA reflects decadal and multidecadal variability, which could reflect detrending methods used in the MXDA that primarily preserve high to medium frequency variability (Stahle et al., 2016).
Independent evidence suggests that this dipole mode reflects a physically meaningful pattern of hydroclimate variability. While historical droughts recorded in the MXDA show a dipole structure strikingly similar to the pattern in Figure 1 (Stahle et al., 2016), this pattern has not been previously been described in lake and speleothem data. However, close inspection of the raw proxy data provides strong evidence of this mode: Laguna X'Caamal, Juxtlahuaca, Yok Balum, Laguna El Gancho, and the Cariaco Basin exhibit a step change to drier conditions near 1400 CE ( Figure S2), while northern sites like Aljojuca, Juanacatlàn, and Rincon de Parangueo exhibit a step change to wetter conditions ( Figure S2). Independent evidence from the instrumental record (e.g., GPCC) shows that recent decadal hydroclimate variations in Mesoamerica show antiphased changes between the southern region and northern Mexico (Figure 1). We next use model simulations of the last millennium to explore the mechanisms underlying this dipole mode of variability.

Dynamical Drivers of Hydroclimate Variability
The dominant mode of Mesoamerican annual precipitation, as well as different seasons' precipitation, in the fully forced CESM LME simulations strongly resembles the dipole mode in proxy data (Figures 2 and S5). The pattern switches from positive to negative loadings at approximately 20°N and accounts for between 36% and 41% of total variance in regional precipitation across the 12 ensemble members. This mode is positively correlated with moisture convergence in southern Mesoamerica and negatively correlated with moisture convergence in northern Mesoamerica (Figure 2b). The dominant contribution to moisture convergence changes associated with the dipole mode comes from changes in the winds (Figure 2d), since this term and total moisture convergence exhibit similar correlation patterns with the dipole mode. Composite plots show that positive (negative) excursions of the dipole mode are associated with weaker (stronger) easterly winds (Figure 2d). The simulated dipole mode's correlations with moisture convergence changes due to thermodynamic effects (e.g., changes in moisture) and transient eddies does not resemble the correlation pattern with the total moisture convergence field, suggesting that these terms are less important (Figures 2 and S6).
During positive excursions of the dipole mode, weaker easterly winds decrease moisture export from southern Mesoamerica, increasing local rainfall, while simultaneously decreasing moisture convergence into dry regions of northern Mesoamerica, decreasing rainfall in these regions (Seager et al., 2009;Wang et al., 2006) ( Figure 2d). Diminished easterly winds across Mesoamerica are ultimately driven by a large-scale pattern of lower pressure over the tropical Atlantic and higher pressure over the equatorial Pacific, reflecting a weaker Atlantic-Pacific pressure gradient (Figure 2e) (Giannini et al., 2001). This anomalous interbasin pressure gradient is also associated with northward displacements of the Atlantic Intertropical Convergence Zone (ITCZ) as well as a weaker southwestern edge of the North Atlantic Subtropical High (NASH) (Kushnir et al., 2010;Wang et al., 2006).
Observational analyses corroborate the dynamical insights provided by the LME. Similar wind anomalies and pressure gradients are associated with the dipole mode in the GPCC, although the pattern is less coherent, possibly reflecting a short, noisy instrumental record ( Figure S8). This Atlantic-Pacific interbasin gradient has also been shown to drive interannual variability over the Caribbean region, especially during the late summer season (Taylor et al., 2002(Taylor et al., , 2011.

Diversity of SST Responses
While the dipole is ubiquitous across simulations with different forcing combinations in the LME, it is associated with different SST patterns ( Figure S6). In volcanic-only ensemble members of the LME, the dipole mode explains between 38% and 42% of regional rainfall variance and exhibits a strong relationship with equatorial Pacific SSTs. These volcanic-only simulations produce weak correlations with tropical Atlantic SST (Figure 3). In contrast, an unforced control simulation produces a dipole mode explaining 38% of precipitation variability, but features stronger correlations with tropical Atlantic variability than volcanic-only simulations. The dipole mode in full-forcing simulations exhibits strong negative correlations with both equatorial Pacific SSTs and weak but significant correlations with northern Atlantic variability (Figure 3).
We next explore how the dynamics underlying the dipole mode last millennium simulations differ across models four other last millennium simulations. The dipole mode is the leading pattern of last millennium hydroclimate variability in all models ( Figure S9). However, while each model features negative associations with Pacific temperatures and positive associations with Atlantic SSTs, the spatial pattern of correlation coefficients and regression coefficients differ between models (Figures 4 and S9). In the Atlantic sector, HadCM3, CESM1.2, and, to a lesser extent, CCSM4 associate the dipole with a horseshoe-like pattern of North Atlantic SST variability, characteristic of Atlantic Multidecadal Variability (AMV) (Ba et al., 2014;Ting et al., 2011). This includes a pattern of interhemispheric temperature contrast resembling the Atlantic "meridional mode" that is typically associated with changes in the position of the Atlantic ITCZ (Xie & Carton, 2004). In contrast, IPSL and MPI produce positive correlations both north and south of the equator in the Atlantic (Figure 4). In the Pacific, MPI and HadCM3 produce weak negative associations centered on the tropical eastern or central Pacific (Figure 4). However, IPSL, CESM1.2, and CCSM4 produce strong negative correlations between the dipole mode and SSTs in the tropical eastern Pacific. In instrumental data, the dipole mode has a much stronger linear relationship with a pattern resembling AMV, and a much weaker relationship with Pacific variability, than is produced in any of the models (Figure 4). This could reflect model bias or could reflect a biased signal from a short, noisy instrumental record.
Different correlation patterns of SST across models suggest the possibility of different mechanisms governing Atlantic-Pacific gradients in each model. In some models, this gradient may be driven by decadal changes in El Niño frequency and/or equatorial Pacific mean state that remotely influence the tropical Atlantic. For instance, delayed warming in the tropical Atlantic both north and south of the equator follows a peak El Niño event (Chang et al., 2006;Chiang & Lintner, 2005) (Figures 4 and S10). CMIP5 models exhibit a large range in the strength of El Niño-Southern Oscillation (ENSO) variability (Bellenger et al., 2014). It is therefore possible that models with stronger decadal variability in the equatorial Pacific will reflect Atlantic-Pacific interbasin gradients dominated by Pacific-centric dynamics and will produce higher correlations between the dipole mode and equatorial Pacific SSTs. HadCM3 and MPI, which produce the weakest correlations between the dipole mode and equatorial Pacific temperature, feature lower decadal variability in the Niño 3 index in their last millennium simulations compared to the other models (Table S4). In other models, Atlantic dynamics may drive changes in the interbasin gradient. Warming in the equatorial and tropical south Atlantic is known to induce cooler conditions in the equatorial and eastern Pacific on decadal timescales, similar to the temperature configuration associated with the dipole mode (Chikamoto et al., 2015;Jia et al., 2019;McGregor et al., 2014). However, low-resolution ESMs with insufficiently resolved Mesoamerican topography, like many simulations used here, may not fully capture the influence of Atlantic processes on equatorial Pacific variability (Xie et al., 2008). Fully assessing the The correlation between the GPCC dipole mode and the COBE SST product. Regions with statistically insignificant correlations have been masked out (Ebisuzaki, 1997).
diversity of Atlantic-Pacific linkages across a suite of models is beyond the scope of this work but is a critical area for future research. These efforts will benefit from high-resolution last millennium temperature records, especially from the tropical Atlantic, where data is limited (Tierney et al., 2015).

Relevance to Future Projections
Intermodel differences in large-scale SST fields likely play a role not just in the simulation of past, but also future, hydroclimate in Mesoamerica. To explore this possibility, we perform an EOF analysis of the intermodel spread in 21st century Mesoamerican rainfall projections across 17 CMIP5 models. The multimodel ensemble mean projects drying across Mesoamerica, with the strongest drying projected to occur in the southern region. Nevertheless, the intermodel standard deviation of future rainfall change is comparable to the mean response projected by models ( Figure S11). Nearly half (48%) of the variance in the rainfall response across models can be explained by a mode featuring a dipole pattern, with wetter conditions in southern Mesoamerica and drier conditions in the north (Figure 5a). The structure of this EOF suggests that a negative (positive) principal component score indicates that the associated model projects more (less) drying in southern Mesoamerica than in northern Mexico, compared to the multimodel mean.
Whether a model has a negative or positive principal component score is strongly related to the change in the interbasin gradient between the Atlantic and the Pacific simulated by that model (Figure 5b). There is a strong positive relationship (r=0.73) between principal component scores for each model and an index of equatorial Atlantic SSTs (4°N to 4°S, 40-20°W following Voldoire et al., 2019) subtracted from Niño 3 temperature (5°N to 5°S, 150-90°W) ( Figure 5). Models with a positive rainfall PC score simulate stronger warming in the equatorial Atlantic and less warming in the equatorial Pacific compared to the multimodel mean, while models with a negative score produce the opposite. This suggests that the eastern Pacific and tropical Atlantic are key centers of action where models disagree about the magnitude and sign of future SST change, although further work is needed to identify the sources of these SST differences (Ba et al., 2014;Ting et al., 2011;Voldoire et al., 2019;Yang et al., 2017). Moreover, this analysis demonstrates that divergent SST projections are a key source of model disagreement in Mesoamerican hydroclimate projections.

Conclusions
We used a compilation of proxy indicators of hydroclimate to show that the leading mode of decadal and longer timescale hydroclimate variability in Mesoamerica features opposite-signed hydroclimate anomalies between southern and northern portions of the region. This pattern appears in tree-ring-based drought reconstructions and in lake and speleothem data. This dipole mode is driven by changes in moisture transport that result from an interbasin gradient of pressure and temperature between the Atlantic and Pacific. This emphasis on zonal asymmetries across the Atlantic and Pacific, which has been noted elsewhere (Lachniet et al., 2017), complements perspectives that focus on meridional migrations of the Atlantic and Pacific ITCZ (Haug et al., 2003;Park et al., 2019), since northward migrations of the ITCZ are typically accompanied by cool eastern equatorial Pacific SSTs and warm conditions in the northern tropical Atlantic. Importantly, climate models appear to disagree on the relative role of Pacific and Atlantic dynamics in driving this variability.
Disagreement about future changes to the interbasin gradient dominate intermodel spread in the magnitude and pattern of future rainfall change in Mesoamerica. It is therefore critical to identify the models that best simulate large-scale changes and interactions between the Atlantic and Pacific. Comparisons of models and data for key paleoclimatic intervals have the potential to serve as a "testing ground" for assessing model skill in simulating changes in past Atlantic-Pacific temperature variability and associated hydroclimate dynamics. This will require new high-resolution temperature records of the last millennium, especially from data-limited regions. Paleoclimatic research is therefore of critical importance for understanding future hydroclimate change in regions projected to experience dramatic departures from preindustrial climate, including Mesoamerica.

Data Availability Statement
All raw proxy data are available for download from NOAA's National Center for Environmental Information Paleoclimatology Database or, in the case of Rincon de Parangueo, available in Park et al. (2019). Matlab code for performing EOF analysis on proxy data and model output will be archived upon publication on the NOAA NCEI Paleoclimatology Database. CMIP5 model simulations and the CESM1.2 are publicly available via the Earth System Grid (https://esgf-node.llnl.gov/projects/esgf-llnl/). Model code for CESM1.2 is publicly available at http://www.cesm.ucar.edu/models/cesm1.2/ website.