Evaluation of Karst Spring Discharge Response Using Time-Scale-Based Methods for a Mediterranean Basin of Northern Algeria

Abstract

Understanding of behavior, variability, and links between hydrological series is a key element for successful long-term water resources planning and management. In this study, various time-scale-based methods such as correlation and spectral analysis (CSA), cross wavelet (XWT), and wavelet coherence transform (WCT) were applied to assess the response of daily rainfall and karst spring discharge for the Sebaou River basin, which is located on Mediterranean basin in northern Algeria. The CSA revealed that the hydrogeological systems under study are characterized by various memory effect (small, poor, reduced, and extensive) with regularization times ranging from 5 to 50 day. XWT between rainfall and discharge time series indicates few marked disruptions in the spectra between the 1980s and 1990s corresponding to the dry period. The annual process is visible, dominant, and more amplified compared to the multi-annual fluctuations that characterize the 1-3- and 3–6-year modes, which explained the multi-annual regulation. The nonlinear relationship of the short-term components seems to be linked to the periods of storage (infiltration). Compared to the WCT components of 2–5, 26, and 52 weeks, there is a strong coherence for 102 weeks, which explains the long-term component, indicating a quasi-linearity of the rainfall-runoff relationship. According to the obtained results, the construction of more water resources structures is recommended to increase the water storage and improve the water supply due to the richness of the hydrographic network. On the other hand, the impacts of human activities on streamflow due to the looting of rocks and sands in the Sebaou River valleys have reached alarmingly high levels that require urgent intervention for the protection of water and ecological resources and their better rational use.

1. Introduction

Karst is a land with unique shapes and drainage. A worn increased porous zone called epikarst is located on the surface, but in the subterranean there are networks of fractures, splits, and conduits of all sizes and shapes, which comprise an extremely and complex active hydrologoelogical system [1,2]. The reverse flow occurs in channels while uniform flows propagate in fractures and fissures [3]. It can also register rapid discharge and base discharge components for karst springs. The tubing network regulates the fast flow component, whereas the fractures and cracks network manage the base flow component of the release [4,5]. Karst lands are very susceptible to anthropogenic factors, and are vulnerable to contamination from many sources; therefore, it is frequently important to measure the amount of rapid flow [6,7,8]. The result of the karst’s internal structure is rapid oscillations of groundwater, often over large amplitudes, causing rapid changes in the velocity and groundwater direction. This implies that floods from one karst spring to another occur depending on the groundwater level, therefore, in many karst locations across the world, adjacent karst springs with overlapping catchments are common [9,10]. The word “overlapping” shows that portions of the catchment are intermittent or permanent in the presence of multiple karst springs. The groundwater exchange system is typically complicated and time-consuming in adjacent karsts [1,11,12,13].

Different runoffs in recent years have received extensive attention in hydrological and climate research for the quantification of consequences of human interactions and climate change [14,15,16,17,18,19,20]. Some human impacts are demographic pressure, traditional management, and random building of storage structures such as dams. In addition, illegal logging by sand smugglers from the main bed of Wadis River and overexploitation of aquifers by industry and agriculture with very advanced pumping are rapidly decreasing the aquifer storage, and this appears to be the reasons affecting the decline of flows and regulatory reserves [21]. The impact assessments of these factors are essential, and a good understanding of the external environment and hydrogeological structure is required. Understanding this system in depth is one of the targets of this current investigation, as based on the results, it may or may not allow the use and manage of groundwater from the aquifer. Therefore, it is very necessary to point out the abrupt moment (time) of change in these factors to assess the individual impact of these the causes on the runoff. It is worth highlighting that changes in river runoff in different rivers throughout the globe are the central focus of a huge literature that seeks to determine the temporal features of these rapid changes before assessing the effect of climate variability and interactions with humans [22].

Numerous processes with spatiotemporal variation happen across the hydrological cycle and are influenced by hydrological time series. By means of correlation functions, the impacts of these processes may be identified. The high dependence or interconnection between time series can significantly increase the correlation values, therefore misunderstandings generally arise about the effects contained in correlation functions, where some partial correlation functions can overcome those inconveniences by eliminating the linear influence. There are relatively new techniques for the arithmetical reliability test in the hydrological arithmetical has been used to analyze the impacts of hydrometeorological fluctuations on the karst streamflow, e.g., Refs. [23,24].

