Stalagmite evidence for Early Holocene multidecadal hydroclimate variability in Ethiopia

Abstract

A multiproxy oxygen and carbon isotope (δ¹³C and δ¹⁸O), growth rate, and trace element stalagmite paleoenvironmental record is presented for the Early Holocene from Ethiopia. The annually laminated stalagmite grew from 10.6 to 10.4 ka and from 9.7 to 9.0 ka with a short hiatus at ~9.25 ka. Statistically significant and coherent spectral frequencies in δ¹³C and δ¹⁸O are observed at 15–25 and 19–23 years, respectively. The observed ~1‰ amplitude variability in stalagmite δ¹⁸O is likely forced by nonequilibrium deposition, due to kinetic effects during the progressive degassing of CO₂ from the water film during stalagmite formation. These frequencies are similar to the periodicity reported for other Holocene stalagmite records from Ethiopia, suggesting that multidecadal variability in stalagmite δ¹⁸O is typical. Several processes can lead to this multidecadal variability and operate in different directions. A hydroclimate forcing is likely the primary control on the extent of the partial evaporation of soil and shallow epikarst water and associated isotopic fractionation. The resulting oxygen isotope composition of percolation water is subsequently modulated by karst hydrology. Further isotopic fractionation is possible in-cave during nonequilibrium stalagmite deposition. Combined with possible recharge biases in drip-water δ¹⁸O, these processes can generate multidecadal δ¹⁸O variability.

INTRODUCTION

A number of major air streams and convergence zones influence the modern climate in Ethiopia and the larger Horn of Africa region (Nicholson, 2017). Rainfall amount and intensity in Ethiopia is determined by the annual migration of the African rain belt, which is associated with the movement of the Intertropical Convergence Zone (ITCZ). The annual migration of the ITCZ determines the onset, duration, and termination of the East African monsoon, leading to a strongly bimodal annual cycle, resulting in two rainy seasons: the “big rains” or summer rains (between June and September), which are dependable and whose maxima migrate with the position of the ITCZ, and a second rainy season, the “small rains” or spring rains, which are less consistent and occur between March and May with maxima in April.

In addition, east–west adjustments in the zonal Walker circulation regulated by the El Niño–Southern Oscillation (ENSO) and the Indian Ocean Dipole (IOD) cause short-term (annual to decadal) fluctuations in the intensity of precipitation in Ethiopia. These are possibly a direct response to sea-surface temperature variations in the Indian and Atlantic Oceans, which are in turn affected by the ENSO and the IOD (Nicholson, 2017; Taye et al., 2021). While the global-scale atmospheric circulation patterns determine the rainy seasons in Ethiopia, local rainfall distribution is modulated by topographic features such as the highland barriers separated by the rift zone (Asrat et al., 2018).

Nearly 80% of the more than 100 million people inhabiting Ethiopia depend on rain-fed agriculture for their subsistence. Both the summer and spring rains in most parts of the country are important for adequate and sustained harvest. However, interannual variability is higher for the spring rains than the summer rains (e.g., Viste et al., 2013), and failure of the spring rains is common (Diro et al., 2008). Failure of the spring crop usually leads to reduced annual productivity (McCann, 1990) and in most cases leads to famine, at least in some worst-hit parts of the country, such as in 1984 and 2009, the two driest years since 1971 (Viste et al., 2013). The southeastern Ethiopian lowlands were affected by failure of the spring rains as recently as the 2013–2014 and 2015–2016 growing seasons.

There has been a general decline in the reliability of the spring rains since 1979 (e.g., Williams and Funk 2011; Viste et al., 2013), and data on the failure of the spring rains for the modern era suggest this occurs at a decadal frequency. For instance, within the 1995–2010 period, Viste et al. (2013) identified a cluster of dry spring seasons nationwide in 1999–2004 (except 2001) and in 2008–2011. The causes for the failure of the spring rains remain unclear. However, some studies (e.g., Segele et al., 2009; Williams and Funk, 2011; Viste et al., 2013) agreed that the failure is usually associated with deflections of the transport of moisture to Ethiopia due to atmospheric circulation anomalies. For instance, the 2009 spring drought was largely attributed to the deflection by southwesterly anomalies of the easterly flow bringing moisture from the northern Indian Ocean and the southeasterly flow bringing moisture from the southern and equatorial Indian Ocean (Viste et al., 2013).

Paleoclimate records provide a useful insight into the processes determining rainfall climate variability (Bar-Matthews et al., 1997; Hu et al., 2008), such as the decadal frequency of failure of the spring rains described earlier. For Ethiopia, annually laminated records such as those widely present in stalagmites have the necessary temporal resolution to investigate past multidecadal climate variability (Asrat et al., 2007, 2018; Baker et al., 2007, 2010). Previous research has shown that the strong seasonality of rainfall leads to the ubiquitous formation of annual growth laminae (Asrat et al., 2008). The warm climate leads to a fast stalagmite annual growth rate of about 100 to 500 μm/yr (Asrat et al., 2008; Baker et al., 2021), permitting high-resolution geochemical analyses. Tectonic activity associated with the East African Rift System adjoining the cave sites leads to discontinuous stalagmite deposition rarely lasting more than 1000 years, with stalagmites often having distinctive cone-shaped morphologies indicative of a drainage of a water source (Asrat, 2012). Two discontinuously forming Early to Middle Holocene stalagmite records from the Mechara caves (Ach-1 and Bero-1 stalagmites) have previously exhibited multidecadal variability in δ¹³C and δ¹⁸O and growth rate (Asrat et al., 2007; Baker et al., 2010). However, multidecadal variability in speleothems can be climatically forced, can derive from the inherent nonlinear properties of karst hydrology, or can arise from a combination of the two; for example, nonlinear karst processes amplifying the signal from extreme climate events (Baker et al., 2012).

