Hybridization of DRASTIC Method to Assess Future GroundWater Vulnerability Scenarios: Case of the Tebessa-Morsott Alluvial Aquifer (Northeastern Algeria)

Abstract

In this study, a new approach integrating a groundwater vulnerability method and a numerical model for predicting groundwater resource sustainability under actual and future conditions of exploitation (2010–2030) is proposed in the semi-arid region of the Tebessa-Morsott alluvial aquifer (northeastern Algeria). The groundwater vulnerability method-based DRASTIC model was used to evaluate and delineate the vulnerable areas using a GIS technique. The MODFLOW code, on the other hand, was used to calculate the dynamics of groundwater level under actual and future conditions of exploitation considering two scenarios. The results of the application of the DRASTIC method to the reference year conditions (year 2010) showed that the high and average vulnerability classes covered a wide zone of the study area, about 97%. These results were validated based on the nitrate concentration values (R² = 0.955). However, the results for predicting future groundwater vulnerability showed that groundwater vulnerability variation over time (period 2010–2030) was closely related to groundwater depth variation caused by the pumping rate, since the decreases in the piezometric level produce a worsening of groundwater vulnerability. To achieve better groundwater management, an experimental site for artificial recharge supplemented by hydro-chemical monitoring of the groundwater could be an effective remediation strategy.

1. Introduction

Water, a rare and vulnerable resource, is the main environmental issue all over the world. The scarcity and continuous deterioration of this resource has triggered researchers to adopt a huge variety of modern investigation techniques for its protection to avoid uncensored overexploitation and its anarchic management [1]. Groundwater is essential for agricultural and industrial developments; however, the rapid extension of irrigated areas, along with population growth and industrial development, has stressed the quantitative and qualitative status of this precious resource [2]. Due to the growing issue related to water resources, people in developed countries have finally begun to realize that it is necessary to change the bad practices surrounding this resource and to implement a coherent policy for rational groundwater management [3], as also remarked on by the SDGs proposed by the United Nations. In particular, the sixth SDG states the importance of ensuring “Clean water and sanitation”, highlighting the primary importance of protecting and restoring water-related ecosystems. The lack of water in arid and semi-arid regions is becoming a major constraint for local and regional sustainable development. Moreover, arid and semi-arid regions are also characterized by significant spatiotemporal imbalances between water demand and water supply, making the implementation of new strategies extremely necessary for the implementation of integrated and sustainable groundwater resources management [4]. Today, vulnerability maps are widely used for groundwater management and protection [5,6]. These maps are generally established based on geological, pedological, hydrological and hydrogeological factors. Moreover, is possible to distinguish two types of vulnerability: (i) intrinsic vulnerability linked to the geological environment and local climatic conditions, and (ii) specific vulnerability linked to the properties of the contaminant/pollutant.

The first groundwater vulnerability maps were established by researchers in the early 1970s [7]. Afterwards, they were then positively tested in other countries around the world. The principle of their elaboration consisted in carrying out the synthesis of a few parameters influencing the aquifer’s vulnerability. In general, groundwater vulnerability maps can be grouped into four major categories: (1) overlay and index methods, (2) statistical methods, (3) modeling-based methods [8,9,10,11], and (4) recently introduced hybrid methodologies.

The latter represent a combination of the above-mentioned categories aimed to increase the reliability of groundwater vulnerability results. They involve the hybridization of rating methodologies with statistical and mathematical models aimed to increase the reliability of results [12,13,14]. However, as a drawback, these new methods become extremely site-specific and data-dependent. Generally, rating and overlay methods are the most widely utilized, and several methods have been proposed and tested in the last 30 years, such as: DRASTIC [15], GOD [16], AVI [17], PRAST [18], ISIS [19], SINTACS [20], EPIK [21], SI [22], RISKE [23], PI [24], and GALDIT [25]. Among these, the DRASTIC method, along with GOD, has been applied in areas with different climates and hydrogeologic characteristics with good results. Moreover, its results are easily understood by water managers. However, these methodologies only offer a static representation of vulnerability, making their hybridization mandatory to ensure proper management as well as enabling short-term vulnerability forecasting.