Correlation spectral analysis (CSA) was developed by Jenkins and Watts (1968), Box and Jenkins (1976), and Mangin (1984) [25,26,27]. It was used for the identification of basic characteristics of karst underground, including research on hydrodynamics of groundwater levels [28,29], transport properties and turbidity dynamics [30], surface and underground flow interactions, cave duct systems, karst underground network connections, storm behavior effects, and response times [31,32,33,34,35,36]. The wavelet method has a different approach because it is directly applied to intervals and can discover latent and hidden aspects in a time series [37]. According to Torrence and Compo (1998) [38], the emergence of wavelet approach represents a qualitative leap in the modern interpretation of results in the time–scale domain. Traditionally, the Short Time Fourier Transform (STFT) is the most used technique to investigate a series in the frequency domain, where it is based only on the assumption of stationarity but, unfortunately, this does not allow a good temporal resolution to the detriment of frequency resolution and vice versa due to the Heisenberg uncertainty principle [39]. One of the main reasons for the success of wavelet transform-based methods is their high spatial precision at small scale. The use of wavelet transform increases significantly year after year in engineering and earth science fields, such as hydroclimatic series analyses [40], large-scale fluctuations [41,42], and analyses of hydrochemical and isotopic characteristics [43]. For example, the relationship between groundwater storage with El Niño, NAO, and AMO variability was investigated by Resende et al. (2019) [44] across Africa. The results indicate that these teleconnection patterns can have a significant impact on groundwater storage and recharge process. In a central-eastern European agricultural region, shallow groundwater (SGW) level fluctuations were assessed according to the climate-driven periodicity. Garamhegyi et al. (2018) [45] revealed a clear relationship between the drought periodicity and the absence of annual cyclic comportment in the SGW level. Meng et al. (2021) [46] analyzed the rainfall effects on karst spring response characteristics in northern China. The results documented that the flow response time to rainfall ranged between 7–16 months and that the karst system in southern China has a low regulating effect compared to the southern part. In Algeria, human effects have also reached record levels in all natural resources, especially in the study area, which requires interventions by the competent authorities to take appropriate deterrent measures and carry out accurate studies that provide realistic solutions for sustainable development and better water and natural resource management. In the absence of studies that discuss streamflows and the nature of karst spring water in Algeria using time-frequency-based methods, this study may provide water resource managers with some pictures for understanding the hydrogeological system of Sebaou River basin. Thus, in this research, CSA, cross wavelet (XWT), and wavelet coherence transform (WCT) were applied on daily hydroclimatic time series of Sebaou River basin (northern Algeria) to:

• Understand the internal structure of aquifers, storage capacity and obtain information about periodic characteristics.
• Classify the karst systems of the basin based on memory effect and regularization time.
• Study the random process of daily hydroclimatic time series using statistical fractals.
• Assess the causal links and linearity between rainfall and runoff for each sub-basin of the study area.
• Extract the significant coherence and covariance through isolated components developed between rainfall and runoff time series at a time-scale domain and identify dry and wet periods as well as anthropological impacts on the daily streamflow of the Sebaou River basin.

2. Study Area and Database

The Sebaou River basin is one of the Mediterranean basins located in the central north of Algeria, located between longitudes E 03°30′ and 04°30′ and latitudes N 36°30′ and 37°00′ with an area of 2500 km². Its elevations reach more than 2030 m above sea level (Figure 1b).

The top of the Sebaou river basin consists of the limestone chain that borders the south of the Kabyle stands (Bouira, Tizi Ouzou, and Bejaia) [47], with flysch in the north and east. In the west, it is the Miocene that lines the entire valley. Argillaceous shales and gray and schistous marls with layers of sandstone occupy 70.5% of the middle Sebaou. On the other hand, the lower Sebaou consists of 43.15% of marl formations, sandstone and conglomerates (Figure 1c). These marls, which constitute 69% of the total area of the basin, crack by dehydration during the dry season and deteriorate by several centimeters, forming a carpet of gravel mixed with clay-limestone dust [48].

According to Djemai (2008) [51], the geology of the study region is marked by the peri-Mediterranean alpine orogeny. It appears in several lithological facies, ranging from the oldest Cambro-Ordovician, which forms the plinth, to the most recent, the Quaternary composite of alluvial deposits. Understanding the overall structure of the region depends on the origin of the flyschs. The soft and brittle tectonics responsible for the (N-70) and (N-140) structures and the current morphology are the products of finite Miocene tectonics. The high permeability area is linked to the limestone formations of Djurdjura and detrital of the alluvial aquifer (Figure 1d) [51].

The limestones of Djurdjura and the alluvial deposits of the Sebaou River and its tributaries are the well-known reservoirs of exploited water (Figure 1). The age of the hydrographic network is after the finished-Miocene phase that led to the formation of longitudinal folds (N-70).

The construction of the dam in the municipality of Beni Aïssi in the Aïssi River started in 1994 and was completed in 2001. It is an earth dam 94 m high and 515 m long. The rainfall distribution over the Sebaou River basin (Figure 2a) indicates that the observed interannual rainfall averages (1972–2010) ranged between 700 and 1200 mm with coefficient of variation between 0.20 and 0.33 (Figure 2b).

The databases used in the analysis were obtained from the Algerian National Agency of Water Resources (ANRH) and the details about the stations are presented in

Table 1. Rainfall and runoff stations, and sub-basins studied in the cross analysis.

Rainfall Station/ANRH ID	Runoff Station/ANRH ID	River Names	Sub-Basin/ANRH Code	Area of the Sub-Basin in ha	Study Periods
Ait Aicha 02-15-09	Boubhir 02-15-13	Boubhir	Acif N’boubhir 02-15	53,830	1987–2011
Fréha 02-16-03	Fréha 02-16-05	Dis	Sebaou Rabta 02-16	43,330	1987–2002
Benni Yenni 02-17-12	RN 30 02-17-15	Aïssi	Aïssi 02-17	47,020	1985–2010
Tizi Ouzou 02-18-10	Belloua 02-18-03	Sebaou	Sebaou Sebt 02-18	30,630	1987–1998
DEM 02-19-02	RN 25 02-19-09	Boughdoura	Boughdoura 02-19	53,450	1973–1993
Baghlia 02-20-02	Baghlia 02-20-01	Sebaou	Sebaou Maritime 02-20	22,890	1964–1998

