Holocene hydroclimatic variability and multi-frequency analyses at Lake Sidi Ali (Morocco)

The North African desert margin is considered one of the most sensitive areas to future climate changes, yet the periodicities, coupling mechanisms and external forcing of Holocene environmental variability remain poorly understood. To investigate millennial- to centennial-scale periodicities in Holocene climatic variability and geomorphological processes, we use a Holocene sediment record from Lake Sidi Ali in the semiarid to sub-humid Middle Atlas with a robust ²¹⁰Pb $/$ ¹³⁷Cs and pollen-concentrates-based ¹⁴C chronology. We use a high-resolution core scanning-XRF record, in order to distinguish between lake-internal (e.g., chemical precipitation) and lake-external (e.g., detrital input) processes. Redfit and Wavelet time series analyses reveal distinct periodicities of millennial to centennial scale. By a correlation analysis of extracted, highly significant, frequency analysis spectra, three XRF-based “Redfit Proxy Groups” (RPGs) which potentially reflect different hydroclimatic forcing mechanisms were derived. Subsequently, we integrated environmental and climatic proxies from the same core (Cedrus pollen abundance, magnetic susceptibility, δ¹⁸O and δ¹³C values of ostracod shells, grain-size endmembers and total organic carbon) and used their wavelet domain to improve the interpretation. Finally, we identified two main periodicity regimes that affected, on the one hand, the hydrological regime and, on the other hand, the lake productivity and catchment erosion dynamics. For RPG 1 (Ca, Sr, $\mathrm{Ca}/\mathrm{Ti}$, $\mathrm{Sr}/\mathrm{Ti}$), we identified 2 and 1 kyr periodicities, which we interpret as precipitation/evaporation related proxies in the context of North Atlantic and solar forcing. For RPG 2 (Fe, Ti, K, $\mathrm{Si}/\mathrm{Ti}$), we observe 3.5 and 1.5 kyr periodicities, which we interpret as driven by lake productivity or detrital input. Overall, our results show that Holocene environmental variability at Lake Sidi Ali was structured by two partly decoupled regimes: hydroclimatic variability linked to precipitation/evaporation changes and catchment-related variability linked to erosion, productivity and terrestrial sediment supply.

1 Introduction

North-western Africa and the Western Mediterranean are known as regions highly sensitive to climate change (Alverson et al., 2003; Giorgi, 2006). The desert margin in North Africa was significantly affected by Holocene climate changes accompanied by landscape transformations (Fletcher and Zielhofer, 2013; Ruddiman, 2014). Especially, large-scale hydroclimatic forcing mechanisms in the context of outer-tropic North Atlantic originated winter rains and the spatiotemporal evolution of the African Humid Period (AHP) are under debate (Dallmeyer et al., 2020; deMenocal et al., 2000; Zielhofer et al., 2019b). The North Atlantic-driven winter-rain variations seem to be a controlling factor of the Holocene Western Mediterranean hydroclimate and respective environmental responses (Zielhofer et al., 2019b). Regarding the AHP, it is suggested that sub-tropical airmass variations, especially enhanced and northward-shifted African monsoon precipitation patterns are the major cause for more humid conditions in the Saharan realm (Braconnot et al., 2007; Shanahan et al., 2015). However, coeval enhanced Atlantic winter rains have been argued to be necessary to reproduce proxy moisture signals in Mediterranean North Africa (Cheddadi et al., 2021). Even if the monsoonal precipitation did not reach Mediterranean Northwest Africa, further hydroclimatic effects of subtropical air masses are possible. In this context, Holocene pollen records from the Moroccan Atlas Mountains indicate variations in summer drought severities that could be caused by subtropical summer temperature variations (Campbell et al., 2017).

However, for a better understanding of large-scale and external hydroclimatic forcing as well as the influences of different air masses, time-series and frequency analyses and the detection of periodicities are of particular importance in paleoenvironmental research (Azuara et al., 2020; Kern et al., 2013; Rösch and Schmidbauer, 2018). Frequency analyses break down the respective time series into individual subharmonics. This allows the identification of long-term and short-term periodicities and periodicity changes through time (Mudelsee, 2014). Additionally, frequency analyses of time-series data can help to develop a robust understanding of proxy behaviour and coupling the extracted periodicities with other archives and climate forcings (Debret et al., 2009; Mudelsee, 2019; Sabatier et al., 2020).

Robust frequency analyses already exist from the Mediterranean Northwest Africa to the southwestern High Atlas from speleothem records. There are two speleothem records from the Middle Atlas Mountains (Grotte du Piste and Chaara caves) that cover the entire Holocene (Ait Brahim et al., 2019; Wassenburg et al., 2016). These records have been interpreted in terms of winter precipitation changes with periodicities at millennial, centennial and sub-centennial scale. Further, Sha et al. (2019) presented a speleothem record from the SW High Atlas Mountains that seems not to be influenced by Atlantic air masses but rather tropical-monsoonal air masses. Other published speleothem records only cover the last ca. 1000 years (Ait Brahim et al., 2017, 2018; Wassenburg et al., 2013) and discuss centennial-scale periodicities as Atlantic air mass variabilities in the context of changes in the North Atlantic Oscillation. In this respect, lake sediments are highly complementary to speleothem records, as these archives reflect external hydroclimatic forcing mechanisms but also dust influx, catchment erosion, and lake-internal dynamics such as chemical precipitation and biomass productivity (Aufgebauer et al., 2012; Boës et al., 2011; Cohen, 2003; Conley and Schelske, 2001; Meyers, 2003; Neff et al., 2008; Zielhofer et al., 2017b). In this context, there is a lack of high-resolution lacustrine records in Mediterranean North Africa that cover the entire Holocene suitable for frequency analyses.

Regarding lacustrine records, high resolution (mm-scale) core scanning XRF (Croudace and Rothwell, 2015; Dunlea et al., 2020) data provide the opportunity for the detection of long-term (millennial) and short-term (sub-decadal to centennial) changes in lake-system behaviour (Lauterbach et al., 2011). Typically, measured elements and associated ratios can reflect a range of lake internal and external processes (Davies et al., 2015; Kylander et al., 2011). Cohen (2003) and Davies et al. (2015) present vast collections of XRF-based proxies and relate them to lake internal or external processes. However, multiple processes can influence the elemental composition such that solitary elemental proxies must be considered alongside other proxies to enhance their interpretation, including e.g., pollen, grain-size distributions, and δ¹⁸O and δ¹³C values on ostracod shells. Overall, multi-proxy approaches can help to enhance the robustness of palaeoclimatological and environmental interpretations (Martin-Puertas et al., 2012; Sanchini et al., 2020; Unkelbach et al., 2019).

