Assessment of Spatial Distribution and Temporal Variations of the Phreatic Groundwater Level Using Geostatistical Modelling: The Case of Oued Souf Valley—Southern East of Algeria

: Since the beginning of the 1980s, several regions in the northern Sahara of Algeria have been confronting the rising groundwater. Among all these regions, Oued Souf Valley represented one of the most acute affected by this phenomenon. Due to the natural topography and the in-sufﬁcient/weakness of water management and miscoordination between different sectors that are represented by intensive exploitation of deep groundwater reservoirs which returns to the shallow aquifer, absence of sewage and drainage network, leakage from drinking water supply system, the groundwater has raised to the surface or near to the surface, affecting the traditional cultural environment and urban areas and degrading all socio-economic aspects of the Oued Souf habi-tants. To preserve the Oued Souf environment, a vertical drainage system has been constructed. Consequently, in this research, an evaluation of the vertical drainage system performance and its impact on groundwater level stabilization has been performed by mapping the water table of the phreatic groundwater level using geostatistical modeling using ordinary kriging (OK) interpolation method, which has been applied to analyze the spatial and temporal structure of groundwater level ﬂuctuation. Meanwhile, hierarchical cluster analysis (HCA) was applied for grouping the wells based on the groundwater ﬂuctuations for 2008, 2009, 2014, 2016, 2018, and 2021. However, the vertical drainage system reﬂected a signiﬁcant decline of groundwater from 2009 to 2018 due to the important drained volumes through it but another rising phenomenon might be threatening the region in the near future and this is what was indicated in the 2021 groundwater level data. Cluster analysis has generated four groups based on their ﬂuctuation means that are increasing from the ﬁrst group to the fourth group ascendingly. The ﬁrst cluster grouped the drains that have a shallow depth (average mean of 5.91 mbgl) and declined over the clusters. The clusters are spatially combined with signiﬁcant separation of the fourth cluster which represents the deepest group (12.89 mbgl). Based on this research, several factors are inﬂuencing the stability of the phreatic groundwater level and even the performance of the drainage system, the most important of which is the overexploitation from deep groundwater reservoirs such as complex terminal and continental intercalary (in drinking and irrigation) and even the illegal use of the phreatic groundwater with important quantities for irrigation and illegal industries.


Introduction
Groundwater is one of the most important elements of the water cycle and is found almost everywhere beneath the Earth's surface (hills, mountains, plains, deserts, coasts, and the ecosystem and the traditional urban character of the city), as well as the lower-lying areas of the city and its suburbs. Then the phenomenon progressed to the point of touching a large part of Souf from Guemar in the north to El Ogla in the south. At the same time, it was reported that other areas in the east of the region, around Z'goum were exposed to the rising groundwater but its effect was limited. Scattered areas around Ourmes, Reguiba in the northwest, Oued Alenda and Mih Ouenssa in the southwest were affected. At that time, serious pollution and health problems have arisen, such as the rise of the polluted water table, which is filling with blackish and nauseating water mixed with all kinds of waste, because of the flooding of septic tanks or cesspools that allows bacteriological pollution to spread in the superficial aquifer and the contamination of the water by nitrates of domestic and agricultural origin. The heavy mineralization of the water is also caused by the heavy evaporation from open water bodies and the dissolution of salts. This is because the Ghouts are also used as wild rubbish dumps. On the other hand, the stagnant water leads to a proliferation of mosquitoes and increases the number of inhabitants affected by waterborne diseases and parasitic diseases, skin diseases, leishmaniasis, malaria, and typhoid. For example, several hundred cases of typhoid fever have been recorded in the town of El Oued, according to the municipal health service [21].
In an attempt to halt this tragedy, Algerian authorities responded in 2005 by launching a megaproject to combat increasing groundwater levels and pollutant spread. The Sewerage Plan, Purification Plan, Evacuation Plan, and Drainage Plan are the four plans that make up this project. Unfortunately, this project failed to combat its main objective due to several problems. However, Oued Souf was and still is an open book for various studies to understand the principal reasons for shallow groundwater upwelling and its forecasting trying to refer this phenomenon to the hydrogeological setting of the region or anthropogenic activities [22][23][24][25][26][27][28].
The fluctuations of the water table and its formation in a region are largely controlled by various conditions such as terrain, slope, topography, lithology, rainfall patterns, surface water bodies, etc. Despite all these factors, we believe that there are influencing factors more related to anthropogenic activities that affect the spatial and temporal fluctuation of the groundwater table, especially in Saharan regions such as Oued Souf valley. Therefore, this research aims to answer the following questions:

•
How is the vertical drainage system performing and what is its impact on fighting upwelling of the phreatic groundwater table and its stability between 2008,2009,2014,2016,2018, and 2021 in terms of spatial distribution and temporal variation of the phreatic groundwater fluctuation? • How can the weakness of water resources management in El Oued contribute to the continuity of the upwelling of groundwater? However, the phreatic groundwater level fluctuations point to the need for a precise and complete investigation to understand how groundwater level fluctuations behave on both temporal and spatial dimensions. Geostatistical modeling is an extremely valuable technique for studying such processes [29]. Artificial Intelligence methodologies, especially statistical modeling techniques that are from the family of machine learning (ML), algorithms, have proven to successfully complement or even replace the classical process-based models, usually requiring considerably less data elaboration and computing time; they are therefore more easily implemented in computes the dynamics of groundwater [30].
Geostatistical techniques with its various interpolation methods such as ordinary Kriging interpolation (OK), simple Kriging interpolation (SK) and universal Kriging interpolation (UK), deterministic interpolation comprise global polynomial interpolation (GPI), inverse distance weighted interpolation (IDW) and planar spline interpolation and local polynomial interpolation (LPI) are commonly used to generate maps of soil properties and have been described in many studies [31][32][33], mapping and predicting groundwater quality and its contamination [34][35][36] and mapping and predicting the groundwater level [37][38][39][40][41][42]. On the other hand, hierarchical cluster analysis (HCA) is one of the most commonly used methods to classify instances of variables according to their variance or similarity. The wide application of this method has recently extended to practical hydrochemistry investigations for clustering the hydrogeochemical processes in groundwater in samples into different groups that are important in the geological and hydrogeological sense [43][44][45][46]. Therefore, for a better understanding of this problem and to answer the questions generated previously, a methodology was used based on geostatistical modeling based on (OK) interpolation method which has been applied to analyze the spatial and temporal structure of groundwater level fluctuation. Meanwhile, (HCA) was applied for grouping the wells based on the groundwater fluctuations.