. In addition, long-term rainfall time series were used in the univariate analysis. Figure 3 represent the daily time series of rainfall and streamflow in the Aïssi River sub-basin (1985–2010), whose code is 02-17, as shown in Figure 2 and Table 1.

3. Materials and Methods

3.1. Correlation and Spectral Analysis

Correlation and spectral analysis (CSA), expressing the variation in input-output covariance for different frequencies, has been successfully applied in regional studies of large karst aquifers in the world [52,53].

3.1.1. Simple Analysis

Correlation analysis corresponds to the processing of a single time series, for example, the input or output. It assesses the dependence of successive measures for accumulative time intervals [27,54]. The values r_k (k = 0, 1, 2, 3, …, m) are the autocorrelation coefficients obtained. According to Box and Jenkins (1976) [26], the choice of truncation (m) is not based on theoretical concepts, but it can be set as m = N/2, m = N/3 or m = 2N/3. The following expression gives r_k:

$$r_k=\frac{C_k(k)}{C_k(0)}$$

(1)

$$C_k(0)=\frac{1}{N}{\displaystyle \sum _{t=1}^N(}{x}_t-\overline{x}{)}^2{C}_k(k)=\frac{1}{N}{\displaystyle \sum _{t=1}^{N-k}(}{x}_t-\overline{x})(x_{t+k}-\overline{x})$$

(2)

where $C_x\left(k\right)$ is the autocovariance and k is the time lag or step.

According to Padilla and Pulido-Bosch (1995) [55], the power spectral density function ${\mathsf{\Gamma }}_x\left(f\right)$ is an unbiased approach of the Fourier transform of the autocorrelation function, where the cyclic phenomena seem like peaks in the graph of ${\mathsf{\Gamma }}_x\left(f\right)$ [27,54], which is defined as:

$${\mathsf{\Gamma }}_x(f)=2\left[1+2{\displaystyle \sum _{x=1}^m{D}_k{r}_k\cos (2\pi fk)}\right]$$

(3)

where f is the frequency and $D_k$ is a weighting function chosen in such a way that the estimated value of the spectrum ${\mathsf{\Gamma }}_x\left(f\right)$ is not biased [55]. The Tukey filters were used in this study.

3.1.2. Cross-Analysis

Cross-analysis simultaneously treats two time series, one considered as the input time series of the system, and the other as the output time series. Therefore, it no longer seeks to directly describe the details of the system, but rather a black box with inputs and outputs. It is, therefore, by studying its operation at the source that the structure of the system and its degree of karstification are deduced [54,56].

• Cross correlograms

Cross-analysis gives the causal relationship between two time series X_t and Y_t of N observations [23,25,57]. When $r_{+k}\ne r_{-k}$, this explains that the cross-correlation function is not symmetric:

$$r_{+k}=r_{xy}(k)=\frac{C_{xy}(k)}{\sqrt{C_x^2(0)}{C}_y^2(0)}{r}_{-k}=r_{yx}(k)=\frac{C_{yx}(k)}{\sqrt{C_x^2(0)}{C}_y^2(0)}$$

(4)

• Cross spectrum

The cross spectrum corresponds to the decomposition of the covariance between inputs and outputs in the frequency domain. A complex number explains the spectral density function that represents the asymmetry of the intercorrelation function given by the following expression:

$${\mathsf{\Gamma }}_{xy}=\left|{\alpha }_{xy}(f)\right|\exp \left[-i\mathsf{\Phi }(f)\right]$$

(5)

From the point of view of its application to the study of hydroclimatic research, the cross-amplitude function ${\alpha }_{xy}\left(f\right)$ expresses the variation of the hydrological input-output covariance for different frequencies.

The phase function ${\mathsf{\Phi }}_{xy}\left(f\right)$ expresses the output delay in relation to the input for each frequency with a variation range of 2π, between -π and + π (Equation (6)):

$${\alpha }_{xy}(f)=\sqrt{{\psi }_{xy}^2(f)+{\mathsf{\Lambda }}_{xy}^2(f)}{\mathsf{\Phi }}_{xy}(f)=\arctan \frac{{\mathsf{\Lambda }}_{xy}(f)}{{\psi }_{xy}(f)}$$

(6)

where ${\mathsf{\Gamma }}_{xy}\left(f\right)$ is the cross spectral density function between ${\mathsf{\Gamma }}_x\left(f\right)$ and ${\mathsf{\Gamma }}_y\left(f\right)$ of the input and output, respectively, i denotes $\sqrt{-1}$, ${\alpha }_{xy}\left(f\right)$ is the amplitude, ${\mathsf{\Phi }}_{xy}\left(f\right)$ is the phase function at the frequency f, ${\psi }_{xy}\left(f\right)$ is the co-spectrum, and ${\mathsf{\Lambda }}_{xy}\left(f\right)$ is the quadrature spectrum.

The coherence function $k_{xy}\left(f\right)$ exhibits the square of the correlation between the cyclical components of the input-output at the corresponding frequency. It provides information about the linearity of the system and is assimilated to an intercorrelation between the events.

