Millennial hydrological variability in the continental northern Neotropics during Marine Isotope Stages (MISs) 3–2 (59–15 cal ka BP) inferred from sediments of Lake Petén Itzá, Guatemala

Lake Petén Itzá (Guatemala) possesses one of the longest lacustrine sediment records in the northern Neotropics, which enabled study of paleoclimate variability in the region during the last ∼400 000 years. We used geochemical (Ti, $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ and $\mathrm{Mn}/\mathrm{Fe}$) and mineralogical (carbonates, gypsum, quartz, clay) data from sediment core PI-2 to infer past changes in runoff, lake evaporation, organic matter sources and redox conditions in the water column, caused by hydrological changes in the northern Neotropics during Marine Isotope Stages (MISs) 3–2. From 59 to 39 cal ka BP climate conditions were relatively wet, and the lake was marked by higher primary productivity and anoxic bottom waters. This wet environment was interrupted for two periods of possible low water level at 52 and 46 cal ka BP, when our data suggest higher evaporation, high terrestrial organic matter input and persistent oxic conditions. Between 39 and 23 cal ka BP, evaporation and input of terrestrial organic matter increased considerably, lake level declined, and lake bottom waters generally became oxic. These conditions reversed during the Last Glacial Maximum (23.5–18.0 cal ka BP), when runoff and lake productivity increased, and rising lake level caused bottom waters to again become anoxic. Comparison of our hydrologic proxy data with sea surface temperature anomalies between the eastern Pacific and the Caribbean suggests that changes in the intensity of the Caribbean Low-Level Jet (CLLJ) may have influenced long-term changes in runoff during MISs 3–2. Higher intensity of the CLLJ during the onset of MIS 3 and the LGM might have led to greater runoff into the lake, whereas the MIS 3–2 transition experienced a weaker CLLJ and consequently less runoff. A refined, high-resolution age–depth model for the PI-2 sediment core enabled us to identify millennial-scale Greenland interstadials (GIs) 14–2, Greenland stadials (GSs) 14–2 and Heinrich stadials (HSs) 5–1. In general, HSs and GSs were characterized by drier conditions. In contrast to GSs and HSs, GIs were characterized by greater runoff and overall wetter conditions, with the most pronounced GI peaks between 40 and 30 cal ka BP. Whereas GSs 9, 8, 7 and 6 began with abrupt increases in evaporation and ended with gradual increases in humidity, GSs 11 and 10 showed reversed patterns. The Lake Petén Itzá paleohydrology record, along with other regional paleoclimate records, led us to conclude that shifts in the position of the Intertropical Convergence Zone (ITCZ) altered moisture delivery to the lake on millennial timescales. During GSs and HSs, high evaporation from Petén Itzá (dry climate conditions) was associated with a more southerly position of the ITCZ, whereas wetter GIs prevailed during a more northerly ITCZ position. Although abrupt millennial-scale shifts in ITCZ and hydroclimate between GSs/HSs and GIs can be linked to instabilities in the Atlantic Meridional Overturning Circulation (AMOC), longer-term changes were additionally influenced by changes in atmospheric convection linked to modulations of the CLLJ in response to ΔSST between the equatorial Pacific and tropical Atlantic.

1 Introduction

The last glacial period (115–11.7 cal ka BP) was marked by strong millennial-scale climate variability in the high-latitude Northern Hemisphere. Marine Isotope Stage (MIS) 3 (59–29 cal ka BP) and the last deglaciation (18–11.7 cal ka BP) were dominated by abrupt climate events (e.g., Heinrich, 1988; Lisiecki and Stern, 2016; Bradley and Diaz, 2021). In the Greenland ice core, these are referred to as Dansgaard–Oeschger (D–O) events (Johnsen et al., 1992) and consist of rapid oscillations between cold Greenland stadials (GSs) and warm Greenland interstadials (GIs). During some of the most extreme GSs, North Atlantic marine sediment cores were found to contain layers of ice-rafted debris (IRD) rich in detrital carbonate derived from Paleozoic bedrock that underlies the Hudson Strait (Heinrich, 1988; Broecker et al., 1992; Hemming, 2004). These so-called “Heinrich events” (Heinrich, 1988; Broecker et al., 1992) were attributed to massive discharges of icebergs from the Laurentide Ice Sheet to the North Atlantic via the Hudson Strait and in some cases to a European origin (e.g., Grousset et al., 2000). GSs that contain Heinrich events are referred to as Heinrich stadials (HSs). Freshwater forcing associated with Heinrich events weakened the Atlantic Meridional Overturning Circulation (AMOC), and the Intertropical Convergence Zone (ITCZ) moved south, resulting in periods of abrupt, substantial rainfall changes in the tropics and subtropics, referred to as “tropical hydroclimate events” (Bradley and Diaz, 2021). Paleoclimate records from the northern Neotropics revealed that during those events, cool and dry conditions prevailed. For instance, such conditions were recorded in a speleothem record from Cuba by a reduction in growth rates and higher δ¹⁸O values (Warken et al., 2019).

Figure 1 Location of Lake Petén Itzá, along with regional climatological and geological information. (a) Location of Guatemala in Central America and geographic position of the lake (white star). Neotropical region (green; shapefile obtained from Löwenberg-Neto, 2014) and ITCZ positions during summer (red) and winter (purple). The extent of the CLLJ in summer (red) and winter (purple) is shown. The dashed lines represent the zonal easterly winds greater than 16 m s⁻¹ at 925 mbar level associated with the CLLJ. Service layer credits: Esri, Garmin, GEBCO, NOAA and NGDC (World Topography and Oceans Maps; Esri Canada Education and Research Group; under CC BY-NC license, 2022; scale not given). (b) Overview of geology in northern Guatemala and location of Lake Petén Itzá (black rectangle) region: Qa – Quaternary alluvium, Tsp – Tertiary upper Oligocene–Pliocene, Tic – Eocene sedimentary rocks, Tpe – Paleocene–Eocene marine sediments, Kts – Cretaceous–Tertiary clastic sediments, Ksd– Cretaceous carbonaceous rocks, CPst – Carboniferous–Permian sand conglomerate (geological information: IGN, 1970). (c) Digital elevation model, principal rivers and streams, and bathymetry of Lake Petén Itzá; location of core sites: PI-2 (red dot) and PI-6 (white dot).

Lake Petén Itzá, northern Guatemala (Fig. 1), is located in a geologically and biologically diverse region that is climatically sensitive (Macario-González et al., 2022), affected by the annual migration of the ITCZ, trade wind intensity and moisture transport from the Caribbean Sea. The lake sediments represent one of the longest (∼400 kyr) paleoclimate archives in the Neotropics (Kutterolf et al., 2016). Few paleoclimate records from the northern Neotropics are older than 100 kyr in age (e.g., Lake Chalco, central Mexico; Martínez-Abarca et al., 2021a). In addition, the relatively high sedimentation rate in Lake Petén Itzá enables high-temporal-resolution sampling, making it possible to investigate abrupt, short-duration climate changes such as GS and GI events during MISs 3–2. This has been possible only for a few paleoclimate records from the Caribbean and the Gulf of Mexico region (e.g., Warken et al., 2020a). For example, the $\mathrm{Mg}/\mathrm{Ca}$ record from a speleothem in Cuba (Santo Domingo Cave) documented relatively humid conditions during GIs and drier climate during GSs (Warken et al., 2019). GSs were identified in earlier studies of Lake Petén Itzá sediments and suggested that mean annual air temperature (MAAT) was 19–20 ^∘C, 6 ^∘C lower than the modern annual atmospheric temperature in the lowland northern Neotropics (Correa-Metrio et al., 2012). The cooler conditions were accompanied by high fire frequency around the lake (Correa-Metrio et al., 2012). Identification of dry conditions during HSs was initially based on the presence of gypsum layers (Hodell et al., 2008; Mueller et al., 2010). Subsequent work suggested an even greater decline in MAAT in the region during HSs, between 6 and 10 ^∘C (Hodell et al., 2012; Grauel et al., 2016). Evidence from an ostracod transfer function suggested Lake Petén Itzá experienced a decline in mean water temperature of 3 ^∘C during HSs (Cohuo et al., 2018; Pérez et al., 2021). In addition, stable-oxygen-isotope measures on ostracod shells revealed a relatively greater evaporation $/$ precipitation ratio during HSs, from which lower lake levels and drier climate conditions were inferred (Escobar et al., 2012). During HS 1 (18–15 cal ka BP), which followed the Last Glacial Maximum (23–18 cal ka BP) and is one of the most studied periods in the Petén Itzá record, the presence of nektobenthic ostracod species such as Cypridopsis vidua and Cypria petenensis (previously called C. okeechobei and Physocypria globula), as well as the appearance of Paracythereis opesta (previously called Limnocythere opesta), suggests a relatively cold lake with high primary productivity and lower water levels (Pérez et al., 2013). These cold and dry conditions promoted the establishment of savanna vegetation dominated by Poaceae and Acacia, and in the lake, littoral “colder-water” ostracod species such as Cytheridella ilosvayi and Paracythereis opesta were abundant (Pérez et al., 2021).

