Quantification of drought during the collapse of the classic Maya civilization

More than two decades ago, a sediment core from Lake Chichancanab (Yucatán Peninsula, Mexico; fig. S1) provided the first physical evidence of a temporal correlation between drought and the sociopolitical transformation of the Classic Maya civilization during the Terminal Classic Period (TCP) (1). The presence of gypsum horizons and a concomitant increase in the oxygen isotope ratio (¹⁸O/¹⁶O) in shells of ostracods and gastropods suggested the TCP was among the driest periods of the Holocene in northern Yucatán. Paleoclimate records produced subsequently provided additional evidence for drought during the TCP (2–9), but the magnitude of hydro-climate change and its influence on Maya agricultural and sociopolitical systems remains controversial (10). The qualitative nature of most climate proxy archives, combined with dating uncertainties, has prevented detailed assessment of the relationship between past climate and cultural changes (10–12).

Recent attempts to quantify estimates of past changes in rainfall amount and assess the impact on ancient Maya agriculture have used isotopes of either oxygen (δ¹⁸O) (6, 13) or hydrogen (δD) (9–11). No study to date has combined the two isotope systems, because the materials used for analysis (i.e., carbonates and leaf waxes) preclude simultaneous measurement of the multiple isotopologs of water. Combined analysis of δ¹⁸O, δ¹⁷O, and δD is a powerful method to estimate past hydrologic changes quantitatively because hydrogen and triple oxygen isotopes each undergo slightly different fractionation during evaporation, leading to changes in the derived d-excess [δD – (8 × δ¹⁸O)] and ¹⁷O-excess [ln(δ¹⁷O + 1) − 0.528 ln(δ¹⁸O + 1)] parameters (14–19). In an effectively closed hydrological basin such as Lake Chichancanab, the primary controls on the isotopic fractionation of lake water during evaporation include the fractional loss of precipitation to evaporation (P/E), normalized relative humidity (RH_n), temperature, and changes in the precipitation source (1–3, 14). The value of d-excess is largely dependent on RH_n and temperature, whereas ¹⁷O-excess is controlled mainly by RH_n (14–19). Because the predicted trends of d-excess and ¹⁷O-excess in evaporating waters display different responses to climate variables, they can be evaluated individually using an iterative model (20).

We took advantage of the benefits of using all isotopologs of water and their derived parameters (d-excess and ¹⁷O-excess) by measuring triple oxygen and hydrogen isotopes in the hydration water of gypsum (CaSO₄·2H₂O) in sediment cores from Lake Chichancanab (fig. S2) (3). Today, the lake water is near saturation for gypsum; during past periods of drier climate, when the lake volume shrank, gypsum precipitated from the lake water and was preserved as distinct layers within the accumulating sediments (1–3). When gypsum forms, water molecules are incorporated directly into its crystalline structure, and this “gypsum hydration water” (GHW) records the isotopic composition of the parent fluid, with known isotopic fractionations (14, 17, 21–26). Unlike oxygen isotope fractionation during formation of carbonate minerals (27, 28), fractionation during gypsum crystallization is practically independent of temperature (24) or biological and kinetic (non-equilibrium) effects (17). Additionally, isotopes of GHW that are measured in the sedimented gypsum inherently record the driest periods, offering a distinct advantage over other traditional climate archives such as speleothems or mollusk shells, which may fail to register peak drought conditions because of growth hiatuses. Absolute differences in δ¹⁸O, δD, ¹⁷O-excess, and d-excess values between modern and paleo–lake water provide an estimate of differences between the lake hydrologic budget during the TCP and today (Fig. 1). Results were evaluated using a numerical isotope mass balance model that must satisfy all isotope variables (20) (fig. S3), and thus provides a more robust constraint on past hydrology than does modeling δ¹⁸O or δD alone.