In this study, we present a high-resolution, palaeolimnological and multi-element core-scanning XRF record supplemented by medium-resolution discrete proxy data from grain size endmembers, Cedrus pollen, δ¹³C and δ¹⁸O values of ostracod shells, magnetic susceptibility and total organic carbon from Lake Sidi Ali in the Moroccan Middle Atlas (Zielhofer et al., 2017, b, 2019b). The objectives of the study are (i) identification of statistically robust millennial- to centennial-scale periodicities in lake-internal and lake-external proxies using complementary stationary (Redfit) and non-stationary (wavelet) frequency analyses; (ii) grouping of XRF-derived proxies according to similarities in their frequency behaviour in order to distinguish dominant process domains controlling sediment composition; and (iii) linking these frequency-based proxy groups to hydrological, catchment-erosion and redox-sensitive processes through integration with independent non-XRF proxies. The overarching aim is to assess the temporal persistence of the identified periodicities throughout the Holocene and to explore their potential links to supra-regional climatic forcing, including North Atlantic variability and solar-paced cycles.

2 Study area

Lake Sidi Ali is located in the Middle Atlas in Morocco (33° 03^′ N, 05° 00^′ W, 2080 m above sea level), a transitional location between the wetter NW and arid SE regions of Morocco (Fig. 1a). The lake is considered to have developed through a combination of tectonic and karstic dynamics in the Pliocene-Quaternary period (Akdim, 2015). The catchment includes mountain ridges up to ca. 2350 m a.s.l. and is characterised by steep slopes with a generally sparse vegetation cover but also including notable stands of cedar (Cedrus atlantica) with individual trees exceeding 800 years in age (Copes-Gerbitz et al., 2019; Zielhofer et al., 2017a). The local climate is influenced by three major air masses (the Atlantic, the Mediterranean and the Saharan) along the NW African desert margin (Knippertz et al., 2003). Nowadays, the Middle Atlas is classified between a Cs and BSk Koeppen-Geiger climate. The regional climate is characterised by summer-drought conditions and humid spring and winter seasons driven by Atlantic and Mediterranean moisture transport. The study area has a mean annual precipitation of 487 mm and a mean annual temperature of 14.5 °C, a mean January temperature with 2.0 °C and a mean July temperature with 19.7 °C (Harris et al., 2020). The vegetation assemblages in the Middle Atlas are typically distributed along an altitudinal gradient, according to changes in temperature and moisture availability (Benabid, 1982; Campbell et al., 2017). The area around Sidi Ali is mainly dominated by Cedrus atlantica, Quercus rotundifolia, Juniperus thurifera and Crategus sp. (Benabid, 1982; Linares et al., 2011). However, the forest vegetation is degraded due to overgrazing and logging, and Cedrus atlantica shows clear signs of die-back (Cheddadi et al., 2022; Rhanem, 2011). Therefore, the current vegetation cover reflects a state of human disturbance (Copes-Gerbitz et al., 2019). Currently, the lake is separated into two sub-basins which are divided by a basalt ridge, which merge into a single, continuous lake at high lake level stands (Barker et al., 1994; Lamb et al., 1999). The main basin reached a water depth of 38 m in September 2012 (Zielhofer et al., 2017a). The closed catchment (ca. 17 km²) shows no significant surface inlet or outlet in the lake (Fig. 1b). Several hydro-limnological measurements show a clear stratification of the water column with anaerobic conditions in the hypolimnion in summer (Dumont et al., 1973; Zielhofer et al., 2017a). The alkaline lake (pH 9.1–9.7) has a specific conductivity between 1.2 and 1.6 mS cm⁻¹ at the lake surface and low Ca⁺ contents, which indicates evaporative conditions (Lamb et al., 1999; Zielhofer et al., 2017a).

Figure 1 Geographical, climatological and topographical context of the study area. (A) Location of Lake Sidi Ali in the Western Mediterranean/Northwest African context. The three major air masses influencing the Middle Atlas hydroclimate are indicated by the coloured arrows. (B) Local topographical context of Lake Sidi Ali based on SRTM (NASA Shuttle Radar Topography Mission (SRTM), 2013). The contour lines have a vertical spacing of 50 m. Isobath intervals are not evenly spaced. Both panels are underlain by Google Earth satellite images (Image data: © NASA 2024, maps data: © Google 2024) (EPSG: 4326).

3 Material and methods

3.1 Data acquisition

In 2012, a 19.63 m sediment record was recovered in 2 m core sections from the deepest position of Lake Sidi Ali using a floating platform and a UWITEC piston corer. The cores show slightly laminated gyttja facies without major sedimentological shifts (Zielhofer et al., 2017a). The cores were split and the working halves were analysed by ITRAX core-scanning XRF (cs-XRF) system (Croudace and Rothwell, 2015; Zielhofer et al., 2017a) (Table 1, Fig. 2a; 1 mm intervals). The cs-XRF proxies comprise 17 863 data points for each element. Subsequently, we sampled the working half in 1 cm intervals. A total of 198 discrete samples were prepared for benchtop XRF (air-drying, 2 mm sieving, homogenisation with the planetary mill Retsch MM 200), and pressed to pellets for 2 min under 20 t of pressure (Vaneox press) (Schmidt et al., 2023b; Zielhofer et al., 2017b). The pellets were measured with a Spectro Xepos X-ray fluorescence device in a helium atmosphere (Fig. 2a). The Bayesian age-depth model for the Holocene sequence of lake Sidi Ali is based on AMS radiocarbon dating of 26 pollen concentrate samples accompanied by ²¹⁰Pb and ¹³⁷Cs dates for the most recent part (Fletcher et al., 2017).

The magnetic susceptibility of the cores was measured on-site using the Bartington MS2C core logging sensor, with measurements taken at 2 cm intervals (Zielhofer et al., 2017a). Samples measuring 1 cm³ in volume were prepared for pollen analysis, spaced at regular intervals of 10 cm. This method yielded a total of 201 samples (Campbell et al., 2017). Grain sizes of 124 core samples were analyzed using a Malvern Mastersizer laser diffraction particle size analyzer. Subsequently, a robust end-member modeling analysis (EMMA), employing eigenspace analysis and scaling procedures, was conducted on all non-zero grain size classes across the 124 samples using the EMMAgeo package in R (Dietze and Dietze, 2019; Zielhofer et al., 2017b). Four to six adult ostracod shells (about 20 µg) of Fabaeformiscandona sp. and Candona sp. were picked for isotope mass spectrometer analyses (δ¹³C and δ¹⁸O values) (Zielhofer et al., 2019b, a). Total organic carbon (TOC) was quantified by subtracting total inorganic carbon (TIC) determined through Scheibler carbonate measurements from total carbon (TC) values from an Elementar CNS analyzer (Zielhofer et al., 2017a).