Previous paleoclimate and paleoenvironmental information about MISs 3–2 was obtained primarily from Site PI-6 (water depth 71 m), one of the seven sites drilled in Lake Petén Itzá. Site PI-2, where the water depth is 54 m, is located 3 km east of Site PI-6 and showed a higher mean long-term sedimentation rate (∼150 cm kyr⁻¹ (Site PI-2) vs. ∼117 cm kyr⁻¹ (Site PI-6); Kutterolf et al., 2016). Therefore, millennial paleoclimate information was recorded at greater temporal resolution at Site PI-2, though the record covers a shorter time span (∼59 kyr). Another advantage of Site PI-2 is that it lies in a flatter area and is thus less affected by episodic slumping that influences sedimentation rate. For this study, we used high-resolution geochemical and mineralogical data from Site PI-2 sediments to explore hydrological variability in the northern Neotropics during MISs 3–2 (59–15 cal ka BP). We compared our results from Lake Petén Itzá with other regional records to better understand potential impacts and drivers of climate changes, particularly mechanisms that caused shifts in moisture transport to the northern Neotropics.

2 Study area

2.1 Location and climate

Lake Petén Itzá is located in the northern lowlands of Petén, northern Guatemala (16^∘59^′39.90^′′ N, 89^∘49^′21.07^′′ W; 110 m a.s.l. – above sea level; Fig. 1). It has a surface area of ∼100 km² and a maximum depth of 165 m (Hodell et al., 2006). The region has a humid tropical climate (“Af” according to the Köppen classification), with mean monthly air temperatures ranging from 22 to 30 ^∘C (INSIVUMEH, 2021). Mean annual precipitation is ∼1600 mm (INSIVUMEH, 2021) but varies intra-annually, with a pronounced dry season from January to May and a wet season that extends from late May through December. The climate of the northern Neotropics, particularly in terms of precipitation, is influenced by different climate forcing mechanisms that are still debated. First, the rainy season is during the summer, when the ITCZ is in its northernmost position. Second, the Caribbean Low-Level Jet (CLLJ), a region of maximum horizontal wind speeds (up to 16 m s⁻¹) at 925 mbar level, is linked to a precipitation maximum along the Caribbean and Central American coasts mainly between June and July (Amador et al., 2000). The CLLJ is a complex system that responds to large-scale atmospheric changes, the Atlantic Warm Pool, and topographic effects and land–sea contrast in addition to variations in the sea surface temperature (SST) gradient between the eastern equatorial Pacific Ocean and Caribbean Sea (Spence et al., 2004; Wang, 2007). When the Pacific has higher temperatures than the Atlantic, the CLLJ intensifies towards the northern Caribbean, and consequently, precipitation increases in Petén (Whyte et al., 2008). Humidity can also be transported by the North American Monsoon (NAM) towards the north of Mexico from the Gulf of Mexico and the Caribbean, mainly during July and September, when atmospheric temperatures reach their maximum, and the monsoon migrates northward (Hu and Dominguez, 2015).

2.2 Geology

Lake Petén Itzá lies in a basin formed by asymmetric faulting (Fig. 1). To the north are marine carbonates of Paleocene–Eocene age, whereas to the south, surface rocks are composed of Cretaceous limestones (IGN, 1970). To the southeast, Quaternary alluvial sediments, rocks and soils dominate, the latter consisting mainly of carbonates and clays (Simmons et al., 1959; IGN, 1970; Mueller, 2009). The karst geology of the region promotes rapid infiltration of rainfall into the groundwater (Hodell et al., 2008), although a slow-flowing ephemeral stream at the southeastern end of Lake Petén Itzá transports water and sediment into the lake during the rainy season.

3 Methods

3.1 International Continental Scientific Drilling Program (ICDP) drilling campaign

During February and March 2006, multiple cores were drilled at seven locations in Lake Petén Itzá (PI-1, PI-2, PI-3, PI-4, PI-6, PI-7 and PI-9) as part of an International Continental Scientific Drilling Program (ICDP) project (Petén Itzá Scientific Drilling Program, PISDP). Maximum drilling depth, 133 m, was achieved at Site PI-7. Sedimentology, magnetic susceptibility and density data for the cores were presented by Hodell et al. (2006, 2008) and Mueller et al. (2010). Here, we focus on the sediment sequence recovered from Site PI-2. The site displayed a relatively high mean sedimentation rate (∼150 cm kyr⁻¹; Kutterolf et al., 2016) and an average core recovery of 86.3 %, which enabled climate and environmental studies at high temporal resolution. The PI-2 sequence was retrieved from a water depth of 54 m in the eastern part of the lake (16^∘59^′58.04^′′ N, 89^∘44^′41.51^′′ W; Fig. 1c). The drill site receives detrital material in the form of colluvium from the steep northern shore. Five holes were drilled at Site PI-2, the longest of which reached a sediment depth of ∼67 m. The composite sequence at Site PI-2 was divided into 11 lithostratigraphic units based on changes in facies composition (Mueller et al., 2010). For this study, we analyzed sediment samples from Units 6 to 2 (67–19 m depth), which span the interval of MISs 3–2 (Correa-Metrio et al., 2012; Escobar et al., 2012).

3.2 Age–depth model

Previous studies used slightly different chronologies for the sediment sequences from Lake Petén Itzá, with a focus on the composite section from Site PI-6. Those studies used linear interpolation between radiocarbon dates (Hodell et al., 2008; Mueller et al., 2010), projection onto the PI-6 sequence of accelerator mass spectrometry (AMS) ¹⁴C dates from other cores (Escobar et al., 2012) and ages of tephra layers (Kutterolf et al., 2016). Mays et al. (2017) constructed a Bayesian age model for Site PI-6. We similarly refined the age–depth model for PI-2 using a Bayesian model (Blaauw and Christen, 2011). In this model (PETEN 02), we used 20 radiocarbon ages from the PI-2 sequence that were reported by Mueller et al. (2010) and Escobar et al. (2012), together with four dated tephra layers identified in the PI-2 core (Kutterolf et al., 2016; Table 1). Radiocarbon ages were calibrated using the IntCal20 curve (Reimer et al., 2020). The model was established using the Bacon package (Blaauw and Christen, 2011) of the RStudio software (v. 4.2). An initial estimated mean sedimentation rate of 1 mm yr⁻¹ (100 cm kyr⁻¹) for 28 segments along the core was determined. Processes that deposit material rapidly at the bottom of the lake, such as turbidite flows and volcanic ashfalls, can confound the chronology by contributing substantial amounts of material in a short time span (Moernaut et al., 2017; Mulder et al., 2019). Hence, we excluded 31 such deposits during age modeling, including carbonaceous turbidites and tephra layers, identified using the descriptions of Mueller et al. (2010; Appendix A).

Table 1 AMS ¹⁴C dates and ages of tephra used for the age–depth model. Data with ^a are measurements of terrestrial organic matter (i.e., woody debris; Mueller et al., 2010). Sample lab IDs with ^b are tephra layers (Kutterolf et al., 2016). All radiocarbon dates were calibrated using a Bayesian model (Blaauw and Christen, 2011) and the IntCal20 curve (Reimer et al., 2020); “avg” in the lab ID column refers to the mean age obtained from two dates on the same sample.

3.3 X-ray fluorescence sediment analysis

X-ray fluorescence (XRF) measurements were carried out at the University of Minnesota, Duluth, using a Cox Analytical Itrax Core Scanner (Cr tube, 30 kV, 55 mA, 15 s exposure) with a resolution of 1 cm. Particular focus was given to titanium (Ti), aluminum (Al), calcium (Ca), iron (Fe) and manganese (Mn) because these elements are reliable indicators of runoff (Ti, Al), evaporation (Ca) and redox conditions (Fe, Mn). Data from XRF, in counts per second (cps), were processed using a centered log ratio (CLR) transformation, following Weltje et al. (2015). CLR-transformed compositional data express elemental quantities in terms of actual concentrations and incorporate uncertainties acquired during core scanning, such as sediment water content, grain size and irregularities of the sediment surface (Dunlea et al., 2020). For each XRF measurement, we divided an element (e.g., Ti_(cps)) by the geometric average of the geochemical set (Ti, Al, Ca, Fe and Mn). The natural logarithm of the ratio was then calculated (e.g., $\ln ({\mathrm{Ti}}_{\left(\mathrm{cps}\right)}/\mathrm{geometric}\mathrm{average})$). Calibration was carried out using Xelerate software (Bloemsma, 2015; Weltje et al., 2015). A principal component analysis (PCA) was carried out with CLR data (Appendix B) to distinguish the processes that control each element. The PCA results are given in Appendix B. We used these results to define the $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratio as an indicator of evaporation and the $\mathrm{Mn}/\mathrm{Fe}$ as a proxy for redox processes (Wersin et al., 1991; Yarincik et al., 2000). Simultaneously, a LOESS smoothing function was applied to the CLR data to dampen the noise signal from measurement errors or anomalous data. LOESS smoothing was carried out in RStudio (v. 4.2) (R Core Team, 2018) with a span value of 3.