The gain function $G_{xy}\left(f\right)$ expresses the variations of the regression coefficient (input variance/output variance), accordingly depending on the frequencies. Therefore, it provides an estimate for the augmentation or reduction of the input signal relative to the output signal.

$$k_{xy}(f)=\frac{{\alpha }_{xy}(f)}{\sqrt{{\mathsf{\Gamma }}_x(f){\mathsf{\Gamma }}_y(f)}}{G}_{xy}(f)=\frac{{\alpha }_{xy}(f)}{\sqrt{{\mathsf{\Gamma }}_x(f)}}$$

(7)

3.2. Cross Wavelet Transform

In this study, the XWT between rainfall (X_n) and runoff (Y_n) is defined by the cross-wavelet power spectrum W^xy = W^x W^y*, where (*) explains the conjugate complex W^xy, and is given as follows [58,59]:

$$D\left(\frac{\left|W_n^X(s)W_n^{X*}(s)\right|}{{\delta }_X{\delta }_Y}<p\right)=\frac{Z_{\nu }(p)}{\nu }\sqrt{P_k^X{P}_k^Y}$$

(8)

where:

$$P_k=\frac{1-{\alpha }^2{}^{}}{{\left|1-\alpha e^{-2i\pi k}\right|}^2}$$

(9)

$P_k$ is the Fourier spectrum with autocorrelation of lag-1.$P_k^x$ and $P_k^y$ are calculated for $X_n$ and $Y_n$ of the variance ${\mathsf{\delta }}_x$ and ${\mathsf{\delta }}_y$, respectively. Z_v (P) is the significance level for the probability (P) density function. For XWT, the user must be aware that a coefficient of XWT can be high because the wavelet power spectrum of the two signals is high [60].

3.3. Wavelet Coherence Transform

According to Torrence and Webster (1998) and Grinsted et al. (2004) [58,59], WTC function is given by the following equation:

$$R_n^2(s)=\frac{{\left|S(s^{-1}{W}_n^{XY}(s))\right|}^2}{S\left({\left|s^{-1}(W_n^X(s))\right|}^2\right).{}^{}S\left({\left|s^{-1}(W_n^Y(s))\right|}^2\right)}$$

(10)

where S is the smoothing operator and resemble the mother-wavelet. According to Torrence and Webster (1998) [58], the most compatible parameter $S$ for Morlet wavelet is given by the following equation:

$$S(W)=S_{scale}\left(S_{time}(W_n(s))\right)$$

(11)

where $S_{time}$ and $S_{scale}$ are smoothing operators in time and scale, respectively. Further information and details on the XWT and WTC theories can be found in Refs. [38,58,59].

4. Results and Discussion

4.1. Overview of the Rainfall Trends

Figure 4a represents some annual rainfall time series of the study area that explain that the rainfall in the Sebaou River basin can reach 1800 mm/year, making it one of the most watered basins in the Mediterranean. The intersection between the progressive and retrograde curves of the sequential Mann-Kendall test (SQMK) (Figure 4b) shows that the first significant negative trend starts in the years 1974 and 1984 until the early 2000s. However, the second non-significant positive trend begins in early 2001. According to Figure 4b, the SQMK test indicates a significant value of U(t) for the 95% confidence interval, approximately equal to −2.5; this tendency was evidenced in all the stations during the year 2001. These climatic trends have been observed in the Mediterranean basin and in northern Algeria, which agree with the presence of negative rainfall trends between the 1980s and 1990s [61,62,63,64,65,66,67].

4.2. Univariate Correlation and Spectrum Analysis

The short-term analysis of rainfall (observation window from 1 to 125 days) is presented in a correlogram, which decreases very rapidly, assuming the value 0.2 at the 3rd day and then becoming zero at 35 days (Figure 5). The rainfall does not have a clear structure and has a quasi-random character in the short-term analysis (<6 months).

The global short-term analysis of daily flow rates shows that the correlograms (Figure 5) decrease rapidly at Boubhir and Freha, a little fast for RN30 before 1999, and slow for Baghlia and RN30 after 1998. This explains why the hydrogeological system of Sebaou River is spatially heterogeneous and complex, or why each sub-basin (karstic system) is characterized by different memory effects each other, according to the time dependence and independence of the events affecting the flows.

Table 2. Classification of six karstic systems of Sebaou River basin and some representative Mediterranean karstic system based on CSA of daily flow.

Author	Region	Karstic System	Memory Effect (Day) (rk = 0.1–0.2)	Spectral Band: Cutoff Frequency	Regularization Time (Day)
Mangin (1984) [27]	Pyrenees (France)	Aliou	Poor (5 days)	Very large (0.3)	10–15
		Baget	Small (10–15 days)	large (0.2)	20–30
		Fontestorbes	Large (50–60 days)	narrow (0,1)	50
		Torcal	Extensive (70 days)	Very narrow (0.05)	70
Bouchaou (1995) [68]	The pleated Middle Atlas (Maroc)	Asserdoune	Extensive (70–80 days)	Very narrow (0.04–0.05)	70–80
Larocque et al. (1997) [69]	Western France	Rochfoucauld	—	—	76
Amraoui et al. (2004) [54]	The tabular Middle Atlas (Maroc)	Bittit	Large (37–45 days)	Very Wide	35
		Ribaa	Extensive (70 days)	large (0.14)	57