Table 1 Data overview with temporal resolution for the Sidi Ali Holocene record parameters and references of the previous publications of the data.

3.2 XRF calibration

Cs-XRF elemental values given as intensities were calibrated with depth-equal benchtop-XRF results given in mg kg⁻¹ using the log-ratio calibration model of the Xelerate software (Bloemsma, 2015; Weltje et al., 2015). This allows the more powerful combination of high resolution cs-XRF data with quantitative values from benchtop-XRF. For the primary selection of elements and ratios (Fig. 2a, Table S1), we used broad compilations from Cohen (2003) and Davies et al. (2015) and previously conducted XRF element proxies from Sidi Ali core (Zielhofer et al., 2017a, b).

Figure 2 Flowchart of the high resolution XRF methodological approach. (A) Acquisition of high resolution XRF data, calibration and forming elements and elemental ratios. ^*Proxy choice was done by common and published palaeolimnological XRF proxies (Cohen, 2003; Davies et al., 2015). (B) Processing of high-resolution data for time series (Redfit and Wavelet) analyses, by detrending and equidistant resampling (precondition for Wavelet analyses). (C) Extraction of significant (line above the critical Level of the AR1 red noise) Redfit power spectra of each proxy. (D) Correlation analysis of the significant Redfit power spectra using the Spearman coefficient. The crosses show the non-significant (0.05 level) correlation values.

3.3 Statistics

3.3.1 Data (pre-) processing

Large-scale cycle-like patterns (multi-millennial) that are not completely represented in the record can have impacts on the results from the Redfit and Wavelet analyses (Borradaile, 2003; Egan et al., 2022; Hochman et al., 2019; Rösch and Schmidbauer, 2018). For this reason, we detrended the data by fitting a 3rd-order polynomial regression to each proxy (Fig. 2b). The order of the model was chosen to ensure that variability with periods longer than approximately half the length of the record is removed prior to spectral analysis. Subsequently, we used the detrended time series for the Redfit analyses, as this method operates with unevenly spaced data (Schulz and Mudelsee, 2002). In contrast, the wavelet analysis requires equidistant spaced time series (Rösch and Schmidbauer, 2018). Debret et al. (2007) used a Spline interpolation for equidistant spacing, but we experienced “overoscillation” with artificial cyclicity for spline interpolations. Consequently, we applied a linear model on each time series using the average temporal resolution of each proxy (Table 1). The detrended and equidistant spaced time series act as input data for the wavelet analysis (Fig. 2b). Further, we smoothed the raw time series data by a kernel function (Härdle and Vieu, 1992) with a 500-year span, in order to extract the major positive and negative peaks in the time series for the comparison of the timing of the cycles (Fig. 2b, red line).

3.3.2 Principal Component Analyses

A principal component analysis (PCA) reduces the dimensionality of multi-proxy datasets and is widely used in the analysis of XRF elemental data (Bertrand et al., 2024; Evans et al., 2019). In this approach, the main trends within the dataset are extracted as continuous axes (Borcard et al., 2011). These axes are calculated based on the eigenvectors of the variables and their Euclidean distances, resulting in orthogonal components that are statistically independent and uncorrelated. Two PCA analyses were conducted. The first analysis PCA_freq used the extracted spectral information from the Redfit analyses as input data. The second analysis PCA_time was based on the raw proxy data in the time domain. The results of both analyses are presented in the Supplement (Figs. S7 and S8).

3.3.3 Redfit analysis

The Redfit algorithm works with a Lomb-Scargle Fourier transformation (Lomb, 1976; Scargle, 1982) to overcome unequal time intervals (Li et al., 2021; Schulz and Mudelsee, 2002). The process uses the Hanning window and three 50 % overlapping segments (Bunn and Korpela, 2021; Schulz and Stattegger, 1997). The raw spectrum is smoothed, segmented, and linearly interpolated to eliminate artificial high frequencies and make it consistent (Ólafsdóttir et al., 2016). Finally, the averaging of all segments results in the periodogram. The significance of the periodogram was tested against the AR1 (first-order autoregressive) process, also called red noise. Redfit uses Monte-Carlo simulations to create a set of AR1 processes and extracts various significance levels. We used critical false alarm level, that is determined by the properties of the data and therefore may be more conservative than 90 %, 95 % or 99 % confidence (Schulz and Mudelsee, 2002; Thomson, 1990). The dplR package (v1.7.1) (Bunn and Korpela, 2021) within the R (+RStudio) environment (R Core Team, 2020) was used.

3.3.4 Wavelet analysis

The wavelet analysis (Torrence and Compo, 1998) also creates a periodogram, but the power of the cycles is linked with the temporal evolution. In contrast to Redfit, wavelet analysis reveals possible non-stationarities of the data (Debret et al., 2007). Periodic components get combined, and cycles and amplitudes can be displayed separately (Torrence and Compo, 1998). This procedure allows the detection of periods occurring only for a part of the time series (Debret et al., 2007). Here, the resulting wavelet power spectra is also tested against the AR1 process to avoid artefacts from the persistence of the data (Rösch and Schmidbauer, 2018). During the Fourier transformation, the length of the time series is increased by adding zero values to avoid false periodic events. However, this creates boundary effects as the wavelet gets closer to the edge. We avoided interpreting data within the area with increased edge effects (cone of influence) due to low accuracy of spectral information (Cazelles et al., 2008). We use the waveletComp package (v1.1) implemented by Rösch and Schmidbauer (2018) within the R environment. The algorithm uses the Morlet mother wavelet, because it keeps its shape through frequency shifts, and therefore, provides a reasonable separation from different frequency band contributors without excessive loss in temporal resolution (Cazelles et al., 2008; Goupillaud et al., 1984). The wavelet function needs equidistant spaced data (Fig. 2b) and we used a significance level of 95 %.

3.3.5 Correlation analysis of significant Redfit spectra