Multi-stalagmite and multiproxy analyses are essential for investigating the reproducibility of paleoclimate records in speleothems (Hellstrom and McCulloch, 2000; Dorale and Liu, 2003). Here, we present a third high-resolution stalagmite paleoclimate record for the Holocene from Achere Cave, southeastern Ethiopia. Stalagmite Ach-3, which formed in the Early Holocene, is dated by U-Th series and annual laminae and analyzed for δ¹³C and δ¹⁸O and trace elements. We use a time-series analysis to investigate the multidecadal geochemical proxy signal in the stalagmite and compare this with other Middle and Late Holocene stalagmite records from the region.

METHODS

Site description

The Achere cave forms part of the bigger Achere–Aynage cave system and has been previously described (Asrat et al., 2007, 2008). The Achere–Aynage cave system developed along numerous NE-SW oriented parallel rifts on the southeastern Ethiopian highlands, close to the Main Ethiopian Rift (MER), indicating their development and modification through time in close association with rift-forming processes (Fig. 1). The maze-like cave network developed within a narrow, 20–25 m vertical zone, parallel to the bedding of Jurassic limestone. A laterally extensive calcareous mudstone/marl horizon within the limestone currently marks the roof of the cave chambers (Brown et al., 1998; Asrat et al., 2007, 2008).

Figure 1. (A) Regional structural setting of Ethiopia showing the location of the Mechara caves. The epicenters of the major earthquakes in the Main Ethiopian Rift and the adjoining highlands are marked (note that earthquake epicenters in the northern Afar depression are not represented). Insets show the mean position of the Intertropical Convergence Zone (ITCZ) in July (summer) and January (winter) over Africa; and the mean monthly rainfall (RF in mm) and mean monthly temperature of the Mechara region, at the Bedesa Meteorological Station (1994–2014 data from the Ethiopian Meteorological Agency). Location of B is marked by a solid rectangle around the location of Mechara; (B) The topography, geology, structure, and drainage system of the Mechara karst area and locations of the entrances to the caves (including Achere–Aynage and Bero). (C) Achere–Aynage cave survey showing the location of stalagmite Ach-3 and a previously published stalagmite, Ach-1. A and B modified from Asrat et al. (2008, 2018, respectively); C modified from Brown et al. (1998).

The aquifer architecture and hydrological flow regimes above the caves are a strong reflection of the tectonic–lithological interaction, which has been changing through time, even within the time frame of a single speleothem growth. Active tectonics in many cases is responsible for developing and continuously modifying the fracture systems that usually refocused groundwater flow paths along newly formed or reactivated fractures and conduits, in many cases leading to the cessation of growth of speleothems, manifested in growth hiatuses (Asrat, 2012). The location of the Mechara caves in close proximity to an active seismic zone of the MER (see Fig. 1) is also manifested in the uniquely short growth phases of stalagmites from Mechara (with median growth duration of 172 years) compared with the median growth duration of 447 years of annually laminated stalagmites globally (Baker et al., 2021).

The limestone terrain in the Mechara area, including the top of the limestone beds forming the Achere–Aynage caves, are overlain by very shallow (generally less than 50 cm deep) soils composed of lime-rich, soft calcareous layers overlain by dark organic-rich humus layers, classified as rendzinas (Bruggeman, 1986). In the wider area, chromic cambisols develop over the sandstones and shales, which form low hills above the limestone sequence. These soils are in most parts strongly eroded (Asrat et al., 2008).

The Mechara area is currently agricultural, with the land above the caves dominated by cultivated fields of teff (Eragrostis tef, a grain native to Ethiopia), maize (Zea mays), and millet (Panicum miliaceum), perennial cash crops like khat (Catha edulis) and coffee (Caffea sp.), and scattered patches of trees and scrub (Blyth et al., 2007). Though no vegetation history of the Mechara area in particular exists, the southeastern Ethiopian highlands were dominated by woody vegetation cover during the Early Holocene (Umer et al., 2007).

The Mechara area, at an altitude of 1500–1800 m above sea level, is characterized by an average annual temperature of 21°C and mean annual rainfall of ~1000 mm (see Fig. 1). Temperature is generally constant, except in the months of November to January, when it is ~2°C lower than the annual average. Precipitation is bimodal and shows strong seasonal variation: the main rainy season extends from June to September (big rains), with an average rainfall of ~160 mm/month, and the small rains fall between March and May, with an average rainfall of ~100 mm/month (Asrat et al., 2008). The small rains typically represent just 25–35% of total annual rainfall, with a total range of 15–43% (data from 20 years of complete data since 1984, Bedessa meteorological station; Fig. 1). The ratio of precipitation to potential evapotranspiration (P/PET; i.e., the aridity index) in the Mechara region is calculated to be 0.86 (FAO New_LocClim) or 0.88 (Wagari, 2005).

The precipitation δ¹⁸O record from the only long-term monitoring station at Addis Ababa shows that there is little seasonal variability in the modern δ¹⁸O (e.g., Baker et al., 2010). The isotopic composition of precipitation in July and August, the peak of the summer (big) rains, has δ¹⁸O, which is more negative than April small rains by ~3‰ (Baker et al., 2010). A recent study on δ¹⁸O and δ²H of precipitation samples collected at daily, weekly, and monthly intervals in different parts of Ethiopia representing local climate regimes confirmed the weak correlation between rainfall amount and δ¹⁸O values of precipitation (Bedaso et al., 2020). The same study further indicated the absence of discernible source region variability among the different stations. The mean moisture back-trajectory paths show the Mechara caves on the southeastern Ethiopian highlands receive most of their moisture from the southwestern and northern Indian Ocean, on southerly and easterly wind trajectories, respectively.

Asrat et al. (2008) reported cave monitoring data indicating that the Ach-3 stalagmite grew in a cave that has a nearly constant within-cave temperature of ~20.5°C. The cave has modern relative humidity of 87.5 ± 11.5% (number of measurements, n = 14) and within-cave pCO₂ content of 745 ± 365 ppm (n = 15). Drip waters in the cave have Ca²⁺ and Mg²⁺ concentrations of 3.13 ± 1.88 mmol/L and 0.66 ± 0.57 140 mmol/L (n = 12), respectively. Compared with the range of drip-water Ca²⁺ concentration (2.63 ± 2.36 mmol/L) in all the monitored caves in Mechara, the Achere cave drip waters have distinctly higher Ca²⁺ concentration, implying “open system” evolution (Baker et al., 2016), wherein the calcareous (limestone, marl, and carbonate-rich mudstone) aquifer readily contributes Ca²⁺ ions to the drip waters, and likely lead to rapid calcite formation that could be out of isotopic equilibrium. Limited cave drip-water oxygen isotope data from Achere cave demonstrate a limited range of δ¹⁸O composition from –1.6‰ to –0.5‰ (n = 10) (Asrat et al., 2008).