Study Area Description
Geographically, the Oued Souf valley lies in the southeast of Algeria in the province of El-Oued, in the middle of a large synclinal basin. It is known as the low Sahara because of its low position between 32 • 30 00 and 34 • 12 00 N as latitude and 6 • 15 00 and 7 • 20 00 E as longitude. Historically and administratively, Oued Souf became a municipality in 1957 and an official province in 1984. Today, the Oued Souf valley covers an area of 11,738 km 2 , divided into 18 municipalities. In 2009, 500,000 inhabitants lived in the 18 municipalities. According to 2015 data from the National Statistics Office, the population of Oued Souf exceeds half a million, mainly concentrated in a group of large cities in the region (41.41 inhabitants/km 2 ) [47]. The Oued Souf region borders the Republic of Tunisia to the east and the province of Tebessa to the northeast, Biskra, and Khenchla to the north, and Biskra to the northwest. At the same time, it is bordered to the west and southwest by Djelfa and Ouargla, respectively, as is shown in Figure 1.
(GPI), inverse distance weighted interpolation (IDW) and planar spline interpolation and local polynomial interpolation (LPI) are commonly used to generate maps of soil properties and have been described in many studies [31][32][33], mapping and predicting groundwater quality and its contamination [34][35][36] and mapping and predicting the groundwater level [37][38][39][40][41][42]. On the other hand, hierarchical cluster analysis (HCA) is one of the most commonly used methods to classify instances of variables according to their variance or similarity. The wide application of this method has recently extended to practical hydrochemistry investigations for clustering the hydrogeochemical processes in groundwater in samples into different groups that are important in the geological and hydrogeological sense [43][44][45][46]. Therefore, for a better understanding of this problem and to answer the questions generated previously, a methodology was used based on geostatistical modeling based on (OK) interpolation method which has been applied to analyze the spatial and temporal structure of groundwater level fluctuation. Meanwhile, (HCA) was applied for grouping the wells based on the groundwater fluctuations.

Study Area Description
Geographically, the Oued Souf valley lies in the southeast of Algeria in the province of El-Oued, in the middle of a large synclinal basin. It is known as the low Sahara because of its low position between 32″30′00″ and 34″12′00″ N as latitude and 6"15′00″ and 7″20′00″ E as longitude. Historically and administratively, Oued Souf became a municipality in 1957 and an official province in 1984. Today, the Oued Souf valley covers an area of 11,738 km 2 , divided into 18 municipalities. In 2009, 500,000 inhabitants lived in the 18 municipalities. According to 2015 data from the National Statistics Office, the population of Oued Souf exceeds half a million, mainly concentrated in a group of large cities in the region (41.41 inhabitants/km 2 ) [47]. The Oued Souf region borders the Republic of Tunisia to the east and the province of Tebessa to the northeast, Biskra, and Khenchla to the north, and Biskra to the northwest. At the same time, it is bordered to the west and southwest by Djelfa and Ouargla, respectively, as is shown in Figure 1.  Agricultural activities in the region, shown on a land use/land cover map published in [48], consist mainly of irrigated agriculture, particularly the cultivation of date palms planted in large hand-dug craters in the dunes. This method of palm cultivation is known in the local vernacular as "Ghout" or "Ghout hollows". Due to the capillarity of the cultivation method, water from the shallow aquifer can rise to the deep roots of the date palms, allowing farmers to irrigate at a low cost. However, in recent years, agriculture in Oued Souf has expanded to include other crops such as potatoes, olives, maize, and other vegetables, all of which are irrigated using groundwater from the aquifer beneath the whole Oued Souf region [49,50]. The topographical nature of the El Oued Souf region, also known as the lower Sahara region because of its low altitude and is located in the Septentrional Sahara. As it is exposed in Figure 2, the study area has an average elevation of 78.54 m above sea level (m.a.s.l), with a notable significant decrease from the southwest to the northeast of the study area, reaching about 66.42 m.a.s.l in drain 13 (D13).
planted in large hand-dug craters in the dunes. This method of palm cultivation is known in the local vernacular as "Ghout" or "Ghout hollows". Due to the capillarity of the cultivation method, water from the shallow aquifer can rise to the deep roots of the date palms, allowing farmers to irrigate at a low cost. However, in recent years, agriculture in Oued Souf has expanded to include other crops such as potatoes, olives, maize, and other vegetables, all of which are irrigated using groundwater from the aquifer beneath the whole Oued Souf region [49,50].
The topographical nature of the El Oued Souf region, also known as the lower Sahara region because of its low altitude and is located in the Septentrional Sahara. As it is exposed in Figure 2, the study area has an average elevation of 78.54 m above sea level (m.a.s.l), with a notable significant decrease from the southwest to the northeast of the study area, reaching about 66.42 m.a.s.l in drain 13 (D13).