So far, cross-spectral analyses are typically bivariate, i.e. performed on pairs of proxies (Mudelsee et al., 2020; Ólafsdóttir et al., 2016). In order to reduce the dimensions (number of variables), we developed a procedure to group the proxies by using the significant spectra of the periodograms of the Redfit analyses. First, we extracted the spectrum of each XRF proxy from the Redfit periodogram. In a second step, we removed the non-significant parts of the spectra (Fig. 2c). The correlation analysis was done with the significant part (the values below the critical level were removed/ set to NA – “Not Available”). Due to the precondition of zero NA values throughout all variables (proxies) for a correlation analysis, we could only focus on periodicities between 1 and 4 kyr. For power values of higher frequencies, at least one variable had NA values. However, the correlation matrix (Fig. 2d) shows the correlation coefficients after Spearman and the configuration of the proxies along the axes was done by similarities of positive or negative rho values. We additionally performed a PCA in the frequency domain (PCA_freq) to evaluate the robustness of the grouping approach. The PCA results (Fig. S7) independently confirm the proxy groups identified by the correlation matrix (Fig. 2d).

4 Results

4.1 Dominant periodicities and inter-proxy similarities derived from Redfit analysis

The results of the Redfit analysis for 18 XRF proxies show major periodicities at ca. 3.5, 2, ca. 1.5, ca. 1 kyr and several in multi-centennial scale (Figs. S1–S3). Two major groups with similar periodicity patterns are visible. As an example, the Redfit periodograms show similar periodicities for Ca and $\mathrm{K}/\mathrm{Ti}$ (2, ca. 1 kyr and some in multi-centennial scale; Fig. 3a and b) and S and Ti (3.5, 1.5 kyr and less in multi-centennial scale, Fig. 3c and d).

Figure 3 Redfit periodograms with significance levels of selected and calibrated cs-XRF proxies (a) Ca (RPG 1), (b) $\mathrm{K}/\mathrm{Ti}$ ratio (RPG 1), (c) S (RPG 2) and (d) Ti (RPG 2). Black dots with numbers show local peaks with periodicity in years that are significant at the highest critical level. Note: x-axis is on logarithmic scale.

The correlation matrix of the significant Redfit spectra (Fig. 2d) shows clearly three Redfit Proxy Groups (RPG) of proxies with very similar spectral behaviour derived from the Redfit analyses. In general, the in-group Spearman coefficients are very high and out-group values very low. The first RPG is dominated by Ca, Si, Sr, $\mathrm{Ca}/\mathrm{Fe}$, $\mathrm{Ca}/\mathrm{Si}$, $\mathrm{K}/\mathrm{Ti}$, $\mathrm{Sr}/\mathrm{Ca}$ and $\mathrm{Sr}/\mathrm{Ti}$. The second RPG consists of Fe, K, S, Ti, $\mathrm{Fe}/\mathrm{Ti}$, $\mathrm{Mn}/\mathrm{Ti}$, $\mathrm{S}/\mathrm{Ti}$ and $\mathrm{Si}/\mathrm{Ti}$. The third RPG consists of Mn and $\mathrm{Fe}/\mathrm{Mn}$.

RPG 1 has periodicities of 2, 1 kyr, ca. 700 years and multiple modes between 600 and 340 years (Fig. 4). All proxies within this group show very similar results. RPG 2 shows dominant modes of periodicity of 3.5, 1.5 kyr, few around 2.2 kyr and between 600 and 340 years. $\mathrm{Fe}/\mathrm{Ti}$ and $\mathrm{Mn}/\mathrm{Ti}$ have additional periodicities around 170 years. RPG 3 has a more diverse pattern. Mn has periodicities of 1.5 kyr and 850 years and multiple modes between 600 and 340 years. $\mathrm{Fe}/\mathrm{Mn}$ has periods of ca. 1.5 and 1 kyr and several modes in the multi-centennial range. The non-XRF proxies show fewer significant periodicities. Bi- to multi-centennial cycles are present in the magnetic susceptibility, Cedrus pollen and δ¹³C values. Millennial or multi-millennial periods are only present in the local endmember values, and δ¹³C and δ¹⁸O values. The TOC (total organic carbon) shows no significant cycle.

Figure 4 Synoptic graph of the positive peaks of the significant Redfit spectra of each XRF and non-XRF proxy. The x-axis is log₁₀ scaled, reflects the period and the red dots show the position of the peaks. The configuration and grouping of the XRF proxies follow the results from the correlation analysis (Redfit Proxy Groups). The blue dotted lines (with labels on top of the graph) show dominant modes of periodicity of the presented proxies.

4.2 Proxy maxima and minima during the Holocene

To facilitate the comparison of long-term proxy behaviour in the time domain, the smoothed time series were reduced to their dominant local minima and maxima (Fig. 5). This simplified representation highlights recurring phases of enhanced or reduced proxy expression and allows comparison of temporal coherency between the different Redfit Proxy Groups (RPGs) and the independent non-XRF proxies. Although Fig. 5 does not represent the complete stratigraphic variability of the records, it provides a condensed overview of the principal long-term trend changes and their temporal alignment throughout the Holocene. The first group of proxies (RPG 1, Fig. 5) shows dominant (absolute and secondary maxima) peaks at around 1.8 and 4.5 ka BP. The absolute and secondary minima peaks are distributed more variably through time but mostly during the Early and Middle Holocene until 6 ka BP. Within the first group the peaks of Si are distributed in opposite directions. The second group (RPG 2, Fig. 5) consists of two sub-groups with counter cyclical behaviour (Fe, K and Ti vs S, $\mathrm{Fe}/\mathrm{Ti}$, $\mathrm{Mn}/\mathrm{Ti}$, $\mathrm{S}/\mathrm{Ti}$ and $\mathrm{Si}/\mathrm{Ti}$). Major minima and maxima of both sub-groups are located around 3.8, 4.5 and 6.6 ka BP. The third group (RPG 3, Fig. 5) has a more diverse peak distribution pattern with only two proxies (Mn and $\mathrm{Fe}/\mathrm{Mn}$).

Figure 5 Temporal distribution of major minima and maxima derived from the smoothed proxy time series throughout the Holocene. Red diamonds indicate local maxima and blue diamonds local minima of the individual proxy records after 500-year kernel smoothing. Encircled symbols mark the dominant absolute and secondary extrema of each proxy series. The arrangement of the proxies follows the Redfit Proxy Groups (RPGs) identified from the correlation analysis of the significant spectral power. Grey-labelled proxies represent independent non-XRF datasets. The figure highlights similarities and contrasts in long-term temporal behaviour between the different proxy groups and serves as a simplified overview of major hydroclimatic, erosional and lake-internal trend changes.

4.3 Non-stationary behaviour of millennial and centennial cycles