Chettih and Mesbah (2010) [70]	Saharan Atlas (Algeria)	Seklafa	2.5	0.4	1.5
		Kerakda			3.5
		Rhouiba			4
Bouanani (2004) [71]	basin Tafna (Western Algeria)	Sebdou	Small	—	5
		Mouilah	Large	0.025	21
		Isser	Extensive	0.018	43
This work	Sebaou River (Algeria)	Boughdoura River	Reduced (18 days)	large (0.21)	20–30
		Aïssi River	Extensive (53–84 days)	Very narrow (0.032)	50
		Acif N’boubhir	Poor (9 days)	Large (0.22)	15
		Sebaou Sebt River	Small (16 days)	Large (0.19)	20–30
		Sebaou Rabta River	Poor (3 days)	Very large (0.44)	5
		Sebaou maritime River	Extensive (66 days)	Very narrow (0.067)	60

shows a summary classification of karst systems of the Sebaou River basin and some representative Mediterranean karstic system based on the memory effect and the regularization time.

Since the rainfall signal is a random input function, the log-spectrum rises dominating periodic structures to f = 0.0028 (hydrologic cycle every 357 days). The amplitudes of the spectra are homogeneously distributed and are subjected to negative trends at all frequencies (Figure 6).

Hydrological time series are often highly random. In order to study the character of the available hydrological time series, an analysis method frequently used in the study of statistical fractals, corresponding to the log-log representation of the variance density spectra, is applied. This method makes it possible to identify the Gaussian, Brownian, or deterministic character of a data series. The slope β of the log-log density spectrum assumes values between +1 and −1 for fractional Gaussian noise and between −1 and −3 for fractional Brownian motion. A zero slope (β = 0) is characteristic for pure Gaussian noise, and a slope β = −2 is characteristic for the pure Brownian domain. Slopes in the range −2 to −3 are characteristic of the persistent Brownian domain, while slopes in the range −1 to −2 are characteristic of the antipersistent Brownian domain.

The spectral analysis of the daily precipitation time series allows us to observe a linear behavior over the scale range, which extends between one day and 15 days (Figure 6a and

Table 3. Statistical fractals of the main hydroclimatic time series of the Sebaou River basin.

Time Series	Stations	Period	Slope (β1)	Scale Invariance Ranges	Slope (β2)	Scale Invariance Ranges
Daily rainfall (mm/day)	Tizi Ouzou	1990–2009	−0.21	14 days–1 year	−0.66	1–13.5 days
	Ait Aicha	1972–1991	−0.15	9 days–1 year	−1.10	1–8.5 days
		1991–2010	−0.32	11 days–1 year	−1.03	10–13 days
	DEM	1967–1988	−0.26	16 days–1 year	−0.82	1–15 days
		1988–2010	−0.002	16 days–1 year	−0.88	1–15 days
	Freha	1972–1991	−0.27	10 days–1 year	−0.89	1–9 days
		1991–2010	0.07	11 days–1 year	−0.88	1–10 days
	Beni Yenni	1972–1991	−0.09	10 days–1 year	−1.10	1–9 days
		1991–2010	−0.10	11 days–1 year	−0.73	1–10 days
Daily runoff (m3/s)	Belloua	1949–1958	−0.26	11days–1 year	−1.25	1–10 days
		1972–1983	−0.22	12 days–1 year	−1.14	1–11 days
		1987–2000	−0.37	12 days–1 year	−2.98	1–11 days
	Baghlia	1963–1985	−0.32	12 days–1 year	−2.85	1–13 days
		1985–1997	−0.01	13 days–1 year	−2.24	1–12 days
	Freha	1986–2001	−0.28	20 days–1 year	−1.60	1–19 days
	Boubhir	1987–2002	−0.13	13 days–1 year	−1.45	1–12.5 days
	RN25	1973–1994	−0.75	14 days–1 year	−2.21	1–15 days
	RN30	1985–1998	−0.48	20 days–1 year	−2.43	1–19 days
		1998–2010	−0.39	30 days–1 year	−1.61	1–29 days

), often encountered in the literature, e.g., [72]. The upper limit of the domain is not very clear. It is always possible to implement, in addition, an automatic detection procedure for linear portions, if the user wishes to make the location of the rupture more objective. The invariance ranges of the analyzed scales are characterized by an exponent of the spectrum less than 1 (−0.002 < β < −1.10).

Short-term noise analysis places the streamflow at Belloua station in the fractional gaussian noise domain with the slope β equal to −0.97 for the 1972–1984 period, and the slope β is strong enough for the high frequencies, corresponding to a fractional Brownian motion, which is −1.40 for the 1987–2000 period (Figure 6b and Table 3). These time series, therefore, represent an unstructured random phenomenon for the first period and typical of a quasi-deterministic phenomenon for the second period.