Climatology of the Study Area
Oued Souf is characterized by a hyper-arid climate characterized by low precipitation, high temperatures, furthermore, intensive evaporation. Based on climatic data from the meteorological station of Geumar Airport over 42 years (1978-2020), the total annual rainfall was 67.63 mm per year. The maximum total rainfall was observed in 2009 at 233.7 mm/year, and the lowest annual total rainfall was remarked in 2020. The origin of precipitation in the Saharan regions is different according to the seasons. During the summer they are due to monsoon depressions, in winter they are due to depressions accompanying the southward migration of the polar fronts. During the intermediate period, these precipitations are due to the Sudano-Saharan depressions crossing the Sahara from the south to the north [51].

Climatology of the Study Area
Oued Souf is characterized by a hyper-arid climate characterized by low precipitation, high temperatures, furthermore, intensive evaporation. Based on climatic data from the meteorological station of Geumar Airport over 42 years (1978-2020), the total annual rainfall was 67.63 mm per year. The maximum total rainfall was observed in 2009 at 233.7 mm/year, and the lowest annual total rainfall was remarked in 2020. The origin of precipitation in the Saharan regions is different according to the seasons. During the summer they are due to monsoon depressions, in winter they are due to depressions accompanying the southward migration of the polar fronts. During the intermediate period, these precipitations are due to the Sudano-Saharan depressions crossing the Sahara from the south to the north [51].
The monthly variability of rainfall and temperature analysis over 42 years (1978-2020) revealed that the maximum rainfall is around 13.89 mm recorded during January, while the minimum is around 0.23 mm recorded during July. The maximum value of temperature recorded between (1978-2020) is 37.43 • C during July. While the minimum observed temperature is 14.76 • C during January. The hydrological year of the study area is branded by the presence of a dry period and the total absence of a wet period over the year, even during January when it was observed lowest temperature and highest precipitation.
The application of the Thronthwaite approach revealed that the potential evapotranspiration (PET) decreases from July to January from 260.13 mm to 62.16 mm, then inversely increases again, the annual amount of PET is 1850.7 mm. After the real evapotranspira-tion calculation (RET), it was observed that the amount of precipitation is lower than the RET value.
The easily utilizable reserve (EUR) [52,53], computation gave negligible value. After the application of the Tixeront-Berkaloff formula, the runoff (R) was also negligible (0.0025 mm) with a runoff coefficient equal to 0.000445. Consequently, the infiltration amount (I) was null due to high evapotranspiration which is more than the precipitation (P). At the same time, according to the obtained results, there is an agricultural deficit (AD) throughout the year, hence the need for irrigation Table 1.

Geomorphology, Geology and Hydrogeology
The valley of Oued Souf is located on the northern side of the Great Oriental Erg and is typified by continental sand dunes formed during the recent Quaternary. The Quaternary formations (sand dunes) are mainly composed of fine-grained, compact, homogeneous, and uniform sand that covers the entire region, especially in the south where the dunes can reach up to 100 m high. Between the sand dunes occur corridors that form deepened plateaus (known as Sahans), which are typically quite extensive, sometimes stony, and covered with Quaternary gypsum deposits. The saline depressions (Sebkha) are located in the northernmost part of the study area and occupy the basal part of the vast lower Saharan basin [54][55][56][57]. The Oued Souf region lies in the northeastern Mesozoic basin of the Sahara (also called the Triassic basin), which is situated to the northeast of the Sahara platform. The region is underlain by sediment formations ranging in age from the Lower Cretaceous to the Quaternary [58]. The base of the sedimentary basin is formed by water-rich marine formations from the Paleozoic era, overlain by a sequence of layers that can extend over 2000 m [59][60][61]. Despite the geology and geomorphology of the region, the prevailing climatic conditions hinder the development of surface water resources [51]. Therefore, groundwater is the only accessible water resource in the Souf region that can be used for various purposes [62,63].
The study area is part of the Northwest Sahara Aquifer System (NWSAS), which is the second-largest reservoir in the world. This aquifer system is shared by three countries: Algeria, Tunisia, and Libya [64]. It has a total area of about one million km 2 , of which 70% is located in Algeria, 6% in Tunisia, and 24% in Libya. The aquifer system of the study area consists of three primary aquifers that differ in depth and physiochemical characteristics [65][66][67]. The phreatic aquifer is the uppermost aquifer in the region [68]. The free water table in this aquifer consists of fine sands interbedded locally with sandy clays and gypsum lenses. This aquifer has a thickness of about 100 m and a depth to groundwater of 1 to 40 m. Underlying this aquifer is an impermeable clay substrate. Thousands of traditionally dug wells extract groundwater from the phreatic aquifer. The total number of traditional wells exceeded 35,000 in 2015. The average permeability of the phreatic aquifer is 10 −4 m/s, while the horizontal transmissivity and storage coefficient is estimated to be 10 −2 m 2 /s and 0.2, respectively. Natural recharge of the phreatic aquifer can be represented by infiltration of precipitation and runoff from the rock formations on the southern margin of the Great Oriental Erg and, simultaneously, by rare torrential local rain events, such as those that fell in the region in April 1947 and May 1967 [69][70][71][72][73][74]. The complex terminal is the term for the next deeper aquifer. This complex aquifer consists of several aquifers in different geological formations. The complex terminal aquifer belongs to the continental formations of the Cretaceous to Miocene. Its depth is between 400 and 600 m, while its thickness is usually around 400 m. This aquifer contains fossil water that is between 20,000 and 30,000 years old. In 2015, about 182 deep wells were drilled in this aquifer, 28 of which are for irrigation and 154 are for municipal purposes and drinking water supply [13,75,76].
The third deepest aquifer is called the continental intercalary aquifer and consists of continental deposits from the Middle Jurassic to the Lower Cretaceous (Barremian and Albian). The lithology of the aquifer consists of sandstones and clayey sandstones. This aquifer is located at varying depths from 1800 to 2200 m and has a thickness of 200 to 400 m. This aquifer contains fossil water that is at least 20,000 to 30,000 years old. Only four deep wells have been drilled in the continental aquifer of the region, all of which are used for drinking water production, as the groundwater from this depth has a temperature of over 70 • C [50,77]. The reservoirs composition of the Northwest Sahara Aquifer System (NWSAS) is revealed in Figure 3. of 1 to 40 m. Underlying this aquifer is an impermeable clay substrate. Thousands of traditionally dug wells extract groundwater from the phreatic aquifer. The total number of traditional wells exceeded 35,000 in 2015. The average permeability of the phreatic aquifer is 10 −4 m/s, while the horizontal transmissivity and storage coefficient is estimated to be 10 −2 m 2 /s and 0.2, respectively. Natural recharge of the phreatic aquifer can be represented by infiltration of precipitation and runoff from the rock formations on the southern margin of the Great Oriental Erg and, simultaneously, by rare torrential local rain events, such as those that fell in the region in April 1947 and May 1967 [69][70][71][72][73][74]. The complex terminal is the term for the next deeper aquifer. This complex aquifer consists of several aquifers in different geological formations. The complex terminal aquifer belongs to the continental formations of the Cretaceous to Miocene. Its depth is between 400 and 600 m, while its thickness is usually around 400 m. This aquifer contains fossil water that is between 20,000 and 30,000 years old. In 2015, about 182 deep wells were drilled in this aquifer, 28 of which are for irrigation and 154 are for municipal purposes and drinking water supply [13,75,76].
The third deepest aquifer is called the continental intercalary aquifer and consists of continental deposits from the Middle Jurassic to the Lower Cretaceous (Barremian and Albian). The lithology of the aquifer consists of sandstones and clayey sandstones. This aquifer is located at varying depths from 1800 to 2200 m and has a thickness of 200 to 400 m. This aquifer contains fossil water that is at least 20,000 to 30,000 years old. Only four deep wells have been drilled in the continental aquifer of the region, all of which are used for drinking water production, as the groundwater from this depth has a temperature of over 70 °C [50,77]. The reservoirs composition of the Northwest Sahara Aquifer System (NWSAS) is revealed in Figure 3.