Regarding the wavelet power spectra, Ca and $\mathrm{K}/\mathrm{Ti}$ show similarities concerning millennial-scale periodicities; significant 2 and 1 kyr cycles for the predominant parts of the Holocene (Fig. 6a and b). In contrast, S and Ti show periodicities of ca. 3.5 kyr and at least partly ca. 1.5 kyr (Ti also ca. 1 kyr) patterns (Fig. 6c and d). All wavelet power spectra of the XRF proxies can be found in the Supplement (Fig. S4–S6).

Figure 6 Exemplary wavelet power spectra of (a) Ca, (b) $\mathrm{K}/\mathrm{Ti}$, (c) S and (d) Ti. Colours show the specific power values. White lines enclose the significant areas and black lines local maxima. Grey facet shows the cone of influence.

To assess the temporal persistence of the periodicities identified by Redfit, wavelet power spectra were analysed for all proxies. There, RPG 1 show almost continuous 2 and 1 kyr cycles throughout the Holocene (Fig. 7). Also, the 600 to 340 years cycles are well represented. RPG 2 shows a clear and consistent 3.5 kyr cycle pattern. In contrast to RPG 1, the 1.5 kyr cycles are more dominant and much less multi-centennial cycles are apparent. RPG 3 shows a more diverse pattern but also clear multi-centennial cycles between 600 and 340 years; mainly in the Early Holocene.

Figure 7 Synoptic graph of extracted significant wavelet power spectra of each proxy. The configuration and grouping of the proxies follow the results from the correlation analysis (Redfit Proxy Groups). The columns of the graph represent the main modes from the Redfit synoptic results (blue, dashed lines in Fig. 4). Each box reflects the whole Holocene time series. Grey bars show the cone of influence, where no significant wavelet periods were calculated due to edge effects. The red-brown colour shows the appearance of the specific cycle (see columns) throughout the Holocene, which occur within the 0.05 significance level and have a power value higher than the 80 % quantile.

5 Discussion

5.1 Numerical robustness and sensitivity of the frequency analyses

Hatvani et al. (2022) developed an online tool to determine the robustness of the spectral information. As the chronological uncertainty of our cs-XRF record, based on a Bayesian age-depth model, ranges from ca. 500 years in the lower half up to ca. 200 years in the upper parts (Fletcher et al., 2017), robust minimum periodicities that can be extracted range from ca. 300 to ca. 120 years. According to that, the 170 years periodicity is less robust and has to be discussed cautiously.

The correlation-based grouping of the proxies should be regarded as a frequency-domain classification of coherent spectral behaviour rather than as a stand-alone process attribution (Mudelsee, 2014; Ólafsdóttir et al., 2016). The applied conservative AR1 red-noise significance threshold emphasizes persistent and statistically robust periodicities, whereas weaker or transient cycles may be underrepresented (Schulz and Mudelsee, 2002). Importantly, the Redfit-based grouping itself is independent of the interpolation procedure, because Redfit operates on the unevenly spaced original datasets. Interpolation was only required for the wavelet analyses used to assess the temporal persistence of the periodicities (Debret et al., 2007; Torrence and Compo, 1998). To further evaluate the robustness of the identified spectral structures, we additionally conducted PCA analyses on both the original proxy dataset (PCA_time) and the extracted spectral information (PCA_freq) (Figs. S7–S9). Although PCA was not used for the primary grouping procedure, the dominant periodicity patterns of the first principal components are comparable to those derived from the original proxy records. This agreement suggests that the identified millennial-scale frequency structures are reproducible across different statistical approaches and are unlikely to result solely from the chosen grouping workflow.

5.2 Element based (XRF) proxy signals

All XRF elements reflect processes during erosion, sedimentation or post-sedimentary alterations, which might be coupled with local or regional triggers that refer to hydroclimatic changes (Bittner et al., 2020; Buggle et al., 2011; Kylander et al., 2011). However, only very rarely are all available elements from cs-XRF shown in palaeoclimate studies (Foerster et al., 2012; Koinig et al., 2003; Schmidt et al., 2019). Often, only selected XRF proxies are shown and used for interpretation (Cohen, 2003; Olsen et al., 2013). In general, our XRF proxies comprise conservative elements that are geochemically fairly stable under different environmental conditions (Boës et al., 2011; Guo et al., 2019) and non-conservative elements that are less stable and affected by biological or lake internal redox mechanisms (Bruland et al., 2014). However, the discussed XRF proxies are highly variable according to the general climatological, geomorphological and geochemical setting of their location (Davies et al., 2015; Foerster et al., 2018).

5.3 Coupling of XRF and non-XRF proxies

Similar spectral behaviour between XRF proxies does not necessarily imply identical environmental forcing, but rather indicates coherent temporal variability that requires environmental interpretation in the context of independent proxy evidence. The identified RPGs therefore represent statistically coherent groups of proxies in the frequency domain. Their environmental interpretation is therefore based on comparison with independent palaeoenvironmental proxies (non-XRF), rather than on spectral similarity alone. The integration of Sidi Ali XRF and Sidi Ali non-XRF proxies is complex in both the time domain and the frequency domain. In the time domain, some similar patterns between Sidi Ali XRF and Sidi Ali non-XRF proxies can be found; e.g. the magnetic susceptibility shows similarities with the peak distribution of Fe and Mn (Fig. 5). In the frequency domain, the temporal resolution of the data generally plays a crucial role in the significance of periodicities (Trauth, 2021). Checking the significance with a high-test level of the red noise of both Redfit and wavelet indicates that the results are statistically robust (Mudelsee, 2014). The number of significant periodicity peaks differs strongly between XRF and non-XRF proxies (Fig. 4). This is probably related to different resolutions (Yu et al., 2016). Hence, a comparison of the extracted periodicities from XRF and non-XRF proxies is also done visually and with less weight on statistical significance.

5.4 Process-based interpretation of the Redfit Proxy Groups (RPGs)

5.4.1 RPG 1 shows spectral characteristics consistent with hydroclimatic variability

