Research ArticleCLIMATOLOGY

Intertropical convergence zone variability in the Neotropics during the Common Era

See allHide authors and affiliations

Science Advances  14 Feb 2020:
Vol. 6, no. 7, eaax3644
DOI: 10.1126/sciadv.aax3644


Large changes in hydroclimate in the Neotropics implied by proxy evidence, such as during the Little Ice Age, have been attributed to meridional shifts of the intertropical convergence zone (ITCZ), although alternative modes of ITCZ variability have also been suggested. Here, we use seasonally resolved stalagmite rainfall proxy data from the modern northern limit of the ITCZ in southern Belize, combined with records from across the Neotropics and subtropics, to fingerprint ITCZ variability during the Common Era. Our data are consistent with models that suggest ITCZ expansion and weakening during globally cold climate intervals and contraction and intensification during global warmth. As a result, regions currently in the margins of the ITCZ in both hemispheres are likely transitioning to more arid and highly variable conditions, aggravating current trends of increased social unrest and mass migration.


The intertropical convergence zone (ITCZ) is the world’s most important rainfall belt, affecting the livelihood of billions of people globally. Understanding its behavior during future warmer climate is thus vitally important. Intergovernmental Panel on Climate Change projections suggest that the Neotropics are particularly vulnerable to warming-induced droughts (1). On millennial time scales, shifts in ITCZ mean position likely reflect changes in the temperature contrast between the two hemispheres, with ITCZ migration toward the warmer hemisphere. However, the exact nature of the ITCZ response to global temperature changes remains unsettled. During the Common Era (CE), substantial evidence that climate change associated with the Little Ice Age (LIA) was global in extent (2, 3) now exists. Proxy records from the Neotropics suggest that the ITCZ shifted to the south during the LIA (46); however, recent model results (7) and limited observations elsewhere (8, 9) suggest that the ITCZ may have expanded in both hemispheres while weakening in the central equatorial core region. The apparent contradiction between paleoclimate proxy data, which seemed to imply wetter Central America with increased warming (10) and model results suggesting drier climate with warming, has been noted (11). Previous observations had several limitations, including combining monsoon and other proxies that are not necessarily reflective of purely ITCZ modulation (5, 8) or records not capturing the full ITCZ annual excursion (4, 6).

Here, we use new monthly-scale speleothem rainfall proxy data from a cave site at the ITCZ’s current northern margin in Central America coupled with published data from the southern margin of the ITCZ to explore its variability throughout the CE (Fig. 1). We then use this information to determine whether the combined data are best explained by (i) meridional ITCZ shifts exclusively or (ii) expansion and contraction of the ITCZ. The results have profound implications for predicting future hydroclimate of the global tropics.

Fig. 1 Locations of Yok Balum Cave (magenta circle) and other records discussed in the text.

The average annual northernmost (July) and southernmost (January) margins are also shown (yellow and blue bands, respectively).


The climate of the study area in southern Belize is tropical. There is a latitudinally controlled north-south and a topographically controlled east-west gradient in annual rainfall, with a total annual rainfall of ~1300 mm in the drier north and west and as high as 4500 mm in the south (12). Yok Balum Cave, where stalagmite YOK-G comes from, is found in southern Belize (16°12′30.780″N, 89°4′24.420″W; fig. S1). The cave has a stable temperature of 22.74° ± 0.09°C (1σ) with high relative humidity approaching 100% throughout the year (fig. S2). The area receives about 400 to 700 mm/month during the rainy season (June to September) and a factor of less than 10 during the dry season (February to April) (12). This contrast in seasonal rainfall is due to annual migration of the ITCZ. The ITCZ reaches its northmost position during the boreal summer, the rainy season in our study area, while the driest season is at a time when the ITCZ is at its southernmost position during boreal winter (13, 14). The ITCZ locally represents a zone of convergence resulting from the coincidence of cross-equatorial winds from the Southern Hemisphere and northeast trade winds blowing from the North Atlantic subtropical high (15). The El Niño–Southern Oscillation (ENSO) modulates Caribbean-sourced moisture through its influence on the Caribbean Low-Level Jet (CLLJ). During warm ENSO events, the CLLJ is stronger and vice versa, although the resultant impact on precipitation can be highly varied (16). The correlation between ENSO and historical precipitation in our study area is weak (17). ENSO influence on Central American precipitation is more seen on the Pacific Coast in which the warm (cold) phase of ENSO results in drier (wetter) conditions (17). Previous studies have shown that the moisture source for tropical Belize is primarily the Caribbean Sea and is carried by the CLLJ (15, 16, 18).