3.4 X-fay diffraction for sediment mineral analysis

To document changes in sediment mineralogy, 100 samples from Site PI-2 cores were used for X-ray diffraction (XRD) analysis. Sampling was carried out at irregular intervals ranging from 0.5 to 1.0 m. Sediment samples were dried at room temperature and ground with a mortar and pestle. Measurements were conducted with a Rigaku MiniFlex 600 (15 mA/40 kV) XRD at the Department of Physical Geography, Freie Universität Berlin. Peaks obtained from 3 to 80^∘ of rotation were identified with X-Pert High Score (Version 1.0b) from Philips Analytical. The peaks were calibrated against the main quartz peak (I100) at d=0.334 nm (3.34 Å). Diffraction potential files of the International Centre for Diffraction Data, USA, were used as references for identification of the different mineral phases. The quartz and gypsum counts, the sum of montmorillonite + vermiculite (clay) and the sum of calcite + dolomite + aragonite + magnesium calcite (carbonate) were added, and the relative percentages were calculated as follows: mineral (%) = (mineral count ⋅ 100)$/\mathrm{\Sigma }$ minerals (Last, 2001).

3.5 Geochemistry of bulk sediment

Analyses of total carbon (TC), total organic carbon (TOC), total inorganic carbon (TIC) and total nitrogen (TN) were carried out on the same sediment samples used for XRD. All measurements were performed at the Department of Geological Sciences, University of Florida, Gainesville. TC and TN were determined with a Carlo Erba NA1500 CNS elemental analyzer. TIC was measured coulometrically (Engleman et al., 1985) with a UIC CO₂ coulometer 5011 (Coulometrics) coupled with an AutoMate preparation device. TOC was calculated as the difference between TC and TIC. $\mathrm{TOC}/\mathrm{TN}$ ratios are reported on a molar basis.

4 Results

4.1 Chronology

Data presented in this study cover the depth interval from 67 to 19 m in the PI-2 record, corresponding to an age range from 59 to 15 cal ka BP (Fig. 2). From the sediment surface to 51 m depth, average age uncertainty is ±0.8 kyr. The uncertainty in the chronology increases in the deeper sediment intervals, especially below 51 m depth (>46 cal ka BP), resulting in ±6.7 kyr at the bottom of the sequence (67 m sediment depth). There were two intervals during which event deposits (mainly carbonaceous turbidites) were abundant: between 67 and 48 m sediment depth (58.5–43.8 cal ka BP) and from 23 to 19 m (16.6–15.0 cal ka BP). Our age–depth model differs from the one of Kutterolf et al. (2016), mainly in the interval between 40 and 34 cal ka BP, where ages calculated by our model are as much as 2 kyr younger. This is, however, within the average uncertainty in our age model at these depths (greater than ±2.3 kyr).

Figure 2 Age–depth model for the PI-2 record: (a) chronology of the PI-2 record. Gray vertical lines indicate the location of the 32 instantaneous deposits that were removed for the modeling. Ages and their uncertainties are indicated in blue. The solid green line shows the chronology obtained by Kutterolf et al. (2016) and used by Cohuo et al. (2018, 2020). The dashed red line represents the average age, whereas the dashed black lines are the minimum and maximum range for each age. (b) Age difference between the model used in this study and the one presented by Kutterolf et al. (2016) in the interval between 59 and 15 cal ka BP. Data were calculated as follows: age difference (ka) = age of this study − age reported by Kutterolf et al. (2016).

4.2 Climate proxies

Results are presented by age and according to the lithostratigraphic units proposed by Mueller et al. (2010; Fig. 3; Appendix B). CLR values for Ti vary between −2.1 (low runoff) and +1.6 (high runoff) throughout the sequence, with an average of $-\mathrm{0.1}\pm \mathrm{0.3}$ (Fig. 3). In Units 6 and 5, Ti is characterized by fluctuating values that range between −1.6 and +0.7. The lowest average values ($\le -\mathrm{2.1}$) were observed in Unit 4. Units 3 and 2 display the highest Ti content, both with maxima of +1.1.

Figure 3 Sediment variables in the PI-2 sequence, between 59 and 15 cal ka BP. Horizontal dashed black lines indicate periods that coincide with unit boundaries described by Mueller et al. (2010). The timing of the MIS 3–2 transition and the Last Glacial Maximum (LGM) is shown on the right (Lisiecki and Raymo, 2005). Runoff and evaporation are indicated by Ti centered log ratio (CLR) values and $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios, respectively. Carbonate (CaCO₃) content is shown as percentage dry weight (% wt). Percentages of quartz, clay (montmorillonite + vermiculite) and gypsum are shown. Total organic carbon is an indicator of primary productivity, ratios of total organic carbon to total nitrogen ($\mathrm{TOC}/\mathrm{TN}$) are indicators of organic matter sources, and the CLR value of $\mathrm{Mn}/\mathrm{Fe}$ is used as a proxy for redox conditions.

CLR values of $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ range from −0.1 (low evaporation) to +6.4 (high evaporation) throughout the sequence, with an average of $+\mathrm{3.0}\pm \mathrm{1.0}$ (Fig. 3). In Units 6 and 5, $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ is variable (+0.7 to +5.9), with a mean of $+\mathrm{2.8}\pm \mathrm{0.8}$ and three pronounced peaks of +5.2 (51.9 cal ka BP, 58.9 m), +5.8 (46.1 cal ka BP, 50.8 m) and +5.4 (40.3 cal ka BP, 45.3 m). The ratio increases in Unit 4 and ranges from +1.5 to +5.9 (mean = $+\mathrm{3.6}\pm \mathrm{0.7}$). Unit 3 displays the lowest $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$, varying between −0.1 and +4.0 (mean = $+\mathrm{1.9}\pm \mathrm{0.5}$). Unit 2 is characterized by highly variable $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$, ranging between −0.1 and +6.4 (mean = $+\mathrm{3.3}\pm \mathrm{1.4}$), and shows the highest values of up to +6.4 at the top of the interval.

CLR values of $\mathrm{Mn}/\mathrm{Fe}$ ratios range between −7.4 (anoxic) and +0.2 (oxic), with an average of $-\mathrm{3.4}\pm \mathrm{0.6}$ throughout the record (Fig. 3). Unit 6 displays the lowest values of the sequence, ranging from −6.6 to −2.7. Unit 5 presents a gradual increase in $\mathrm{Mn}/\mathrm{Fe}$ ratios between 52.7 and 46.1 cal ka BP (60.0–50.8 m), reaching as high as −1.1. Units 5 and 4 show similar values, which vary from −7.4 to +0.2, with maxima of +0.2 (46.2 cal ka BP, 50.8 m), −1.6 (40.2 cal ka BP, 45.3 m), −1.7 (37.1 cal ka BP, 40.5 m) and −0.8 (29.6 cal ka BP, 34.3 m). In Units 3 and 2 values fluctuate between −7.3 and −2.6 (mean = $-\mathrm{3.8}\pm \mathrm{0.5}$).

The mineralogical composition of the Lake Petén Itzá record, determined by XRD, is dominated by carbonates, including calcite, dolomite, aragonite and magnesium calcite. About 65 % of the sediment samples contained more than 50 % carbonate, with calcite (mean 45.9±2.7 %) and dolomite (mean 13.4±1.3 %) being most common. Quartz shows the highest average values in Units 6 (3.0±1.5 %), 5 (3.2±2.0 %) and 3 (5.1±2.5 %). The clay content (montmorillonite + vermiculite) is highly variable throughout the record, ranging from 0 % to 13.7 % (mean 3.6±3.0 %). It shows a gradual decrease from Unit 6 (4.5 %) to Unit 4 (2.1 %). Unit 3 contains the highest clay content (mean 8.5±2.9 %) but decreases upward to an average of 1.9 % in Unit 2. The gypsum content is highly variable and ranges from 0 % to 100 %, with an average of 32.6±3.6 % (Fig. 3). Units 6, 5 and 3 show the lowest gypsum content of the record, with average values of 5.1 %, 14.9 % and 6.2 %, respectively. In Unit 5, maxima of 96.8 % and 77.7 % occur at 50.9 cal ka BP (56.8 m) and 46.3 cal ka BP (51.1 m), respectively. In Units 4 and 2, a large increase in gypsum content is observed, with average values of 53.9 % and 78.2 %, respectively.

CaCO₃ values obtained by XRD vary from 4.2 % to 75.0 %, with an average of 33.9±15.3 %. Units 6 and 5 show similar values of CaCO₃, oscillating between 8.7 % and 73.4 % (mean 40.7±12.4 %). A reduction in the CaCO₃ content characterizes Unit 4, with an average of 28.9±14.9 %. In Unit 3, concentration values vary from 26.3 % to 58.7 % (mean = 34.7±7.5 %). Unit 2, with an average of 20.9±17.9 %, is characterized by a slight decrease in CaCO₃ content but shows a maximum of 75.0 % at 15.5 cal ka BP.

