Northeast Yucatan hurricane activity during the Maya Classic and Postclassic periods

The collapse of the Maya civilization in the late 1st/early 2nd millennium CE has been attributed to multiple internal and external causes including overpopulation, increased warfare, and environmental deterioration. Yet the role hurricanes may have played in the fracturing of Maya socio-political networks, site abandonment, and cultural reconfiguration remains unexplored. Here we present a 2200 yearlong hurricane record developed from sediment recovered from a flooded cenote on the northeastern Yucatan peninsula. The sediment archive contains fine grain autogenic carbonate interspersed with anomalous deposits of coarse carbonate material that we interpret as evidence of local hurricane activity. This interpretation is supported by the correlation between the multi-decadal distribution of recent coarse beds and the temporal distribution of modern regional landfalling storms. In total, this record allows us to reconstruct the variable hurricane conditions impacting the northern lowland Maya during the Late Preclassic, Classic, and Postclassic Periods. Strikingly, persistent above-average hurricane frequency between ~ 700 and 1450 CE encompasses the Maya Terminal Classic Phase, the declines of Chichén Itza, Cobá, and subsequent rise and fall of the Mayapán Confederacy. This suggests that hurricanes may have posed an additional environmental stressor necessary of consideration when examining the Postclassic transformation of northern Maya polities.

3 C, weighed, and then combusted for 4.5 hours at 500º C to remove any remaining organic matter. Post-ignition masses are used to quantify downcore changes in coarse grain particle deposition, and are expressed in mass per unit of volume (D>63 μm mg/cm 3 ). Ratios of total organic carbon to total nitrogen (C:N) were measured on downcore bulk sediment samples to characterize any past changes in organic matter source [6]. Sample spacing was approximately 30 to 60 cm, however sediment was sampled at higher spatial resolutions around coarse deposits or sedimentary features. For further discussion see Supplemental Text S4.

S1.3. Radiocarbon Dating and Age Model Development
Age-depth models were developed from plant fragments (e.g., seeds and leaves, n = 13) that were radiocarbon dated at the National Ocean Sciences Accelerator Mass Spectrometry (NOSAMS) facility at the Woods Hole Oceanographic Institution. Terrestrial plant macrofossils provide the best age estimate for a given stratigraphic interval since plants from terrestrial environments incorporate radiocarbon exclusively from the atmospheric carbon pool. Conventional radiocarbon results arising from plant macrofossils were calibrated into calendar years using IntCal20 [7]. Downcore age models for each core were developed using a Bayesian statistical approach within the Bacon v2.5 software package operated in the R computing environment v4.0.3 [8]. Twenty-four radiocarbon dates were initially extracted from the two cores, however 11 were excluded from the age model since they did not depict a coherent age-depth relationship consistent with their stratigraphic order.
All radiocarbon measurements on organic matter fragments provided calibrated age results that are less than ~2200 years old, suggesting the recovered successions accumulated during the Common Era (CE). However, while 24 radiocarbon dates were extracted from the 4 two cores (15 from C1, and 9 from C2; Figure 2), only 13 followed a stratigraphically consistent chronology. The remaining 11 were all older than could be reasonably expected for the depth from which they were recovered. Anomalously old ages can sometimes be attributed to reworked older material (e.g., the 760 ± 15 BCE date from 1050 cm in C2), or material that is impacted by hardwater effects. While care was taken to ensure that only terrestrial plant remains were selected for AMS dating, not all the material extracted for this purpose could be positively identified as derived from the terrestrial surface (e.g., seeds, leaves). As such, it is possible that unknown plant fragments extracted from the upper part of the cores came from submerged aquatic vegetation and were impacted by a hardwater effect. This is supported by the more enriched δ 13 Corg values of radiocarbon results near the coretops ( Figure S13), which may be more indicative of an aquatic carbon source [9]. Six of the upper inconsistent dates fall between 100 ± 30 CE and 600 ± 40 BCE, which makes them ~2000 years older than the estimated ages for their depths based on the position of younger age results.
It is important to note that ten of the 11 inconsistent dates were extracted from within Phase 1 of the two cores (seven from C1, and three from C2). While these dates may have been older reworked material, it also possible that these upper dates were impacted by hardwater effects from in situ produced organic material. To ensure chronological parity of the records, stratigraphic tie-points were selected by tracing clear sedimentary contacts between cores 1 and 2 which were assumed to be contemporaneous. Two tie-points were identified within C1 and C2 bracketing a prominent sedimentary feature clearly discernable in each core (Figure 2). This feature presents as a slump in C1 (identified as Sl-6) and a corresponding hiatus in C2 (Figure 2, S6). The uppermost tie point (H1; Figure 2) was extracted from a coarse bed directly above the slump in C1 (304 cm depth in core). Omitting the tie points produced an asynchronous chronology between C1 and C2 for features that are stratigraphically synchronous.