Sample description

The Achere–Aynage cave system contains abundant speleothems. The Ach-3 stalagmite was sampled in Achere cave in April 2004 from a narrow chamber leading to the bigger Moenco Chamber (where Ach-1 was sampled; Asrat et al., 2007), about 200 m from the cave entrance. Ach-3 developed on a low, narrow ledge 2 m beneath a roof marked by a mudstone layer. The chamber was dry, and the speleothem was inactive at the time of sampling, though some soda straw stalactites in the vicinity of the chamber indicate recent seasonal dripping.

Ach-3 is a 420-mm-long, slender stalagmite, narrowing from the bottom (120 mm diameter) to the top (60 mm diameter; Fig. 2). The stalagmite was sectioned into two halves, and one half was polished and scanned at high resolution and used to conduct lamina counting in triplicate with Image analysis software (Image-Pro 5, Media Cybernetics). The laminae show similarity to calcite layers in other speleothems in the region such as Bero-1 and GM-1 (Baker et al., 2010; Asrat et al., 2018). Continuous laminae of calcite were visible throughout the sample, marked by changes in calcite fabric, alternating between brownish dense and white porous calcite layers (Fig. 2). Some slight shifts in the growth axis mark the position of one of the growth hiatuses. The other half of Ach-3 was continuously milled down its long profile using a handheld dental drill (drill bit diameter = 500 μm) for δ¹³C and δ¹⁸O analysis at ~0.51 mm resolution (825 samples) and trace element analysis at ~4.6 mm resolution (91 samples). Additional samples for δ¹³C and δ¹⁸O were also drilled following some individual growth layers in order to perform the Hendy test. The fast growth rate of individual lamina of Ach-3 (with lamina width ranging between 200 μm and 1300 μm and average width of 490 μm), allows drilling of individual growth layers even at the flanks of the stalagmite. Seven samples for U-Th dating were similarly drilled using a dental drill, with samples located at the top and base of the stalagmite, on either side of possible growth hiatuses, and regularly spaced within growth phases (Fig. 2).

Figure 2. Ach-3 hand-section in both scanned image (left) and sketch (right), showing the four growth phases; locations of the major and minor growth hiatuses; and sampling for isotopes, trace elements, U-Th analyses, and U-Th ages. The central inset is a sample of a high-resolution scan (not to scale) along the central growth axis showing the annual laminae of Ach-3.

Geochemical analyses

Our methods follow those previously published in Asrat et al. (2007, 2018) and Baker et al. (2010). δ¹³C and δ¹⁸O were analyzed at the National Environmental Isotope Facility at Keyworth, UK. The calcite samples were reacted with phosphoric acid and cryogenically purified before mass spectrometry using an Isoprime plus multiprep dual inlet mass spectrometer. The Hendy test samples were analyzed at the University of New South Wales (UNSW Sydney) Analytical Centre using a MAT 253 mass spectrometer with a Kiel carbonate device. By comparison with laboratory marble standards KCM (Keyworth) and IAEA603 (UNSW), the sample ¹⁸O/¹⁶O and ¹³C/¹²C ratios are reported as δ¹³C and δ¹⁸O values in per mil (‰) versus VPDB. Analytical precisions are 0.07‰ for δ¹⁸O and 0.04‰ for δ¹³C on the standard marble (KCM) and 0.05‰ for δ¹⁸O and δ¹³C (IAEA603).

Trace elements were analyzed from 91 powders at UNSW, Sydney. Samples of approximately 0.05 g were dissolved in 1:1 hydrochloric acid, diluted, and analyzed for Ca and Mg using a PerkinElmer Optima 7300DV ICP-OES instrument. Ba, Sr, Al, Cu, Fe, K, Na, P, Pb, S, Zn, and U were analyzed using a PerkinElmer NexION 300D ICP-MS instrument.

Seven U-Th analyses were performed in the Uranium Series Chronology Laboratory, Institute of Geology and Geophysics, Chinese Academy of Sciences. The powdered subsamples of approximately 0.1 g were totally dissolved and spiked with a mixed ²²⁹Th-²³³U-²³⁶U. Uranium and thorium fractions were separated on 2 ml anion-exchange columns following standard techniques (Edwards et al., 1987). Then, the separated uranium and thorium solutions were measured on a multi-collector inductively coupled plasma mass spectrometer (Neptune Plus MC-ICP-MS instrument). The procedures followed those described in Cheng et al. (2013).

Time-series analysis

Variogram analysis on the annual growth rate time series was undertaken to determine the flickering parameter (f), information content (IC), and range (r) (Mariethoz et al., 2012). Flickering quantifies the growth rate acceleration from one year to the next through the lag-one autocorrelation of the detrended growth rate series, where 0 is no flickering (monotonous increases or decreases in growth rate) and −0.5 is the signal obtained from white noise (Baker et al., 2021). The observed flickering parameter (Mariethoz et al., 2012; Asrat et al., 2018) typically ranges between −0.5 and 0, the more negative f values indicating stronger flickering, interpreted as large changes in growth rate from year to year, indicative of a karst store filling and draining. To enable such a large interannual variability while maintaining continuous deposition over hundreds of years, a sufficiently large-volume karst store is hypothesized. Other statistical measures of information contained in the growth rate data are the variogram properties IC and r. IC quantifies the proportion of correlated signal in the time series as opposed to noise and varies between 0% (pure noise) and 100% (noiseless correlated signal). Range is the autocorrelated part of the signal, that is, the minimum timestep for which reliable variability might be observed from growth rate time series.