TOC values vary between 0.2 % and 6.1 %, with an average of 2.1±1.1 % throughout the sequence (Fig. 3). In Unit 6, TOC values range from 1.1 % to 4.5 %. Unit 5 displays TOC values that range from 0.5 % to 6.1 % and decrease substantially around 47.2 cal ka BP (52.3 m). From 37.8 to 30.9 cal ka BP (42.0 to 35.0 m) in Unit 4, there are three distinct maxima, with TOC values of 3.7 % (37.3 cal ka BP, 41.1 m), 5.2 % (35.4 cal ka BP, 38.4 m) and 3.4 % (33.0 cal ka BP, 36.5 m). The TOC content in Unit 3 increases to a maximum of 3.7 %, with an average of 2.2±0.6 %. Unit 2 is characterized by lower TOC values, which vary from 0.2 % to 2.5 % (mean 0.7±0.4 %).

Molar $\mathrm{TOC}/\mathrm{TN}$ ratios range from 4 (high aquatic-organic-matter content) to 38 (high terrestrial-organic-matter content) throughout the record (mean 16.8±6.3). Unit 6 shows $\mathrm{TOC}/\mathrm{TN}$ ratios that average 10.7±3.2, with a gradual increase to 18.9 at the top of the unit (Fig. 3). Units 5 and 4 display similar $\mathrm{TOC}/\mathrm{TN}$ ratios, ranging between 4.4 and 38. Four maxima of 31.3 (51.8 cal ka BP, 58.6 m), 35.9 (46.4 cal ka BP, 51.0 m), 27.4 (40.4 cal ka BP, 45.3 m) and 38 (31.9 cal ka BP, 35.7 m) were identified, as were minima of ∼6.8 between 40.3 and 39.2 cal ka BP (45.2–43.9 m depth). At the top of Unit 4 (31.0–23.4 cal ka BP, 35.0–29.0 m depth), $\mathrm{TOC}/\mathrm{TN}$ ratios decrease gradually from 26.1 to 10.8. Units 3 and 2 have the lowest values of the entire sequence, ranging between 4.6 and 17.5, with a slight increase to 20.9 at the top of Unit 2.

5 Discussion

Discussion of the Lake Petén Itzá paleoenvironment focuses on three key objectives: (1) hydrological variability during MISs 3–2 and millennial oscillations (GSs, GIs and HSs) recorded at Site PI-2, (2) comparison of our Site PI-2 data with previously reported data from Site PI-6, and (3) identification of climate forcing mechanisms that may have provided humidity to the Petén area through transregional proxy comparison. Our mineralogical and geochemical data are presented and interpreted according to the five periods that correspond to Lithological Units 6–2 proposed by Mueller et al. (2010) (Appendix C): (a) onset of MIS 3 (Unit 6: 59–52.7 cal ka BP), (b) mid-MIS 3 (Unit 5: 52.7–39.3 cal ka BP), (c) MIS 3–2 transition (Unit 4: 39.3–23.5 cal ka BP), (d) Last Glacial Maximum (Unit 3: 23.5–18.0 cal ka BP) and (e) deglaciation (Unit 2: 18.0–15.0 cal ka BP). The lithological units were based on sedimentology and stratigraphy and verified with core-logging data (density, magnetic susceptibility) and with seismic reflection profiles throughout Lake Petén Itzá basin.

5.1 Hydrological and environmental responses during MISs 3–2

5.1.1 Onset of MIS 3 (59–52.7 cal ka BP): high runoff and low evaporation

The onset of MIS 3 corresponds to Unit 6. Ti is widely used as an indicator for runoff (Mason and Moore, 1982; Davies et al., 2015). High Ti (up to 0.7) in this unit suggests high runoff, particularly around 56 cal ka BP (Fig. 3). This is supported by concomitantly high clay and quartz content, with averages of 4.5 % and 3.0 %, respectively. Clay minerals, particularly montmorillonite and vermiculite, are generated during chemical weathering of unstable siliciclastic minerals such as plagioclase, especially in tropical regions where atmospheric humidity is high, resulting in the enrichment of quartz (Weltje, 1994; Boggs, 2006; Van De Kamp, 2010). Quaternary alluvial sediments could be the source of siliciclastic minerals (Simmons et al., 1959; Mueller, 2009). We suggest that the high clay and quartz content in Unit 6 is indicative of wet conditions, which provided abundant runoff to the lake.

Ca provides information on supersaturation and precipitation of ions in solution, triggered by high evaporation (Engstrom and Wright, 1984; Boyle, 2001). Ca in lake sediments, however, can also originate from detrital inputs in karst regions such as in Petén and from biogenic sources, e.g., shell formation by mollusks and ostracods and CO₂ withdrawal by plants for photosynthesis. Therefore, we normalized Ca values to indicators of detrital input such as Ti and Fe (Appendix B) and used the $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ as an indicator for past variability in evaporation (Mason and Moore, 1982; Yarincik et al., 2000). Low $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios (mean 2.3±0.6) and gypsum content (mean 5.0±5.1 %) at the onset of MIS 3 indicate low evaporation of lake water and more humid conditions, with only a few intermittent dry periods, at ∼55.1 and 53.8 cal ka BP (Fig. 3).

The $\mathrm{TOC}/\mathrm{TN}$ ratio has long been used as an indicator for the source of organic matter in lacustrine sediment sequences (Talbot and Johannessen, 1992; Meyers and Ishiwatari, 1995). In the interval between 59 and 53 cal ka BP, low $\mathrm{TOC}/\mathrm{TN}$ ratios (mean 10.7±3.2) indicate that the sediment organic matter is predominantly of aquatic origin (Meyers, 2003). The $\mathrm{Mn}/\mathrm{Fe}$ ratio is used to infer variations in redox conditions in lakes (Wersin et al., 1991; Naeher et al., 2013; Friedrich et al., 2014; Ortega-Guerrero et al., 2020). Lower values of the $\mathrm{Mn}/\mathrm{Fe}$ ratio generally indicate low concentration of O₂ in the water column, which is explained by more rapid reduction in Mn compared to Fe under anoxic conditions, favoring preferential Mn deposition. In contrast, higher values of the $\mathrm{Mn}/\mathrm{Fe}$ ratio suggest high levels of O₂, as Fe oxidizes faster than Mn in oxic environments. Low $\mathrm{Mn}/\mathrm{Fe}$ ratios (mean $-\mathrm{3.7}\pm \mathrm{0.3}$) indicate anoxic conditions in the lake's bottom waters. Anoxia may have been associated with a high water level (as is observed today), which likely prevented wind-driven mixing of the entire water column.

High values of magnetic susceptibility measured at Site PI-6 (Fig. 4) support our interpretation of a wet climate, compared to the less humid environment inferred for the end of MIS 3 (48–23 cal ka BP; Hodell et al., 2008). Previously reported ostracod data from Sites PI-2 (Cohuo et al., 2020) and PI-6 (Pérez et al., 2021) are consistent with substantial precipitation in the study area during the deposition of Unit 6 and are associated with high water levels, reflected by the dominance of Cypria petenensis. Moreover, the low content of terrestrial organic matter at Site PI-2 is consistent with the pollen-based inference for extensive grasslands around Lake Petén Itzá, indicating relatively little terrestrial plant biomass in the watershed (Bush et al., 2009; Fig. 4).

Figure 4 Sediment variables in cores from Sites PI-2 and PI-6 (3 km to the west; Fig. 1). Centered log ratio (CLR) values of Ti, $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ and $\mathrm{Mn}/\mathrm{Fe}$. $\mathrm{TOC}/\mathrm{TN}$ and magnetic susceptibility (Hodell et al., 2008). Pinus/Poaceae pollen and charcoal datasets (Correa-Metrio et al., 2012) as well as δ¹³C values from ostracod valves (Escobar et al., 2012). Dashed brown lines indicate unit boundaries described by Mueller et al. (2010).

Predominantly wet conditions during the onset of MIS 3, like those inferred from Lake Petén Itzá sediments, have also been reported from other sites in the northern Neotropics (Fig. 5) using geochemical records. They include lakes Chalco, central Mexico (Martínez-Abarca et al., 2021b) and Babicora, northern Mexico (Roy et al., 2013). In contrast, pollen and sedimentological information from Lake Fúquene (Colombia) indicates moderately dry conditions and a reduction in lake surface area associated with a decline in precipitation in northern South America (Groot et al., 2013). Differences in precipitation between the northern and southern limits of the American tropics were likely associated with the latitudinal position of the ITCZ (Jaeschke et al., 2007; Gibson and Peterson, 2014; Martínez-Abarca et al., 2021b; Lozano-García et al., 2022).

