Evidence of large water-level variations found in deltaic sediments of a tropical deep lake in the karst mountains of the Lacandon forest, Mexico

Lake Tzibaná is one of the largest (1.27 km2) and deepest (Zmax = 52 m) karstic lakes in the UNESCO’s Biosphere Reserve “Nahá-Metzabok” and in the Lacandon Forest, southeastern Mexico. It archives sediments from multiple sources and the inflowing Nahá River forms deltaic deposits. In 2019, the water level in Lake Tzibaná declined by ~ 15 m, persisting for 4 months and exposing the Nahá River Delta. A geophysical profile on the exposed delta revealed an accumulation of ~ 20 m of such deposits. Three sediment outcrops from an inactive channel in the Nahá River Delta, which ranged in height from 0.6 to 1.43 m, were sampled and a multi-proxy analysis of biological remains and geochemical variables was conducted. Four facies were observed: (1) massive-coarse sand, (2) fine sand, (3) dark leaf litter and (4) massive silty clay, each characterized by specific microcrustacean, testate amoebae and diatom taxa. Six dark leaf litter horizons were radiocarbon dated and revealed a complex depositional history including inverted ages making the establishment of an age model difficult. Nevertheless, past lake-level changes and the formation of the four facies match three characteristic water-level stages, which can also be observed on recent satellite images: (1) Massive-coarse sand deposits, with compositional and sedimentological characteristics of a shoreline environment and fluvial lateral banks, were formed during large-magnitude reductions in the lake level, similar to the one in 2019, (2) Interbedded layers of fine sand and dark leaf litter, currently found in low-energy fluvial environments, were formed during shorter and less pronounced decrease intervals, and (3) Massive silty clay, with distinctive microorganisms from low-energy lacustrine environments, is deposited during high water-level stages, when the delta is covered by water. Our findings illustrate how hydrological changes alter sedimentary dynamics in deltaic areas of lakes. Despite the complexity of their depositional processes, deltaic records can serve as a complementary source of paleolimnological information to records from distal zones due to their sensitivity to variations in water level, especially during extreme and prolonged desiccation events. Future research should attempt to combine evidences from deltaic and sediment sequences from deeper zones of Lake Tzibaná to reconstruct water-level variations during the entire Holocene. Understanding past lake-level reductions is not only relevant for the local indigenous communities but also crucial for the conservation of this ecosystem of international importance.