Stable isotope and annual growth rate time-series data were analyzed for their spectral properties. Spectral analysis was performed using the Spectrum software for unevenly spaced paleoclimate time series (Schulz and Statteger, 1997). Lomb-Scargle Fourier transforms were conducted, with five windows used (Bartlett, Hanning, Rectangular, Welsh, and Triangular) in order to undertake the spectral analysis of oxygen, carbon, and growth rate time series and the coherency between isotopic time series. The autocorrelation of the stable isotope time series was investigated by determining the autocorrelation function.

RESULTS AND INTERPRETATION

Chronology

Stalagmite Ach-3 is laminated throughout with 925 laminae. The triplicate lamina counts show insignificant lateral thickness variation. In hand-section, likely growth hiatuses with dissolution features were visually identified at lamina number 675 (growth hiatus 1, which is also marked by a slight shift in the growth axis) and 243 (growth hiatus 2) from the top, separating the sample into three growth phases (Fig. 2): growth phase 1 (laminae 925–676); growth phase 2 (laminae 675–244), and growth phase 3 (laminae 243–1). A third possible growth hiatus was identified near the end of the third growth phase (growth phase 3a: laminae 243–28; growth phase 3b: laminae 27–1).

The results of seven U-Th analyses are provided in Table 1. An age–depth model confirming the three growth phases is given in Figure 3. A basal date of 10,643 ± 82 years was used as an anchor point to constrain the initial growth phase containing 243 laminae. The six other U-Th ages occurred in stratigraphic order from 9850 ± 68 years to 9045 ± 55 years. The stratigraphically youngest three ages are all very similar to one another, despite the presence of possible hiatuses, suggesting that the growth phase 3 was of short duration. Following the method of Liu et al. (2013), for growth phases 2 and 3, the age–depth profile based on the laminae was aligned with that of the U-Th ages using two criteria: (1) for each growth phase, the mean age deviation between the two age–depth models was minimized; and (2) the age–depth models for growth phases 2 and 3 allowed for the observed hiatus between growth phases. The close agreement between the duration of stalagmite formation after hiatus 1 as determined by U-Th (the difference between the corrected U-Th ages ACH3-1 and ACH3-6 of 805 ± 93 years, 1σ) and the number of laminae (675 laminae) is indicative that the laminae of Ach-3 are annual in nature. This would agree with the widespread observation of annual laminae in other Ethiopian speleothems, which are due to the strong seasonality of rainfall with a distinct dry season (Asrat et al., 2007, 2018; Baker et al., 2007, 2010). Ach-3 laminae thickness has an average of 450 μm, and this is equivalent to the annual accumulation rate observed in Holocene and last interglacial Ethiopian stalagmites: Ach-1 (530 μm/yr); Bero-1 (450 μm/yr), Merc-1 (290 μm/yr); Asfa-3 (320 μm/yr) and GM-1 (440 μm/yr) (Asrat et al., 2007, 2018; Baker et al., 2007, 2010). We are therefore confident that the laminae are annual in nature.

Figure 3. An age–depth model for Ach-3. Depth measured as distance (mm) from the top of the speleothem. Locations of ages and hiatuses are marked.

Table 1. ²³⁰Th dating of stalagmite Ach-3 (2σ error).

^a δ²³⁴U = ([²³⁴U/²³⁸U]_activity – 1) × 1000.

^b Corrected ²³⁰Th ages assume the initial ²³⁰Th/²³²Th atomic ratio of 4.4 ± 2.2 ×10⁻⁶. Those are the values for a material at secular equilibrium, with the bulk earth ²³²Th/²³⁸U value of 3.8. The errors are arbitrarily assumed to be 50%. Bold type indicates dates used in constructing the age model. All years are since the year of analysis (2018).

^c δ²³⁴U_initial was calculated based on ²³⁰Th age (T), i.e., δ²³⁴U_initial = δ²³⁴U_measured × e ^δ234×T.

Geochemical proxies

The 825 δ¹³C and δ¹⁸O analyses are presented in Figure 4A as scatter plots of oxygen versus carbon isotopes down the growth axis, as well as for analyses made along seven growth laminae (Fig. 4B) equivalent to the classic Hendy test (Hendy, 1971). Figure 4A shows that the two isotopes are positively correlated along the growth axis in all growth phases, except for growth phase 2, and Figure 4B shows that the two isotopes are positively correlated along all 260 sampled growth laminae, including those in growth phase 2. This correlation between δ¹³C and δ¹⁸O is similar to results reported for other Ethiopian stalagmites (Asrat et al., 2007, 2018; Baker et al., 2010) and demonstrates that deposition is not in isotopic equilibrium (Fantadis and Ehhalt, 1970; Mickler et al., 2006; Wiedner et al., 2008). The δ¹³C/δ¹⁸O gradient is between 3.0 and 3.5 along growth laminae, and for growth phases 1, 3a, and 3b is 2.1, 1.0, and 0.5, respectively, with no correlation between δ¹³C and δ¹⁸O in growth phase 2. These gradients observed in stalagmite Ach-3 are similar to the mean value of the gradient of δ¹³C/δ¹⁸O of 3.8 observed along vertical transects and 3.9 observed spatially across calcite deposited on glass plates by Mickler et al. (2006). These were attributed to kinetic fractionation during calcite deposition out of isotopic equilibrium due to ¹⁸O and ¹³C Rayleigh-distillation enrichment in the HCO₃⁻ reservoir during progressive CO₂ degassing and calcite precipitation. They are also similar to the gradient of δ¹³C/δ¹⁸O of 1.4 ± 0.6 for the fast degassing of CO₂ in carbonate precipitation experiments (Wiedner et al., 2008). Though the classic Hendy test might not be conclusive in predicting the equilibrium or nonequilibrium deposition of calcite (e.g., Dorale and Liu, 2003), our cave monitoring and modern speleothem records from the Mechara caves further confirm that calcite deposition out of isotopic equilibrium is likely for Ach-3. The lowest values of the predicted equilibrium calcite δ¹⁸O variations from measured modern drip-water δ¹⁸O data in various caves in the region are not observed in speleothem δ¹⁸O records, indicating calcite deposition out of isotopic equilibrium (Baker et al., 2007; Asrat et al., 2008). However, in Ach-3, we note a trend over time in the δ¹³C/δ¹⁸O gradient and extent of nonequilibrium deposition. In phase 1, the gradient is 2.1, and in the last years of deposition (phase 3), the gradient is 1.0 (phase 3a) and 0.5 (phase 3b), which could indicate a change in the extent or type of isotopic fractionation; for example, additional evaporative fractionation due to slower drip rates and/or increased kinetic fractionation due to increased drip-water pCO₂.