S1.4. Event-bed Identification
Given the lack of classic riverine inputs to the study site, abrupt increases in coarse grained deposition are most likely related to abrupt re-suspension and initial settling-out of coarse sedimentary particles following an increase in wind and wave energy from an intense storm event (e.g., Brown, Reinhardt [2]), such as a hurricane or winter storm. Sediment samples from the shallow fringing margin of the lagoon were predominately sand sized (> 63μm diameter) particles indicating a likely sediment source during storm passage.
Hurricane force winds traveling from the northeast along the lagoon's long axis would generate significant internal waves. Coarse sediment entrained by these waves would be deposited within the cenote. Observed coarse beds are not shell beds (either ostracodes, gastropods, or bivalves), which precludes an invertebrate-based ecologic explanation. In fact, very little shell material was recovered at all in the stratigraphy. To interpret these eventdriven coarse beds they must be considered separately from the background autogenic carbonate-driven sedimentation. Identifying the likely background sedimentary signal required first collapsing all contiguous anomalous deposits of assumed instantaneous deposition greater than 5 cm in thickness (Figure 2, S6, S7) to single 1 cm thick beds. Next, a five-point running mean was used as a first pass filter to identify and remove obvious peaks from coarse fraction data. Peaks were omitted from the record if they exceeded the mean value of that five-point running average. Resulting gaps were then linearly interpolated before the record was again smoothed using an 11-point moving mean to identify the background sediment signal and minimize decadal-scale variability in autogenic carbonate 7 deposition. This background signal was subtracted from the raw grain size data to generate the coarse-anomaly record.
While the most significant coarse deposits likely resulted from anomalously high depositional energies, it is unlikely that all coarse beds exceeding background sedimentation represent clear evidence of near-to-direct striking storm events. Minor variations in coarse particle deposition may result from internally driven basin processes (i.e. biogenous sediment such as calcite tubules that may develop around aquatic vegetation). To isolate coarse deposits indicative of near-to-direct striking storm events we first calculated the standard deviations of the coarse-only signal using an 11-point moving window. Deposits were deemed significant event beds if their coarse fraction value surpassed the 95 th percentile for that interval and if that value exceeded the mean coarse fraction of the two adjacent intervals by 50%. This second step was required to avoid overcounting in intervals with highly variable background sedimentation that could not be completely accounted for (such as through large scale sedimentary shifts). The record was then manually adjusted in accordance with the following criteria. If two or three successive 1 cm intervals were categorized as significant, the depth of the maximum value within that couplet or triplet was assigned to the event and the other interval(s) were grouped with that event. For continuous coarse deposits exceeding 3 cm in width, peaks separated by at least 2 cm of decreasing values were assumed to be independent events. This method of distinguishing event beds differs from those used in other coastal karst environments [11,12]. The long-term decrease in mean event bed density within Cenote Muyil meant that those methods were inappropriate as they assume a degree of consistency among all event beds within a sediment archive.
8 Event bed count were summed over a 100-year moving window to quantify centennial-scale event distribution. The separate frequency counts were each scaled to their standard z-score which reports how many standard deviations each observation is from the mean. Finally, the separate standardized records were combined into a single composite record for the cenote. Individual coarse beds may not be recorded in both cores due to issues of internal basin geometry and/or the energy and direction of the depositional event. Deposits present in one core but not the other may indicate under-sampling within an individual core or over-sensitivity (and thus over-sampling) within the other. To resolve this, the standardized records for the overlapping portions of the two cores were averaged into a single record showing standard deviations around the mean return interval.

S.2. Slump-like Deposits
Of the 15 (11) slumps present in C1 (C2), nine (five) display clear fining upward sequences. These thicker deposits may have been produced by slow moving tropical storms driving continuous mobilization of the grains fringing the edge of the cenote. However, gravity driven sediment flow resulting from slope failure on the sides of the hole cannot be ruled out. Kenter [13] found that slopes dominated by fine-grained carbonate sediments have angles of repose rarely surpassing 15°. The northwest slope of Cenote Muyil is currently ~10° ( Figure S2) and may not require an external mechanism to initiate a slope failure. indicate an initial passing storm followed by prolonged basin-wide sloshing and sediment resuspension and transport [16].
Seiches are oscillatory waves that may form within semi-enclosed basins [17,18]. A persistent standing wave oriented along the transverse axis of the lagoon could potentially generate the sustained energy needed to produce these thicker deposits. Seiches can result from earthquakes or tangential wind-stress and may mobilize or resuspend bottom sediments for deposition in the cenote. Additionally, lake seiches may be produced by hurricane driven winds [19,20], though the development of the seiche may depend on the orientation of the passing storm and not all storms may generate seiches within the lagoon. First order wave height calculations [21] show that in a fetch-limited basin such as Chunyaxché (7 km along the long axis) significant wave heights may exceed 3.6 m during wind conditions comparable to a Category 3 hurricane strike (sustained winds > 96 kt). Wave heights may reach 6.0 m when exposed to sustained Category 5 wind speeds such as those generated by Hurricane Gilbert in 1988.
The disparity in slump thickness and shape between the two cores may be the result of core location within the hole and the fining of the deposits across the basin. C1 is closer to the center of the cenote and is proximal to the gentler slope (~6°) of the southeastern sides ( Figure S2) and C2 is closer to the steeper northwest side of the hole (~10°). Thicker deposits in C1 may be indicative of sediment transported from the east or south of the hole that fine toward the northwest. The possible density driven flow that produced Sl-6 may have originated on the steeper northwest slope closer to C2. Slope failure on the steeper side of the hole may have been more likely to scour the C2 location and transport the material towards C1.