RPG 1 consists mainly of proxies related to soluble elements which are known to be sensitive to authigenic precipitation. Ca, Sr and related ratios have frequently been interpreted as indicators for a changing precipitation $/$ evaporation ratio (Davies et al., 2015; Foerster et al., 2012; Kylander et al., 2011). Although Si has the same periodicity, the peaks run in exactly opposite directions (Fig. 5). It is a well-known phenomenon that in alkaline-rich lakes, Ca and Si abundances are opposed (Czymzik et al., 2016). Typically, this is interpreted as the dominance of internal or external processes, since Si is either formed by diatoms or is delivered as detrital input (Regattieri et al., 2016). Zielhofer et al. (2017a) argued for the Sidi Ali record that changes in the water balance might affect the carbonate content, due to an increased relative enrichment of carbonates during low lake levels. The periodic signature of RPG 1 is dominated by 2 and 1 kyr cycles. The latter is in accordance with the 1 kyr periodicity of the hydroclimatic δ¹⁸O signal from the same core (Zielhofer et al., 2019b) (Fig. 4). Increased δ¹⁸O values are in accordance with low precipitation $/$ evaporation ratios. Specifically, multiple millennial-scale peaks in δ¹⁸O during the Early and Middle Holocene (e.g., 5.0, 6.2, 8.2, 9.3 and 10.2 ka BP) reveal short-term decreases in the precipitation $/$ evaporation ratio (Zielhofer et al., 2019b) which can be linked to RPG 1 (Fig. 7) and North Atlantic cooling (Bond events) (Bond et al., 2001). This may indicate a Holocene teleconnection between the North Atlantic climate and Western Mediterranean winter rain. In addition, δ¹³C values (Fig. 5) are interpreted as a groundwater inflow-related proxy (Zielhofer et al., 2017a). This subterranean process is more dominant when less winter rain feeds the lake. We observe a strong ca. 2 kyr periodicity during the Early and Middle Holocene (until 4 ka BP) and link δ¹³C to RPG 1 (Fig. 8f).

The Cedrus pollen show for the Early to Middle Holocene a 1 kyr periodicity (Fig. 8a). After ca. 4 ka BP, the 1 kyr pattern disappears and an ascending (around 1.5 kyr) periodicity becomes visible. In addition, the Early Holocene 2 kyr pattern prolongs until 6 ka BP. Therefore, we associate the Cedrus pollen abundance to RPG 1, which reveal significant 2 and 1 kyr cycles. As discussed in Campbell et al. (2017), Cedrus atlantica is sensitive to prolonged summer drought. For the Early Holocene, Zielhofer et al. (2019b) discussed that Cedrus responded very sensitively to the orbitally-forced summer insolation maximum. Therefore, the absence of Cedar could act as summer drought proxy in the Early Holocene. Furthermore, there is a bi-millennial-scale coincidence of short-term Cedrus peaks at 8.2 and 10.2 ka BP and decreases in Early Holocene subtropical summer sea surface temperatures (deMenocal et al., 2000) which point to a subtropical atmospheric impact at Sidi Ali during the summer season. The TOC shows a clear ca. 2 kyr periodicity and some solitary 1 kyr patterns (Fig. 8c). This fits to RPG 1. This is also supported by a high negative correlation (−0.86) between carbonates (Ca) and TOC (Zielhofer et al., 2017a). They argue that Ca and TOC are connected via aerobic/anaerobic lake-bottom conditions (organic matter preservation) due to lake-level changes.

Also belonging to RPG 1 is the $\mathrm{K}/\mathrm{Ti}$ ratio, which can be used as a proxy for Saharan dust inputs. Zielhofer et al. (2017b) were able to couple the $\mathrm{K}/\mathrm{Ti}$ signal from the Sidi Ali core with remote dust from the Saharan direction using provenance analyses. A peak in the $\mathrm{K}/\mathrm{Ti}$ ratio at 4.2 ka BP corresponds with a hyper arid phase in the central Sahara (van der Meeren et al., 2022). The Saharan dust signal is in accordance with periodicities that are known from North Atlantic Ocean-atmosphere pattern and from monsoonal variability (Cruz et al., 2021; Gasse and Van Campo, 1994). Specifically, the Early to Middle Holocene 2 kyr periodicity is known from the Saharan and Sahelian domain (Gasse and Van Campo, 1994; Tjallingii et al., 2008).

5.4.2 RPG 2 – interpreted catchment erosion and terrestrial input signals

This periodicity of RPG 2 is dominated by distinct 3.5 and 1.5 kyr patterns (Figs. 5 and 7). However, the 3.5 kyr cycles are present throughout the Holocene, whereas the 1.5 kyr pattern occurs more explicitly during the Middle and the second half of the Late Holocene (Fig. 7). Although the 3.5 kyr cycle is somewhat less trustworthy due to the short length of the total record, the dominance and grouping in the dataset seems clear.