Figure 5 Comparison of paleoclimate and runoff/precipitation records spanning MISs 3–2. (a) Centered-log-ratio Ti record from Lake Petén Itzá (this study). (b) Reflectance data from core MD03-2621, Cariaco Basin (Deplazes et al., 2013b). (c) Standardized record of Ti obtained from ODP Site 1239 in the equatorial Pacific Ocean (Rincón-Martínez et al., 2010). (d) δ¹⁸O data from a stalagmite from Santo Tomas Cave, Cuba (Warken et al., 2019). (e) Standardized record of Ti in two sediment sequences from Lake Chalco, central Mexico (Site CHAVII obtained from Lozano-García et al., 2015; CHA08 site obtained from Martínez-Abarca et al., 2021b). (g) Standardized TiO₂ data obtained from Lake Babicora, northern Mexico (Roy et al., 2013). Original data were standardized by subtracting the mean of the data and dividing by the standard deviation.

Reflectance data from marine sediments collected in the Cariaco Basin, Venezuela (Deplazes et al., 2013a), and δ¹⁸O data from marine record MD02-2529 retrieved from the eastern equatorial Pacific (Leduc et al., 2007a) indicate a more northerly position of the ITCZ during the onset of MIS 3. Hodell et al. (2008) suggested that increases in runoff in Petén Itzá co-occurred with increases in SST in the subtropical northeast Atlantic (Site SU8118, Portugal; Bard et al., 2000). SSTs between sites in the region, however, have not yet been considered. Therefore, we compared the Ti data from Site PI-2 with the SST difference (ΔSST) calculated between Site MD02-2529 in the eastern Pacific (Leduc et al., 2007a) and ODP Site 999A in the Caribbean (Schmidt et al., 2006a; Fig. 6). It should be noted though that ΔSST yields considerable uncertainties associated with the prediction error of the quantitative estimation of SST in both records, their age–depth models and the calculation of Δ itself because of the use of different temperature proxies ($\mathrm{Mg}/\mathrm{Ca}$ planktonic foraminifera for ODP-999A and alkenone U${{}^{\prime }}^K$₃₇ for MD02-2529) and temporal resolutions. Even with these limitations, the ΔSST record, used as a CLLJ proxy for variations in the strength of the CLLJ, shows an apparent coupling of the ΔSST with the low, frequent variations recorded in our Ti data (Fig. 6). Decreases in ΔSST indicate that the Caribbean SST was lower than the value in the Pacific. This could have promoted an increase in the intensity of the CLLJ, since a cold Caribbean Sea favors expansion of the North Atlantic Subtropical High, displacing the CLLJ to the west and thereby increasing humidity in Petén Itzá (Wang, 2007; Lachniet et al., 2009; Fig. 7). During early MIS 3, ΔSST decreased, and it became drier until ∼45 cal ka BP. After that, ΔSST increased, and wetter conditions occurred around the lake. Although the strength of the CLLJ does not depend solely on the ΔSST between the eastern Pacific and the Caribbean, but also on the regional topography, the position of the ITCZ and the expansion/contraction of the North Atlantic Subtropical High (NASH; Wang, 2007), we find clear evidence for a close relationship between ΔSST and Ti and hence runoff to Lake Petén Itzá regarding low, frequent variations. This implies that in addition to changes in the ITCZ, spatiotemporal changes are linked to variations in CLLJ that steered the influx of moisture into Petén during MISs 3–2. Additionally the onset of MIS 3, a relatively warm period (Bradley and Diaz, 2021), might have been characterized by a strong NAM, which could have promoted moisture transport from the Caribbean towards northern Mexico. This was shown in records from northern Mexico (e.g., lakes Babicora and Sayulita), where wet interstadials were inferred (Metcalfe et al., 2013; Roy et al., 2013).

Figure 6 Low-frequency co-variability between runoff (Ti) and changes in the intensity of the Caribbean Low-Level Jet (CLLJ). We use the Pacific–Caribbean differences in sea surface temperature (ΔSST) as a proxy for the CLLJ. A stronger ΔSST, as a consequence of cooling in the Caribbean and/or warming in the Pacific, causes a stronger CLLJ, with more precipitation and runoff at Lake Petén Itzá. Red line: ΔSST calculated by subtracting annual SST values inferred from the marine sediment record at Site MD02-2529 in the eastern Pacific (Costa Rica; Leduc et al., 2007b) from values at ODP Site 999A, located in the Caribbean (Schmidt et al., 2006b). Site 999A SST data were obtained using $\mathrm{Mg}/\mathrm{Ca}$ ratios of the planktonic foraminifera Globigerinoides ruber, whereas the MD02-2529 dataset is based on the U${}_{\mathrm{37}}^{K^{\prime }}$. As the two records possess different temporal resolutions, we approximated the closest ages between both records and subtracted the available SST data. Gray line: Ti from Site PI-2 in Lake Petén Itzá (this study). Ti data were averaged using a LOESS smoothing function (factor = 10). The different time periods discussed in the text are linked to the stratigraphic units defined by Mueller et al. (2010).

Figure 7 Maps showing regional paleoclimate records and limnological/climate characteristics inferred for (a) the onset of MIS 3, (b) mid-MIS 3, (c) the Last Glacial Maximum (LGM) and (d) deglaciation. Blue refers to mainly humid conditions, whereas red indicates primarily dry conditions. Possible mechanisms providing moisture to the northern Neotropics are indicated. Abbreviations: OM aqua (aquatic organic matter), OM terr (terrestrial organic matter), OM mix (mixed organic matter), Mod (moderate), SST (sea surface temperature), AMOC (Atlantic Meridional Overturning Circulation), ITCZ (Intertropical Convergence Zone), CLLJ (Caribbean Low-Level Jet). Records: (1) Lake Babicora (Roy et al., 2013), (2) marine record MD02-2575 (Ziegler et al., 2008), (3) Lake Tulane (Grimm et al., 2006), (4) Bahamian speleothem record AB-DC (Arienzo et al., 2017), (5) Lake Pátzcuaro (Bradbury, 2000), (6) Lake Chalco (Martínez-Abarca et al., 2021b), (7) Lake Tecocomulco (Caballero et al., 1999), (8) Piscina de Yuriria (Holmes et al., 2016), (9) Lake Petén Itzá (this study, marked with a star), (10) Blanchard Cave (Royer et al., 2017), (11) Larga Cave (Warken et al., 2020b), (12) marine record MD02-2529 (Leduc et al., 2007a), (13) Cariaco Basin record MD03-2621 (Deplazes et al., 2013a), (14) Lake Fúquene (Groot et al., 2011), (15) marine record ODP Site 1239 (Rincón-Martínez et al., 2010).

5.1.2 Mid-MIS 3 (52.7–39.3 cal ka BP): shift to higher evaporation and lower lake levels

The mid-section of MIS 3 corresponds to Unit 5. High Ti (mean $-\mathrm{0.1}\pm \mathrm{0.3}$) and stable quartz contents (mean 3.2±2.0 %) indicate that precipitation and runoff to Lake Petén Itzá remained high during the deposition of Unit 5 but were lower than in Unit 6 (Fig. 3). An increase in TOC content (mean 2.8±1.0 %) at 51 cal ka BP indicates increased productivity, and a pronounced increase in $\mathrm{TOC}/\mathrm{TN}$ ratios (mean 18.7±5.9) indicates a relatively greater delivery of terrestrial organic matter to the lake (Meyers, 2003). The TOC content reaches maxima in two intervals, around 46.2 and 40.2 cal ka BP, coinciding with increased evaporation ($\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios of up to 5.9). These two periods correlate with peaks in oxic conditions, as suggested by maxima of the $\mathrm{Mn}/\mathrm{Fe}$ ratio (0.2 and −1.6; Fig. 3). Our data thus suggest more oxygenated bottom waters, possibly associated with lower lake levels, which resulted in greater wind-induced mixing of the water column (Yu et al., 1984; Gale et al., 2006). A short-lived increase in runoff terminates Unit 5, a result of more humid conditions that may have delivered greater nutrient loads to the lake. The latter would have favored the input of aquatic organic matter to the sediments between 40.4 and 39.3 cal ka BP. Wetter conditions and greater autochthonous production are supported by high Ti (up to 0.7) and low $\mathrm{TOC}/\mathrm{TN}$ ratios (6.7), respectively.

A decrease in runoff and lake level lowering, which interrupted the general wet climate during mid-MIS 3, were also inferred by Pérez et al. (2021) from Site PI-6 based on ostracod assemblages that reflect higher conductivity. This was explained by an increase in ion concentrations in the water column associated with higher evaporation $/$ precipitation ratios. In addition, the presence of the high-conductivity-preferring diatom Cyclotella petenensis was reported from Site PI-6 sediments between 40 and 39 cal ka BP (Paillès et al., 2018), suggesting high evaporation rates. Enhanced evaporation was associated with an increased contribution of terrestrial organic matter. This may at first appear contradictory. However, it is consistent with pollen-based support from Site PI-6 for an increase in relative abundance of non-arboreal taxa, e.g., Poaceae (grasses) during Unit 5 (Fig. 4; Correa-Metrio et al., 2012). We interpret the increase in terrestrial organic matter as not reflecting changes in the terrestrial vegetation, but rather a shift in the proximity of the plant community (grasslands) to Sites PI-2 and PI-6 because of shrinkage in the lake surface area. This was observed at Lake Tzibaná (southern Mexico) as a consequence of a 20 m decline in water level, during which vegetation colonized subaerially exposed lake bottom, and large amounts of organic matter entered the remnant basin during sporadic periods of rain (Martínez-Abarca et al., 2023).