Several modifications have been proposed in the last years: Busico et al. [14] achieved the optimization of rating methods combining statistical analysis with parameter objectivation; Khosravi et al. [26] enhanced DRASTIC performance using artificial intelligence; and Duc-Vu et al. [27] combined DRASTIC and MODFLOW to increase vulnerability performance. A complete review of all methodologies and modifications is available in Machiwal et al. [1] and in Taghavi et al. [28].

On the other hand, the hydrogeological modeling approach remains an important tool to grasp the objectives of groundwater management in different sectors: irrigation, drinking water supply and industry. Numerical modeling is a very practical tool and is often used by researchers for the study of groundwater flows [29] and contaminant flows in aquifers [30]. It has been used in groundwater vulnerability assessment but is usually confined to a small area since it requires an over-parametrization of the aquifer system. However, it is a powerful tool allowing the simulation of different flows and groundwater contamination scenarios (actual and future states of the model) [31]. However, the implementation of a model requires an understanding of the aquifer’s behavior and a prior geometric description of the hydrogeological system [32]. In recent years, Algeria has considered a new policy for the management of water resources.

This policy aims to protect water resources in their deposits, before they are contaminated, through legislative tools that respond to socioeconomic and demographic constraints. This research concerns the Tebessa-Morsott alluvial aquifer located in a semi-arid region, a stressed aquifer that, today, supplies several urban centers and cities with a population of around 150,000 inhabitants. Currently, the aquifer is threatened by many sources of pollution, in particular, discharges from domestic water and small industries in the region, which discharge into the aquifer without prior treatment and without any measure to protect the water resource.

The agricultural expansion experienced by the region after the application of the National Program for the Development of Agriculture (NPDA) in 2000 and that of the National Fund for the Regulation of Agricultural Development (NFRAD) in 2001 has led to an increase in water demand, exposing this resource to both groundwater depletion and quality degradation. To address these threats, the present study aims to investigate the hydrodynamic behavior of the aquifer by combining hydrogeological modeling and groundwater vulnerability to propose measures and ensure groundwater management of this resource for the near future (2020–2030). To assess aquifer vulnerability in actual and future conditions, we propose a new approach by coupling the groundwater flow model using the MODFLOW model with the vulnerability assessment method (DRASTIC model) using GIS technique. The groundwater flow model will be used to calculate the dynamics of groundwater level under actual and future climate and exploitation conditions in order to delineate the actual and predict the future status of groundwater vulnerability in the aquifer.

2. Materials and Methods

2.1. Description of the Study Area

The Tebessa-Morsott plain is located in the northeastern part of Algeria, 650 kilometers southeast of the capital city, Algiers. It occupies a part of the large Medjerda watershed—classified as basin number 12 according to the orohydrographic delimitation proposed by the National Water Resources Agency (ANRH)—and the entire Ksob wadi sub-basin. It corresponds to a large, closed depression extending along the NW–SE direction [33] and covering a surface of 600 km² [34]. It lies between longitudes of 7°64′ E and 8°18′ E and between latitudes of 35°22′ N and 35°46′ N (Figure 1).