Stalagmite YOK-G is a very clean aragonitic speleothem with high uranium (U) and low detrital thorium (Th) concentrations that grew from 400 to 2006 CE (the year when it was collected). We measured 52 U-series ages (table S1) using minimal amounts of powder to reduce errors related to drill-hole sizes and were able to obtain mean age uncertainties of 7 years (2σ) across the entire record (Fig. 2A). Drip monitoring data from the site reported in a previous study (12) demonstrates that large rain events recorded by a weather station directly above the cave site affect the drip hydrology within a few hours. In addition, in the same study, detailed high-resolution analysis of radiocarbon (14C) of the younger 50 years was performed to address fluid residence time directly. We found that the first fraction with 14C enrichment corresponds to CE 1954 to 1956 peak in the Northern Hemisphere bomb 14C compilations (19), in good agreement with the U-series age of CE 1955 for the same layer, demonstrating that there is no discernable lag between surface hydrology and the drip in the cave. Details of our chronology and age model construction are discussed in Materials and Methods. Two previous studies (12, 20) reported on the upper portion of YOK-G (since 1550 CE). Here, we present new YOK-G δ13C and trace element datasets that extend the previously published YOK-G δ13C record by 1150 years. The complete YOK-G isotope record now consists of 7151 δ13C data points at ~0.22-year mean resolution overall, with monthly resolution across key intervals. The application of δ13C as a rainfall proxy in YOK-G was previously established (12, 20, 21). Several factors can influence the δ13C values in speleothems. The amount of organic matter in the soil above a cave can vary, leading to higher (lower) organic contribution and lower (higher) δ13C values. Degassing due to prior calcite precipitation (PCP) increases the δ13C value of the solution from which a stalagmite forms. Drier conditions lead to greater PCP due to longer residence time in the epikarst and increased ionic strength of the fluid. Changes in vegetation, from C3 to C4 plants, reflect more arid conditions. This point is addressed below using U concentration data. Last, δ13C values in speleothems can be affected by temperature-dependent fractionation, although this is not a likely factor due to the stability of temperature in the cave (fig. S2). That may not have been the case during past drier climate but that also moves δ13C values in the same direction—higher values with increased warming. Using high-resolution sampling from YOK-G, we previously demonstrated that there is an excellent match between the monthly rainfall amount from Punta Gorda rain gauge station, the closest station to our site, and the monthly YOK-G speleothem δ13C data (r = 0.99, P < 0.001) (12). There is also very good correlation between historical rainfall data from the region (within about 125-km radius of our site) and YOK-G δ13C data (fig. S4). Thus, moving forward, we will discuss δ13C variability synonymously with rainfall variability.

Fig. 2 Age model and proxy comparison for stalagmite YOK-G.

(A) U-series–based age model for speleothem YOK-G. The model was constructed using the COPRA algorithm (37), and the 2σ uncertainty limits are shown as green (top) and blue (bottom) curves. Individual sample ages (52 ages) and their errors are denoted as black bars. (B) δ13C and U concentration data for YOK-G. Note the inverted axis for U concentrations (brown line) compared with δ13C (black line). ppm, parts per million.

In addition, here, we use U concentrations to confirm δ13C as rainfall proxy. Because U has partition coefficient in speleothem aragonite that is greater than one (3.74 ± 1.13) (21), the amount of prior aragonite precipitation occurring above stalagmite YOK-G largely modulates U concentrations in the stalagmite. The excellent match between δ13C and U concentrations across the study interval (r = −0.74, P < 0.001) (Fig. 2B), two elements with different exogenic cycles, demonstrates that both proxies robustly reflect hydrologic variability, in contrast to nonhydrologic variables, such as change in type of vegetation, for example. While the correlation between δ18O and U concentrations is not as strong (r = −0.44, P < 0.001) (fig. S6) as that between δ13C and U concentrations, δ18O time series does capture the long-term changes in hydroclimate. Similarly, the δ18O annual signal is less clear than the δ13C annual signal, consistent with previous interpretations of YOK-G δ18O as reflecting both the rainfall amount and the rainfall δ18O value, which is linked to tropical cyclones (20).