Figure 4. Scatter plots of δ¹⁸O vs. δ¹³C for: (A) each growth phase (numbers shown are slopes of best-fit lines) and (B) Hendy tests along growth laminae in stalagmite Ach-3. Note that similar nonequilibrium deposition was observed in Ach-1 and Bero-1 (Asrat et al., 2007; Baker et al., 2010).

Trace element data for the 91 samples are presented in Supplementary Table 1. Elements were normalized to calcium and analyzed using principal component analysis (PCA) (Supplementary Figure 1). Three components explained 80% of the variability in the data. PC1 (36% of the variance explained) correlated with the elements P, Na, K, and Zn; PC2 (22% of the variance explained) correlated with Mg, Sr, and U; and PC3 (22% of the variance explained) correlated with Fe, Al, Ba, and Pb. We interpret PC1 as soil- or cave sediment–derived elements, given the presence of nutrients and organic-associated metals (Borsato et al., 2007; Hartland et al., 2012). PC2 is interpreted as bedrock-derived dissolution elements, and PC3 as elements derived from sediment, colloidal, and particulate material (Borsato et al., 2007). Time series of the three principal components show that all three components have high scores at the start of growth and decline over the first growth phase (Fig. 5). PC2 then has a long-term decrease over the rest of the period of deposition, indicative of a decrease in bedrock-derived metals over time (Fig. 5). PC1 increases to its highest value, and PC2 increases by a lesser amount, over the last years of deposition, while at the same time PC3 decreases to its lowest score.

Figure 5. Time series of the first three principal components (PC1–PC3).

The time series for δ¹³C and δ¹⁸O are presented in Figure 6, together with representative trace element data for PC1 (P/Ca) and PC2 (Sr/Ca, Mg/Ca) and annual growth rates. The 825 isotopic analyses represent an approximately annually resolved record. In the first deposition phase, from ~10.6 to ~10.4 ka, there is a trend toward lower ratios in Mg/Ca and Sr/Ca and more negative δ¹⁸O, indicative of generally increasingly wetter conditions or a shorter vadose zone water residence time. Higher concentrations of elements derived from soil or cave sediment (e.g., P), soluble elements, and detrital material in the lowermost growth laminae suggest the flushing of these materials into the cave at the beginning of deposition.

Figure 6. Time series of growth rate and geochemical proxies in Ach-3: (A) annual growth rate, (B) δ¹³C, (C) δ¹⁸O, (D) Sr/Ca, (E) Mg/Ca, and (F) P/Ca.

Stalagmite deposition from ~9.7 to ~9.0 ka in growth phases 2 and 3a has a long-term trend to more negative δ¹³C and lower Sr/Ca and Mg/Ca. This could be indicative of the continuation of the trend to increasingly wetter conditions or a shorter vadose zone water residence time and decreasing prior calcite precipitation along the flow path over this period (Fairchild et al., 2000). Growth rates and oxygen isotope composition exhibit no long-term trend and instead have multidecadal variability.

Over the possible short-duration growth phase 3b at ~9.3 ka, that is, the last 28 years of deposition, geochemical trends reverse, with increasing PC1 (soil- or sediment-derived elements) and PC2 (bedrock-derived elements) and decreasing PC3 (colloidally transported elements) (Fig. 5), increases in δ¹⁸O, and an increase in growth rate. Taken as a whole, these are indicative of a change in hydrology. Similar changes in geochemical, growth rate, and isotopic trends have been observed previously at the end of stalagmite deposition during Middle and Late Holocene (Asrat et al., 2007, 2018) and interpreted as a change in hydrological regime as the hydroclimate dries; for example, disconnection from the soil-water store or decrease in fracture flow component. In these records, the role of active tectonics in controlling speleothem growth duration by changing the flow regimes has been common.

The mean δ¹⁸O composition of the Early Holocene Ach-3 (−5.86 ± 0.42‰) is more negative compared with all other modern (Merc-1: −1.22 ± 0.31‰; Asfa-3: −1.37 ± 0.37‰; Baker et al., 2007) and Middle to Late Holocene (Bero-1: −3.42 ± 1.45‰, Baker et al. 2010; Ach-1: −3.20 ± 0.35‰, Asrat et al., 2007) samples from the region. All the published stalagmite records have evidence of calcite deposition out of isotopic equilibrium. Assuming a similar extent of calcite deposition out of isotopic equilibrium in all the stalagmites, including Ach-3, this indicates that drip water was ~2‰ more negative in the Early Holocene (Ach-3: −5.86 ± 0.42‰) compared with the Middle Holocene (Ach-1: −3.20 ± 0.35‰; Bero-1: −3.42 ± 1.45‰).

Time-series analysis

A summary of the results of spectral analysis on both stable isotopes and growth rate time series and variogram analysis and flickering of growth rate time series is presented in Table 2. Full spectral analysis results are presented in Supplementary Table 2 and Supplementary Figure 2, and autocorrelation plots in Figure 7.

Figure 7. Autocorrelation functions for δ¹⁸O and δ¹³C.

Table 2. (A) variogram analysis for stalagmite Ach-3 for the three growth phases (phase 1: oldest; phase 3: youngest): range (r), information content (IC), and flickering (f); and (B) summary of geostatistical properties for Ach-3, Bero-1, and Ach-1: univariate spectral analysis, showing the dominant and statistically significant (from red noise) periodicities in the oxygen isotope, carbon isotope, and growth rate time series.^a