In general, the log-spectral analysis of the daily streamflow time series allows the classification of the annual spectra into two different groups according to the average slope value β. The results show a composite character with a steep slope β (−1.25 < β < −2.98) for high frequencies, where the series is non-conservative, whereas for low frequencies, the slope is greater than −1 (−0.01 < β < −0.48). For this group, two scales with linear behavior are highlighted, with a clear break in the spectrum around 8.5 to 29 days approximately, which is clearly demonstrated in the rainfall spectra (Figure 6b). Two ranges of scale invariance separated by a break were retained, the first between 1 day and 3 weeks and the second between 3 weeks and 1 year, i.e., two distinct types of variability regimes. In addition, this break is identified as the synoptic maximum, defined as the maximum time scale for which synoptic events persist. According to Rossi et al. (2011) [73], the values of the β spectral exponents provide information about the variability of each signal, depending on the time scale considered: a low rate of β is linked to a high variability in the data, while a high rate of β is related to low variations. According to the results of Tessier et al. (1993) and De Lima and Grasman (1999) [74,75], all obtained values of β were less than 1 and highlighted the conservative character of the rainfall time series, which is very much in agreement with the results of this study. Therefore, the two diverse slopes ß1 and ß2 characterize two categories of variability regimes, respectively related to short- and long-term patterns. The results of the log-log spectral analysis of the available rainfall time series are presented in Table 3.

4.3. Cross Analyses

The system’s response to a rainfall impulse occurs in two modes according to Marsaud (1997) [76]:

-

A transmissive function corresponding to the peak explains a well-developed drainage.

-

A slower decrease explains the capacitive effect.

The cross-correlograms in Figure 7 consist of well-individualized responses for the Acsif N’Boubhir and Sebaou Rabta systems and, less obviously, for the Boughdoura system, which explains the well-developed drainage. At Sebaou Rabta, the 5-day peak represents a rapid flow, and the 36-day peak represents a delayed flow (Figure 7). However, the Aïssi River system does not show an increase in the correlogram, which indicates that the response to rainfall impulses is probably very low, explaining that the system performed a filtering of the rainfall signal. In addition, the human impact and the construction and ordering of the Taksebt dam since 1998 may have influenced the rainfall-runoff relationship (Figure 7).

The cross-amplitude function (Figure 8a) shows that the cause-and-effect relationship is strong for low frequencies (f equal to 0.032, which corresponds to 31 days), especially for the Boughdoura, Sebaou Rabta, and Acsif N’Boubhir systems. From this cutoff frequency of 0.032 day⁻¹, the covariance becomes negligible, with the exception of some frequencies lower than 0.2, which have some characteristic peaks. According to El Hakim (2005) [77], the high covariance indicates that these periods correspond to a good cause-and-effect relationship. The amplitude function shows that the impulse response is valid for low and medium frequencies, especially for the A’N boubhir and Rabta sub-basins. This confirms the low inertia behavior of these systems (Figure 8a).

The lag representative of the rainfall-runoff relationship is calculated for a cutoff frequency (i.e., 0.032 = 32 days) and gives a delay value of about 10 days (Figure 8b). This value is relatively high, which means that the rise in flood duration is a little long and can only be noticed at the outlet after 7 to 10 days.

The gain function, which reflects the way in which the rainfall is put into reserve, shows poor attenuation at high frequencies (Figure 8c). There are attenuations for the three frequency modes of the different systems. This attenuation of rainfall, on the one hand, may correspond to a short-term effect associated with a reserve in flood events, except for some insignificant low frequencies peaks in the Acsif N’Boubhir and Sebaou Rabta systems. This amplification is due to the low filtering of the system.

The coherence function represents an indicator of linearity, which shows highly oscillating values varying between 0.014 and 0.810 with a negative trend of overall average coefficients equal to 0.430. The coherence function (Figure 8d) shows a probably poor linearity of the rainfall-runoff relationship, which is related to the replenishment of reserves during some days of the year.

The streamflow variation depends on the groundwater-rainfall system response and specifies the conductivity of the fractures, which expresses the porosity. This latter contributes to the propagation of rainfall in the form of infiltration in heterogeneous aquifers. This is perfectly evident in the case of flow: sources and rivers, in the case of heavily karstified areas [78].

The Morlet XWT (Figure 9b,c) of daily data of the Sebaou River basin system was calculated to highlight the temporal change of the rainfall-runoff relationship by a common power spectrum, in which the continuous wavelet coefficients of precipitation are multiplied by the other one of the discharges [53]. The results are mapped beforehand into the time-frequency domain.

The XWT spectra (Figure 9b,c) reveal common significant power at 95% confidence interval. The cross spectra do not highlight the seasonal components, except for few insignificant localized events (Figure 9a,b).