Based on field observations, topographic map data and satellite images (downloaded from https://earthexplorer.usgs.gov/ (accessed on 12 August 2022) Digital Elevation, SRTM 1 Arc-Second Global; the model domain was processed using ArcGIS 10.8), this region is surrounded by mountains ranging in altitude from 1713 m to 780 m a.s.l (Figure 1). We cited Zitouna (1324 m) and Edyr (1471 m) mountains in the north, Doukkane (1713 m) and Anoual (1545 m) mountains in the south, Belkefif (1338 m) mountain in the west and Djebissa (1120 m) mountain in the southeast. In contrast, the plain is characterized by a relatively flat surface with elevations ranging from 950 m to 700 m a.s.l at the south, east and west borders from the center and the north of the study area, forming a slight slope (average 0.25%) that favors runoff towards the depressions and the wadis. This morphology allows the formation of an important alluvial aquifer (Figure 1). The plain is drained by an important permanent wadi named El Ksob wadi, which follows a SE–NW direction, and its two tributaries, El Kebir wadi and Chabro wadi, which cross it for a length of about 60 km [35]. The study area is characterized by a semi-arid climate with wet winter months (November to May) and dry summer months (June to August). The annual average precipitation is about 380 mm, and the annual average temperature is 17 °C. The evaporation rate is about 650 mm/year. Other climatic parameters, such as wind speed and evaporation, show remarkable seasonal variation. Specifically, between winter (December) and summer (August), the wind speed varies from 3.2 m/s to 1.9 m/s, and the humidity rate varies from 18% to 76%.

2.2. Geological and Geophysical Properties

Geologically, the study area is very heterogeneous (Figure S1). It is part of the autochthonous North Auresian structure of the Saharan atlas [33,36]. It consists essentially of the following formations: (i) The Triassic outcrops around the Bekkaria region (Djebissa mountain) in the southeast, the Bir Dheb region (Belkefif mountain) in the west, and the north part of the Morsott region; this formation contributes to the high salinity of the aquifer [35]; (ii) the Cretaceous carbonate formations represented by the significant calcareous–marly and marl formations, which are very remarkable on the edges of the plain.

Finally, (iii) important Mio-Plio-Quaternary alluvial deposits (conglomerates, gravels, sand and sandstones) occupy the rest of the plain. These Mio-Plio-Quaternary formations constitute the main reservoir of the plain [37]. In addition, the borehole logs carried out in this plain indicate a thick permeable sedimentation (aquifer) that varies between 50 m and 150 m [38,39]. To determine the aquifer’s geometry and depth and the nature and form of the substratum, an analysis of a geophysical survey (Figure 2), carried out by the General Company of Geophysics [40], was necessary. This analysis shows that the aquifer presents in the form of an elongated gutter (synclinal form) in a SE–NW direction. However, there are: (i) a relatively resistant level (alluvial deposits) of Mio-Plio-Quaternary, with resistivity values varying between 10 ohm.m and 100 ohm.m, and (ii) marl formations (substratum, conductor) of Dano-Montian (Paleogene age), with resistivity values less than 10 ohm.m (varying between 2 ohm.m and 8 ohm.m). All these formations rest on Cretaceous limestone. The thickness of the alluvial deposits is very important, reaching 300 m in the central part of the plain (Ain Chabro region) but not exceeding 50 m in the northern part of the plain (Morsott region). Geological and geophysical studies revealed the existence of the Mio-Plio-Quaternary alluvial aquifer (unconfined aquifer) formed by alluvial deposits (conglomerates, gravels, sand and sandstones). Its substratum is made of marls of the Paleogene age (Dano-Montian). The thickness of the aquifer varies between 150 and 300 m, constituting an aquifer favorable for drilling.

2.3. Hydrogeological Setting

The Tebessa-Morsott groundwater aquifer is mainly recharged by direct precipitation, by later inflow from mountains around the plain and from leakage through faults in the deep aquifer (Cretaceous formations) [35]. The general flow of groundwater converges to the north, parallel to the current course of the wadis (Figures S2 and S3). The mountains of the south and southwest regions that are in contact with this aquifer constitute the permeable limits, which play a role in the recharge of the aquifer. The pumping tests of several water wells (boreholes) have been investigated to calculate the hydrodynamic parameters of the aquifer [38]. The hydraulic conductivity value was estimated at about 3 × 10⁻⁴ m/s, and the transmissivity values vary from 10⁻³ to 5 × 10⁻² m²/s. In 2010 in the Ain Chabro region, the measured water discharge in 22 water wells was 4700 m³/day, confirming the importance of the aquifer in terms of productivity.

2.4. The Conceptual Framework

The aim of this study is to propose a new framework to predict actual and future groundwater vulnerability under the climate conditions of semi-arid regions using the overlay–index method coupled with the groundwater flow simulation model. The work has been structured following a well-defined procedure to assess actual and future groundwater vulnerability. It has been divided into three steps: (1) The first step will provide actual groundwater vulnerability mapping through application of the DRASTIC method; (2) the second step will be the realization of the groundwater flow simulation model of the study area, under a steady state and a transient regime; (3) the final step will consider the modeling of future groundwater vulnerability using the DRASTIC approach coupled with the groundwater flow simulation model. Figure 3 shows the procedure adopted for the groundwater vulnerability mapping in this study.

2.5. Groundwater Vulnerability Using DRASTIC Method

The DRASTIC method was proposed by the EPA (Environmental Protection Agency) in the United States in 1985 and was developed by Aller et al. [15]. The methodology assesses the vertical vulnerability based on seven parameters: depth of water (D), net recharge (R), aquifer media (A), soil media (S), topography (T), the impact of the vadose zone (I) and hydraulic conductivity (C). A weight was assigned to each parameter, the values of which vary between 1 and 5 depending on the degree of influence. The highest weight (5) was attributed to the depth of water and the impact of the vadose zone. The lowest weights (2 and 1) were attributed, respectively, to the soil media and the topography. Each parameter was mapped with an index, called a “rating”, typically ranging from 1 to 10, where the lowest rate represents the conditions of a lower vulnerability to contamination. The DRASTIC vulnerability index was calculated as follows:

DRASTIC Index = Dw Dr + Rw Rr + Aw Ar + Sw Sr + Tw Tr + Iw Ir + Cw Cr

(1)

where: D, R, A, S, T, I and C are the abbreviations of the seven parameters, and “r” and “w” are the corresponding rating and weights, respectively.

Table 1. Evaluation classes of vulnerability degree using the DRASTIC model.

DRASTIC Vulnerability	Vulnerability Index
Low	23–100
Average	101–140
High	141–200
Very high	200–226

shows the final DRASTIC classifications according to [41], while Table S1 shows the ratings and the weighting values assigned to each parameter in the DRASTIC model. The obtained map presents the relative vulnerability degree of a sector (of the study area). The potential pollution increases in the same direction as the DRASTIC vulnerability index, which represents the evolution of the contamination risk of an aquifer formation.

2.6. DRASTIC Validation

To assess the reliability of the DRASTIC method, nitrate concentrations were used for validation. The presence of high NO₃ concentrations in water can cause health problems, such as stomach cancer for adults and methemoglobinemia for children [42]. NO₃ concentrations below 50 mg/L are usually acceptable to consumers [43,44]. According to [45,46], the NO₃ concentrations of the study area were divided into two classes. Specifically, all water samples with NO₃ concentrations below 50 mg/L were considered suitable for drinking, while those with a concentration higher than 50 mg/L were labeled as unsuitable for potable use.

2.7. Groundwater Modelling

In this work, the MODFLOW [47] code was chosen to simulate the groundwater flow of the Tebessa-Morsott alluvial aquifer system. The MODFLOW code allows the calculation of the hydraulic head on each mesh of a three-dimensional domain by the finite difference method [48], using seeking to solve the diffusivity equation with partial derivatives of flows by a combination of the continuity equation and Darcy’s law. The representative elementary volume adapted to the flow in the aquifer must consider the flow domain over the entire wetted height between the coasts z1 and z2. Level z1 represents the impermeable substratum of the aquifer, and level z2 represents either the impermeable cover of a confined aquifer or the piezometric surface of an unconfined aquifer (as in our case), where the coast is identified as the head h. The governing equations for 2D or 3D groundwater flow are based on water mass balance and Darcy’s law [49]:

$$\left[\frac{\partial }{\partial x}\left(K_x\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_y\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_z\frac{\partial h}{\partial z}\right)\right]=S_s\frac{\partial h}{\partial t}+q_s$$

(2)

By integrating and taking into account: ∂h/∂z = 0 (Dupuit’s hypothesis):

$${{\displaystyle \int }}_{z1}^{z2}\frac{\partial }{\partial x}(K_x\frac{\partial h}{\partial x})dz+{{\displaystyle \int }}_{z1}^{z2}\frac{\partial }{\partial y}\left(K_y\frac{\partial h}{\partial y}\right)dz={{\displaystyle \int }}_{z1}^{z2}{S}_s\frac{\partial h}{\partial t}dz+{{\displaystyle \int }}_{z1}^{z2}qdz$$

(3)

Again, assuming that z1 and z2 vary in relation to x and y:

$$\frac{\partial }{\partial x}(({{\displaystyle \int }}_{z1}^{z2}{K}_{X}dx)\frac{\partial }{\partial x})+\frac{\partial }{\partial y}(({{\displaystyle \int }}_{z1}^{z2}{K}_{y}dy)\frac{\partial }{\partial y})=({{\displaystyle \int }}_{z1}^{z2}{s}_{s}dz)+({{\displaystyle \int }}_{z1}^{z2}qdz)$$

(4)

In this work, the modeling steps were:

1. Mesh and geometric criteria: the mesh is an important step that should be adapted to the problem. The size of the modeled domain depends on several factors: (i) The desired precision in the calculation; (ii) the number and distance of water points (boreholes and wells); and (iii) the shape of the piezometric surface in relation to the limit of the aquifer (modeled domain).

2. Boundary conditions: the different types of boundary conditions of the model are inspired by the geology of the region. They make it possible to delimit the domain within which the model will be responsible for calculating the hydraulic head of each mesh from the initial loads introduced as input. It is important to ensure that the boundaries imposed on the model create the same effect as its natural boundaries: (i) The model is limited vertically by a substratum—this is a zero-flow limit; (ii) the lateral limits are also considered as zero-flow (impermeable) limits; and (iii) at the surface, the water table is free, and its load is expressed by an infiltration rate—this is a condition of distributed flow over the surface of the aquifer roof.

3. Results and Discussions

3.1. Model Geometry and Boundary Conditions

The Tebessa-Morsott alluvial aquifer system was conceptualized to develop a flow model in a steady state and in a transient regime using two different development scenarios.

The aquifer was discretized into 70 × 82 square meshes with 500 m sides. The model was made up of 1200 active meshes, i.e., an area of 600 km² (Figure 4).

The piezometric map from May 2010 [33] was used as the reference piezometric state in the numerical model (Figure S2). This piezometric map was utilized as a top layer in the numerical model and for the calibration of the steady-state regime. However, the bedrock and the piezometric maps were exploited to calculate the saturated aquifer thickness map and to develop the aquifer geometry model.

Using the Thornthwaite’s method, the recharge rate was evaluated between 2 cm/year and 18 cm/year [33].

Based on the pumping tests in twenty-two boreholes capturing the aquifer, the permeability values varied between 2 × 10⁻⁵ m/s and 3 × 10⁻⁴ m/s [36]. Moreover, the total porosity ranged between 28% and 80%.

Hydrogeological conceptual models are collections of hypotheses describing the understanding of groundwater systems, and they are considered one of the major sources of uncertainty in groundwater flow and transport modeling. According to basic hydrogeological information, the following boundary conditions were set in the steady-state groundwater model. The aquifer was considered as unconfined with a heterogeneous layer. The entire study area received recharge by rainfall (flow boundary). The wadis of the region (El Ksob wadi, El kebir wadi and Chabro wadi) were defined as the second flow boundary. In some lateral boundaries, the boundaries were closed and considered watertight (no-flow boundary). Conversely, some boundaries were simulated as open flow boundaries, where hydraulic gradients allow flow across these boundaries (fractured limestones of the Maastrichtian and Turonian). The boundary conditions are given in Figure 4, and all the aquifer’s parameters are shown in

Table 2. The aquifer parameters used for model setup and calibration.

Parameter	Value	Unit
Groundwater level	1.0–78	m
Number of layers	1	U
Aquifer thickness	150–300	m
Hydraulic conductivity	2 × 10−5–3 × 10−4	m/s
Specific storage	8–22	%
Recharge	84,354	m3/day
Withdrawal (wells)	36,986	m3/day
Stream leakage	195,320	m3/day
Boundary limits	147,953	m3/day

.

Finally, the simulation was developed over a period of 20 years (period 2010–2030).

3.2. Flow Simulation Model

3.2.1. Steady-State Model

Figure 5a shows the calibrated map in the steady-state regime. The value of the simulated groundwater head was comparable with observed data with a difference less than or equal to 1 meter, except in the zone located at the outlet, where the difference was more than 1 meter. However, the aquifer budget balanced (Figure 5b).

Figure 6 shows a scatter diagram of the observed heads and calculated heads in 17 monitoring wells. Clearly, in the steady state, the measured groundwater levels of the observation wells matched very closely with the calculated values (with a coefficient of determination R² = 0.96).

3.2.2. Transient Model

The purpose of this part was to use the results of this simulation to predict future groundwater vulnerability in the Tebessa-Morsott alluvial aquifer located in a semi-arid region. The results of [50] show that the climatic conditions and the net recharge would not highly contribute to a change in groundwater vulnerability in the future, especially in semi-arid and arid regions. Accordingly, two scenarios were simulated maintaining the same recharge rate but modifying the exploitation rate over a period of 20 years.

(a) First scenario: With the existing water points and their current flows, a period of 20 years was simulated. The piezometric map obtained with this scenario (Figure 7a) shows that the aquifer essentially keeps the same flow directions as in the case of the reference piezometry. Moreover, the aquifer budget is quite balanced (Figure 7b).

(b) Second scenario: In this case, the operating rate of the most productive boreholes was doubled in four different sectors of the study area, namely, El Hammamet, Ain Chabro, Tebessa and Bekkaria. The piezometric map obtained for this second scenario (Figure 8) shows a considerable drawdown of up to 40 m in the El Hammamet region. This situation could lead to aquifer pollution since this region is known for its agricultural activities. Moreover, the aquifer budget is unbalanced, marking a worrying situation.

3.3. Data Preparation and Analysis of Vulnerability Factors

Seven thematic maps were established for the calculation of the DRASTIC index:

• Depth to water

The depth to water is an important parameter in the DRASTIC model. For the Tebessa-Morsott aquifer, this parameter was determined during the piezometric surveys from May 2010. The map of this parameter (Figure 9a) shows the presence of six (6) classes. The lowest values are located near the west (north of the Hammamet city) and southwestern borders (Bekkaria city) and the central zone of the study area, while the higher values are located in the southern (north of the Tebassa city) and northern zones of the study area.

• Net Recharge

Recharge plays a key role in the transfer of water from the ground surface to the aquifer. To estimate this parameter in the Tebessa-Morsott alluvial aquifer, the results of the mathematical model established in the region was used [38]. The map of the “Net recharge” parameter (Figure 9b) shows three classes: (i) The first class covers the major central part of the plain, where, due to a lithological composition mainly of clay-loamy alluvium, the recharge does not exceed 50 mm/year; (ii) a second class with a recharge of 50 to 100 mm/year occupies the banks of the wadis (El kebir, El Ksob and Chabro wadis); (iii) a third class occupies the piedmonts of the Maestrichtian and Turonian limestone borders with a recharge of 100 to 180 mm/year.

• Aquifer media

The map of this parameter was produced using a geological map covering the region, the geological sections, the borehole logs for the aquifer and the geoelectrical cross-sections [40]. The analysis of the map relating to this parameter (Figure 9c) shows the presence of three zones: (i) Zone 1, characterized by the predominance of silt and clay formations with intercalations of gravel or sand; (ii) zone 2, characterized by clay and gravel formations of medium permeability (less vulnerable); (iii) zone 3, characterized by very permeable formations (gravel, pebbles and limestone).

• Soil media

The soil media influence the penetration of water-borne pollutants into the aquifer. The map was produced using the information collected from the soil maps and the borehole logs. The examination of the map relating to the soil media parameter (Figure 9d) revealed three soil classes: (i) A gravelly texture, which presents a high risk of vulnerability, located along the wadis; (ii) a clayey-gravelly texture, which presents a lower risk of vulnerability; (iii) a silty-clayey texture with a low risk of vulnerability that occupies the areas of the western, northeastern and southwestern borders.

• Topography

In the regions where the slope is steep, the runoff is high, and therefore, the contamination of groundwater is low. In the study area, the slope values were estimated by referring to the 1/25,000 topographic maps of Tebessa, El Hammamet and Morsott. The map of the slope parameter (Figure 9e) shows that more than 75% of the study area has a very low slope (0–6%), which favors the migration of the pollutant towards the aquifer. The least vulnerable areas are located at the edges.

• Impact of the vadose zone

The role of the vadose zone is very important in the DRASTIC method, since the nature and thickness of the zone located above the piezometric level controls the aquifer’s vulnerability. The geophysical study (geoelectrical cross-sections) and the borehole logs of the aquifer were used to establish the map of this parameter. This map (Figure 9f) shows the presence of three classes: (i) Class (1), characterized by large gravel formation or even very permeable limestone pebbles (high permeability), is present in the east part and the El Hammamet locality in the southeastern area of the plain, as well as the Morsott cities in the northern part; (ii) class (2), characterized by the clay and gravel formations (average permeability), is located mainly in the Bekkaria zone in the south of the plain and the Tebessa-Ain Chabro localities and the extreme east of the study area; (iii) class (3), characterized by the predominance of clay and silt formations (low permeability), occupies the entire central part of the plain.

• Hydraulic Conductivity

The hydraulic conductivity of the aquifer provides information on the propagation rate of pollutants in the aquifer. The pumping tests in twenty-two boreholes capturing the aquifer [38] were used to establish the map of this parameter. According to the established map (Figure 9g), three classes of conductivity are distinguished: (i) Class (1) presents a value of hydraulic conductivity oscillating from 14.7 × 10⁻⁵ to 32.9 × 10⁻⁵ m/s—this class is located in the Ain Chabro locality (south of the plain); (ii) class (2) presents a values of hydraulic conductivity varying between 4.7 × 10⁻⁵ and 14.7 × 10⁻⁵ m/s—this class is located in an area limited to Tebessa and El Hammamet cities in the south and Morsott city in the north; (iii) class (3) has very low values of hydraulic conductivity that vary between 4.7 × 10⁻⁷ and 4.7 × 10⁻⁵ m/s—this class is distributed over the entire plain.

3.4. Groundwater Vulnerability Using DRASTIC Index

Figure 10 shows the calculation of the DRASTIC index in the Tebessa-Morsott alluvial aquifer, classified in three classes of vulnerability: low, medium and high. The classes of average and high vulnerability cover more than 90% of the area. The aquifer is characterized by an average–high vulnerability, which is also negatively influenced by agricultural activities and urban discharges (Morsott, Hammamet and Bekkaria cities). However, in these localities were implanted the majority of boreholes for drinking-water supply and agricultural needs. These areas must now be subject to rigorous control to define the appropriate protective measures. The increase in the degree of vulnerability is mostly related to the pedological characteristics, vadose zone media and aquifer material (permeable).

3.5. Drastic Model Validation

Worldwide, nitrate is considered the main anthropogenic pollutant. It is linked to the leaching of store rocks and to human activity (agricultural and industrial activities). To validate the results of vulnerability, a simple linear correlation was made between the DRASTIC index and NO₃ concentrations. The choice of this method was motivated following its successful use in several studies [14,45,51,52]. In the study area, a value of R² = 0.955 was calculated (Figure S4); this is a very acceptable value for this work.

3.6. Future Vulnerability Assessment

In this section, the work aims to evaluate the future DRASTIC classification in the Tebessa-Morsott alluvial aquifer according to the two simulated scenarios for the transient model previously proposed in this study in response to a groundwater level predicted in the aquifer for the year 2030.

3.6.1. First Scenario

Assuming the conditions in the first scenario (current operating rate), Figure 11a represents the future distribution of the vulnerability classes according to the DRASTIC method for the year 2030. The results show a decrease in the average–high vulnerability classes compared to the reference groundwater vulnerability map (2010) in the study area. Conversely, the zone of the low vulnerability class will increase in the southern part (west of Bekkaria city). This class increases by 2.5% compared to the map based on the data in 2010.

3.6.2. Second Scenario

The results in Figure 11b indicate the future distribution of the vulnerability in relationship to the second scenario (increase the operating rate) for the year 2030 according to the DRASTIC method. A very important and remarkable displacement of the low vulnerability class over a period of 30 years (2030) compared to the reference groundwater vulnerability map (2010) due to overexploitation was observed. The map shows the presence of two important zones for this class located in the east (north of Hammamet city) and in the southern part (west of Bekkaria city). The increase in the area has a rate of 5.3% and 3.1% compared to the reference and the first scenario maps, respectively. In this study and based on the modeling conditions applied to the MODFLOW model, the variations in climate conditions (very low recharge rate) do not significantly influence the groundwater vulnerability using the DRASTIC method in the Tebessa-Morsott alluvial aquifer; conversely, overexploitation in this area can decrease the risk of vulnerability. Moreover, the piezometric level changes have a strong relationship and directly influence groundwater vulnerability.

The new proposed approach considering the hybridization between a numerical model (MODFLOW model) and a groundwater vulnerability method (DRASTC) to predict actual and future groundwater vulnerability represents a step forward from the canonical concept of static representation of vulnerability. The possibility to forecast short- and long-term situations according to the most suitable climatic and socioeconomical scenarios is a valuable tool for groundwater management. At the regional scale, this proposed concept can be easily adapted for different aquifer types (karst, fissured or porous) and climatic regions. On the other hand, the concept’s suitability for regional-scale assessment is currently the main conceptual drawback, since the numerical simulation requires detailed aquifer information that is not available everywhere, especially on the national scale. Moreover, it is of paramount importance to include in the assessment the pollutants’ transport dynamic, nitrogen formation and removal processes in both vadose and saturated zones [53], along with the role of fracture and sinkholes in nitrogen leaching [54], to ensure proper management of land and water resources.

4. Conclusions

Today, in many regions worldwide, groundwater is the main source for domestic, agricultural and industrial water supply. However, the risk of pollution phenomena has become a serious problem for the quality of groundwater reserves. This problem is widely known in countries with semi-arid and arid climates due to low rainfall and high overexploitation of aquifers following the increase in demand for drinking-water supply and water for the agricultural and industrial sectors. This high demand for water causes a rapid drawdown in aquifers. In this work, the concept of groundwater vulnerability was used in the Tebessa-Morsott alluvial aquifer (northeastern Algeria). A new proposed approach for assessing groundwater vulnerability was considered that integrates a numerical model (MODFLOW model) and an index–overlay method (DRASTIC method) to predict actual and future groundwater vulnerability and the groundwater resource sustainability of the Tebessa-Morsott alluvial aquifer.

The model was calibrated and used to predict the piezometric level fluctuations for the period from 2010 to 2030 under two scenarios of water demand, where the exploitation rate was the important parameter in the simulation. The results of this calibration show that the simulated piezometric level fluctuations continue to decline over a period of 20 years, with the drawdown varying between 12 and 40 meters. The obtained results of the DRASTIC method application based on the data for the reference year (2010) indicates that the risk of pollution is significant since the aquifer is characterized mainly by average and high vulnerability (about 97%). Furthermore, nitrate concentration data confirmed this critical situation.

The results of the new proposed approach, which considers the integration of the MODFLOW model and the DRASTIC method to predict actual and future groundwater vulnerability for the year 2030 in the Tebessa-Morsott alluvial aquifer, indicate:

• The study area is classified at a average–high pollution risk of groundwater and environmental deterioration.
• The increase or decrease in the pollution risk for groundwater (high or low vulnerability) is closely related to the piezometric level variation (groundwater depth) caused by the pumping rate (overexploitation or pumping reduction).
• Overexploitation ensures the protection of the water table (deep piezometric level) from all pollution types observed on the ground surface.
• Conversely, overexploitation has a negative effect on the hydrodynamic state of the aquifer: significant drawdown, depletion of the water table, etc.
• Some measurements should be proposed to protect this groundwater system.

For better groundwater management of the experimental site, artificial recharge by treated wastewater from the study area supplemented by a hydrochemical control is recommended. This technique has been successfully tested in arid and semi-arid regions by several authors [55,56,57,58]. Finally, the results of this new proposed approach are satisfactory and can be applied to predict actual and future groundwater vulnerability in many regions with semi-arid and arid climates.