^a Oxygen and carbon isotope time series have coherent periodicities at 15–16 years (Ach-3), 16–17 and 25 years (Bero-1), and 16–17 years (Ach-1).

Variogram analyses revealed short periods of autocorrelations in the growth rate data, which means that periodicities on decadal timescales can yield meaningful climate information. These are the range, r = 28 years in growth phase 2, and a much shorter range of r = 12–13 years in growth phases 1 and 3 (Table 2A). The r values are low compared with a global analysis of the growth rates of laminated stalagmites in Mariethoz et al. (2012) and Baker et al. (2021), but similar to other Ethiopian samples. The information content, IC, in the growth rate time series ranges from 50% to 67%, highest and relatively similar in growth phases 2 and 3. An IC >50% means that the stalagmite growth rate data contain significant useful signal. An IC >50% and r < 150 years classifies Ach-3 as a “Type A” stalagmite of Mariethoz et al. (2012), which is likely to be suitable for interpreting multidecadal information, with the higher IC in phases 2 and 3 suggesting that these are less noisy. The presence of flickering, f, of −0.26 (phase 1), −0.37 (phase 2), and −0.34 (phase 3), is indicative of a water-filled store supplying the stalagmite with sufficient volume to maintain continuous deposition for at least several decades, with hydrologically controlled year-by-year variations in water level controlling interannual growth rate variations. Phase 1 of deposition has a lower IC and relatively short range, which suggests that the first growth phase contains the least climate information.

Inspection of the autocorrelation of δ¹³C and δ¹⁸O time series for each growth phase (Fig. 7) shows that the autocorrelation for both stable isotopes is similar for growth phases 1 and 3. Between growth phases, there is a slight decrease in autocorrelation from growth phase 1 to growth phase 3 and a slight decoupling of the δ¹³C and δ¹⁸O autocorrelation functions in growth phase 2. If soil processes were the dominant control on speleothem δ¹³C, the slow decomposition of soil carbon over years to centuries (Carlson et al., 2019; Markowska et al., 2019) would lead to a relatively constant soil carbon isotope composition, and the resulting speleothem would be expected to lead to a stronger autocorrelation in δ¹³C compared with δ¹⁸O. This is not observed in Ach-3. The lower autocorrelation of δ¹³C compared with δ¹⁸O in growth phase 2 agrees with the observed lack of correlation between δ¹³C and δ¹⁸O through time in growth phase 2 and a possible decrease in the extent or a change in the type of isotopic fractionation in this phase. Overall, the similarity in the autocorrelation functions of δ¹³C and δ¹⁸O, combined with the evidence of isotopic fractionation from the correlation between δ¹³C and δ¹⁸O over time and along growth layers, suggests the dominant control of in-cave isotopic fractionation processes on the composition of both δ¹³C and δ¹⁸O, strongest in growth phases 1 and 3.

Spectral analysis on the δ¹³C, δ¹⁸O, and growth rate time series is presented in Table 2B and Supplementary Figure 2. There are similar and consistent periodic components in the δ¹³C and δ¹⁸O time series at around 15–25 years and 19–25 years in all three growth phases. Bivariate analysis of δ¹³C and δ¹⁸O demonstrates a coherency at 16–17 years. In growth phases 1 and 3, these periodic components in the stable isotope time series occur at time periods greater than the value of r obtained from the growth rate data, suggesting an independent forcing mechanism is dominant. Evidence that isotopic fractionation is occurring during deposition and that this is likely to be from within-cave fractionation processes suggest that within-cave isotopic fractionation processes are the dominant driver of the observed multidecadal periodicity in the stable isotope time series. We consider this further in the “Discussion.” These within-cave isotopic fractionation processes can be climatically forced, and we cautiously interpret these spectral frequencies as representative of an indirect hydroclimatic forcing affecting in-cave isotopic fractionation processes. Spectral analysis on the growth rate time series demonstrates that there are no periodic signals shorter than the range, r, for all growth phases (Table 2B). Table 2B also presents the results of previously published spectral analyses on Holocene Ethiopian stalagmites, demonstrating a consistent multidecadal periodic signal in δ¹⁸O time series between different time periods and different caves.

DISCUSSION

Conceptual model of stalagmite deposition

We present a conceptual model of the hydrogeochemistry and associated stalagmite growth in Figure 8. Stable isotope and trace element geochemical data and time-series analyses, combined with our hydrogeological understanding of the unsaturated zone properties of the limestone (Asrat et al., 2007), suggest that stalagmite Ach-3 is fed by a mixture of diffuse flow, through porous limestone and calcareous mudstone, as well as solutionally enlarged fractures. The latter are relatively small in volume and more important than diffuse flow contributions, as indicated by the 28 year range in growth rate time series, as explained in the previous section. This is indicative of a relatively small water store that controls growth rate variability through limits on the extent of prior calcite precipitation in the fracture and can determine drip rate. Considering the whole period of stalagmite formation, trace element data identify an initial sediment- or soil-derived elemental signal, potentially indicative of an initial flush of trace elements from the soil or interactions with cave sediments, and a loss of this elemental signal in the last decade of deposition. The duration of this last growth phase is the same as the range in the variogram analysis of growth rate and consistent with the inferred small water volume of the karst fracture. δ¹³C and δ¹⁸O have very similar autocorrelation functions and coherent, periodic signals in the time series, and there is strong correlation between δ¹³C and δ¹⁸O along growth laminae and within growth phases. This indicates a common control on both isotopes of within-cave isotopic fractionation.

Figure 8. Conceptual model for the deposition of stalagmite Ach-3: (A) growth phase 1: initiation and flushing from soil dominating the flow; (B) growth phase 2: wet and continuous growth from full storage, with multidecadal variability due to within-cave processes (such as drip rate or water saturation); (C) growth phase 3a: similar flow conditions to those of phase 2 but with less water storage; and (D) major tectonic process leading to the redirecting of flow regimes and relocation of drip sources resulting in rapid shutoff and growth cessation. Cartoons modified from Asrat et al. (2007, 2018).