The annual scale, which is represented by the 256–512-day band, and the multi-annual scale, which is represented by spectral bands greater than 512 days, indicate a fairly continuous component during the study periods with strong coefficients for all the spectra (Figure 9a–c and

Table 4. Extracted significant correlations between daily rainfall and streamflow time series using XWT and WCT.

Spectral Bands	Years with Significant Correlations Based on cross Wavelet Spectra (XWT) Analysis
	Acif N’boubhir River	Aïssi River	Boughdoura River	Sebaou Rabta River	Sebaou Sebt River	Sebaou Maritime River (Outlet)
6–8 month	1998, 2000 2006, 2007	1997, 2007	1974, 1976	1993, 1995, 1998, 2000	1994–1995	1974, 1978, 1983, 1987 1989
1 year	1987–2010	1986–2002 2002–2010	1974–1994	1989–2001	1990–1997	1974–1998
1–3 year	1996–2004	1996–2004	1983–1991	1990–1999	1993 1994–1997	1975–1981 1986–1993
3–6 year	1996–2005	1997–2005	—	1995–1999	1994–1997	1979–1983
6–8 year	1996–2004	1993–2004	1981	—	—	1984–1987
8–12 year	1996–2004	1996–2001	—	—	—	—
Spectral Bands	Years with Significant Correlations Based on Wavelet Coherence Spectra
6–8 month	1988–2000	1986	1977, 1979, 1981, 1984, 1987, 1990	1988–1990 1992–1994 1998–2001	1991, 1997	1973, 1976 1983, 1985
1 year	1988–2010	1991–1999	1974–1976 1979–1994	1988–2001	1990–2000	1973–1975 1979–1993
1–3 year	1988–2010	2003–2010	—	1988–2001	—	1987–1993
3–6 year	1995–1999	—	1978–1990	1992–1999	—	—
6–8 year	1996–2004	1999–2004	—	—	—	1984–1987
8–12 year	—	—	—	—	—	—

), with the exception of a few marked discontinuities between the 1980s and 1990s, which are related to the dry periods. The amplifications of the annual components are visible, dominant, and more amplified compared to the multi-annual fluctuations, characterizing the 1–3- and 3–6-year modes (Figure 9a–c and Table 4). This can be explained by the multi-annual regulation shown by the CSA. The attenuations of short-term components seem to be linked to the periods of storage (infiltration).

As a first observation on the WCT spectra (Figure 9d–f), we can notice a decrease in the correlation surfaces (less important dark red color) compared to the XWT spectra. Thus, the wavelet coherence transform really indicates the degree of dependence between the input-output systems.

At small scales (32–128 day), dispersed structures are marked due to successive passages of floods particularly visible in the system of Acsif N’Boubhir and Sebaou Rabta rivers. At an annual scale (256–512 day), the most significant consistency between rainfall and streamflow is clearly evident during the study period, with a transition to the multiannual modes at 512–1024-day band (Acsif N’Boubhir and Sebaou Rabta rivers) (Figure 9d–f), but sometimes there is intermittence. The interruption in the annual component at Aïssi River between the 1999–2002 period (Figure 9d–f and Table 4) was previously explicated by the construction and ordering of the Taksebt Dam (Figure 9c,d). According to Chinarro et al. (2012) [53], high coherence values explain the linear response of the system during wet periods, whereas the response to the rainfall event is not linear when the period is dry. Therefore, it can be concluded that two main modes are the wet mode with linear behavior that corresponds to the annual process, while the dry mode with non-linear behavior corresponds to successive drought events and anthropological activities, which affect the annual flow and regulatory reserves. According to Hadjou (2008) [79] and the report of the department of Water Resources of Tizi Ouzou, the boreholes in the downstream part of the Aïssi bridge are currently only producing a quarter of their initial capacity.

This drop in water production is due to the lack of water resources after the dam was filled. The piezometric level statement for the year 2001 show that the water table of the Aïssi River is drying up, leaving only a meter of wet alluvium in place. Khelifa et al. (2021) [80] also observe this in the Seybouse River (Algeria), where further water is exploited and used for agriculture using advanced water pump technology. This decrease in water amount and uninterrupted industrial toxic discharges decrease water quality and inflict damage on the aqua-fauna system. In the Medjerda River (Algero-Tunisian basin), Kadir et al. (2020) [81] observed overexploited groundwater with intensive water extraction for agriculture during the 2000−2012 period (0.004 Mm³/year).

From the XWT and WCT spectra, the evolution of some components was isolated and analyzed using the theory of the cross wavelet and wavelet coherence [82], i.e., the isolated cross wavelet spectrum components (Figure 10a), the real squared part (Figure 10b), and the phase components (Figure 10c) of the wavelet coherence function. For high frequency processes, 2- and 4-week components were chosen. However, the 26- and 52-week components were chosen for the medium-term and the 102-week component was selected as the long-term components (Figure 10a–c). Very poor consistencies were observed for the 2–4-week components, since the real part of the consistency is very low, which agreed with the phase that is diverse from zero. This indicates that the rainfall-runoff relationship was not linear (Figure 10a–c), except for some characteristic peaks. Between the 26-week seasonal and the 52-week annual components, there are some highly consistent values just for the Sebaou Rabta and Sebaou Maritime river system, but the phase is always non-zero. This also supports the rainfall-runoff nonlinearity relationship (Figure 10a–c). There is a strong low frequency coherence of the 102-week component, as the coherence is somewhat strong in the real part, while the phase oscillates close to zero for some time interval. This approves the pseudo-linearity of the rainfall-runoff relationship for this component, in particular for the Sebaou Maritime River (outlet) and the Boughdoura River (Figure 10a–c). According to Chettih and Mesbah (2010) [82], the highest magnitudes of coherence are related to successive events of floods and rainfalls; this shows how quickly the systems respond. During wet periods, a “good” match is observed, while during dry periods no match is evidenced.