RPG 2 consists of proxies that are commonly interpreted as indicator of detrital influx (Fe, Ti, K, and related ratios; used and described for Sidi Ali in Zielhofer et al. (2017a) or bioproductivity ($\mathrm{Si}/\mathrm{Ti}$, S, $\mathrm{S}/\mathrm{Ti}$) processes (Davies et al., 2015; Kylander et al., 2011). Additionally, the $\mathrm{Fe}/\mathrm{Ti}$ and $\mathrm{Mn}/\mathrm{Ti}$ ratios show also multi-centennial (340–600 years) periodicity patterns (Figs. 5 and 7). This corresponds partly with patterns of RPG 1 and stronger with RPG 3. The general temporal pattern of the included proxies (Fig. 5) shows that this group can be differentiated in two subgroups (Fe, K, Ti vs. $\mathrm{Fe}/\mathrm{Ti}$, $\mathrm{Mn}/\mathrm{Ti}$, S, $\mathrm{S}/\mathrm{Ti}$, $\mathrm{Si}/\mathrm{Ti}$). The observed opposing behaviour may be a result of varying detrital sedimentation rates versus a biogenetic sedimentation rate. The XRF-based interpretation is supported by the magnetic susceptibility as a local erosion proxy (Zielhofer et al., 2017a). The magnetic susceptibility shows a distinct and persistent 3.5 kyr and a Middle Holocene 1.5 kyr periodicity (Fig. 8b). Furthermore, multi-centennial periodicities are less pronounced. This spectral pattern is similar to that of RPG 2. Zielhofer et al. (2017a) reported a high correlation between Fe and magnetic susceptibility due to terrestrial input. The coarse-grained local endmember derived from grain-size EMMA is interpreted to reflect local catchment erosion and short-distance sediment supply to Lake Sidi Ali (Zielhofer et al., 2017b). However, Redfit and wavelet analyses reveal no clear or persistent expression of the ∼ 1.5 and ∼ 3.5 kyr cycles in the frequency domain of this proxy (Figs. 4 and 8d), indicating that local erosion processes might not only be paced by these millennial-scale periodicities but also by human land use and disturbance in the catchment (Bourchachen et al., 2025; Campbell et al., 2017; Cheddadi et al., 2015).

5.4.3 RPG 3 – redox-sensitivity

RPG 3 consists only of two XRF-derived proxies, Mn and $\mathrm{Fe}/\mathrm{Mn}$. Both are known as redox sensitive (Davies et al., 2015; Davison, 1993). However, Mn is more affected by redox processes compared to Fe (Boyle, 2002). Fe corresponds with RPG 2 and seems to be mainly reflecting detrital input. The effect of redox processes is strongly dominated by lake-level variations and subsequent strengths and duration of lake stratification and circulation (mictic) properties (Davies et al., 2015). The lake level of Sidi Ali is controlled by winter-rain variations and its hydroclimatic consequences (Zielhofer et al., 2019b). Therefore, a possible interpretation for RPG 3 is that it reflects combined influences of processes related to both RPG 1 and RPG 2. The periodicity pattern of RPG 3 supports this explanation, as it shows patterns of both, RPG 1 and 2. The mix of ca. 2, 1.5, 1 kyr and several multi-centennial cycles is clearly in contrast to those of RPG 1 and 2 (Figs. 5 and 7). Fe and Mn are associated with detrital input which is also driven by catchment related processes. The subsequent redox-related alteration rather corresponds with periodicities that are associated with RPG 1. Moreover, multi-centennial periodicities occur also in redox-sensitive proxies, like the magnetic susceptibility record (Figs. 5 and 7).

Figure 8 Wavelet power spectra of (a) Cedrus pollen proportion (Campbell et al., 2017), (b) Magnetic susceptibility scan and (c) TOC (Zielhofer et al., 2017a), (d) Grain Size Endmember Analysis (EMMA) ratio (local provenance) (Zielhofer et al., 2017b), (e) δ¹⁸O values and (f) δ¹³C values of ostracod shells (Zielhofer et al., 2019b). Colours show the specific power values. White lines enclose the significant areas. Grey facet shows the cone of influence.

5.5 Solar cycles and their forcing on hydroclimatic variability in the Western Mediterranean

5.5.1 Hydroclimatic pattern of the bimillennial Hallstatt cycle

The bimillennial cycle is a well-known periodicity in Holocene palaeoclimate studies (Azuara et al., 2020; Goslin et al., 2018). Generally, it is considered to be triggered by solar forcing via the Hallstatt cycle (Steinhilber et al., 2012; Usoskin et al., 2016). In the North Atlantic realm, Stewart et al. (2017) presented a Bromine record from a peat bog in Scotland that shows a distinct 2 kyr cycle pattern (Fig. 9a). This Bromine record measured by cs-XRF reflects the storminess of the adjacent Atlantic Ocean. Stewart et al. (2017) reconstruct past variations in North Atlantic storminess and atmospheric circulation, reflecting changes in the NAO. As such, it represents an important archive of North Atlantic climate variability, which is directly relevant for interpreting hydroclimatic changes at Sidi Ali due to the strong influence of NAO-like westerly variation on moisture transport into the Mediterranean region.

Regarding the Western Mediterranean, the composite speleothem record from the Chaara Cave located in the Middle Atlas covers the entire Holocene. Ait Brahim et al. (2019) presented there a δ¹⁸O value-based Western Mediterranean rainfall index. The wavelet power spectrum shows high power values for a 2 kyr cycle since the Middle Holocene (Fig. 9b). We suggest that the Western Mediterranean precipitation $/$ evaporation ratio might be influenced by a bi-millennial periodicity. However, the same frequency analysis for the δ¹⁸O speleothem record of the nearby Grotte du Piste Cave (Wassenburg et al., 2016) did not reveal any bi-millennial periodicity (Fig. 9c). Hence, the climatic imprint of bi-millennial periodicities in Mediterranean hydroclimatic records remains still under debate. In contrast, a speleothem record located in the SW High Atlas Mountains (Wintimdouine Cave) is discussed as the airmass trajectory which is driven by the West African Monsoon system and consequently by the tropical realm (Sha et al., 2019). Wavelet analyses show a ca. 2 kyr cycle pattern (Fig. 9d) which is in accordance with the Saharan and sub-Saharan realm (Gasse, 2000; Gasse and Van Campo, 1994; van der Meeren et al., 2022).

Figure 9 Wavelet power spectra of (A) Br record (Stewart et al., 2017), (B) Rainfall index (Ait Brahim et al., 2019), (C) δ¹⁸O record from Grotte de Piste cave (Wassenburg et al., 2016) and (D) δ¹⁸O record from Chaara cave (Sha et al., 2019), (E) NAO index (Olsen et al., 2012). Colours show the specific power values. White lines enclose the significant areas. Grey facet shows the cone of influence. Panel F shows a map of the geographic locations of the datasets presented in this figure (EPSG: 4326).

5.5.2 Hydroclimatic forcing of the millennial Eddy cycle

The eddy cycle is discussed as solar forcing mechanism behind this millennial cycle pattern (Debret et al., 2007; Steinhilber et al., 2012) and seems to affect the Atlantic precipitation regime (Zielhofer et al., 2017a). Wavelet analysis of the total solar irradiance (TSI) shows a millennial cycle which is present during the last 9.5 ka BP, although its amplitude varies through time (Steinhilber et al., 2012). The 1 kyr millennial cycle is well-known in many palaeoclimatic records (Debret et al., 2009; Soon et al., 2014; Zhao et al., 2021). During the Holocene, the Scottish Shebster peat record (Stewart et al., 2017) shows a significant 1 kyr cycle of North Atlantic storminess until 4 ka BP (Fig. 9a).

In the Western Mediterranean, the 1 kyr pattern in the rainfall index from the Chaara Cave (Ait Brahim et al., 2019) seems present from the onset of the Holocene until ca. 8 ka BP and again from ca. 4 ka BP (Fig. 9b). The Grotte du Piste speleothem record (Wassenburg et al., 2016) shows only a minor, but not highly significant 1 kyr cycle from ca. 10 to ca. 5 ka BP (Fig. 9c). Regarding the potential impact of sub-Saharan and Saharan air masses, the δ¹⁸O record from Sha et al. (2019) does not show a clear 1 kyr pattern (Fig. 9d). Kuhlmann et al. (2004) argues that there is a distinct border between the tropical monsoonal impact and the North Atlantic impact between 27 and 30° N. They used the 1 kyr cycle pattern for the differentiation between the two impacts. In this context, they described that there is no 1 kyr periodicity in the African tropics in the Holocene.

5.5.3 Hydroclimatic forcing of multi-centennial cycles

The multi-centennial periodicity patterns cluster around 700 and between 600 to 340 years (Fig. 4). These patterns are also related to solar forcing, with unnamed cycles at ca. 340, 500 and 710 years (Steinhilber et al., 2012). Furthermore, Soon et al. (2014) demonstrated that the multi-centennial cycles are a result of Atlantic Ocean dynamic responses of the solar forcing. The multi-centennial periodicities are also reflected in the Moroccan speleothems from the Middle Atlas (Ait Brahim et al., 2019; Wassenburg et al., 2016; Fig. 9b and c), suggesting periodic winter-rain variability. The NAO record from Olsen et al. (2012) which is related to winter rain (over Greenland), also shows periodicities for the 700 and 340 years cycles (Fig. 9e). Surprisingly, the Atlantic storminess record from Stewart et al. (2017) includes almost no visible multi-centennial periodicities (Fig. 9a). As we find the majority of multi-centennial cycles in the hydroclimatic proxy group (RPG 1), we suggest that it is dominated by Atlantic forcing mechanisms which is in accordance with Zielhofer et al. (2019b).

5.5.4 Hydroclimatic forcing of the de Vries (∼ 200 years) cycle

The observed cycle length around 200 years is known as de Vries cycle (Steinhilber et al., 2012). The cycle has been proposed to be modulated in amplitude by the longer Hallstatt cycle, implying that its climatic imprint may vary through time rather than being continuously expressed (Komitov, 2024; Usoskin et al., 2016). Periodicities of similar lengths are found in many proxy-based studies (Ojala et al., 2015; Soon et al., 2014). Two speleothem records from North and South Morocco that cover the last 1000 years, included ca. 200 years periodicities and linked them to NAO-like variations (Ait Brahim et al., 2017, 2018). They argue that the de Vries cycle varies between dry and humid conditions. This would also affect the lake level of Sidi Ali, thus altering the redox processes at the lake floor and within the lake sediments. Accordingly, these lake-level–driven changes are reflected in the detrital proxies of RPG 2 and in the redox-sensitive proxies of RPG 3. The latter additionally exhibit a ∼ 170 years periodicity, which may tentatively be linked to the de Vries cycle, although the statistical robustness of this specific periodicity is limited and its interpretation therefore remains cautious (see Sect. 5.1).

5.5.5 Indications for a 1.5 kyr cycle?

The 1.5 kyr cycle in the Holocene, on the other hand, is subject of many studies (Debret et al., 2007; Fletcher et al., 2013; Soon et al., 2014). Nevertheless, there is an open debate about the frequency framing (distinct separation between 1.8 to 1.5 kyr often not achievable, due to dating uncertainties) and the fundamental forcing of this pattern (Dima and Lohmann, 2009; Soon et al., 2014). The derived NAO index of Olsen et al. (2012) shows, even though not highly significant, a ca. 1.5 kyr pattern (Fig. 9e). Due to the time span of the NAO record, the pattern is only visible from ca. 5 to ca. 2 ka BP. Other studies also suggest a 1.5 kyr periodicity that is stable since the Middle Holocene (Azuara et al., 2020; Smith et al., 2016). The marine record of Bond et al. (1997) shows a similar pattern. There is a transition from 1 kyr periodicity to 1.5 kyr periodicity in the Middle Holocene. However, Witt and Schumann (2005) detected a significant ca. 1.5 kyr cycle in Greenland ice core data, which is valid for the entire Holocene. Furthermore, Soon et al. (2014) identified a persistent 1.5 kyr cycle in a nitrate record in an Antarctic ice core. One can argue that the persistent Holocene-wide cycling is due to solar forcing, whereas the Middle Holocene shift from 1 to 1.5 kyr seems to be the dynamic response of the Atlantic Ocean circulations. As we see the onset of the 1.5 kyr cycle in RPG 2 in the Middle Holocene (Fig. 7), we assume an external forcing that is specifically affecting catchment processes and terrestrial supply. The clear mechanisms remain still unknown. The speleothem records in the Middle Atlas do not show a distinct 1.5 kyr pattern (Fig. 9b and 9c). The speleothem record of Sha et al. (2019) from the SW High Atlas shows a ca. 1.5 kyr pattern between ca. 7.5 and 2 ka BP, although statistically not highly significant (Fig. 8d). The Saharan-related processes, thus, seem to follow this pattern. Similar to the 1.5 kyr pattern, the 3.5 kyr periodicity is difficult to attribute to a specific forcing mechanism. It seems that catchment-related processes and local terrestrial supply differ from supra- regional North Atlantic or sub-tropical hydroclimatic forcing.

6 Conclusions

This study demonstrates that the combination of calibrated high-resolution core-scanning XRF data and multi-frequency time-series analyses (Redfit and wavelet) provides a robust framework to investigate Holocene hydroclimatic variability at the North African desert margin. The calibration approach allows quantitative interpretation of mm-scale cs-XRF data while reducing compositional effects, and the combined stationary and non-stationary spectral analyses enable both the detection and temporal evaluation of significant periodicities under conservative red-noise assumptions.

A central outcome is the identification of three Redfit Proxy Groups (RPGs) derived from correlations of significant Redfit spectra. These groups differ not only in proxy composition but also in their characteristic frequency patterns, indicating that distinct environmental processes dominated lake-system behaviour during the Holocene. RPG 1 is characterised by persistent ∼ 2 and ∼ 1 kyr periodicities and additional multi-centennial modes. These proxies primarily reflect changes in the precipitation/evaporation balance, lake level and authigenic carbonate formation. Their frequency coherence with independent proxies from the same core (δ¹⁸O, δ¹³C, TOC and Cedrus pollen) indicates a dominant control by hydroclimatic variability, likely linked to North Atlantic–Western Mediterranean winter-rain dynamics. RPG 2 is dominated by ∼ 3.5 and ∼ 1.5 kyr cycles and shows fewer multi-centennial periodicities. This group is interpreted to reflect catchment erosion, terrestrial input and productivity-related processes, supported by its close correspondence with the magnetic susceptibility record. The emergence of the ∼ 1.5 kyr mode mainly after the Middle Holocene suggests a temporal shift in the sensitivity of catchment processes to external forcing. RPG 3 displays mixed millennial and multi-centennial periodicities and is interpreted to record redox-sensitive lake-internal processes, integrating both hydrological and detrital influences.

Overall, the results indicate two main Holocene periodicity regimes at Lake Sidi Ali: a hydroclimatic regime dominated by ∼ 2 and ∼ 1 kyr cycles, and a catchment-related regime characterised by ∼ 3.5 and ∼ 1.5 kyr variability. These findings highlight that hydroclimate, productivity and erosion respond to different, only partly interacting forcing mechanisms at millennial to multi-centennial scales in Mediterranean Northwest Africa.