In growth phase 1, there is an initial input of soil- or sediment-derived material. There is low information content in the growth rate time series in this growth phase, indicating a relatively noisy signal due to the combination of growth rate controls from the initial flush of soil-derived material as well as a hydrological control. The periodic signal in the growth rate time series and range are identical at ~12 years, suggesting relatively limited water storage to the stalagmite during this growth phase (indicated by an empty reservoir in Fig. 8A). In growth phase 2, the best information content and largest range is observed, which we interpret as the karst store relatively full of water (full storage reservoir in Fig. 8B) compared with other growth phases. Decreasing Sr/Ca and Mg/Ca ratios over this growth phase further indicate increasing water availability. In this growth phase, the δ¹³C and δ¹⁸O data show some evidence that isotopic fractionation processes have less dominant control on isotopic composition than in the other phases. In growth phase 3a, the range in the growth rate time-series analysis decreases, but all other proxies are identical to phase 2 and indicative of persisting high water availability (half storage reservoir in Fig. 8C). Throughout these growth phases, there is a consistent multidecadal variability in δ¹³C and δ¹⁸O, which is interpreted as being forced by nonequilibrium deposition processes. Finally, in phase 3b, we have a 28 year period of deposition in which trace element data indicate a decrease or loss of soil connectivity. This results in an increase in growth rate until growth cessation (an empty reservoir in Fig. 8D). Given the preceding growth indicated progressive increases in water availability, we infer that tectonic activity disrupted the water flow path to the stalagmite between growth phases 3a and 3b.

Multidecadal variability in Ethiopian stalagmite δ¹⁸O

Multidecadal variability in δ¹⁸O, combined with the similarity in the autocorrelation functions of δ¹³C and δ¹⁸O and the correlation between δ¹³C and δ¹⁸O over time and along growth layers, suggests that the multidecadal variability in stable isotopes is due to changes in the extent of isotopic fractionation, through nonequilibrium fractionation processes, such as changes in drip rate or drip-water calcite saturation that control the extent of ¹⁸O and ¹³C enrichment in the HCO₃⁻ water film during progressive CO₂ degassing and stalagmite precipitation (Mickler et al., 2006; Scholz et al., 2009).

Spectral analyses on δ¹⁸O for the three stalagmites, Ach-3 and the previously published Bero-1 and Ach-1, demonstrate a multidecadal variability throughout the Holocene (Table 2; Supplementary Fig. 3). The amplitude of this variability is ~1‰. The dominant statistically significant frequencies are between 13 and 30 years. We observe spectral frequencies in this range in stalagmites from different sites with different hydrogeology and flow paths (Asrat et al., 2007; Baker et al., 2010). Depending on flow path, this multidecadal variability in Ethiopian stalagmite δ¹⁸O can derive either from water isotope fractionation processes or from a direct signature of the δ¹⁸O of precipitation. The former include the partial evaporation of soil and shallow epikarst water that may increase with drier conditions (drier = more positive δ¹⁸O), and within-cave fractionation due to changes in the extent of isotopic nonequilibrium during stalagmite formation (increased drip-water pCO2 = more positive δ¹⁸O isotopic composition). A direct signature of the δ¹⁸O of precipitation is also possible in cases with limited water mixing and a fast-flow component to the hydrology (wetter = more negative). Figure 9 quantifies these processes for the specific example of Ethiopian stalagmites.

1. 1. Precipitation δ¹⁸O (Fig. 9, process A): The summer (big) rains have more negative δ¹⁸O (by ~3‰) than the small rains. Low rainfall amounts during the small rains could lead to more negative recharge water δ¹⁸O, but as the small rains represent only about one-third of the total annual rainfall, any effect is expected to be less than ~1.2‰ in annual weighted mean isotopic composition of precipitation.

2. 2. Mixing in the karst (Fig. 9, process B): Recharge waters will likely mix with water of different ages, depending on the flow path and the presence and volume of any subsurface karst water stores, such as solutionally enhanced fractures. Where well-mixed water from a single store is the source of drip water and no soil or epikarst evaporation is significant, there will be a more negative δ¹⁸O signal deriving from the precipitation δ¹⁸O. Any changes in the annual mean δ¹⁸O of precipitation due to changes in the relative proportion of small and big rains (see point 1) will be decreased in amplitude due to the mixing of waters toward the long-term weighted mean δ¹⁸O of precipitation.

3. 3. Selective recharge (Fig. 9, process C): A single mixed store is a simplification of actual karst hydrology, where multiple water flow paths are more common (Tooth and Fairchild, 2003; Fairchild et al., 2006; Hartmann and Baker, 2017); for example, an additional fracture or bypass flow that allows a fast-flow, less mixed flow component. In these instances, a recharge bias in the δ¹⁸O signal may be preserved in the drip-water δ¹⁸O. In the global meta-analysis of drip-water δ¹⁸O, Baker et al. (2019) demonstrated drip waters that were up to 2‰ more negative than the annual mean of precipitation, most commonly observed in regions with very distinct wet seasons in otherwise water-limited environments. Considering the relatively high P/PET ratio (~0.86) of the region, cave drip waters in the Mechara area might be expected to be up to 1‰ more negative than the annual mean of precipitation due to selective recharge. In-cave fractionation processes could operate in the opposite direction to this effect (see point 5).

4. 4. Partial evaporation of water (Fig. 9, process D): Precipitation that contributes to the soil-water store, and in some cases the shallow epikarst water, can undergo evaporation, leading to the remaining water δ¹⁸O becoming increasingly isotopically positive (Cuthbert et al., 2014). Partially evaporated water having a more positive δ¹⁸O than the original precipitation may be subsequently recharged to the cave. In a global meta-analysis, Baker et al. (2019) identified the presence of drip water that was exceptionally up to + 2.8‰ compared with weighted mean precipitation δ¹⁸O, and for water-limited environments with P/PET similar to the Mechara region, up to + 1.7‰. Partially evaporated δ¹⁸O has previously been hypothesized as forming part of the δ¹⁸O in an Ethiopian stalagmite (Baker et al., 2010), where forward modeling for the modern growth phase of the Bero-1 stalagmite identified a positive isotopic offset of 2.0‰ to 2.5‰, attributed to evaporative fractionation processes between rainfall and the stalagmite. However, the effect of possible changes in the relative proportion of small and big rains on the partial evaporation of soil or epikarst waters is unclear. For example, if the small rains led to the recharge of more partially evaporated water than the big rains, due to relatively low rainfall amounts in the former, then relatively dry small rain seasons could lead to more negative drip-water δ¹⁸O.