The highly resolved YOK-G δ13C record enables the detection of seasonal- to millennial-scale climate variability (Fig. 2B). The early part of the record, ca. 400 to 900 CE, is characterized by large shifts in hydroclimate, including the most arid interval between 625 and 825 CE. A somewhat more humid climate was present from 900 to 1400 CE, during the Medieval Climate Anomaly, followed by an extended wet interval during the LIA and lastly by progressive drying since the late 1800s CE (Fig. 1B). The transition to the wettest interval from ~1400 CE to the late 1800s corresponds with a departure from the strong coherent variability with sites north of our study area, such as the annual-scale δ18O moisture proxy (22) from Juxtlahuaca (JUX) Cave (Fig. 3B) in southwestern Mexico (17.44°N, 99.16°W) (Fig. 1). The JUX-MX record and YOK-G δ13C data show coherent variability before 1400 CE (r = 0.67, P < 0.001) (Fig. 3A). In contrast, after the initiation of the LIA at ~1400 CE, the two records diverge markedly. The LIA in southwestern Mexico (Fig. 3, A and B) is characterized by increased aridity, whereas in the YOK-G record, this interval is characterized by high rainfall. Similarly, before 1400 CE, a very good correlation (r = 0.61, P < 0.001) exists between YOK-G δ13C and the Chaac stalagmite from Tzabnah Cave (20.73°N, 89.47°W) located near the Yucatán Peninsula’s northwestern tip, north of Belize (Figs. 1 and 3B) (23), a region currently receiving ~75% less annual rainfall than southern Belize (24) with little ITCZ influence.

Fig. 3 Comparison between YOK-G δ13C and other rainfall proxy records at the northern and southern extents of the ITCZ.

From north to south (A to C): (A) The JUX Cave (stalagmite JUX-MX) δ18O time series (light brown) from southwest Mexico (22). (B) Tzabnah Cave (Chaac stalagmite) δ18O time series (green line) from the northern Yucatan, Mexico (23). (C) Huagapo Cave δ18O time series (blue line) from the central Peruvian Andes (19). (D) Cariaco Basin Ti concentration record (purple) (4). Note that the Cariaco Ti axis is showing the opposite moisture trend (wet top). The best match for the full record is achieved between the two regions (C) that are located at the northern and southern margins of the ITCZ. Correlation statistics in each panel are for the time intervals indicated.

We also compare the YOK-G data to records from the Southern Hemisphere ITCZ domain (Fig. 3B) represented by Huagapo Cave and the Cariaco Basin rainfall proxy datasets (Fig. 3C). The Huagapo Cave speleothem record (stalagmites P00-H1 and P09-H2) (25) from the central Peruvian Andes is one of the best-dated high-resolution records of South American monsoon dynamics. The large variability observed in the record was previously partially attributed to latitudinal shifts in Atlantic ITCZ position (26). A significant positive correlation (r = 0.63, P < 0.001) exists between the YOK-G δ13C and the Huagapo Cave δ18O records throughout the CE (Fig. 3C). This is remarkable given their large separation across the equator (Fig. 1). The peak LIA increase in South American monsoon strength implied by the Haugapo Cave record has been observed elsewhere, such as in the paired microalgal and mangrove biomarker 2H/1H record from the Galápagos (Fig. 1) (6), which was attributed to a southward shift of the ITCZ (6). Similarly, variability in the Cariaco Basin sediment titanium (Ti) concentration was previously attributed to changes in hydroclimate of the onshore watershed that drains into the basin, driven by meridional shifts in ITCZ position (4), where high Ti concentrations correspond to wetter than normal climate and low Ti concentrations correspond to drier than normal climate conditions. The Cariaco Ti concentration data are strongly correlated with the YOK-G δ13C data (r = 0.72, P < 0.001, 42-year lag) (Fig. 3D) after ~1300 CE. The lead-lag issue is likely related to the larger age errors of the untuned Ti time series. The drying trend during the LIA in the Cariaco Basin watershed was attributed to southward migration of the ITCZ (4).