Vertical Drainage System
Oued Souf was and still is among the areas that knew an aggressive rising of polluted groundwater over the last years, Figure 4 addresses the main effects that exacerbated the phenomena. For this reason, a sanitation project (mega project) was constructed to counteract the rising water level in the phreatic aquifer via an additional drainage system (vertical drainage system). It consists of draining excess water from agriculture through a network of 58 drainages drilled at a depth of 21 to 40 m and equipped with submersible pumps capable of pumping 6 L per second. The drainages are spaced about 500 m apart and connected to 37 km of pipelines to equalize the water level in Oued Souf City. To stabilize the water table in the flooded neighborhoods and ensure a minimum depth of 1.5 m, as well as to reverse the phenomenon of upwelling in the rest of the city and prevent new areas from being flooded, a general drawdown of 5 to 10 m deep will be ensured, allowing autonomous sanitation of the areas that cannot be connected. The drainage system is currently managed by the El-Oued National Sanitation Office (l'Office National d'Assainissement d'El-Oued ONA). The principle of this system is to collect all the infiltration water under the agglomeration of El-Oued by pumping it into 58 drilled wells and to use as much of this water as possible on-site for the irrigation of green spaces, while the rest of this water is collected in the pumping station and piped over a distance of about 4200 m to the collector of the treated wastewater of the treatment plant to be treated by aerated lagoons then discharged in the discharging site (Chott Halloufa) [79].
tical drainage system). It consists of draining excess water from agriculture through a network of 58 drainages drilled at a depth of 21 to 40 m and equipped with submersible pumps capable of pumping 6 L per second. The drainages are spaced about 500 m apart and connected to 37 kilometers of pipelines to equalize the water level in Oued Souf City. To stabilize the water table in the flooded neighborhoods and ensure a minimum depth of 1.5 m, as well as to reverse the phenomenon of upwelling in the rest of the city and prevent new areas from being flooded, a general drawdown of 5 to 10 m deep will be ensured, allowing autonomous sanitation of the areas that cannot be connected. The drainage system is currently managed by the El-Oued National Sanitation Office (l'Office National d'Assainissement d'El-Oued ONA). The principle of this system is to collect all the infiltration water under the agglomeration of El-Oued by pumping it into 58 drilled wells and to use as much of this water as possible on-site for the irrigation of green spaces, while the rest of this water is collected in the pumping station and piped over a distance of about 4200 m to the collector of the treated wastewater of the treatment plant to be treated by aerated lagoons then discharged in the discharging site (Chott Halloufa) [79].

Data Collection
The research covered 58 monitoring wells that are part of the vertical drainage system and are located in Oued Souf valley. A total of 55 drains are located in El Oued, while three drains are located northwest of Bayadha. Groundwater levels were measured in 2008,2009,2014,2016,2018, and 2021 by ANRH (Agence National des Roussources Hydrauliques-National Agency for Water Resources) and ONA (Office National de l'Assainissement-National Sanitation Office). The measurements of the groundwater level were carried out in May of each year to avoid seasonal fluctuations in the groundwater level.

Data Collection
The research covered 58 monitoring wells that are part of the vertical drainage system and are located in Oued Souf valley. A total of 55 drains are located in El Oued, while three drains are located northwest of Bayadha. Groundwater levels were measured in 2008,2009,2014,2016,2018, and 2021 by ANRH (Agence National des Roussources Hydrauliques-National Agency for Water Resources) and ONA (Office National de l'Assainissement-National Sanitation Office). The measurements of the groundwater level were carried out in May of each year to avoid seasonal fluctuations in the groundwater level.
The groundwater depth measurements were taken from the land surface using level probes and piezometers installed in each drain of the entire system. As all drains are equipped with pumps, the measurement was carried out before or after the pumping phase, leaving enough time for the groundwater level to return to its static state. In the meantime, the water flow and volumes were measured.
The geographical coordinates of these monitoring wells were determined using a Globe Positioning System (GPS). Climate data for the period 1978 to 2020 were obtained from the meteorological station of El Oued-Guemar airport. These data include monthly rainfall and monthly temperature data. The meteorological station is located in the coordinates Latitude: 33 • 30 N and Longitude 06 • 47 E with an altitude of 63 m.a.s.l.
On the other hand, the deep groundwater reservoirs data including the extracted volumes and their usage from the Complex terminal and Continental intercalary aquifers were obtained from the ANRH company.

Geostatistical Modelling
Geographical information systems (GIS) have increasingly facilitated environmental protection and resource planning recently has extended its utilizations in different fields such as geology, hydrology, environmental monitoring by integrating spatial data with other information for the generation of the spatial distribution of groundwater parameters and their evaluation using geostatistics [80][81][82][83][84][85][86]. Kriging is a stochastic interpolation approach for spatial surface prediction. Because it incorporates statistical models, it is adaptable and provides for the analysis of data spatial autocorrelation. The data are assumed to come from a stationary stochastic process in Kriging, and some methods necessitate that the data be regularly distributed. The process of kriging is separated into two parts: measuring the spatial structure of the data and producing a predicted surface. Kriging will utilize the fitted model from the variogram, the spatial data configuration, and the values of the measured sample points surrounding the prediction location to predict an unknown value for a specific location [36,87,88]. Since ordinary kriging is widely employed as a dependable estimate method among the many kinds of kriging, it was introduced in this research for producing spatial prediction maps [89][90][91].
The selection of the best-fitted semi-variogram models was done based on the mean error (ME) and root mean square standardized error (RMSSE) values. Furthermore, each one of the models is evaluated according to the precision of the predictions when RMSSE is close to unity and ME is minimized.
Several classes were accustomed to illustrating the structural spatial dependency and this is relying on nuggets' variance/sill ratio. Strong when the ratio is less than 25% when the ratio is between 25% and 75% it is moderate, and weak when the ratio is larger than 75%. All the geostatistical analyses were made using ArcGIS 10.4.1 software.

Clustering Analysis
Hierarchical clustering analysis (HCA) is an unsupervised grouping method that distinguishes variables based on their similarity. Euclidean distances were used to classify the parameters into initial clusters, while Ward's agglomeration method was used to link the resulting initial clusters [44,[92][93][94][95]. In this research, this technique was used to group the groundwater depths of the phreatic groundwater aquifer for elaborating the differences according to their similarity using Originpro 2021 software.

Hierarchical Clustering of Groundwater Levels
Hierarchical clustering analysis (Q-mode) has been applied to classify the monitoring wells (drains) from 2008, 2009, 2014, 2016, 2018, and 2021 according to their fluctuation variability using Ward's linkage technique and Euclidean distance. It is an important statistical tool that divides each case into a distinct group or cluster and links each cluster until only one is remaining. The application of HCA has generated a dendrogram of spatial hierarchical clustering analysis. However, four groups have been created ( Figure 5). The mean of each group was the distinguishing factor that contributes to the classification of these groups (it is increasing from the first group to the fourth group ascendingly), ( Table 2).
The Spatio-temporal variation in the wells of different clusters in Figure 6 revealed that 22 wells in cluster 1 showed similar fluctuation variability, with the mean fluctuation ranging between 2.64 in 2009 to 11.44 m below ground level (mbgl) (2018), and the mean fluctuation of this group represent the shallowest level (5.91 mbgl). The clustered wells of the first group are concentrated in the center and expand to the southeast of the study area, while only one well is located in the northeast of the stud area. The second cluster consists of eighteen wells whose mean fluctuation varied from 4.33 in 2009 to 14.30 mbgl (2018). This group has a mean fluctuation that can be considered deeper than the first cluster. This clustered group is spread over the study area (the southern east, the center to the northwest, and in one well located in the northeast of the study area). The Spatio-temporal variation in the wells of different clusters in Figure 6 revealed that 22 wells in cluster 1 showed similar fluctuation variability, with the mean fluctuation ranging between 2.64 in 2009 to 11.44 m below ground level (mbgl) (2018), and the mean fluctuation of this group represent the shallowest level (5.91 mbgl). The clustered wells of the first group are concentrated in the center and expand to the southeast of the study area, while only one well is located in the northeast of the stud area. The second cluster consists of eighteen wells whose mean fluctuation varied from 4.33 in 2009 to 14.30 mbgl (2018). This group has a mean fluctuation that can be considered deeper than the first cluster. This clustered group is spread over the study area (the southern east, the center to the northwest, and in one well located in the northeast of the study area).   The third cluster entails the group of wells that have a depth level that fluctuates from 6.90 (2009) to 10.96 mbgl (2018), while the total mean fluctuation of this cluster was 9.31 mbgl. The wells of this group are graduating spatially from the center to the northest of the study area, except for one well located in the northeast. The deepest fluctuation means were represented in the fourth cluster, whereas the mean fluctuations were ranging from 11.64 (2009) to 15.61 mbgl (2018). On the other hand, the total mean fluctuation of this group was 12.98 mbgl. The wells of this cluster are located separately in the southwest of the central part to the southwest of the study area. The third and fourth groups consist of the same number of wells as the second group.

Geostatistical Modelling and Groundwater Level Fluctuations
The first step in statistical analysis entails creating frequency tables, histograms, and basic statistical calculations. Summary statistics, including the mean, coefficient of variation, minimum, maximum, skewness, and kurtosis values, are given in Table 3. Based on the histograms, kurtosis, and skewness values, 2008, 2014, 2016, and 2021 can be considered approximately normal distributed Figure 7. A logarithmic transformation has been applied to the data from 2009 and 2018 before performing a geostatistical analysis. The anisotropy of the phreatic aquifer was verified using semivariance analysis and ordinary kriging to understand more about the upwelling of groundwater in the study area. For The third cluster entails the group of wells that have a depth level that fluctuates from 6.90 (2009) to 10.96 mbgl (2018), while the total mean fluctuation of this cluster was 9.31 mbgl. The wells of this group are graduating spatially from the center to the northest of the study area, except for one well located in the northeast. The deepest fluctuation means were represented in the fourth cluster, whereas the mean fluctuations were ranging from 11.64 (2009) to 15.61 mbgl (2018). On the other hand, the total mean fluctuation of this group was 12.98 mbgl. The wells of this cluster are located separately in the southwest of the central part to the southwest of the study area. The third and fourth groups consist of the same number of wells as the second group.

Geostatistical Modelling and Groundwater Level Fluctuations
The first step in statistical analysis entails creating frequency tables, histograms, and basic statistical calculations. Summary statistics, including the mean, coefficient of variation, minimum, maximum, skewness, and kurtosis values, are given in Table 3    There is a slight degree of anisotropy, which means that the created experiment semivariograms for the different directions have some differences in their sills and infl ence ranges. As a result, the discovered anisotropies were overlooked, and isotropic mo els were used instead. At the same time, mean error (ME) and root mean square standar ized error (RMSSE) values have been considered for the selection of the best fitted sem variogram models. Whereas, the evaluation of the selected model was conducted bas There is a slight degree of anisotropy, which means that the created experimental semivariograms for the different directions have some differences in their sills and influence ranges. As a result, the discovered anisotropies were overlooked, and isotropic models were used instead. At the same time, mean error (ME) and root mean square standardized error (RMSSE) values have been considered for the selection of the best fitted semi-variogram models. Whereas, the evaluation of the selected model was conducted based on the predictions, especially, when RMSSE is close to unity and ME is very low. The adopted model parameters for the examined annual groundwater levels are summarized in Tables 4 and 5. Hypothetical models were then fitted to the more typical non-directional semivariograms as illustrated in Figure 8. on the predictions, especially, when RMSSE is close to unity and ME is very low. The adopted model parameters for the examined annual groundwater levels are summarized in Tables 4 and 5. Hypothetical models were then fitted to the more typical non-directional semivariograms as illustrated in Figure 8.   Tables 4 and 5. The range of influence and the sill values of the adopted models reveals some significant differences. The range of influence is the distance within which the groundwater level values are spatially dependent. The minimum spatial dependence (702.55 m) was encountered in 2018, and the maximum one (5061.97 m) was marked in 2021. The presence of the nugget effect indicates inherited variability or variability across shorter distances than the spacing between observation wells, as well as other unaccounted for recording errors. Based on [96] structural spatial dependence, the aquifer structure exhibits strong spatial dependence that remained almost constant over the years except in 2009 and 2021 when the spatial structures were moderate. Very low nugget effects show that the groundwater level fluctuation is severely time-correlated and depicts a strong temporal structure.
As observed by the generated prediction maps of ordinary kriging for the water table status in 2008, several patterns were detected in the study area. Variations extend from the northwest to the southeast of the study area, where groundwater depths are almost 3-4 m and upwelling near the surface (1-3 m) till up to the surface (−0.3 in D20), mainly in the center El-Oued municipality. In the southwest of the research region, where the depth to groundwater was more than 6 mbgl and dropping (going deeper). Another observation was that the depth was decreasing by around 2 m near two drains, D05 (northeast) and D45 (south). In 2008 the beginning of the vertical drainage system functioning occurred, the groundwater levels were very shallow and, in some drains, the groundwater level was up to the surface such as D20 with −0.3 m below the ground level (mbgl). Besides this drain, other drains recorded very low levels that are near to the surface such as D13, D17, D18,32, and D33 where the levels ranged from 0.43 to 1.88 m below the ground surface. On average, the depth of groundwater fluctuated by about 5.42 mbgl during 2008. Similarly, during 2009 the groundwater level was continuing its upwelling from the northwest part to the southeast, where the shallowest level was recorded in D20 with −0.4 mbgl, while the deepest level was noted in the same drain as 2008, which is D52 with 15.17 mbgl. D13, D17, D18, 32, and D33 had the same trend as in 2008, where their levels varied from 0.41 to 1.75 mbgl. The mean groundwater level in 2009 was 5.06 mbgl which was less than in 2008 indicating increasing in the upwelling average in the study area during 2009 as is exposed in Figure 9A,B.
The groundwater level started its significant decline from 2014 and 2016, respectively. However, the shallowest groundwater levels were observed as usual in a gradual upwelling state from the central northwest part to the southeast of the study area to reach the minimum value that was marked in D13 with 0.6 mbgl in 2014 and 1 mbgl in 2016. The deepest levels were noted in the southwest of the study area, where the levels were more than 8 mbgl and increased deeply (significant decline) to 18.23 and 19.53 mbgl marked in D52 in 2014 and 2016, respectively, as is illustrated in Figure 9C,D. Fluctuation rates ranged from 7.97 mbgl in 2014 and 7.81 mbgl in 2016, Table 3. Moreover, the minimum groundwater level in 2018 has shifted to D32 with 5.1 mbgl and is located in the southeast part of the city, while the deepest also has shifted to D37 with 28.  The groundwater level started its significant decline from 2014 and 2016, respectively. However, the shallowest groundwater levels were observed as usual in a gradual upwelling state from the central northwest part to the southeast of the study area to reach the minimum value that was marked in D13 with 0.6 mbgl in 2014 and 1 mbgl in 2016. The deepest levels were noted in the southwest of the study area, where the levels were more than 8 mbgl and increased deeply (significant decline) to 18.23 and 19.53 mbgl marked in D52 in 2014 and 2016, respectively, as is illustrated in Figure 9     Suddenly, the groundwater level status in 2021 revealed a considerable change that happened in the study area in terms of groundwater level fluctuation, decline, and rising rates. The deepest groundwater level can be highlighted in the center of the study area where the deepest level was reached to a maximum in D34 (15 mbgl), and in the northwest of the study area where the deepest level was registered in D10 (15 mbgl). On the other hand, other drains recorded deep values such as D52, D53, and D54 in Figure 11F and Table 4. Meanwhile, the minimum groundwater level appeared in the central east of the study area registered in D32 (2.1 mbgl). Although a significant decline during the previous years, the depth of groundwater fluctuated by around 8.87 mbgl in 2021, which indicates a significant rise compared especially to 2018 with 3.9 m on average as is shown in Figure  10.
Since the water table in the study area is determined by several factors. Therefore, there are several explanations for the decline of groundwater and its fluctuations during the study period. However, the independence of groundwater use in the study area from the phreatic aquifer has promoted the decline in water returning to the phreatic aquifer. Similarly, 94.7% of the groundwater withdrawn from the entire aquifer system is used for various purposes. In addition, 93.5% of the water comes from the complex terminal aquifer (72.2% from the Mioleocene, 21.2% from the Pontian, 0.1% from the Eocene), and 1.2% from the continental intercalary (Barrimian). The rest (5.3%) is drained by the vertical drainage system, which drains the rest of the phreatic groundwater (Quaternary). On the other hand, a large part of groundwater resources is used for drinking water supply (40.4%) and irrigation (51.3%), and the percentage may be higher, while 8.3% is drained Suddenly, the groundwater level status in 2021 revealed a considerable change that happened in the study area in terms of groundwater level fluctuation, decline, and rising rates. The deepest groundwater level can be highlighted in the center of the study area where the deepest level was reached to a maximum in D34 (15 mbgl), and in the northwest of the study area where the deepest level was registered in D10 (15 mbgl). On the other hand, other drains recorded deep values such as D52, D53, and D54 in Figure 11F and Table 4. Meanwhile, the minimum groundwater level appeared in the central east of the study area registered in D32 (2.1 mbgl). Although a significant decline during the previous years, the depth of groundwater fluctuated by around 8.87 mbgl in 2021, which indicates a significant rise compared especially to 2018 with 3.9 m on average as is shown in Figure 10.
Since the water table in the study area is determined by several factors. Therefore, there are several explanations for the decline of groundwater and its fluctuations during the study period. However, the independence of groundwater use in the study area from the phreatic aquifer has promoted the decline in water returning to the phreatic aquifer. Similarly, 94.7% of the groundwater withdrawn from the entire aquifer system is used for various purposes. In addition, 93.5% of the water comes from the complex terminal aquifer (72.2% from the Mioleocene, 21.2% from the Pontian, 0.1% from the Eocene), and 1.2% from the continental intercalary (Barrimian). The rest (5.3%) is drained by the vertical drainage system, which drains the rest of the phreatic groundwater (Quaternary). On the other hand, a large part of groundwater resources is used for drinking water supply (40.4%) and irrigation (51.3%), and the percentage may be higher, while 8.3% is drained by the drainage system ( Figure 12A,B). The rapid growth of agriculture and its development around the city of El Oued Souf also contribute to the rapid decline of the phreatic aquifer.
Water 2022, 14,1415 19 of 25 by the drainage system ( Figure 12A,B). The rapid growth of agriculture and its development around the city of El Oued Souf also contribute to the rapid decline of the phreatic aquifer. At the same time, the primary results of the drainage vertical system have indicated bad results (upwelled phreatic groundwater level with a rising rate of −0.36 m) between 2008-2009 then these bad results were shifted to be better over the study period. In terms of groundwater level, due to the changing of the sewage system (from septic tanks to vertical drainage system) of El Oued the groundwater started to decline over time and reduced the contamination level significantly in the phreatic groundwater aquifer. Also connecting most of the agglomerations to the sewage network and stopping wastewater disposal into the environment with huge quantities which were returning to the phreatic aquifer previously. About, 27.5 million cubic meters per year are discharged through the drainage system to the treatment station as schematized in Figure 13, with 7.74 l/s as discharge flow. However, the results of 2021 can generate several questions about the groundwater level and can increase fears about returning to the groundwater upwelling catastrophe again. Numerous influences might be leading the study area to another rising groundwater level such as the inefficient drinking water supply network, which is struggling with considerable water leakages, for example, the unaccounted withdrawals from connections (lack of meters or failure of the measurement system), and the illegal withdrawals (such as illegal connections, manipulation of meters). At the same time, the primary results of the drainage vertical system have indicated bad results (upwelled phreatic groundwater level with a rising rate of −0.36 m) between 2008-2009 then these bad results were shifted to be better over the study period. In terms of groundwater level, due to the changing of the sewage system (from septic tanks to vertical drainage system) of El Oued the groundwater started to decline over time and reduced the contamination level significantly in the phreatic groundwater aquifer. Also connecting most of the agglomerations to the sewage network and stopping wastewater disposal into the environment with huge quantities which were returning to the phreatic aquifer previously. About, 27.5 million cubic meters per year are discharged through the drainage system to the treatment station as schematized in Figure 13, with 7.74 L/s as discharge flow. However, the results of 2021 can generate several questions about the groundwater level and can increase fears about returning to the groundwater upwelling catastrophe again. Numerous influences might be leading the study area to another rising groundwater level such as the inefficient drinking water supply network, which is struggling with considerable water leakages, for example, the unaccounted withdrawals from connections (lack of meters or failure of the measurement system), and the illegal withdrawals (such as illegal connections, manipulation of meters).
In terms of the drinking water supply network losses, we can highlight that the network of Oued Souf is struggling with several problems such as defective joints of pipe connections, and conduits (steel, asbestos, concrete, polyvinyl chloride) and valves, leaking tanks, etc. Similarly, the bad connections with other supply systems and finally other nonremunerated withdrawals such as taking for firefighting, withdrawals for inspection work, maintenance of the network, rinsing of pipes, etc. In addition to this, the existence some industries that consume a considerable amount of water such as the gypsum production industry. On the other hand, some problems in the drainage system can generate a rerising in the groundwater level such as the power outage break which affects the pumping process, and also the maintenance of the system. During the measurement of 2021, several drains had not functioned since 2018. Another explanation can be mentioned here but we cannot adopt it since it is out of the scope of this research, which is the non-homogeneity of the unconfined aquifer (phreatic aquifer) and its intercalation by clay lenses at a shallow depth that act as a substrate to support groundwater to rise [25]. Furthermore, the responsible purification plant for purifying the wastewater collected from the study area is suffering from different problems such as the breakdown of some pumps in different purification stages, especially in the Desander and this is mainly because of the high level of wastewater in the grit trap causes the two suction pumps to break down and the third insufficient, which creates considerable pressure on the plant and consequently reducing the required purified of the drained wastewater volumes. In terms of the drinking water supply network losses, we can highlight that the network of Oued Souf is struggling with several problems such as defective joints of pipe connections, and conduits (steel, asbestos, concrete, polyvinyl chloride) and valves, leaking tanks, etc. Similarly, the bad connections with other supply systems and finally other non-remunerated withdrawals such as taking for firefighting, withdrawals for inspection work, maintenance of the network, rinsing of pipes, etc. In addition to this, the existence some industries that consume a considerable amount of water such as the gypsum production industry. On the other hand, some problems in the drainage system can generate a rerising in the groundwater level such as the power outage break which affects the pumping process, and also the maintenance of the system. During the measurement of 2021, several drains had not functioned since 2018. Another explanation can be mentioned here but we cannot adopt it since it is out of the scope of this research, which is the nonhomogeneity of the unconfined aquifer (phreatic aquifer) and its intercalation by clay lenses at a shallow depth that act as a substrate to support groundwater to rise [25]. Furthermore, the responsible purification plant for purifying the wastewater collected from the study area is suffering from different problems such as the breakdown of some pumps in different purification stages, especially in the Desander and this is mainly because of the high level of wastewater in the grit trap causes the two suction pumps to break down and the third insufficient, which creates considerable pressure on the plant and consequently reducing the required purified of the drained wastewater volumes.

Conclusions
Under the arid and semi-arid climate, groundwater is the only reachable resource, especially in Oued Souf valley. However, the decline and rise in groundwater levels in the study area were attributed primarily to human activities.
The availability of huge groundwater reservoirs in Oued Souf turned out to be a disadvantage instead of an advantage for the region due to the failure of draining large quantities that have been originally pumped from deep aquifers and consequently returned to the unconfined aquifer (phreatic aquifer). This is affecting the non-stability of the phreatic

Conclusions
Under the arid and semi-arid climate, groundwater is the only reachable resource, especially in Oued Souf valley. However, the decline and rise in groundwater levels in the study area were attributed primarily to human activities.
The availability of huge groundwater reservoirs in Oued Souf turned out to be a disadvantage instead of an advantage for the region due to the failure of draining large quantities that have been originally pumped from deep aquifers and consequently returned to the unconfined aquifer (phreatic aquifer). This is affecting the non-stability of the phreatic groundwater level and this is reflecting directly on the leaving environment in the Oued Souf region qualitatively and quantitatively.
The used methodology for studying the non-stability of the phreatic groundwater aquifer indicated that the anthropogenic activities are the direct reasons for the rising of groundwater. Despite the vertical drainage system reflecting a significant decline in groundwater from 2009 to 2018, other rerising phenomena might be threatening the region in the near future and this is what was indicated in the 2021 groundwater level data.
Cluster analysis has indicated that the most vulnerable area is represented by group 1 which has a shallow depth (average mean of 5.91 mbgl). The non-vulnerable area is totally separated from the other clusters and represented by the fourth cluster with a deep depth (12.89 mbgl) that signified no threat of rising groundwater levels.
Geostatistical analysis performed by kriging revealed that a certain heterogeneity appeared in the 2009 and 2021 data and were represented by moderate spatial dependence as opposed to 2008, 2014, 2016, and 2018 which have a strong spatial dependence. The used interpolation technique indicated the most vulnerable area over the years located and extended from the center to the east of the study area.
Several factors affect the stability of the phreatic groundwater level and the vertical drainage system performance. However, all the most crucial factors are of anthropogenic origins such as the overexploitation from deep groundwater reservoirs from the complex terminal and continental intercalary (in drinking and irrigation purposes) and even the illegal use of the phreatic groundwater with important quantities for irrigation and illegal industries although its contamination and unsuitability for the utilization.
The results of this research can be highly important to draw the attention and increase the awareness of all water experts locally and globally about the consequences of inadequate water management and negligence in the implementation of the integrated management of water resources.
For future research work, deep and machine learning techniques (Support Vector Machines (SVMs), Gradient Boosting Machine, and M5 and M5-cubist Convolutional Neural Networks) will be used, as well as several modified DRASTIC methods, to estimate the degree of pollution and place the areas most vulnerable to pollution due to the heavy use of underground water for irrigation and industry, which makes the first step towards developing a model for some desert areas.