5. 5. Nonequilibrium deposition (Fig. 9, process E): All stalagmites analyzed in Ethiopia to date demonstrate conclusive evidence of calcite deposition out of isotopic equilibrium. In Ach-3, there is strong correlation between δ¹³C and δ¹⁸O along growth laminae and over time, with δ¹³C/δ¹⁸O gradients <3. Bero-1 and Ach-1 also had δ¹³C/δ¹⁸O gradients <3. The similar range in δ¹³C/δ¹⁸O gradients of the three stalagmites with laboratory experiments (Wiedner et al., 2008) and the meta-analysis and field observations of Mickler et al. (2006), combined with the strong correlations between δ¹³C and δ¹⁸O for each stalagmite, and similar and coherent multidecadal spectral frequencies between δ¹³C and δ¹⁸O, suggests a dominant in-cave control. One such mechanism is a change in drip rate, which controls nonequilibrium isotopic fractionation during the progressive degassing of CO₂ from the water film during stalagmite formation. All three stalagmites have similar amplitude in multidecadal signal (up to ~1‰). The iSOLUTION model of oxygen and carbon isotope composition of stalagmite calcite (Scholz et al., 2009; Deininger and Scholz, 2019) models nonequilibrium isotopic fractionation processes and produces this magnitude of oxygen isotope fractionation for high pCO₂ drip waters and relatively slow drip rates. Kinetic isotopic fractionation due to rapid degassing from high pCO₂ drip waters could also lead to this magnitude of isotopic fractionation for faster drip rates (Mickler et al., 2006; Wiedner et al., 2008) and would be considered likely given the fast growth rates of Ethiopian stalagmites.

Figure 9. Isotopic composition conceptual diagram. The changes in oxygen isotope composition are based on observed Addis Ababa IAEA monthly δ¹⁸O precipitation (process A); observed global range of epikarst and soil evaporative fractionation (open arrow) and range for ratio of precipitation to potential evapotranspiration (P/PET) = 0.9 (filled arrow) (Baker et al., 2019) (process B); well-mixed drip-water δ¹⁸O (process C); observed global range of recharge bias (open arrow) and range for P/PET = 0.9 (filled arrow) (Baker et al., 2019) (process D); and modeled nonequilibrium fractionation factors (Scholz et al., 2009) (process E). WRT, with respect to.

We provide multiple lines of evidence that the multidecadal variability in stalagmite δ¹⁸O in Ethiopian stalagmites is likely due to a complex set of drivers such as the interannual variability in the relative amounts of small and big rains, karst hydrological processes on water mixing, evaporative fractionation of water in the soil, shallow vadose zone, or in the cave, preferential recharge, and isotopic fractionation processes operating with opposite signs in δ¹⁸O from the same climate forcing, as visualized in Figure 9. In years of decreased recharge, decreased drip rate to the stalagmites leads to the potential of increased isotopic fractionation due to calcite deposition out of isotopic equilibrium. Decreases in drip rate do not necessarily have a linear relationship with surface hydroclimate forcing, due to the nonlinear nature of karst hydrology and mixing of waters in karst stores and fractures. A recent global study of speleothem δ¹⁸O demonstrated that within-cave speleothem and drip-water δ¹⁸O variability are driven by karst hydrology due to the influence of fractures on flow paths (Treble et al, 2022). Our observation of multidecadal spectral frequency in δ¹⁸O is therefore likely to be due to individual extremes of dry years, which determine the volume of recharge to these karst stores, and in turn the drip rate from the store, including both the mean annual drip rate and/or the duration of dripping in one year. With drier conditions, in-cave isotopic fractionation and evaporative fractionation effects operate with the same sign, increasing drip-water δ¹⁸O due to increased evaporation at the same time as nonequilibrium deposition increased with lower drip rates. However, for some samples with a fast-flow or bypass-flow component, preferential recharge could be significant in controlling drip-water δ¹⁸O, and this signal could dominate over fractionation processes and generate a multidecadal signal with the opposite sign. Superimposed on all flow types is the possibility of kinetic isotopic fractionation due to high drip-water pCO₂, which is likely given the very fast growth rates of Ethiopian stalagmites.

CONCLUSIONS

We use trace element, growth rate, δ¹⁸O, and δ¹³C of Early Holocene stalagmite Ach-3 to understand the processes occurring during its deposition. The trace element composition identifies an initial growth period with a flush of soil-derived material and a final growth period in which there is a change in hydrology, indicative of drying conditions. We observe a multidecadal δ¹⁸O variability in the Early Holocene Ach-3 and the other two Middle and Late Holocene Ethiopian stalagmites of ~1‰ amplitude. Covariation of δ¹⁸O and δ¹³C demonstrates that all three stalagmites are dominated by isotopic fractionation, likely due to nonequilibrium effects during the progressive degassing of CO₂ from drip waters with a high pCO₂ during stalagmite formation. The amplitude of multidecadal variability in δ¹⁸O is similar to that modeled due to changes in drip rate. Rapid growth rates, fast drip rates, and isotopic fractionation effects are likely the primary controls on the isotopic geochemistry, while active tectonics has played an important role in determining the growth duration of the three Ethiopian stalagmites, with additional influences possible from evaporative fractionation and, for samples with very short water residence time, a small primary precipitation seasonality signal. Despite the extent of calcite deposition out of isotopic equilibrium, differences in mean stalagmite δ¹⁸O through the Holocene are larger in magnitude than the multidecadal variability. Thus long-term (centennial and longer) trends in stalagmite δ¹⁸O are likely to be good proxies for climate, as they record long-term climatic forcing on precipitation δ¹⁸O and drip-water δ¹⁸O.