It is not possible to reconcile these very high-quality hydroclimate records from both sides of the equator by invoking meridional ITCZ shifts exclusively. If the LIA was characterized by a simple southward ITCZ shift, it should have led to drying in our study area, which is contrary to the observed shift to much wetter conditions. Our data, considered on their own, would suggest a northward shift of the ITCZ. All the seemingly contradictory cross-equatorial datasets are reconciled by invoking an alternative model whereby the ITCZ expands meridionally but weakens in the central equatorial core region during globally cold climates and contracts meridionally but strengthens in the core region during warm climates. This explanation is supported by recent modeling results (7) and modern rainfall data that reveal a narrowing of the Pacific ITCZ under modern warming (27).

A wavelet power spectra (28, 29) analysis of our YOK-G δ13C data resampled at a constant 0.25-year resolution (Fig. 4A) illustrates that the pre–1400 CE portion of the spectra is characterized by weak multidecadal variability but exhibits no annual-scale signal, whereas the post–1400 CE data show significant (95% confidence interval outlined in black) annual-scale variability, consistent with the dominant role of the modern annual ITCZ cycle. Broadly, we refer to the time after 1400 CE as the “ITCZ period,” during which the seasonal ITCZ excursion dominates hydroclimate at our site. Although ENSO exerts a strong influence on tropical hydroclimate, our wavelet analysis (Fig. 4A) does not show ENSO type (interannual) periodicities, especially during the ITCZ period. Previous studies have shown that the impact of ENSO on Central American rainfall on the Caribbean side can be spatially variable (17). This observation is supported by modern rainfall data; there is no statistically significant correlation between ENSO and rainfall in our study area (fig. S7). Our monthly-scale δ13C rainfall proxy data permit the reconstruction of paleoseasonality. We compare two decadal intervals, one from the pre-ITCZ period (842 to 852 CE) (Fig. 4, B1) and another from the ITCZ period (1742 to 1842 CE) (Fig. 4, B2), both sampled at about 1.3-month resolution. The pre-ITCZ period data (Fig. 4, B1) do not contain annual δ13C cycles, while the ITCZ period data (Fig. 4, B2) contain a well-developed annual signal. Strong seasonal rainfall cycles are expected in a regime dominated by the annual ITCZ cycle, as is currently the case. In contrast, the large interannual changes and very muted (or entirely absent) seasonal changes before 1400 CE suggest that the ITCZ did not reach the site of Yok Balum Cave prior ca. 1400 CE.

Fig. 4 Wavelet spectral analysis of high-resolution YOK-G δ13C data and detailed look at seasonal rainfall variability.

(A) Wavelet spectral output of the high-resolution (0.25-year resolution) YOK-G δ13C data. (B) Difference in rainfall variability between the pre-ITCZ period and ITCZ period. Both segments were samples at the same, ~1.3-month, temporal resolution (small filled circles).

The contraction and expansion of the ITCZ do not preclude relative hemispherical shifts (the mean position of the ITCZ) in response to differential hemispherical warming and cooling, as we showed previously in response to aerosol injections (12). Both modeling and hemispherical-scale rainfall data from both hemispheres support a meridional shift in the ITCZ starting in the mid-20th Century (30, 31). However, this effect is small compared to the pronounced drying at both peripheries of the tropics (Fig. 3) in response to the rapid warming in both hemispheres starting in the late 1800s (32).

A salient feature of our high-resolution data is that the Central America region is rapidly emerging from what we have called the ITCZ period that began in ~1400 CE (Fig. 2B). The bi-hemispherical relationships discussed here suggest that additional future warming is likely to lead to very dry conditions, last seen around 600 to 800 CE (Fig. 3). Our results imply that the ITCZ is not an intrinsic part of the Central American climate system but is instead a transient feature, which most recently began affecting the region in ~1400 CE and whose influence has gradually waned since the mid-1800s. It is conceivable that the observed modern rainfall decrease could be the precursor to the ITCZ’s complete abandonment of Central America in the future. Our results reconcile what seemed like contradiction between modeling results showing drier Central America with future warming (11) and paleoclimate data from the southern margin of the ITCZ, which implies wetter climate with warming (6, 10). Central American countries currently within the ITCZ’s northern margin are experiencing growing food insecurity driven partially by climate change (33), leading to social upheaval and mass migration (34). Our results would seem to suggest that these conditions will be aggravated further with future warming.


Uranium-series methods