S.3. Sedimentary Phases
The mechanisms behind the distinct changes in mean coarse bed density occurring within the three sedimentary phases are unknown. However, they may be related to baseline changes in depositional energy (i.e., stronger storms) or changes to sediment supply within the lagoon system. While we cannot discount changes to past storm intensity, we may consider that changing environmental conditions over the centuries could have altered the adjacent mangrove wetland surrounding Laguna Chunyaxché or promoted an increase in fine grain marl sedimentation directly through internal basin process. Considering the oligotrophic and phosphorus-limited nature of the ground water within the Sian Ka'an reserve [22], it is unlikely that the wetland underwent periods of continued expansion, but persistent aridity during the Early, Classic, and Postclassic periods may have encouraged growth of the more drought tolerant fringing mangrove species C. erectus [23] which favors the supratidal zone. A return to generally wetter conditions in the late Postclassic [24] would have encouraged the growth and dominance of R. mangle, which is common in the wetland today. Continued dominance of subtidal mangrove species would have increased sediment retention along the lagoon edge, inhibit the re-entrainment of settled particles, and limited coarse particle mobilization [23,[25][26][27][28]. Additionally, mangroves have a baffling effect on local flow which may inhibit the transport of coarser grains during periods of increased hydrodynamic activity. Mangrove driven siltation was also identified as the probable cause for the narrowing and infilling of lagoonal waterways at the coastal Mayan site of Vista Alegre in the northeast Yucatan [29,30].
Reductions in precipitation would not have driven decreases in lake levels through evaporative processes since surface water elevation is mediated by local sea level [31,32] and the lagoon is replenished through the subsurface aquifer. Regional sea levels have been relatively stable during the Common Era and exhibited only gradual rise [33], though even modest increases in water level may alter the sensitivity of a basin and promote local changes in sedimentation. Increased mixing of thermally distinct groundwater masses (perhaps driven by minimal sea level change) may have resulted in CO2 degassing of the alkaline and calcium rich waters. This process would have continually encouraged fine-grain lacustrine marl precipitation within the basin. If this process was driven by gradual sea-level rise, then we would expect a continuous progression of finer background sediment over the Common Era such as that observed in the cores. It is likely that a combination of these factors (increased mixing of water masses, changes to shoreline vegetation, and increased sediment baffling) dampened coarse grain transport over time, with the most pronounced change occurring ~1450 CE.
Further evidence for increased lagoon siltation was observed in the archaeological record [34]. Core 1 and Core 2 each exhibit changes in the ratio of atomic carbon to nitrogen over time that are broadly consistent with the three distinct sedimentary phases ( Figure S14).
Measurements taken from C2-Phase 3 (n = 6) have a peak C:N value of 13.35 and a mean of 12.30 (Table S2) These results indicate that organic material in the lagoon sediments was primarily aquatic in origin with limited terrestrial input. Figure S1. Comparison of modern rainfall totals and major hurricane passage within 100 km of Cenote Muyil. Multi-decadal (20-year) relationship between wet season rainfall (blue) and hurricane frequency (fuchsia) shows a low correlation. This indicates that time-averaged records of rainfall are insufficient proxies for tropical cyclone frequency. Monthly total rainfall data was extracted from 11 0.5° x 0.5° Global Precipitation Climatology Centre v2018 [35] grid cells that fell within 100 km of Cenote Muyil. Grid cells were only included if at least 50% of the land covered by the cell fell within the 100 km radius. Storm counts came from NOAA's best track historical database [36].          However, the 2020 record shows that 900 CE averaged less than six hurricanes a century. The period of least activity depicted in the 2014 reconstruction occurs ~1800 CE with fewer than 5 storms per century recorded. This coincides with one of the most active intervals in the 2020 record with ~16 storms per century recorded.   Tables   Table S1. Mean sediment density (mg/cm 3 ) per phase for beds that satisfy the criteria for event deposits.