The phenomenon of looting the rocks and sands of the valleys of Sebaou River exists, as the looting of sand is a source of profit for thousands of families and sand barons, who take advantage of the nocturnal conditions and the lack of inspection on some places of the valleys to proceed with these illegal mining activities. These activities have succeeded in draining large amounts of sand, especially at the level of the Bouqandora Valley, threatening to cause a serious disaster in agricultural land adjacent to those valleys, as well as in eastern regions such as Ait Issa Maimon and Aknoun, which are still suffering from the emergence crises in shallow drinking waters [83].

The water reserves of this valley, which extend from the municipality of Tadmaet (Tizi Ouzou) and flow into the Mediterranean Sea in the region of Benchod (Boumerdes) for a length of 20 km, are threatened by pollution due to their mixing with the seawater. An unusual percentage of salinity appeared in the chemical elements in water, which are chlorine, calcium, and bicarbonate. The increase in water salinity is due to the low water level of the valley and the theft of sand, especially during the years of low rainfall, as well as the excessive exploitation of wells for agriculture and the construction of dams, which contributed significantly to the aggravation of this problem. According to the Algérie Presse Service [84], there was a 30% decrease in agricultural production in the region due to the contamination of water wells used for irrigation by seawater.

Large violations of the riverbed by some careers and excavators turned up the riverbed and extracted large amounts of sand and rock, despite pronouncements by the control authority to prevent these crimes to natural resources (Figure 11). According to Akanwa (2021) [85], the sands extracted from the river banks and riverbeds are the most desired sands for construction purposes due to their capability to bond to concrete. All over the world, many researchers have been dedicated and widely debated the problems of human effects on streamflow and rivers. It was concluded by [86] that sand mining in the Progo River (Indonesia) provides employment opportunities for the population of the area, which creates a good economic level, but the impact of sand mining on the stability of the riverbed is witnessing a deteriorating trend of −0.43 m/year. In the lower Gangetic Plain in India, the impact of human intervention on downstream channel behavior was evaluated using Landsat images and GIS, and an intensive field survey was carried out to analyze the river cross-section, indicating a modification in the morphological characteristics due to human-nature interface, which induced the upliftment of the riverbed by continuous siltation as well [87]. In the Yellow River basin, it was found that anthropological activities contributed about 80% to the decrease in runoff, while the other 20% was attributed to rainfall [88]. In the Yangtze River, significant changes have been observed since the 1950s in the runoff and sediment discharge dynamics, which has attracted great social attention. In addition, it was found that the anthropological activity impact on sediment discharge was greater than the runoff from 1953 to 2010 [89]. In southwestern China at the dammed Lancang River, impacts of the climate variability and anthropological activities on the flow regime were assessed by Han et al. (2019) [15]. The analysis showed that the construction of the reservoir was the anthropological activity that most severely affected the streamflow, reducing the streamflow during the rainy season.

5. Conclusions

The Sebaou River basin is experiencing significant demographic and socioeconomic development. Groundwater in this basin is one of the major resources for this growing population, as it supports agricultural and industrial activities and ensures the supply of drinking water. However, it is worth highlighting that the ground and surface waters are currently overexploited. In this study, several time-scale-based methods were used to assess and understand the karst spring discharge response for better planning of these water resources.

The CSA revealed that the hydrogeological systems in the study area are characterized by various memory effects, which can be classified as small, poor, reduced, and extensive, with regularization time ranging between 5 and 50 days. The results indicate that the aquifer systems of the Sebaou River basin have a very complex hydrodynamic behavior.

Amplifications of the annual components are visible, dominant, and more amplified compared to the multi-annual fluctuations and characterized by the 1–3- and 3–6-year modes. This can be explained by the multi-annual regulation shown by the CSA. Attenuations of short-term components (i.e., 2–5-, 26-, and 52-week) seem to be linked to periods of storage (infiltration) with two main modes; i.e., wet mode with linear behavior corresponding to the annual process and dry mode with non-linear behavior corresponding to successive drought events in the 1980s and 1990s and anthropological activities, which affected the annual flow and regulatory reserves.

The impacts of human activities on the streamflow due to the looting of rocks and sands from the valleys of the Sebaou River reached alarming levels, which requires urgent intervention to protect water and ecological resources and their better rational use. Therefore, in-depth field studies in the lower course of the river must be initiated and urgent interventions and practical short- and medium-term strategies must be proposed, which include high monitoring of water flows. In addition, the publication of serious statements and the imposition of heavy financial fines to curb the activities of thieves of sand, gravel, and rocks of the river, especially where invasive aquatic species have reached a catastrophic situation, requiring the intervention of specialists in the field in order to revive aquatic ecosystems. Coordinated actions between land use planning, agriculture, energy, and industry managers are essential to present a better solution and achieve the sustainable development of water and other natural resources at a basin scale.

Finally, the main strength of this study is the use of wavelet-based techniques to provide additional information in a time-scale domain, which allows us to understand and make appropriate decisions of the behavior of the hydrogeological system for better water storage and supply. However, it is worth noting that the main weakness of this study was the use of short time series, which prevented us from dividing them into sub-series to detect the change in behavior and some characteristics of the studied time series, especially because the wavelet approach gives good results on stationary and non-stationary phenomena when applied to long time series.