Whereas our Lake Petén Itzá record – together with other lake and speleothem archives from the circum-Caribbean and Gulf of Mexico, such as Lake Tulane, Florida (Grimm et al., 2006); Abaco Islands, Bahamas (Arienzo et al., 2017); Santo Tomas Cave, Cuba (Warken et al., 2019); and Larga Cave, Puerto Rico (Warken et al., 2020a) – reveals moderately wet conditions during the midst of MIS 3, several paleoclimate records from northern and central Mexico, including those from Lakes Babicora (Roy et al., 2013), Tecocomulco (Caballero et al., 1999), Pátzcuaro (Bradbury, 2000) and Chalco (Ortega-Guerrero et al., 2020), indicate reduced rainfall (Fig. 5). This spatial climate contrast may have been associated with the beginning of the southward migration of the ITCZ. The ΔSST between the eastern Pacific and the Caribbean shows a slight decrease, suggesting that the CLLJ decreased in intensity (Fig. 6), thereby promoting low runoff into Lake Petén Itzá. Nevertheless, climate at Petén Itzá remained relatively humid, possibly attributed to moisture supply associated with a stronger AMOC. Indeed, it has been shown that strengthening of the AMOC may promote enhanced northward heat transport, lower cross-equatorial ΔSST and consequently wetter conditions in the tropical Northern Hemisphere (Fig. 7; Clark et al., 2001; Waelbroeck et al., 2018).

5.1.3 MIS 3–2 transition (39.3–23.5 cal ka BP): great hydrological instability

The final stage of MIS 3 and the beginning of MIS 2 correspond to Unit 4. This unit is characterized by high variability in runoff and evaporation. The decrease in Ti (mean $-\mathrm{0.3}\pm \mathrm{0.3}$) indicates a decline in runoff. Moreover, variable but higher $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios (mean 3.6±0.7) indicate an increase in evaporation relative to Unit 5 (Fig. 3). In addition, the carbonate content shows a trend similar to that of the Ti dataset. This coupling may suggest an input of detrital carbonate from the karstic bedrock in the catchment. Carbonate in lake sediments is commonly related to (1) inorganic calcite precipitation and consequently to periods of high productivity and/or evaporation and (2) high nutrient input promoting eutrophication, an increase in pH and consequently enhanced biogenic calcite formation (Dean and Megard, 1993; Dean, 1999). Mineralogical data display higher values of gypsum (mean 53.9±3.4 %), suggesting drier conditions and hence lower chemical weathering in the catchment. High and stable values of the $\mathrm{Mn}/\mathrm{Fe}$ ratio (mean $-\mathrm{3.3}\pm \mathrm{0.5}$) indicate oxic bottom waters within this unit, possibly associated with a decline in lake level and greater wind mixing. $\mathrm{TOC}/\mathrm{TN}$ ratios rapidly increase to values > 20, suggesting a relative decline in the contribution of aquatic organic matter input relative to that from terrestrial sources. After 32 cal ka BP, however, $\mathrm{TOC}/\mathrm{TN}$ ratios decreased gradually to values of ∼10, indicating a shift once again to dominance of aquatic organic matter. In addition, lower TOC values (mean 1.9±1.0 %) in Unit 4 may suggest poor preservation of sediment organic matter as a result of oxidation favored by greater water-column mixing at Site PI-2, as a result of the decline in lake level.

Lower lake levels during the MIS 3–2 transition were inferred previously using the Site PI-6 core (Cohuo et al., 2020; Pérez et al., 2021) from an increase in the relative abundance of the endemic ostracod Paracythereis opesta, a benthic species associated with shallow-water environments (<40 m depth) (Cohuo et al., 2017; Echeverría Galindo et al., 2019). In addition, Paillès et al. (2018) found abundant Cyclotella petenensis and C. cassandrae in deposits from that unit in PI-6. The two diatom species are characteristic of high-conductivity waters, which agrees with our interpretation. Lower water levels, which enabled oxygenation of the water column, are also supported by high δ¹³C values measured in the ostracod Paracythereis opesta (previously called Lymnocythere opesta) at Site PI-6 (Fig. 4; Escobar et al., 2012). The gradual decrease in runoff, associated with dry conditions inferred from Site PI-2, correlates with a decrease in magnetic susceptibility at Site PI-6 (Fig. 4). On the other hand, the high proportions of terrestrial organic matter prior to 32 cal ka BP coincides with an increase in sediment δ¹³C_TOC values in the PI-6 record, which was interpreted to reflect a greater proportion of C₄ plants in the catchment, indicative of drier climate (Mays et al., 2017). The contribution of allochthonous organic matter may also have declined because of the spread of C₄ savanna vegetation in the lake catchment, with high relative abundances of Cyperaceae and Poaceae (Fig. 4; Correa-Metrio et al., 2012).

Paleoclimate conditions different from those around Lake Petén Itzá were inferred from lakes of central and northern Mexico, where a wetter climate and less saline conditions prevailed (e.g., Lake Babicora, Chávez-Lara et al., 2012; Lake Chalco, Caballero et al., 2019; Fig. 5). This inter-regional comparison of paleoclimate records indicates that the Petén region and northern Gulf of Mexico were dominated by dry conditions during the MIS 3–2 transition. We suggest that during the transition (39.3–23.5 cal ka BP), the ITCZ was mainly located farther south, accounting for the lower precipitation and consequent reduction in runoff to Petén Itzá. The gradual trend towards drier conditions observed after 29 cal ka BP, however, resembles the trend of reduced SST in the eastern Pacific that may have lessened the intensity of the CLLJ, resulting in decreased moisture transport to the Petén region. Therefore, it is possible that much of the runoff to Lake Petén Itzá depended on the position of the CLLJ during the transition.

5.1.4 Last Glacial Maximum (23.5–18.0 cal ka BP): a sudden large increase in runoff

The Last Glacial Maximum (LGM) corresponds to Unit 3, which was studied previously using sediments from Site PI-6 (Hodell et al., 2008; Bush et al., 2009; Pérez et al., 2011; Mays et al., 2017). Ti is highest in Unit 3 and varies between −0.6 and 1.1, reflecting an abrupt increase in humidity and runoff (Fig. 3). This is also supported by the increases in clay and quartz contents (means of 8.5±2.9 % and 5.1±2.5 %, respectively). $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios were lowest in this unit (mean 1.9±0.5), indicating substantially reduced rates of evaporation (Fig. 3). $\mathrm{TOC}/\mathrm{TN}$ ratios < 10 suggest high sedimentation of aquatically derived organic matter. Lower $\mathrm{Mn}/\mathrm{Fe}$ ratios (mean $-\mathrm{3.7}\pm \mathrm{0.3}$) and generally higher TOC content (2.2±0.6 %) suggest that the deep lake was characterized by persistent bottom-water anoxia that facilitated the preservation of organic matter.

More humid conditions and higher lake levels during the deposition of Unit 3 were also inferred from sediments at Site PI-6. High magnetic susceptibility and density values suggest high detrital input (Fig. 4; Hodell et al., 2006, 2008), whereas more negative δ¹⁸O and δ¹³C values of ostracod shells indicate higher lake levels and establishment of an anoxic hypolimnion at Site PI-6 (Escobar et al., 2012; Pérez et al., 2013). Ostracod assemblage analysis revealed the presence of deep-water species (>40 m depth) such as Cypria petenensis (Pérez et al., 2021), while the presence of the diatom Discotella gabinii indicates low-conductivity and alkaline water during the LGM (Paillès et al., 2018).