The modern climate around Lake Chichancanab is characterized by a mean annual precipitation of ~1200 mm, a mean annual surface water temperature of ~26°C, and a net annual water deficit of 300 to 400 mm/year (3, 22). Large changes in precipitation and RH_n occur between the dry season (November to May) and the rainy season (June to October) (13, 29). Measured δ¹⁸O and δD of precipitation and groundwater samples from the Yucatán Peninsula, collected from 1994 to 2010, define a local meteoric water line (LMWL) with a slope of 7.7 (Fig. 2). Evaporation enriches the lake in the heavier isotopes of oxygen and hydrogen in water [2.6 per mil (‰) < δ¹⁸O < 3.8‰ and 10.1‰ < δD < 17.2‰], evolving along an evaporative line defined by δD = (5.1 × δ¹⁸O) – 3.1. This evaporation line intersects the LMWL at δ¹⁸O = –4.7(±1.2)‰ and δD = –27.5(±10.7)‰, which is within error of the mean oxygen and hydrogen isotope values recorded in local rivers and groundwater from the International Atomic Energy Agency’s regional Global Network of Isotopes in Precipitation stations (δ¹⁸O = –4.1‰, δD = –24.3‰) (29) and this study (δ¹⁸O = –4.0‰, δD = –23.5‰).

The gypsum deposited during the droughts of the Terminal Classic and early Postclassic periods was used to calculate δ¹⁸O, δ¹⁷O, and δD values of the paleo–lake water, which ranged from 3.6‰ to 4.9‰ for δ¹⁸O, 1.9‰ to 2.5‰ for δ¹⁷O, and 13.7‰ to 18.8‰ for δD (Fig. 1). Mean values of the paleo–lake waters (δ¹⁸O = 4.2‰, δ¹⁷O = 2.2‰, δD = 16.4‰) during drought episodes are significantly greater than modern lake values (δ¹⁸O = 3.1‰, δ¹⁷O = 1.6 ‰, δD = 12.7‰). Age uncertainty associated with the lake record and with periods of gypsum precipitation was calculated using Bayesian age-depth analysis of radiocarbon ages obtained from the sediment cores (3) (Fig. 1). We found high probabilities of drought occurring specifically during the onset (~750 to ~850 CE) and the end (~950 to ~1050 CE) of the TCP (P > 0.85 and P > 0.95, respectively) (20). Multiple proxy climate records across the Maya Lowlands also provide evidence of drought synchronicity, with only slight temporal variations across the region (10).

To estimate quantitatively the magnitude of drought during the TCP, we used a transient model that explicitly simulates the evolution of the isotopic and chemical composition of the lake water, including the gypsum flux to the lake sediments (fig. S3). The modeled gypsum flux can be compared to observed variations in the gypsum content of the sediments, as expressed by variations in sediment bulk density (3). Changes in lake surface area/volume ratio were obtained from the lake bathymetry (fig. S4). The model was run at submonthly resolution in a series of millennial-duration experiments, forced with North American Regional Reanalysis (NARR) data for local precipitation and RH_n. We first tested the model using the climate forcing across the modern sampling period from 1994 to 2010 (fig. S5). It successfully reproduced the mean of modern isotope data, with insignificant gypsum precipitation. This time interval, which was fortuitously one of the driest of recent decades, was then used as the baseline for comparison to paleo-simulations.

To provide scenarios that are directly comparable to the GHW data, we performed long transient simulations in which rainfall and RH_n were reduced by variable amounts to simulate a series of multidecadal-scale droughts. The use of a model allows us to compare, directly and quantitatively, climate conditions that affect the modern lake with conditions that would plausibly lead to drought. First, only the intervals over which the model produced gypsum deposition (modeled sediment density >1.1 g/cm³) were selected. The periods of modeled gypsum accumulation were then aggregated into drought conditions for a given scenario via two pathways: (i) All model variables were averaged across all of the droughts, and (ii) probability density functions (PDFs) were constructed incorporating the variability within and between each decadal-length drought. Data-consistent scenarios were then selected by excluding those model runs that fell outside the 1σ range of the isotope data and where, on average, the model failed to produce significant gypsum accumulation (cutoff of average density <1.2 g/cm³ based on 1σ range; fig. S6). Two possible scenarios were tested subsequently: (i) a reduction in precipitation with accompanying shifts in the isotopic composition of rainwater (i.e., the amount effect), and (ii) a reduction in precipitation with accompanying decreases in RH_n (Fig. 3).