The 234U-230Th (uranium-series) chronology was performed at the Radiogenic Isotope Laboratory, the University of New Mexico. The U-Th separation chemistry is described by Asmerom et al. (35). Part of this method was described by Ridley et al. (12) for the top 456 years. Because YOK-G is composed of aragonite, it has high uranium concentration [6.5 parts per million (ppm), on average]; thus, we were able to use small samples, between 20 and 120 mg for obtaining a precise and accurate uranium-series age control per distance in stalagmite YOK-G. An aliquot of 233U-236U-229Th high-purity mixed spiked was added to each sample. U and Th were separated using 200- to 400-mesh chloride form anion exchange resin. The U and Th isotopes were measured on a Thermo Neptune multicollector inductively coupled plasma mass spectrometer. 233U, 236U, and 232Th were measured on Faraday cups with 1012 ῼ resistors, while 235U and 238U were also measured on Faraday cups with 1011 ῼ and 1010 ῼ resistors, respectively. 234U and 230Th were measured on secondary electron multiplier (SEM) that sits behind a high-abundance filter. All measurements were performed on static mode. Standards NBL-112 and an in-house 230Th-229Th solution were analyzed several times during the run sessions to accurately establish the gain between the SEM and Faraday cups. Initial 230Th/232Th ratio was estimated from the analysis of drip water and two carbonate powders obtained from the top of YOK-G during 2004 and 2006. We were able to derive an empirical relationship between 230Th/232Th atomic ratio, such that 230Th/232Thatomic ratio(ppm) = 322.67 × 232Thppb-0.269. Half-life values from (36) were used for calculation of the dates (table S1). The accuracy and precision of the uranium-series ages were validated by Ridley et al. (12) for the top 456 years using the 14C bomb signal and annual counting of δ13C seasonal oscillations; the seasonal δ13C cycles were used as the main chronology for that publication. For the longer record discussed here, we constructed a robust age model using the computer algorithm COPRA (COnstructing Proxy Records from Age models) (37). The COPRA model ages and the annually resolved ages agree within the age uncertainties of the model ages. We have used the U-Th chronology for the entirety of the record here, including the top 456 years, to provide a continuity of method, but the U-Th chronology is extremely similar to the seasonal δ13C cycle counting chronology, although the latter is probably somewhat more precise and was used in previous publications regarding YOK-G (12, 20).

Uranium concentration data

Uranium concentration data were obtained at the Radiogenic Isotope Laboratory at the University of New Mexico. Samples, between 3 and 10 mg, were dissolved in 3% nitric acid, diluted, and spiked with 115In. The samples and standards were analyzed using a Thermo X-series II inductively coupled plasma mass spectrometer. Analytical precision (2σ) ranged between 0.2 to 6%, depending on concentration.

Stable isotope analysis

A total of 7151 stable isotope analyses were measured along the growth axis of YOK-G covering a period between 2006 and 404 CE, with an average resolution of 0.22 years, but with higher than monthly resolution across certain intervals. To get this unprecedented resolution, YOK-G was continuously milled at 0.1 mm across most of the sample and at 0.2 mm across some intervals. Stable isotope analysis for the first 456 years was reported by Ridley et al. (12), and here, we discuss the record back to 404 CE. Most of the stable isotope analyses were conducted at the Durham University using a Thermo MAT 253 isotope ratio mass spectrometer with a GasBench II (external precision of about 0.05 to 0.10‰) on-line gas preparation and introduction system. Each sample (220 to 250 μg in size) was reacted with 10 drops of orthophosphoric acid (H3PO4) under a helium (grade 5) atmosphere. The solution was left to digest at 50°C for 2 hours. In addition, 14 external and in-house standards were run with each batch (50) of samples. In total, 5644 stable isotope samples were run at Durham, and the others were obtained at the Yale University (1132 samples) and the University of New Mexico (375 samples) following similar procedures as at Durham. Interlab comparisons were facilitated using replicate samples, alternating samples, and standards. Additional comparison analyses were performed on a Kiel IV automated carbonate preparation device via phosphoric acid digestion connected to a Thermo Delta V Plus mass spectrometer at the Las Vegas Isotope Science Laboratory at the University of Nevada, Las Vegas. The replicate δ13C results agreed within their respective internal precisions (see fig. S3).

Growth rate inference