Wet conditions during the LGM were also inferred from δ¹⁸O values in speleothems from the Caribbean region – including those from Blanchard Cave, Guadeloupe (Royer et al., 2017); Larga Cave, Puerto Rico (Warken et al., 2020a); and Santo Tomas Cave, Cuba (Warken et al., 2019) – as well as southwestern (Juxtlahuaca Cave, Guerrero: Lachniet et al., 2014) and northwestern Mexico (Lake Babicora: Metcalfe et al., 2002). In contrast, dry conditions were inferred from other paleoclimate records from central Mexico (e.g., La Piscina de Yuriria: Holmes et al., 2016; Lake Chalco: Lozano-García et al., 2015), the northern Gulf of Mexico (e.g., Lake Tulane, Florida: Grimm et al., 2006; Donders et al., 2011), the Colombia Basin (Lake Fúquene: Groot et al., 2011, 2013; Vriend et al., 2012) and the eastern equatorial Pacific ODP-1239 record (Rincón-Martínez et al., 2010) (Fig. 5). This is consistent with the multi-record analyses of Ramírez-Barahona and Eguiarte (2013), which showed dry conditions for western central Mexico and concomitantly wet conditions in eastern and southern Mexico. Our regional comparison using records from Mexico, the Gulf of Mexico and the Caribbean indicates that the ITCZ may have been positioned mainly over Central America during the LGM (Fig. 7). It is also possible that the decrease in atmospheric temperatures reduced the amount of humidity that entered Mexico from the Caribbean (Hu and Dominguez, 2015). Previous studies showed that during the LGM, an analog to the negative phase of the El Niño Southern–Oscillation (ENSO) persisted, which resulted in a weakening of the NAM, explaining the droughts in central and northern Mexico (Lachniet et al., 2013). The ΔSST suggests that water temperatures in the eastern Pacific were ∼1 ^∘C higher than those in the Caribbean during the LGM (Fig. 6), which may have promoted a more intense CLLJ and thus increased moisture transport to the Caribbean, including to the Petén region. This is consistent with estimates from fossil records (e.g., Trend-Staid and Prell, 2002; Kucera et al., 2005) and numerical simulations (e.g., Kitoh and Murakami, 2002; Otto-Bliesner et al., 2009).

5.1.5 Deglaciation (18.0–15.0 cal ka BP): low runoff and enhanced evaporation

At the onset of Unit 2 (18.0–16.3 cal ka BP), a decrease in Ti to 0.1 and an increase in $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios and gypsum content (averages of 3.3±1.4 % and 78.2±32.8 %, respectively) indicate reduced runoff, drier conditions and an increase in evaporation, which would have resulted in lower lake level (Fig. 3). As in Unit 4, carbonate and Ti values are coupled, indicating that input of carbonates into the lake may depend partly on runoff and subsequent precipitation. This is consistent with ostracod-based inferences from Site PI-6 that indicate lower water levels and greater conductivity, with a dominance of littoral ostracod species such as Cypridopsis vidua, Heterocypris putei and Paracythereis opesta (Díaz et al., 2017; Cohuo et al., 2018; Pérez et al., 2021). Likewise, high δ¹⁸O values in ostracod shells and gypsum hydration water indicate high evaporation (Escobar et al., 2012; Hodell et al., 2012; Grauel et al., 2016). An estimated 56 m lake level decline was reported by Anselmetti et al. (2006) and Hodell et al. (2006), based on the occurrence of a paleo-shoreline in seismic profiles. Comparatively low $\mathrm{TOC}/\mathrm{TN}$ ratios at Site PI-2 (mean 11.5±4.8) indicate a decrease in the relative contribution of organic matter from terrestrial sources. Furthermore, n-alkane distributions in the PI-6 record show reduced loading of terrestrial plant material between 18 and 16 cal ka BP (Mays et al., 2017). This may be explained by the predominance of grassland in the vicinity of the lake (Fig. 4; Bush et al., 2009), which may have contributed little to the sedimentary organic matter.

Dry climate conditions, similar to those at Petén Itzá, were also inferred from sediments in lakes of central Mexico and the Caribbean (Bradbury, 2000; Lozano-García et al., 2015; Royer et al., 2017; Caballero et al., 2019; Martínez-Abarca et al., 2019; Warken et al., 2020a), whereas records from the northern Gulf of Mexico suggest wet conditions that originated from extra-tropical fronts that moved south during the early phase of deglaciation (Ziegler et al., 2008; Donders et al., 2011; Roy et al., 2013; Arienzo et al., 2017). This can be explained by the latitudinal migration of the ITCZ to the south after the LGM and the collapse of the AMOC during HS 1, which promoted dry environments throughout central Mexico and the Caribbean (Fig. 7).

Unit 2 terminates with a sudden increase in Ti (up to 0.8) between 16.3 and 15.0 cal ka BP, indicating an increase in precipitation and runoff. A decrease in $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ (up to −0.1) suggests evaporation remained low, and low $\mathrm{TOC}/\mathrm{TN}$ ratios (<10.9) indicate a predominance of aquatic organic matter. Although dry conditions may have prevailed at the beginning of the deglaciation corresponding to the first stage of HS 1, the second part (17–15 cal ka BP) was characterized by an increase in runoff and precipitation (Hodell et al., 2012; Pérez et al., 2013). This may have been linked to a slight AMOC recovery that preceded the second stage of HS 1, indicated by lower δ¹⁸O values in foraminifera tests from the North Atlantic (McManus et al., 2004).

5.2 Millennial-scale climate variability

Whereas long-term transient climate changes have been inferred from multiple records in the region, our high-resolution geochemical records from Site PI-2 also reveal short-term variations, including abrupt shifts in hydrological proxies between 59 and 15 cal ka BP, similar to those observed in the GICC05 ice-core record from Greenland (Svensson et al., 2008). Low $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios are associated with high Ti during GIs. This suggests wet conditions and enhanced runoff. Maxima of $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ and increased gypsum content at Site PI-2 are largely contemporaneous with increases in the L^* reflectance values in marine core MD03-2621 from the Cariaco Basin, offshore Venezuela (Fig. 8), particularly during GSs 11–2 (Peterson et al., 2000; Deplazes et al., 2013a). Transport of moisture to the Petén region between 42 and 24 cal ka BP seems to have been dominated by changes in SSTs between the eastern Pacific and the Caribbean and in turn by the intensity of the CLLJ. However, variability in runoff and evaporation, suggested by our $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios during GSs 11–2 and in the Cariaco L^* reflectance record, may have been largely influenced by millennial-scale changes in the position of the ITCZ.

Figure 8 Temporal variation in climate proxies in the Petén Itzá PI-2 sediment record compared to other paleoclimate records covering the time span from 60 to 14 cal ka BP. (a) North Atlantic ice-rafted debris (IRD) stack derived from 15 individual sediment cores (Lisiecki and Stern, 2016). Oxygen isotopes (δ¹⁸O) from the GICC05 ice core from Greenland (Svensson et al., 2008). More negative values represent colder intervals used to define Greenland stadials (GSs); more positive values represent warmer periods corresponding to Greenland interstadials (GIs). (b) Reflectance L^* from Cariaco Basin sediments, reported by Deplazes et al. (2013b). Higher values indicate a southward migration of the Intertropical Convergence Zone (ITCZ). (c) $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios, (d) gypsum content, (e) uncertainty in the age–depth model. Gray shaded areas represent HSs.

Maxima in $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ and hence more humid conditions are more frequent during GSs 11–2. Consequently, the intensity of evaporation (dry conditions) was more pronounced between 42 and 32 cal ka BP, suggesting that the Lake Petén Itzá basin was more susceptible to changes in precipitation, evaporation and drainage during the middle and end of MIS 3. No apparent increases in $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ were observed during GSs 14–12. However, the Petén Itzá chronology has a ±2 kyr uncertainty for ages > 45 cal ka BP, and we therefore cannot exclude the possibility that the $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ increases at 48.4, 46.2 and 43.8 cal ka BP correspond to low values of L^* reflectance in the Cariaco Basin record at 49.3 (GS 14), 48.1 (GS 13) and 43.7 (GS 12) cal ka BP, respectively.

The high resolution of the geochemical data, together with the high sedimentation rates at Site PI-2 compared to PI-6, enabled us to study the impact of several GSs on the northern Neotropics (Fig. 8). GSs 11 and 10 show a gradual increase in evaporation over ∼1200 years, followed by a rapid decrease over ∼400 years. In contrast, GSs 9–6 started with a rapid increase in evaporation over ∼500 years that was followed by a gradual reduction in evaporation over ∼1200 years. This difference is likely explained by changes in the intensity of the CLLJ. During GSs 11 and 10, the CLLJ was stronger and provided more moisture to the region (Fig. 6). Consequently, millennial-scale oscillations in the position of the ITCZ and particularly its displacement to the south were possibly muted by the CLLJ. With respect to GSs 9–6, the intensity of the CLLJ and its capacity to transport moisture decreased, so any millennial-scale change in the ITCZ produced abrupt changes in lake evaporation and catchment runoff. $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ maxima associated with GSs 5–3 are less pronounced, which may have been a consequence of a weaker CLLJ during the onset of MIS 2, as high values of ΔSST between the Caribbean and eastern Pacific reveal between 34 and 24 cal ka BP. Evidence of GS 2 is missing and may be linked to the onset of the LGM, which was characterized by high humidity in Petén (see Sect. 5.1.4) and may have suppressed the GS 2 signal in the sediment record.

The response to GIs and GSs in Lake Petén Itzá has been studied in sediments from Site PI-6. For example, Correa-Metrio et al. (2012) provided a high-resolution record of pollen and carbonized material (∼200-year resolution), which revealed an increase in forest fires in this region during GSs. This is consistent with high evaporation in response to dry conditions in the northern Neotropics. The latter is more apparent for GS 13, GS 9 and GS 5, whose presence in our record is synchronous with the North Atlantic IRD deposits associated with HSs 5–3, respectively. The enhanced $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios during GS 13, GS 9 and GS 5 in Lake Petén Itzá may have been a response to dry conditions that dominated during HSs. In addition, higher abundances of gypsum were found during GS 13 and GS 9, suggesting high evaporation of lake water (Fig. 8). Similar increases in $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ and gypsum content are observed during GS 8, GS 7, GS 5, GS 4 and GS 3. We did not notice an increase in gypsum content during GS 11, GS 10 and GS 6, which may be related to the low sampling resolution in this interval. The higher $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratios, however, point to an increase in evaporation.