Introduction
Past water-level variations are typically reconstructed from sedimentary records collected from profundal or littoral sites of lakes that are not directly influenced by the sedimentation from inflowing streams (Martínez-Abarca et al. 2021a;Ortega-Guerrero et al. 2021). However, sedimentation at these sites may eventually be interrupted during extreme desiccation events when the absence of a water column and the resulting subaerial exposure lead to a loss of sediment and consequently a discontinuity in the paleoenvironmental record (Lozano and Ortega 1998;Caballero et al. 1999;Israde-Alcántara et al. 2002;Wood and Scholz 2017). In such cases, delta deposits might serve as complementary paleolimnological record. Several studies have shown that shifts in the deltaic sediment composition provide information about the paleo shoreline and large water-level variations that occur during extreme and prolonged desiccation events (Hyne et al. 1979;Machlus et al. 2000;Sáez et al. 2018;Betancurth and Cañón 2021;Zhang et al. 2021). However, deltaic sediments are rarely used for paleolimnological inferences because processes such as erosion, re-suspension, and posterior re-deposition of old material on top of more recent sediments make their interpretation difficult (Xin et al. 2019).
The Lacandon Forest is one of the most important national protected areas in the tropical mountains of Mexico, because of its biodiversity, ecological processes and natural resources (Medellín 1994;SEMARNAT 2018). The Area of Protection for Flora and Fauna (APFF) of Nahá and Metzabok, a 7200ha large zone in the northern Lacandon Forest, is categorized as a Ramsar site of international importance and, in 2010, was included in the UNESCO World Network of Biosphere Reserves (RAMSAR 2004;Echeverría Galindo et al. 2019). With an area of only 0.4% of the country's surface, the Lacandon Forest contains 48% of the bird's species, 33% of bats, 11% of the reptiles and 25% of the mammals of Mexico, making it the most important region for biodiversity in North America (Vásquez and Ramos, 1992;Medellín 2000). Moreover, it supports vulnerable, endangered or critically endangered species and is a fundamental provider of natural resources for the native communities of the region (Medellín 1994;Hernandez-Ruedas et al. 2014). The region is an important water resource for the surrounding human populations, since it hosts a large number of relatively large and deep karst lakes. However, only in recent years few lakes have been studied in more detail (Charqueño-Celis et al. 2020;Bücker et al. 2021;Hoppenbrock et al. 2021). For instance, lakes Nahá, Amarillo, Lacandón, Metzabok and Tzibaná have been the focus of ecological and paleolimnological studies, because their sediment records have proven to represent rich archives of environmental and climatic change in the region (Domínguez-Vázquez and Islebe 2008; Díaz et al. 2017). Several cores were retrieved from the deep zones of lakes Ocotalito, Lacandon and Tzibaná (Vázquez-Molina et al. 2016;Cisneros-García 2017;Díaz et al. 2017). Bioindicators (pollen, cladocerans, ostracods, among others) and geochemical data revealed past water-level fluctuations in these lakes during the last 7 kyr. An important decrease in water levels occurred during the mid-Holocene (5-4 kyr BP), associated with a decrease in runoff (Díaz et al. 2017). This period has been linked with the southern migration of the Intertropical Convergence Zone, a region of low atmospheric pressure that promotes humidity in the tropics. Despite of the international importance of the aquatic ecosystems in the tropics and their known sensitivity to climatic and hydrological variations during the Holocene, no attempts have been undertaken so far to use deltaic sediments to better understand extreme desiccation events.
Between February and July 2019, Lacandon Lakes Tzibaná and Metzabok experienced dramatic declines in water level (Jiménez 2019;Romero 2019). Deltaic deposits of the Nahá River are usually submerged in Lake Tzibaná but were exposed between August and November 2019 as a consequence of the rapid drop in water level. Until early 2020, the lake level rose and the deltaic sediments were again covered by water. The desiccation event provided a unique opportunity to sample outcrops on the river delta during a field campaign in October 2019. The analysis of these deltaic sediments, a previously unexplored environmental archive in the Lacandon Forest, was expected to provide evidence of earlier large water-level variations similar to the 2019 event. This study aims at (1) describing and analyzing for the first time the uppermost deltaic sediments from the Nahá River Delta, (2) understanding sedimentological processes taking place across the river delta and their relation to water-level variations, as well as (3) estimating the thickness and the internal morphology of the deposits of the Nahá River Delta using electrical resistivity tomography. Fig. 1 Study area. A Location and topographic map of the state of Chiapas, southern Mexico; the dark blue polygon indicates the limits of the Lacandon Forest and the white square the location of Lake Tzibaná. B Regional topographic map of the study area including the location of the water bodies with superficial connections to Lake Tzibaná (topographic data from INEGI, 2013 andSGM, 2008). C Satellite image of the Metzabok-Tzibaná lake system. The three main lakes of the system are shown, as is the Nahá River. D Images of three recent stages in Lake Tzibaná: high-water stand (March, 2018), moderate lake water-level reduction (February, 2019), and lowwater stand (August, 2019). The yellow points indicate the location of the study zone on the Nahá River Delta. Satellite images obtained from Google Earth Pro (2020) CNES/Airbus images, Image date: 29/12/2010 and 14/02/2019, and Planet Explorer (2017; https:// api. planet. com) a stream that originates from Lake O'ton K'ak to the east, (c) a stream from Lake Metzabok to the west and (d) a stream to the north that connects Lakes Tzibaná and Metzabok, when water levels are high. Sediments transported by the Nahá River formed a delta of approximately 0.11 km 2 in Lake Tzibaná.

Geology and climate
Lake Tzibaná is located in the municipality of Ocosingo in the northern part of the Sierra de Chiapas mountain range. The rocks of the region comprise dolomite and limestone deposits of uppermost Cretaceous (Maastrichtian) age in the north, and Paleogene in the south (Quezada Muñeton 1987). Lake Tzibaná lies in limestone bedrock  prone to develop karst features. The most important tectonic structures are NW-SE oriented faults and anticlines (SGM 2008). The region has a tropical climate (Aw) according to the Köppen climate classification. Mean annual precipitation amounts to 1600 mm yr −1 and exhibits strong seasonality, with the dry season extending from November to April and the rainy season from May to October (SMN 2020). Mean atmospheric temperature oscillates around 28 °C, reaching the minimum (15 °C) during dry season and the maximum (36 °C) in the rainy season (SMN 2020).

The sudden water-level decline in 2019
In August 2019, several media reported a sudden, significant reduction in the water level of Lake Tzibaná by approximately 15 m (Jiménez 2019;Romero 2019). This reduced the lake surface area to 0.71 km 2 , i.e., to only 66% of its surface area during high stands (Fig. 1D). Whereas seasonal variations in the level of Lake Tzibaná of approximately 10 m are known to occur (Altmann 2020), major desiccation events, such as the one observed in August 2019, are rare. Riparian Lacandon communities interviewed during the field campaign in October 2019 remembered a similar event to have occurred ~ 70 years ago approximately during the month of April. Interestingly, lake levels seem to have dropped much faster in 2019 (2-3 days) compared to several months for the event ~ 70 years ago. The specific cause of the sudden desiccation event in 2019 remains open to discussion.
Meteorological causes cannot be ruled out completely, but are relatively unlikely: Although the mean annual precipitation reached a minimum of ~ 1500 mm yr -1 in 2019, similar reductions in precipitation observed during the past decades have not led to similar, drastic water-level reductions (Altmann 2020). Furthermore, other lakes in the study area, such as Lakes Nahá and Ocotalito did not present evidence for comparable water-level changes in 2019.
Geophysical surveys at Lakes Tzibaná and Metzabok carried out in 2018 and 2019 report a strong interconnectivity between surface water and the underlying karst aquifer . During the field campaign in October 2019, rapid water-level variations in small remnant water bodies of Lake Metzabok (various meters within a few hours) were observed, which were clearly fed and drained through a series of crevices in the underlying limestone bedrock. Based on these observations, the drainage of Lakes Tzibaná and Metzabok through underground karst conduits becomes the most likely explanation.
High magnitude desiccation may also have occurred during ancient Maya occupations of Lake Tzibaná. A petroglyph found on a limestone cliff located at the shore of the lake has been dated between CE 1250 and 1400 (Palka et al. 2020). The petroglyph, which is situated below the water line during normal lake-level stands, shows the god Quetzalcoatl, who in Mayan mythology is the god of rain. This petroglyph and ceramics found in its surroundings during the 2019 water-level drop indicate an important water-level reduction during its creation (Lozada and Palka 2020;Zorich 2020).
Irrespective of its particular causes, the water-level drop in 2019 exposed a large part of the normally water-covered delta deposits in the southern area of Lake Tzibaná. Field observations and satellite images (Fig. 2) show large sand and gravel banks covering the exposed delta. The temporal variation of an inactive channel (hereafter referred to as "Cin-A") was observed (ESM1). This channel is highly sensitive to lake-level variations: (1) it is activated during moderate water-level reductions (no more than ~ 10 m), (2) it becomes inactive during drastic water-level reductions (similar to August 2019), and (3) it is completely submerged during high water-level stands.

Fieldwork and initial outcrop descriptions
Fieldwork was carried out in October 2019 in the southern part of Lake Tzibaná. We analyzed two sedimentary profiles in the inactive channel Cin-A. This site was selected due to its susceptibility to variations in water level and accessibility of outcrops. Moreover, Cin-A and the surrounding sand banks have been stable topographic/bathymetric features since at least 36 years as can be seen from satellite imagery (Fig. 2, ESM1). Profile 1 (0.5 m above water level in October 2019), with a total thickness of 0.92 m was located 2 m south of the shoreline in October 2019. Profile 2 (0.5 m above water level in October 2019) was 0.6 m thick and located 2 m south of Profile 1. Because of the proximity of profiles 1 and 2, we were able to use visible marker horizons along the outcrops to correlate and combine both profiles into one stratigraphic sequence, representative of the superficial sediments of the inactive channel Cin-A (Fig. 3). Profile 3, 3 m above the water level in October 2019 and 1.43 m thick, was located inside the sand bank "Ba1" (Fig. 2) 70 m southwest of Profiles 1 and 2, and 10 m south of the shoreline in October 2019. During fieldwork, we made preliminary stratigraphic descriptions of sediment color (following the Munsell soil color chart), thickness, texture, structures, composition and boundary type. We sampled approximately 0.5 kg of sediment in strata where sedimentological differences were visually detected. Stratigraphic columns for each outcrop were drawn using the free software SedLog 3.1 (Zervas et al. 2009).

Texture, composition and geochemistry
Grain-size distribution and sorting of sands were determined by sieving samples (70-160 g) collected from five sandy strata of the outcrops (Fig. 3) in order to obtain hints on the respective depositional processes. These samples were previously dried at 40 °C, weighed and subsequently sieved at 1 phi increments from 0 to 4 ϕ on the Krumbein scale. Each size fraction was weighed and its percentage of the initial dry sample weight was calculated (wt%). We calculated the grain-size standard deviation (SD) using the ϕ percentage data. The SD enables quantification of sorting, using the scheme of Folk (1980). Smear slides were prepared from six silt-clay samples (Fig. 3), following the procedure of Marsaglia et al. (2013). We made semi-quantitative estimates of the major sediment components in each sample by counting 200 points, following the method of Gazzi-Dickinson (Ingersoll et al. 1984) and the scheme proposed by Schnurrenberger et al. (2003).
Total Carbon (TC) and Total Nitrogen (TN) were measured in 10 samples that had been dried at 40 °C and ground to powder with a mortar and pestle. Samples were analyzed with a Carlo-Erba NA1500 CNS Analyzer. Total Inorganic Carbon (TIC) of powders was analyzed using a UIC 5017 CO 2 Coulometer, coupled with an AutoMate automated carbonate preparation device. Total Organic Carbon (TOC) was computed as TOC = TC-TIC. The TOC/N mass ratio was calculated and used to estimate the relative contributions of terrestrial and aquatic organic matter input to the sediment using the approach of Meyers and Teranes (2001).
Biological remains of ostracodes, gastropods, bivalves and testate amoebae in 1 g of sediment were identified, and dominant taxa reported. Samples were sieved with a 63-µm mesh sieve using distillated water. Biological remains were observed Fig. 3 Stratigraphy of the profiles 1, 2, and 3, representing the four facies, types of contacts, and grain size of each unit. Letters indicate sections recognized during fieldwork. Cin-A sequence was constructed based on the correlation between Profiles 1 and 2, and in particular, the indicated marker horizons. Triangles mark sample locations under a ×40 magnification stereoscopic microscope. Bivalves and gastropods were identified to genus level using Rubio-Sandoval (2019), ostracode identification followed Echeverría Galindo et al. (2019), and testate amoebae were identified using Echeverría Galindo (2017) and Charqueño-Celis et al. (2020). Poorly preserved adult and juvenile bivalves and gastropods could not be identified, and therefore only their presence/absence is reported. Diatom sample preparation followed Krahn et al. (2021). Mixed and dried sediment samples (~ 0.1 g) were cleaned by adding hydrochloric acid (37%) and hydrogen peroxide (30%), and then heated up to 70 °C for 20 min. Subsequently, samples were centrifuged and washed until the suspension became pH neutral. For determination of diatom concentrations, known quantities of Divinylbenzene microspheres were added to the diatom suspension. Concentrated drops therefrom were placed on cover glasses, dried and embedded in Naphrax®. If possible, at least 400 diatom valves were counted in each sample using an Axio Imager.M2 (Zeiss) research microscope equipped with a plan-apochromatic oilimmersion objective (× 100/NA1.4). Diatom species were mainly identified based on literature such as Krammer and Lange-Bertalot (2010) and Spaulding et al. (2020). Subfossil Cladocera analysis were carried out according to the standard methodology (Frey 1986). Fresh sediment samples (approximately 2 cm 3 ) were treated with 10% HCl to remove carbonate and then heated up in 10% KOH solution. The residue was mixed by a magnetic stirrer (approximately 20 min) to eliminate humic matter. Subsequently, samples were sieved with a 38-μm mesh sieve, transferred into a test tube and filled up to 5 cm 3 with distilled water. Samples were stained with safranin. Depending on the number of remains, two to eleven slides (0.1 cm 3 liquid suspension) were examined per sample. Obtained results were normalized to 1 cm 3 . The taxonomic identification followed Wojewódka et al. (2016;2020a, b) (ESM2).

Chronology
Six leaf samples from organic-rich strata mainly composed of leaf litter were Accelerator Mass Spectrometry (AMS) radiocarbon dated by Beta Analytic Inc., Miami (Table 1). Samples were selected based on their stratigraphic positions representing leaf litter layers clearly distributed along the profiles. Results were calibrated using IntCal20 (Reimer et al. 2020) and converted to calendar ages (AD) with the free software Calib Rev 8.1.0 (Stuiver and Reimer 1993). Stable carbon isotope (δ 13 C) measurements were also obtained for all six samples. Table 1 Radiocarbon dates for six leaf samples from the dark leaf litter (DLL) facies. Samples within each sequence are listed in stratigraphic order. We report ages in the post-bomb calibration and calendar ages (AD). All mean AD data are in the common era (CE) Internal structure of the delta deposits

Geophysical measurements
Electrical resistivity tomography data were collected on the exposed delta along a NW-SE trending transect (Fig. 2), crossing the sedimentary deposits of the inactive channel Cin-A and the bank Ba1. The geophysical transect provided information about the internal morphology of the delta deposits. We used a Syscal Pro Switch 48 (IRIS Instruments, France) with 48 stainless steel electrodes installed at 5-m spacing in the central part of the transect (electrodes 13-36), and at 10-m spacing at the ends of the transect (electrodes 1-12 and 37-48), resulting in a total profile length of 345 m. Following the approach by Flores Orozco et al. (2020), we deployed a dipole-dipole configuration, which combines measurements with dipole lengths of 1, 2, and 4 times the electrode separation and 10 potential readings for each current injection. Data processing included (1) the removal of erroneous measurements defined as those readings with current injections below 1 mA or negative apparent resistivity values, and (2) the visual inspection of pseudosections to remove further obvious outliers. Filtered data were then used to compute resistivity images, using a smoothness-constraint, least-squares inversion algorithm (Kemna 2000).

Composition of deltaic sediments
Both the inactive channel Cin-A and the bank Ba1 contain four main sedimentary facies based on their composition and texture: (1) massive-coarse sand, (2) fine sand, (3) dark leaf litter, and (4) massive silty clay (Table 2, Fig. 4). The sequence Cin-A has a total thickness of 1.18 m, whereas the sequence Ba1 has a total thickness of 1.43 m. Strata of Ba1 are inclined 30° from the horizontal towards the lake and show no evidence of deformation or diagenesis (e.g. decomposition of organic matter or oxidation).

Massive-coarse sand facies
This facies corresponds to massive medium-tocoarse sand and gravel of light brown color (2.5Y 4/3 Munsell). Massive-coarse sand forms the basal layers in the Cin-A and Ba1 sequences. The observed deposits are 30-40 cm thick and have superior erosive contacts. The sand and gravel fractions are mainly composed of rounded limestone grains (80%) with small fragments of bivalves and gastropods between 1 and 10 cm (10%), as well as cladocerans, ostracodes, testate amoebae, as well as leaf and wood fragments < 5 cm (10%). Gastropods of the genus Pachychilus were identified and may belong to the species P. indiorum and/or P. glaphyrus. Testate amoebae (Arcella, Centropyxis) and ostracodes (Cypridopsis) were scarce, while cladocerans were more abundant (Alonella, Chydorus, Karualona). Grain-size analysis (Fig. 5) indicates a moderately sorted sediment (SD = 0.75). TOC/N values in massive-coarse sand samples of Cin-A range between 20 and 22, but display values of up to 60 in Ba1. TIC values in massive-coarse sand samples from Cin-A had values ~ 6%, whereas TIC values in Ba1 were slightly higher around 8.7%. Rounded silt intraclasts were found in the massive-coarse sand deposit of sequence Cin-A, with no spatial organization, such as imbrication or alignment. The mean diameter of these intraclasts varied between 3 and 21 cm (Fig. 4). Smear-slide analyses suggest that the intraclasts are composed of 55% detrital and 45% biogenic material. The detrital material is composed mainly of rhombic carbonate (30%), corroded carbonate (20%), and feldspar (5%). The biogenic fraction is composed of diatoms, ostracodes, cladocerans, and testate amoebae (30%), phytoliths (5%), amorphous organic matter (5%), and charcoal (5%). Diatoms encountered mainly belong to the genera Cyclotella, Gomphonema, and Eunotia, and possess frustules showing evidence of dissolution. Ostracodes belong to the genera Cypridopsis, Chlamydotheca, and Darwinula, cladocerans to the genera Alona, Chydorus, Disparalona, Karualona, and Liederobosmina, and testate amoebae belong to the genera Arcella, Difflugia, and Centropyxis.
The TOC/N ratio ranges between 27 and 30 in the inactive channel Cin-A and between 32 and 36 in the bank Ba1. TIC in dark leaf litter samples (Fig. 5B) varies between 4.5 and 5.8% in the inactive channel Cin-A and between 6 and 7% in the bank Ba1, values that are lower than those measured in all other facies. δ 13 C values in the six dark leaf litter facies samples range from −31.0 to −28.0‰.

Massive silty clay facies
This facies is composed of massive strata of brownish silty sediments (2.5Y 5/4), with a thickness between 5 and 40 cm. In some cases, horizons show subtle laminations and intercalation with fine sand facies. Massive silty clay layers exhibit sharp erosive boundaries with the underlying strata. Smear-slide analyses indicate that massive silty clay samples are composed mainly of 90% detrital and 10% biogenic material. The detrital fraction contains rhombic carbonate (60%), corroded carbonate (28%), and feldspar (2%). The biogenic fraction is composed of charcoal (5%), ostracodes, diatoms, testate amoebae (2%), amorphous organic matter (3%), and a small fraction of gastropod, ostracode, and testate amoebae remains. Gastropods have mean sizes that range from 3 to 5 mm. The dominant genus of testate amoeba is Centropyxis. Cypridopsis and Darwinula are the dominant ostracode genera. Diatom frustules, mostly Cyclotella species, display a high degree of dissolution. The TOC/N ratio in samples from the inactive channel Cin-A shows values around 29, whereas samples from the bank Ba1 display values between 28 and 37. TIC varies between 6 and 6.5% in Cin-A and between 7.5 and 8% in Ba1.

Chronology of superficial deposits
AMS radiocarbon measurements on samples from the inactive channel Cin-A and the bank Ba1 provided dates that range from CE 1058 to 1940 (Table 1) with a gap between CE 1250 and 1400. It is important to highlight that radiocarbon dates are not in stratigraphic order. In sequence Cin-A, only two superficial dates are stratigraphically congruent (CE 1433 and 1169), whereas the bottom ages are younger (CE 1870 and 1750). The two dates of the sequence Ba1 are not in stratigraphic order either (CE 1490 and1644). Given this lack of chronological control, the low number of age determinations, as well as the fact that most of the facies probably represent rapid or near-instantaneous deposition, we did not attempt to develop age-depth models for the sediment sequences (Fig. 6).

Internal morphology of the Nahá River Delta deposits
The electrical resistivity image along transect TZI19-1 (Fig. 7) sheds light on the internal morphology of the delta deposits. High resistivity values (> 30 Ωm) close to the surface are associated with the mostly sandy and gravely deposits, which were sampled and analyzed by direct methods. These coarse resistive materials can be clearly distinguished from a deeper conductive horizon, which we interpret as clayey and silty lacustrine sediments. Profile 3 spans only the uppermost part of Ba1. According to the resistivity section, Ba1 could reach a thickness of up to 20 m (Fig. 7). The inactive channel Cin-A of the Nahá River is also underlain by resistive (> 30 Ωm) sand deposits. However, in contrast to the much thicker bank deposits, the thickness of the sandygravely cover reduces to only ~ 5 m below Cin-A. Elevated electrical resistivity values (> 60 Ωm) observed within a second bank ("Ba2") at the center of the transect suggest higher content of coarse material (sand) than in bank Ba1. Southeast of Ba2, the underlying conductive unit reaches the surface. In October 2019, this area was part of a shallow remnant water body (Fig. 2). Limestone bedrocks, which clearly stand out

Sedimentary facies and past lake water levels
The abrupt water-level decline in Lake Tzibaná, which occurred in August 2019, exposed the usually submerged deltaic deposits of the inflowing Nahá River. Whereas small seasonal variations in the water level of Lake Tzibaná can be observed on satellite Fig. 6 Composite stratigraphy of outcrops. The sequence Cin-A is composed of Profiles 1 and 2 from the inactive channel zone Cin-A (Fig. 2). Ba1 is composed of Profile 3, located on bank Ba1, south of Cin-A images (Altmann 2020; Fig. 2, ESM1), major desiccation events, such as the one in 2019, have never been reported before in scientific literature. This unique opportunity to obtain information about the internal structure and sedimentary facies of the Nahá River delta provides a first record of possible variations in lake water level. Satellite images reveal the inactive channel Cin-A as susceptible area to lake water-level variations during 2017 to mid-2019 (ESM1). In addition, the geophysical transect shows high values of electrical resistivity in the sediments subjacent to Cin-A (> 30 Ωm), indicating a channel bed composed of sand in the first 5 m depth. Both morphology and sedimentological characteristics of the Cin-A outcrops make them points to infer possible water-level changes in the past.
The sedimentological characteristics of the massive-coarse sand, particularly the roundness of grains, the moderate sorting, the comparably low abundance or even absence of diatoms, testate amoebae, cladocerans, and ostracodes, as well as the presence of Pachychilus gastropods, suggest a possible deposition under low water-level conditions (Lam-Gordillo et al. 2012;Sigala et al. 2016;Charqueño-Celis et al. 2020). This interpretation is further supported by the higher abundance of cladocerans belonging to the genus Karualona that encompasses species living near the shoreline in shallow water (Wojewódka et al. 2016(Wojewódka et al. , 2020b, over sandy or rocky substrata, or within macrophytes (Fuentes-Reines and Elmoor-Loureiro 2015).
The first possible process that may explain the massive-coarse sand deposition involves the formation of sand banks resulting from lateral sedimentation within the river bed of the Nahá River. Horizons that sit atop the massive-coarse sand facies in the bank Ba1 lie at an angle of 30° and show no evidence of deformation, which suggests a syn-sedimentary deposition (Fielding 1984). We interpret these deposits as having been associated with an older, now abandoned channel of the Nahá River. This interpretation is corroborated by the high TOC/N ratio (60) of the unique massive-coarse sand deposit in Ba1, which indicates high input of terrestrial organic matter whose value greatly exceeds the TOC/N ratio (20-22) of the only massive-coarse sand deposit in the sequence Cin-A (Fig. 8B). These observations indicate the formation of massive-coarse sand facies of Ba1 from lateral bank sedimentation within the Nahá River. In several studies conducted in the Valdamo Basin (Italy), the Yangtze River Delta (China), and the Hainan Delta (China), coarse sand and gravel deposits with textural characteristics similar to those of the massive-coarse sand facies were described and linked to lateral deposition in a fluvial stream, which further supports our interpretation (Billi et al. 1991;Cheng et al. 2021;Zhang et al. 2021).
The second possible explanation interprets the bank as a beach deposit, where wave action has from the inversion of electrical resistivity tomography data from the exposed delta. Black circles along the surface indicate electrode positions and the dashed lines show the inferred limits of the main geological units reworked and sorted sand and gravel along the lake shoreline. This process is more likely to have formed the massive-coarse sand deposit in the inactive channel Cin-A. The development of similar massive sand banks in lacustrine environments has been linked to wave-reworked sediments in similar studies (Miall 2014). Such intraclasts are often linked to clayey or silty sediments having been reworked by wave action near the shoreline (Tye and Coleman 1989;Pratt 2002). Furthermore, Li et al. (2016) proposed that mud intraclasts embedded in sandy banks along deltaic lake deposits are associated with a sudden ingress of sediment along the delta. In this case, intraclasts may have formed at the delta front and deposited in shallow water (e.g. medial delta area with low slope).
Recent measurements on water samples of the Nahá River yielded a relatively low conductivity (180-300 µS cm −1 ) and a neutral to slightly acidic pH (7-6), in contrast to higher conductivities (316-328 µS cm −1 ) and a slightly alkaline pH of around 8 in lake water (Martinez-Abarca, unpublished data, January 2022). The presence of the diatom species Eunotia in the intraclasts, a diatom indicating nutrient-poor, acidic waters and low conductivities (Mora et al. 2017), suggests an important inlet of fluvial waters from the Nahá River. The cladoceran Liederobosmina sp. was especially abundant in the intraclast sample as well. Specimens of this taxon are generally found in the pelagic zone of deeper lakes, however, Wojewódka et al. (2016Wojewódka et al. ( , 2020a reports this genus also in much shallower waterbodies (max. depth of Fig. 8 Photographs of the study area and deposited sediments. A View to the northwest of the lake. Filled circles indicate the locations of the outcrops (red: sequence Cin-A; blue: sequence Ba1). Dashed lines indicate deposits and superficial contacts between facies. The inactive channel Cin-A of the Nahá River crosses the picture from left to right. B View to the south of the lake. The sand bank Ba2 delimits Cin-A to the east. The blue arrow indicates the former direction of Cin-A. The remnant water body (Re) is situated behind the sand bank. Examples for massivecoarse sand (MCS) and fine sand (FS) facies. C Photograph of the outcrop corresponding to Profile 1. D Photograph of the outcrop corresponding to Profile 3. E Photograph of the current mouth of the Nahá River; the dashed rectangle shows the location of image F. F Laminated deposits of the dark leaf litter (DLL) and fine sand (FS) facies along the shore of the Nahá River approximately 2 m). Similarly, the periphytic diatom Gomphonema sp. found in the intraclasts is mostly found in littoral zones of lakes (Krammer and Lange-Bertalot 2010). Additionally, the cladoceran assemblage in the massive-coarse sand facies was characterized by the highest abundance of littoral species in comparison to the other facies. These evidences support the idea of decreasing lake water-levels and shallow water conditions (littoral zone). The good preservation of biological remains suggests that intraclasts were covered by the massive-coarse sand facies in a short time.
Although we cannot say conclusively which of the two processes (i.e., lateral deposition in a stream or wave reworking along the shoreline) formed the massive-coarse facies, both sedimentary processes indicate times of low lake water levels, similar to the desiccation event during August-October 2019 when we also observed the ongoing formation of sand banks at other parts of the delta (Fig. 9A). The downward continuation of the massive-coarse sand deposits below the outcrop was not visible in the field. However, the electrical resistivity results suggest a possible total thickness of these delta deposits of up to 15 m below the Ba-1 outcrop and up to 5 m below the Cin-A outcrop.
The fine sand and dark leaf litter facies consist of interbedded strata, the thickness of which ranges from 0.15 to 0.58 m in the inactive channel Cin-A and from 0.53 to 1.05 m in the bank Ba1 (Fig. 6). The high percentage of well sorted (SD = 0.17) sediment in the fine sand facies, along with the subrounded angular limestone fragments, indicates low-energy transport and a relatively long transport distance. Calcite concretions represent 34% of the fine sand facies. Such concretions typically occur in organic-rich clastic Fig. 9 Possible sedimentation processes associated with facies found in the inactive channel Cin-A (location indicated by red circles) and the bank Ba1 (blue circles). The cross sections correspond to the NW-SE transect A'-A'' shown in the satellite images on the left. Stratigraphic columns of both sequences are interpreted in terms of three different lake stages: High level (HL), reduction and filling of Cin-A, and low level. The geometry of the limestone bedrock and the layer of fine sediment (possibly lacustrine) is based on the internal morphology of these units revealed by the geophysical section (Fig. 7). A During low-level periods the massive-coarse sand facies is formed by lateral sedimentation in deltaic channels or by wave action along a shore line. In the sequence Cin-A, intraclasts caused by waves (indicated by arrows) suggest a formation process along the shoreline. This scenario is analogous to that occurring during desiccation of August 2019. B Moderate reduction of water level, during which the inactive channel Cin-A was activated (indicated with a black arrow). During this period, the dark leaf litter and fine sand facies were deposited. During the same period, a paleo-channel may have been activated on sand bank Ba1. The scenario is analogous to the beginning of the water-level reduction in February 2019. C High lake-level stage leading to the deposition of the massive silty clay facies. This scenario is analogous to the high water level observed in March 2018 sediments of non-active deltas and form during early diagenesis by compaction of the sediments (Curtis and Coleman 1986). Therefore, the calcite concretions in the fine sand facies suggest the onset of sediment compaction.
Plant remains of the dark leaf litter facies show no evidence of burning, and thus were not produced by forest fires (Fig. 4). Good preservation of these remains, with relatively little decomposition, indicates that dark leaf litter deposits were covered rapidly by the overlying layers. Diatoms and cladocerans were almost absent in dark leaf litter layers suggesting flowing waters (Sa-ardrit and Beamish 2005) as well as rapid covering by overlying layers. The δ 13 C values from the six dark leaf litter samples (−31 to -28‰) indicate a C3 vegetation source (Meyers and Teranes 2001). Nowadays, the vegetation of the region is composed of montane rainforest taxa with evergreen species such as Brosimum alicastrum and Hirtella americana (Domínguez-Vázquez and Islebe 2008). As these plants use the C3 photosynthetic pathway (Kohn 2010), the plant remains of the dark leaf litter facies correspond fairly well with the current vegetation.
Plane-parallel layers in sediments similar to fine sand and dark leaf litter facies (Fig. 9D) were described by Fielding (1984) and Gomez et al. (2001), and were associated with low-energy fluvial environments. Such features are usually found in the mouth of the principal channel of a river, in a lake or floodplain (Miall 2014;Cheng et al. 2021). Moreover, in deep-lake drilling projects (Martinez-Abarca et al. 2021b), such laminations were linked to fluvial regimes. The intercalation of fine sand and dark leaf litter layers consequently indicates a water-level reduction, during which the inactive channels Cin-A and the one inferred from Ba1 were flooded and located close to the mouth of the Nahá River. During low lake stages, large quantities of leaves and wood are transported away from the part of the Nahá River directly upstream from the delta, where they accumulate during high lake-level stages, when the water moves very slowly. High abundance and similarity of ostracodes, cladocerans, diatoms, and testate amoebae species between fine sand and dark leaf litter layers support the inference of a low lake stage. Moreover, the co-existence of both ostracodes genera Cypridopsis and Cytheridella in dark leaf litter layers suggests warm waters with an electrical conductivity below 1000 µS/cm (Cohuo et al. 2020;Pérez et al. 2021). A similar situation was observed in February 2019 (Fig. 9B), when the water level lowered about 10 m and the surface reduced to 0.1 km 2 (Altmann 2020).
Although the dark leaf litter facies in the delta channels might indicate individual depositional events, the complexity of sedimentation dynamics of a river delta, which include successive deposition, erosion and re-deposition processes, might result in the re-deposition of various older strata of vegetation remains during one single depositional event. Radiocarbon dates from the dark leaf litter horizons in Cin-A and Ba1 range from CE 1158 to 1940 and due to the complexity of the stratigraphic order of these dates, the establishment of a reliable age model is not feasible. Inverted ages along the stratigraphic sequence may be explained by the repeated deposition, erosion and re-deposition of the dated vegetation remains. Nevertheless, the deltaic sediment record indicates multiple episodes of plant-remain input into the river-lacustrine system over recent centuries. Moreover, the interruption of intercalated leaf litter and fine-sand facies by massive silty clay deposits still evidences a sequence of low and high lake stands recorded in the delta deposits.
Layers of massive silty clay in Cin-A and Ba1 exhibit sharp erosive contacts with the underlying layers, suggesting subaerial exposure of the underlying layer prior to massive silty clay deposition (Billi et al. 1991). The diatom Cyclotella, the ostracodes Cypridopsis and Darwinula, as well as the testate amoebae Centropyxis in the massive silty clay sediments confirm the sites were submerged during sediment deposition (Sigala et al. 2016;Charqueño-Celis et al. 2020;Pérez et al. 2021). The cladocerans belonging to the genus Liederobosmina include species living in the pelagic zone, in the open water column (Flössner 2000). In Central America, specimens of this genus were present mainly in deeper lakes and predominate in conditions of higher conductivity (Wojewódka et al. 2016). Diatoms display a high degree of frustule dissolution. Poor diatom preservation was probably a result of high lake-water pH, a consequence of the local limestone geology (Hurd et al. 1981;Ryves et al. 2001). Moderately high TOC/N ratios (28-37) indicate relatively little aquatic organic matter (algal) input during the deposition of the massive silty clay sediments (Meyers and Teranes 2001). Thus, we propose that the massive silty clay facies are associated with maximum high-water stands, as observed in March 2018 (Fig. 9C.), when Lake Tzibaná covered an area of around 1.27 km 2 (Altmann 2020). However, the increase in water level should not have been greater than 15 m, as suggested by the presence of Darwinula, a benthic ostracode that lives not deeper than 15 m (Pérez et al. 2015;Díaz et al. 2017). In addition, the TIC shows intermediate values (6-8%) in comparison with the other facies. The surrounding karst bedrock is the ultimate source of carbonate, which is derived from detrital sources, biogenic processes (photosynthetic withdrawal of CO 2 and shell formation), and evaporation-driven precipitation. Smear-slide analyses indicate that detrital material was the major carbonate source in the massive silty clay facies in both sequences suggesting detrital input from runoff. The higher TIC values in the bank Ba1 compared to those of the inactive channel Cin-A (Fig. 5B) indicate that the Ba1 outcrop received more carbonate input via detrital deposition or other processes.

Potential and limitations of deltaic records for paleolimnological studies
Our sedimentological analysis of outcrops on the Nahá River delta in Lake Tzibaná illustrates both potential and limitations of investigating delta sediments for paleolimnological and paleohydrological studies. As recognized earlier, delta deposits may provide relevant general information of sedimentary processes influencing lakes (Machlus et al. 2000). Although deltaic sediments do not necessarily provide a continuous history of the lake evolution, they may help understanding the effect of phenomena, such as extreme desiccation events, avulsion of riverain tributaries or changes in the delta morphology, on the sediment record at distal locations (Blair and McPherson 2008).
Such distal sediment cores in the center of lakes are usually used to study droughts and desiccation of lakes (Marshall et al. 2011). As argued above, sedimentation at distal sites may eventually be interrupted during extreme desiccation events when the absence of a water column and subaerial exposure result in a discontinuity in the paleoenvironmental record. Several records in Central Mexico are examples of loss of sediment during the Last Glacial Maximum (20 kyr cal BP) as a response to low sedimentation rates and high evaporation (Lozano and Ortega 1998;Wood and Scholz 2017). In this case, the deltaic sediments might be able to provide information that lacks in the distal sedimentary sequences.
The effect of avulsion (switch of location of a river channel) in lacustrine deltas has received few attention in scientific literature (Li et al. 2022). For marine deltas, the deposition of multiple delta lobes on coastal plains during periods of relative sea level decline has been described as a consequence of main tributary avulsion . However, there is little information on the response of the sedimentation in distal zones as function of the avulsion processes that occur in the deltaic zone (Morozova and Smith 2000). Multiple avulsions in the same delta favor the formation of new lobes providing different sources of sediment input to the lacustrine system (Jerolmack 2009). A change in the sediment input zone to the lake can alter the rate of sedimentation at distal locations. Furthermore, the rate at which the lake level declines may affect the avulsion response in the deltaic zone. A rapid drop in the water level favors stronger and faster erosion of the walls of the active channels in the delta, consequently, there will be higher sedimentation rates. In a distal lacustrine record, an increase of the sedimentation rate would usually be interpreted as an increase in runoff typically associated with wet periods and, thus, high water levels; an interpretation that contradicts the actual process. In addition, changes in the delta morphology (e.g. slope, active subaqueous channels or direction of the deposition of a delta lobe) can be the result of a short event (e.g. earthquakes). The effect of sudden events on the sedimentation rate in the center of the lake has been studied in Lake Geneva (Switzerland/France), where the link between sediment cores and deltaic deposits allowed to analyze the respose of both archives to flooding events and earthquakes (Kremer et al. 2015). There are several surficial opened basins worldwide where deltaic sediments could provide additional hydrological information on past water-level variability such as Lake Tzibaná and other like Lake Petén Itzá (Hodell et al. 2008) or Lake Izabal, the largest freshwater ecosystem in Guatemala (Obrist-Farner et al. 2019;. Our analysis illustrates that a large amount of organic matter from terrestrial sources (e.g., leaves, branches, and trunks) is retained in both the active and inactive channels of the river delta. Inverted radiocarbon dates of the dark leaf litter horizons imply that old organic material is eventually eroded and then re-deposited above younger sediments, mainly when the channels of the delta are temporarily activated during the onset of a desiccation event. There is no reason to assume that the re-deposition of old organic matter is limited to the delta. Instead, it is likely that part of this organic material is also transported to and deposited in the deep, distal zones. Consequently, radiocarbon dates from organic remains found in distal sediment cores should also be interpreted with caution, especially when delta deposits indicate the action of strong erosion and re-deposition processes e.g. during large desiccation events as the one in 2019. A respective analysis of radiocarbon dates from distal sediment cores from Lake Tzibaná should explore this possibility in future research.
Instead of radiocarbon dating of organic material, OSL (Optically Stimulated Luminescence) dating of mineral grains could further improve the control of the timing of depositional processes. OSL measures the time elapsed since the last exposure of the mineral to sunlight. Several studies have used OSL in estuarine and delta sediments showing that this method provides generally reliable ages (Madsen et al. 2007;Shen and Mauz 2012;Gao et al. 2017;Xu et al. 2020). However, there are some limitation that should be considered including: 1) the incomplete bleaching of the OSL signal caused by the attenuation of the solar spectrum in the water environment and 2) different grain sizes had a different bleaching history being finer grains better bleached.

Conclusions
Sediment records from the Nahá River delta, located in the southern part of Lake Tzibaná, Lacandon Forest, Mexico, provide evidence of recurring variations in the lake water levels, and display four distinguishable facies: (1) massive-coarse sand, (2) fine sand, (3) dark leaf litter and (4) massive silty clay. We propose that the massive-coarse sand facies were produced by two processes: (a) reworking of sediments along the shoreline and a sudden entrance of sediments from the delta front or (b) by lateral deposition in former channels of the Nahá River. Both processes are linked to prolonged reductions in the lake water level similar to that observed in August 2019. Alternating layers of fine sand and dark leaf litter facies have likely been produced in low-energy environments during flooding of currently inactive channels cutting across the delta, also indicating low lake-water levels, similar to that observed in February 2019. Massive silty clay facies are deposited in low-energy environments and thus linked to high lake-water levels. Radiocarbon data do not follow a stratigraphic order and range from CE 1158 to 1940, indicating multiple, protracted waterlevel variations in Lake Tzibaná that produces an activation of deltaic channels and the repeated deposition, erosion and redeposition of leaf litter material. Electrical resistivity data indicates a maximum thickness of the sandy delta deposits of up to 20 m, of which the present study only assessed the uppermost 1.5 m. Field observations right after the desiccation event in 2019 suggest that lake water drained through a series of crevices in the underlying limestone bedrock. Although further hydrogeological studies need to be conducted to reveal the specific causes of modern and past lake water-level variations, this study reveals that similar desiccations of Lake Tzibaná have occurred repeatedly in the past, and are likely to recur in the future. The analysis of the outcrops of the Nahá River delta is an example of how deltaic sediments may produce relevant information to the comprehensive understanding of the distal sedimentary records, particularly analyzing the response of the lake to particular events.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.