In addition to the isotopic and elemental proxies, stalagmite growth rate can be an indicator of hydroclimate (12, 38), although the response depends on the moisture regime. In moisture-limited regimes, growth thickness often correlates positively with rainfall amount (38, 39), whereas in very humid regimes, the correlation can be inverse (40). Consequently, the growth rate response to rainfall can vary considerably on short time scales, depending on drip water ionic strength changes. We were able to infer rough growth rate using the distances between the model age points obtained using COPRA (37), as shown in fig. S6. In the case of YOK-G, significant negative correlations exists between δ13C and growth thickness (r = −0.64, P < 0.001) (fig. S6), suggesting that more rainfall produces thicker laminae. Remarkable correlations between our δ13C time series and other regional and global climate proxies discussed here provide additional strong support for their use as hydroclimate proxies.

Age model

A robust age model with fully propagated errors was constructed using the program COPRA (37) that uses a Monte Carlo simulation and a translation procedure that allows for calculation of proxy time series age uncertainties from radiometric date uncertainties. Annual geochemical layer counting of the top 365 mm (1550 to 1983 CE) shows excellent agreement with the uranium-series ages (12). Although the chronology used for the earlier shorter studies (12, 20, 21) was based on the geochemical layer count, here, we used exclusively a uranium-series chronology to maintain the same dating technique throughout the entire extended record.

Data treatment and statistics

Smoothing, interpolation, and correlation were performed in MATLAB. The YOK-G δ13C data have a mean resolution of 0.22 years, whereas the uranium concentration data have a mean resolution of ~3 years. A number of other records have different resolutions. For the correlation work, the data were interpolated to the lower resolutions of two records. The significance of the r values were determined first by getting a new effective degree of freedom (EDOF) for each dataset using the method of Bretherton et al. (41), as followsEDOF=N(1r1r2)(1+r1r2)(1)where, N is the length of the time series and r1 and r2 are the lag-one autocorrelation. Using the new reduced degree of freedom, a set of t values was calculated. The P value for each correlation coefficient shown is less than the stated threshold value for two-tailed t values. Thus, a P value between 0.005 and 0.001 is shown as P < 0.005. The values that are lower than 0.001, regardless of how low, are simply shown as P < 0.001.

Unlike YOK-G and JUX-MX, the Tzabnah Cave (Chaac), northern Yucatan, Mexico (23) is composed of low-uranium calcite, and consequently, the ages have large errors associated with uncertainties in the initial 230Th/232Th ratios. The ages for Tzabnah Cave were tuned, within the age uncertainties, using the accurate and precise YOK-G age model. The uranium-series ages for the three speleothems compared in Fig. 3 (A and B) (YOK-G, JUX-MX, and Tzabnah Cave) were performed in the first author’s lab. The data were tuned within the age uncertainties using the accurate and precise YOK-G age model (shown in Fig. 3B). Wavelet analysis was performed using a MATLAB program by Grinsted et al. (28), following the initial formulation by Torrence and Compo (29).


Supplementary material for this article is available at

Fig. S1. Location map.

Fig. S2. Temperature and relative humidity data.

Fig. S3. Results from duplicate analyses using Kiel and GasBench inlet systems.

Fig. S4. δ13C versus rainy season (June to October) rainfall data.

Fig. S5. Comparison between isotope data and growth rate.

Fig. S6. δ18O and U concentration data for YOK-G.

Fig. S7. Annual Oceanic Nino Index versus rainy season monthly average rainfall data.

Table S1. Uranium-series chronology data.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank the Belize Institute of Archeology for research permit (2006-2008) to K.M.P. and access to the Cave by the land owner, Santa Cruz village and the UKAA. Funding: This research was supported by grants from the NSF: NSF-BCS 0620445 (to K.M.P.), NSF-HSD 0827305 (to K.M.P. and Y.A.), NSF-EAR-0326902 (to Y.A.), NSF-HSD 0827312 (to D.J.K.); and the Alphawood Foundation (2009-2014 to K.M.P.); and the European Research Council (ERC240167 to J.U.L.B.). Meteorological records provided courtesy of HYDROMET, Government of Belize. Author contributions: Y.A. wrote the draft; V.J.P., Y.A., and V.V.A. were responsible for the uranium-series chronology; Y.A. and V.J.P. were responsible for the trace element work; and J.U.L.B., H.E.R., and C.G.M. were responsible for the stable isotope analysis. Y.A. did the statistical and spectral analysis. D.J.K., K.M.P., Y.A., and S.F.M.B. conceived the original project. K.M.P., S.F.M.B., J.U.L.B., and H.E.R. were responsible for the field work and logistics. All authors contributed to refining the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article