The new age–depth model established here shows that gypsum sand layers were deposited during cold stages HSs 4–1, in agreement with the results of Hodell et al. (2008), Mueller et al. (2010) and Escobar et al. (2012) and with the HS chronology from multiple glacio-marine records across the North Atlantic (Labeyrie et al., 1995; Jullien et al., 2006, 2007; Hodell et al., 2010; Appendix D). With regard to the age uncertainty, HS 5 is likely found between GIs 13 and 12 (49.3–46.0 cal ka BP), also corresponding to a period characterized by gypsum deposition. Lake Petén Itzá sediments associated with HSs 5–1 show Ti below average for the entire sequence, revealing high evaporation and low runoff, associated with dry conditions during their deposition.

Similarities between records from Sites PI-6 and PI-2 confirm the overall suitability of Lake Petén Itzá sediments as a hydroclimate archive. Ostracod data, for example, suggest that HSs 3–1 were characterized by high lake-water conductivity, associated with an increase in evaporation and a reduction in runoff (Cohuo et al., 2018; Pérez et al., 2021). In addition, pollen records show the establishment of savanna vegetation, associated with dry conditions and a drop in MAAT of as much as 6 ^∘C (Correa-Metrio et al., 2012; Hodell et al., 2012). Moreover, δ¹³C and δ¹⁸O values in ostracods indicate a lake level lowering such that Site PI-6 was in such shallow water that the overlying water column did not stratify thermally (Escobar et al., 2012). This is in agreement with paleoclimate records from the eastern equatorial Pacific, the Caribbean and central Mexico (Leduc et al., 2007a; Arienzo et al., 2015, 2017; Medina-Elizalde et al., 2017; Hodell et al., 2017; Caballero et al., 2019), which suggests that HSs 4–1 were generally dry and cold. The drying response of Lake Petén Itzá during HSs correlates with the globally recognized “tropical hydroclimatic events” (THEs), during which extreme regional anomalies in rainfall occurred (Bradley and Diaz, 2021). The dominance of dry conditions in the Caribbean region during THEs was favored by meltwater input to the North Atlantic that, in turn, reduced the AMOC. As a result, the mean position of the ITCZ moved approximately 1^∘ to the south, promoting droughts not only in the northern Neotropics but also in Africa and the Arabian Peninsula (Tjallingii et al., 2008; Zarriess et al., 2011).

A key finding of our high-resolution record is confirmation of the abruptness of millennial-scale hydroclimate shifts, which enabled us to link them to AMOC shifts, rather than to gradual longer-term changes caused by orbital forcing, which underlie the short-term, abrupt climate excursions. The spatiotemporal paleoclimate changes (Fig. 7) identified from multiple records across the Neotropics imply that the causes of transient vs. rapid shifts in the region are more complex than mere latitudinal shifts in the ITCZ in response to AMOC. In particular, meridional inter-basin changes in the SST gradient between the Pacific and Atlantic appear to have influenced shifts in the strength and location of deep convection linked to the CLLJ, which steered moisture into some parts of the Neotropics while, other areas were dry. Our results are conceptually consistent with findings from recent climate modeling of the SE Asian monsoon region during the late glacial, which show that regional changes in atmospheric deep convection in response to SST changes can cause strong regional hydroclimate contrasts linked to meridional shifts in low-pressure centers, rather than only latitudinal shifts in the position of the ITCZ (Hällberg et al., 2022).

6 Conclusions

We used the sediment record from Site PI-2 in Lake Petén Itzá to infer changes in runoff, evaporation intensity, organic matter provenance and redox conditions during MISs 3–2 (59–15 cal ka BP). The onset and middle of MIS 3 (59–39.2 cal ka BP) was characterized by periods of high precipitation and runoff into Lake Petén Itzá. At that time, lake levels were high, the lake was very productive, and bottom-water anoxia persisted. The MIS 3–2 transition (39–24 cal ka BP) was characterized by a gradual decline in runoff and an increase in evaporation but was accompanied by high climate variability indicative of regional climate instability. This regime shift was accompanied by input of mixed aquatic and terrestrial organic matter and more oxygenated bottom waters, related to lower lake levels near the end of MIS 3. MIS 2 was a dry period with low runoff, except for the LGM, during which wet conditions prevailed. During the LGM, lake sediments received greater inputs of aquatic organic matter, and the lake hypolimnion was persistently anoxic, presumably associated with incomplete water-column circulation.

Both GSs and GIs were detected in the sediment record at Site PI-2. Compared to previous studies, the record from Site PI-2 provides greater temporal resolution and demonstrates that GSs 9, 8, 7 and 6 began with abrupt increases in evaporation and ended with gradual increases in humidity, whereas GSs 11 and 10 showed the reverse pattern. This can be explained by the lack of moisture transported by the CLLJ during GSs 9–6, which made the local climate susceptible to changes in the latitudinal position of the ITCZ. The high resolution of the PI-2 record highlights the abruptness of climate shifts, which are tied to abrupt changes in the strength of AMOC.

Comparison of our climate proxy records with other terrestrial and marine paleoclimate data from the northern Neotropics, the Caribbean, Florida and northern Mexico indicates that variability in precipitation and, in turn, runoff during MISs 3–2 was strongly influenced by latitudinal shifts in position of the ITCZ in response to changes in the strength of AMOC. Changes in the intensity of the CLLJ and the transport of moisture by the NAM to northern Mexico also influenced fluctuations in runoff. Short-term changes in Ti and $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ indicate that latitudinal ITCZ variations during MISs 3–2 modified runoff on short (millennial) timescales in response to AMOC shifts, whereas the CLLJ, whose strength/weakness depends on the ΔSST between the Pacific and the Caribbean, among other factors such as the extension of the NASH and interaction with the ITCZ, additionally influenced the ingress of moisture into the region on longer timescales. Whereas our study suggests the importance of inter-basin SST differences in defining regional moisture changes via the strength of the CLLJ on longer timescales, in addition to shifts in the ITCZ, more records and additional modeling studies are required to disentangle the importance of inter-hemispheric changes (ITCZ) and inter-basin changes (CLLJ) on different timescales.

Appendix A

Table A1 A total of 31 strata were identified as “slumps” during the development of the age model. Their sedimentology suggests they were deposited rapidly (days to weeks), and consequently their thickness was omitted for age–depth modeling. The omitted deposits include carbonaceous turbidites and volcanic tephras.

Figure A1 Photographs of two cores from the PI-2 record, showing deposits from short-duration events: (a) carbonate turbidite (PET06-2C-6H-1 core) and (b) tephra C5 (Kutterolf et al., 2016) (PET06-2C-9H-1 core).

Appendix B

A principal component analysis (PCA) was carried out to discern the geochemical variables that are related to each other and that potentially respond to the same hydrological variables.

X-ray fluorescence (XRF) data were used. Previously they were weight-averaged with a LOESS smoothing (factor = 3) to remove noise from the signal and were subsequently converted to centered log ratio (CLR). Titanium, aluminum, calcium, iron and manganese were included in the statistical treatment as they are potential indicators of runoff, evaporation and redox processes.

The PCA (Fig. B1) explained more than 91 % of the original variation in the datasets using only the first two principal components. The first principal component (PC1) was mostly associated with detrital elements (Ti, Fe) and those from authigenic mineral formation (Ca) that reflect changes in runoff and evaporation, archived in the sediments. For this reason, we suggest the use of the $\mathrm{Ca}/(\mathrm{Ti}+\mathrm{Fe})$ ratio as an indicator of evaporation and Ti as a proxy for runoff. PC2 is related to detrital elements (Ti, Fe) and elements related to changing redox conditions at the lake bottom (Mn). As Mn deposition is preferentially favored in oxic environments, we propose the use of $\mathrm{Mn}/\mathrm{Fe}$ as an indicator of oxic/anoxic conditions at the lake bottom.

Figure B1 Ranking obtained from the PCA for the geochemical variables measured in the PI-2 record. The ordination analysis is related to changes in runoff and evaporation (PC1) and variations in the oxygen content in bottom waters (PC2).

Appendix C

Table C1 Main characteristics of Lithological Units 6 to 2 discussed in this study and previously defined by Mueller et al. (2010). We present the depth and age intervals spanned by each unit; the latter values were calculated in this study.

Appendix D

Table D1 Chronology of Heinrich stadials (HSs) in the PI-2 record compared to glacio-marine records. The ^a refers to the beginning of each HS, and the ^b refers to the time of maximal IRD deposition. MCD: master composite depth.