In the first scenario, precipitation δ¹⁸O was reduced with an increase in rainfall according to the amount-effect relationship (i.e., δ¹⁸O_{precipitation}/Δprecipitation_volume = –0.0121‰/mm; fig. S7) with associated changes in δD and δ¹⁷O that track the global MWL (i.e., no changes in d-excess or ¹⁷O-excess). No scenarios with these assumptions are able to reproduce the relationship among δ¹⁸O, d-excess, and ¹⁷O-excess observed in the data. If the constraints provided by d-excess and ¹⁷O-excess are removed and only δ¹⁸O and gypsum precipitation are used, our model permits reductions in precipitation that average 50% over all drought intervals (Fig. 3, blue lines). This estimate is in broad agreement with previous work that relied on carbonate δ¹⁸O-derived precipitation estimates (using the local amount effect), which predicted reductions of up to 40% (6, 13). Our greater estimate of 50% is in part a consequence of the peak drought δ¹⁸O values recorded by gypsum, as well as the integration of simulated gypsum formation and true lake bathymetry in the model. Crucially, however, the added information from the d-excess and ¹⁷O-excess data suggests that multidecadal shifts in the δ¹⁸O of precipitation (caused by the amount effect) were not the dominant factor that affected the isotope budget of Lake Chichancanab during the TCP.

In the second scenario, we reduced precipitation without changes in the δ¹⁸O of precipitation, but instead with concurrent changes in RH_n. In this case, we observed excellent agreement in the modeled evolution of all isotopic data, with increases in δ¹⁸O accompanied by decreases in d-excess and ¹⁷O-excess (Fig. 3, red lines). This analysis yielded plausible scenarios of precipitation reduction that average 47% across all droughts (with a 1σ level of 41 to 54%) accompanied by RH_n reductions of 4% (1σ level of 2 to 7%). This result provides a robust, quantitative estimate of the mean annual hydrological conditions of the combined drought periods during the TCP at Lake Chichancanab.

Although the time evolution of our model is not a direct reconstruction of climate conditions, the model permits heterogeneity within and between each decade-long drought. The ±1σ range determined from the PDFs indicates that the precipitation reduction could vary from 20 to 70% throughout the modeled droughts (Fig. 3). This variability represents the transition into and out of drought phases and demonstrates that the severity of the droughts could be intense (up to a 70% reduction in precipitation) while maintaining the isotope balance and without desiccating the lake. Although variability in the seasonal delivery of rainfall (or lack thereof) is difficult to constrain because the residence time of the lake water is greater than an annual cycle, our results provide quantitative estimates for the total annual reduction in the water available to the ancient Maya for agricultural and domestic use. Note that recorded Colonial-period accounts of later droughts (e.g., 1535–1560 and 1765–1773), during which high mortality, famines, and population displacement were reported (30), are not manifest as intervals of gypsum precipitation in Lake Chichancanab. The lack of gypsum formation is likely a result of shorter duration and/or lower severity of these droughts, providing further evidence that the TCP was an unusually dry period for the Holocene on the Yucatán Peninsula.

Using triple oxygen and hydrogen isotope data to independently deconvolve climate variables of precipitation, RH_n, and the amount effect, we constrained the changing hydrological conditions at Lake Chichancanab. This approach provides a substantial advance over previous attempts to estimate the magnitude of rainfall reduction during the TCP droughts [e.g., (6, 13)]. Furthermore, these quantitative estimates of past rainfall and RH_n can serve as input variables in crop models, thereby clarifying how drought affected agriculture (e.g., maize production) in the northern Maya Lowlands during the TCP (12).