Comprehensive Assessment of the Jebel Zaghouan Karst Aquifer (Northeastern Tunisia): Availability, Quality, and Vulnerability, in the Context of Overexploitation and Global Change

Abstract

Karst aquifers in the Mediterranean region are crucial for water supply and agriculture but are increasingly threatened by climate change and overexploitation. The Jebel Zaghouan aquifer, historically significant for supplying Carthage and Tunis, serves as the focus of this study, which aims to evaluate its availability, quality, and vulnerability to ensure its long-term sustainability. To achieve this, various methods were employed, including APLIS and COP for recharge assessment and vulnerability mapping, SPEI and SGI drought indices, and stable and radioactive isotope analysis. The findings revealed severe groundwater depletion, primarily caused by overexploitation linked to urban expansion. Minimal recharge was observed, even during wet periods. APLIS analysis indicated moderate infiltration rates, consistent with prior reservoir models and the MEDKAM map. Isotopic analysis highlighted recharge from the Atlantic and mixed rainfall, while Tritium and Carbon-14 dating showed a mix of ancient and recent water, emphasizing the aquifer’s complex hydrodynamics. COP mapping classified 80% of the area as moderately vulnerable. Monitoring of nitrate levels indicated fluctuations, with peaks during wet years at Sidi Medien Spring, necessitating control measures to safeguard water quality amid agricultural activities. This study provides valuable insights into the aquifer’s dynamics, guiding sustainable management and preservation efforts.

1. Introduction

Karst formation refers to a special type of land where soluble carbonate rocks such as dolomite and limestone have dissolved. High porosity and permeability are often associated with karst terrain, allowing water to move quickly via underground networks forming aquifers that are, highly productive but extremely vulnerable to contamination [1,2,3]. Various Mediterranean countries exhibit a significant prevalence of carbonate rock outcrops, with 15% of the Mediterranean catchment area characterized by such formations (outcrops) [4,5], show that the Mediterranean region has long been relevant for the research of karst aquifers. Essential for water supply, karst aquifers contribute at least 25% to the region’s domestic water requirements, excluding withdrawals for industrial, agricultural, and tourism purposes. It is worth noting that Karst Groundwater-Dependent Ecosystems (KGDEs) play a crucial role in providing ecosystem services and fostering biodiversity in the Mediterranean region. Furthermore, the Mediterranean region experiences strong climate seasonality, featuring prolonged, dry, and hot summers, which poses a threat to its karst springs. These valuable water resources face increasing pressure from human activities and are susceptible to climate change and anticipated anthropogenic pressures, including heightened water extraction and alterations in land cover and land use [6]. Among the various threats, human disturbances emerge as the predominant type, accounting for 48% of the identified pressures [7]. The overarching concern is to ensure a sustainable freshwater supply from these karst formations, considering potential challenges and impacts on their availability due to environmental changes and human activities [8].

Particularly, in the southern Mediterranean, situated within semi-arid zones, the evaluation of karstic aquifers is of paramount importance due to several interconnected factors that amplify the challenges associated with water resources management. These regions often face a complex set of issues regarding water resources management due to rising water demand induced accelerated population growth, agricultural and urban expansion, and the impact of climate change [9,10]. Thus, karstic aquifers, susceptible to contamination and overexploitation, face intense competition for water resources. Comprehensive assessment becomes imperative for equitable resource allocation and sustainable water supply [11,12].

The scarcity of comprehensive hydrogeological data in these semi-arid areas is a pervasive issue, further intensified when dealing with the complexities of karstic aquifers. This limited data availability hinders a profound understanding of aquifer behavior and vulnerabilities, making the assessment of these systems critical for informed decision-making [13].

Moreover, the vulnerability of semi-arid regions to alterations in precipitation patterns as a result of climate change introduces uncertainties regarding water availability. Karstic aquifers, known for their responsiveness to changes in recharge, are especially vulnerable to shifts in the frequency and intensity of precipitation. It is imperative to comprehend the repercussions of climate variability and change on these aquifers to ensure effective management of water resources [10,14].

Drought events, a common occurrence in semi-arid regions, lead to fluctuations in groundwater levels. Karstic aquifers, with their rapid drainage and recharge, may respond differently to prolonged droughts, impacting water availability [15]. A thorough assessment of karstic aquifers during drought conditions is essential for developing effective water resource management to meet growing water demand and preserve ecosystems [7,16,17]. Particularly, the integration of urban expansion necessitates comprehensive strategies for balancing water needs in the face of rapid urbanization [18,19]. The assessment of karstic aquifers in semi-arid regions is multidimensional and multidisciplinary. It addresses data limitations, verifies the suitability of conventional investigation methods while studying the impacts of climate change and drought, and ensures sustainable water management practices.

In karst hydrogeology, various methods are utilized to study aquifer systems. The choice of methods and their sequence depends on research questions, prior knowledge, and resource availability. Below is a summary in

Table 1. Summary of Methods and Applications in Karst Hydrogeology.

Method	Application
Geological	Elucidates the aquifer’s geometry and flow paths.
Geophysical	Identifies geological structures and fractures.
Speleological	Maps conduit networks, providing direct insights.
Hydrological and Hydraulic	Establish water balances and determine hydraulic properties.
Hydrochemical and Isotopic	Trace water origins, movements, and interactions.
Artificial Tracers	Delineate connections and flow velocities [20]

:

Each method has specific applications, advantages, and limitations, contributing uniquely to understanding karst aquifer dynamics [21,22].

In Tunisia, situated in the southern Mediterranean basin, karst resources constitute more than 13% of the groundwater reservoirs, with the karst hydrogeology influenced by carbonate outcrops [23,24,25]. These karsts are a vital source of water, prioritized first for domestic supply and then for agricultural use. However, nowadays they are facing a decline in water availability and quality and are subject to conflicts over the use, inducing in some cases, social tensions. The evaluation of the percentage of karst formations in Tunisia, subject to commentary and refinement, particularly in the North, is of paramount importance especially during drought conditions. The lack of continuous and historical monitoring further accentuates the need for testing and implementing suitable methods [26,27,28]. This imperative extends not only to Tunisia but also to developing countries in Africa, where effective strategies are essential to address the specific challenges posed by drought, limited monitoring resources (poorly investigated sites) and local conflicts [27,28].

The Zaghouan aqueduct in Tunisia is an example of an infrastructure built to supply water from a karst aquifer to the city of Tunis, specifically for the city of Carthage [29]. The Zaghouan karst has often been studied from a geological and hydrogeochemical point of view (i.e., [25,30]). Availability, quality and vulnerability aspects remain little explored.

This study conducted a comprehensive assessment of the Jebel Zaghouan aquifer, a critical water source in Northeast Tunisia, addressing both groundwater availability and quality. The research further seeks to map the aquifer’s vulnerability to pollution.

The assessment of aquifer availability, quality, and vulnerability is essential as these factors are interconnected and define groundwater sustainability. Availability refers to the aquifer’s capacity to meet current and future water demands, quality ensures the water’s usability for domestic, agricultural, and industrial purposes, and vulnerability evaluates the risk of contamination or depletion due to natural or anthropogenic factors. By considering these dimensions together, the aim is to provide an integrated understanding of the aquifer’s long-term resilience and its ability to support sustainable socio-economic development.

Through the examination of data from nine boreholes and galleries, the study analyzed historical trends and assessed the sustainability of the aquifer. Using the APLIS method (altitude (A), slope (P), lithology (L), infiltration (I), and soil type (S); cf. Section 3.3), spatial recharge rates in the Jebel Zaghouan aquifer will be evaluated. The research will also focus on physicochemical analyses conducted by the Gdhir El Golla (Tunis, Tunisa) water laboratory of SONEDE (National Water Distribution Utility, Tunis, Tunisia), spanning the years 2004 to 2021, to understand variations in groundwater quality and identify potential influencing factors. Stable and radioactive isotopes, tritium and Carbon-14 will be analyzed to understand the hydrodynamic system and functioning of the aquifer. The vulnerability mapping will be based on the COP (Overlying layers (O) factor, Concentration of flow (C), precipitation (P); cf. Section 3.6) method. Emphasizing the importance of ongoing monitoring, these investigations should provide us with information on trends in the use of this aquifer and provide decision-makers with tools to improve the management of the karstic system to avoid its depletion, ensuring its long-term sustainability.

2. Study Area

Jebel Zaghouan, located roughly fifty kilometers from Tunis (Tunisia), constitutes a pivotal Jurassic formation within the Zaghouan massif, it spans an area of approximately 19.6 km² (Figure 1). The Zaghouan region belongs to an upper semi-arid to subhumid climate, with an average annual rainfall of 467 mm characterized by heterogeneous spatial distribution and significant temporal fluctuations (ranging from 245 to 625 mm). The average annual temperature hovers around 17.7 °C.

2.1. Geological Context, Aquifer Geometry and Springs

The Zaghouan anticline is mainly constituted by Jurassic limestone. It is limited by the rock-fall and the cretaceous formations. Figure 2 shows that the geology of the Jebel is characterized by the presence of southern and transverse faults that have created individualized blocks. These faults which allow an infiltration of meteoric waters are between jumps: the Kef El Orma fault, the Great Peak fault and the Achilles fault. The Zaghouan karst aquifer’s eastern section is conducive to seepage water storage, while the western part, influenced by substantial marl deposits, displays a diminished storage capacity [31].

The Jurassic limestone block is a trapezoidal, cavernous, and fissured limestone with a longitudinal dimension of about 8 km along a north to 40° East direction, and an average of 2.4 km in the transverse direction, along a north to 45° West direction. It has a surface area of l9 km² at an altitude of 300 m NGT. The massif is surrounded by marly soil acting as a watertight barrier and can be subdivided into three compartments running from north to south:

• Small Zaghouan which gives birth to Ain Haroun.
• Transmission station massifs, Kef El Orma, Kef El Blidah and Jebel Stâa; they are the most extensive compartment, which give rise to the most important springs including Water temple (Nymphée), Aïn Ayed and Aïn Oued El Guelb.
• The great peak massif which gives birth to the source of Sidi Medina.

The general dip of the limestone layers and the topographical configuration towards the north-west explain the presence and importance of the springs of both slope and east massif. The massif contains 14 springs, the most important springs flowing from the karst, they are on massif north-western slope among them: Nymphée (then captured by Galerie 44), Ain Ayed, Ain El Guelb and Galerie 47 and Ain Haroun (Figure 3).

2.2. Water Resources

The geological composition of Jebel Zaghouan has positioned it as a vital water reserve, harnessed since the roman era, the karst springs of Zaghouan, supplied drinking water for the local cities and for the capital (Carthage then Tunis) through a Roman aqueduct (120 A.D.), of 132 km, still visible along the road from Tunis to Zaghouan. Due to the drought, three galleries have been built in front of the Nymphée and still exist today. The first was built in 1928 (Ain Ayed). Following the drought of 1941–1944, a second gallery was built in 1944 at an elevation of 274.41 m. Pressure for water led to the construction of another gallery in 1947, at an elevation of 265 m. The two galleries (Galerie 44 and Galerie 47) are approximately 300 m apart. They are all equipped with control valves that allow consumers to be served according to their needs. A series of boreholes were also installed from the nineties. Currently, the aquifer is exploited by mainly nine boreholes and galleries intended for the drinking water supply of the city of Zaghouan and the surrounding rural agglomerations. Three of these wells used as commercialized mineral water (Cristaline, Aqualine and Prestine). Since Galerie 44 and Galerie 47 are dry and to cope with the water shortage that Zaghouan city suffers from, two other boreholes (Water temple and Ain Haroun 3bis) were drilled in 2017 and 2018.

3. Materials and Methods

3.1. Historical Data Processing

3.1.1. Discharge Flows

The available flow data corresponding to the natural flow period was recorded from 1915 to 1927. Figure 4 presents the Zaghouan springs production before the digging of the galleries and Zaghouan springs production with exploitation by the galleries respectively. The natural flow period was marked by heavy rainfall of the 1920–1921 and a low rainfall during the 1926–1927 hydrological years, which resulted in high spring flow 6.5 and a very low flow of 1.9 million cubic meters, respectively. These observations are in conformity with the natural flow of the resurgences during this period.

Errors and uncertainties are principally due to methods and time step measurement as well as data processing. Indeed, discharge series were available in graphical form (Figure A1) from 1915 to 1927 at irregular time scales. Discharges were then obtained by linear interpolation on a daily scale. Thus, several sources of uncertainties are issued from the rough data and its interpolation. In fact, measurements were taken on a weekly, twice-weekly, or once-monthly basis rather than daily (Figure A1). In addition, there are the uncertainties associated with the measurement techniques used at the beginning of the 20th century, based on weirs and their calibration curves.

Besides, and since the installation of the valves to control the flows supplied by the galleries, the system is no longer natural. The flows observed from 1958 to 1962 and from 1971 to 1995 are very random and indicate that they are highly dependent on the openings of the gates. These openings depend on several contingencies, in particular SONEDE’s daily demands to meet its users’ water needs.

3.1.2. Precipitations

The rainfall data comes from seven stations located in the study area.

Table 2. Meteorological stations characteristics.

Meteorological Station	Observation Period of Data Series	Location (WGS 84)		Mean Annual Precipitation (mm)	Altitude (m)
		X	Y
Zaghouan PF	1961 to 2017	10.129	36.425	470	130
Mograne CSA SM	1971 to 2020	10.090	36.432	485	151
Zaghouan SM	1909 to 2016	10.153	36.403	481	165
Zaghouan DRE	1976 to 2017	10.144	36.403	449	238
Zaghouan Contrôle	1915 to 1951	10.149	36.395	497	241
Zaghouan Sidi Bou Gabrine	1991 to 2017	10.116	36.375	492	677
Zaghouan Poste Optique	1912 to 2008	10.140	36.441	516	945

summarizes location, precipitation series length, altitude and mean annual precipitation of these stations. Figure 5 shows the spatial variation of precipitation over the karst. The dataset comprises a daily rainfall time series spanning from 1909 to 2020, exhibiting varying lengths and occasional gaps. Only complete daily series were included in the analyses. However, for stations situated at higher altitudes, where data collection is challenging, cross-validation was conducted to verify annual rainfall totals.

3.1.3. Evapotranspiration

Due to the lack of daily potential evapotranspiration data, it was calculated on a monthly scale using Thornthwaite’s method based on monthly mean temperatures and the massif’s latitude. The average monthly temperatures were taken at the Mograne CSA SM station.

3.2. Geochemical and Isotopic Analysis and Methods

The Gdhir El Golla Water Laboratory of SONEDE conducted physicochemical water analysis at active springs, galleries, and boreholes (Figure 3) each year from 2004 to 2021. The analysis was performed twice a year (Generally in January and June or July, rarely in August). The historical data from these analyses has been collected, added to a database, and linked to Carthage coordinates using the QGIS software (QGIS Development Team. QGIS Geographic Information System. Open Source Geospatial Foundation Project. Version 3.28, 2023. Available at: https://qgis.org. (accessed on 26 August 2024). Nitrates issued from this database were used to analyze the impact of rainfall intensity on groundwater vulnerability.

A weekly campaign of groundwater sampling started in April 2021 with the support of the CRDA (Regional Agricultural Development Commission, Zaghouan, Tunisia) at Jebel Zaghouan (Temple) and Ain Ayed 3 bis boreholes. pH varies between 6.12 and 7.93. Electrical conductivity ranges from 1043 μS/cm to 1147 μS/cm. Temperature varies between 15 °C and 21.1 °C (Figure 6). A sampling campaign of the Jebel Zaghouan groundwater aquifer was conducted in September 2022. The water samples were analyzed at the central laboratory of the Ghdir El Golla water complex. Major anions and cations were measured using liquid phase chromatography (Metrohm ion chromatography). Alkalinity was determined via the potentiometric method. The ion balance of the samples ranged from 0 to 5%. A protocol of weekly sampling and local isotopic measurements of Jebel Zaghouan groundwater, was set up as follows:

• From mid-April 2021 to mid-July 2021 at the Temple borehole
• From early October 2021 to October 2022 (at the Ain Ayed 3bis borehole) for the latter have been analyzed for the moment only samples until mid-April 2022.

Samples were analyzed at the Lama laboratory of the UMR 5151 Hydrosciences (University of Montpellier, France). In September 2022 exploratory investigations of radioactives were made, 2 measurements of tritium (at the Temple and Ecomusée boreholes) and 1 of carbon 14 (Ecomusée) have completed the protocol.

3.3. Recharge Estimation: APLIS Method

The APLIS method [32] was developed to estimate the average annual recharge of carbonate aquifers under Mediterranean conditions, expressed as a percentage of precipitation. It is GIS-based method that contains inherent variables that affect recharge: altitude (A), slope (P), lithology (L), infiltration (I), and soil type (S).

The APLIS method can help determine the typical rate of annual recharge in carbonate aquifers by utilizing intrinsic characteristics of the aquifers. It can also be used to create a map of the distribution of recharge rates based on the intrinsic properties of the aquifers. This information is crucial for effectively managing and preserving groundwater resources in terms of both amount and quality. It also allows for an estimation of the average annual resources.

For each of these layers of information, it attributes a score system to create a map of each variable. Scores range from 1 to 10, following an arithmetic progression with a step width of 1, with the aim to easily equate to aquifer recharge rates. The value 1 indicates minimal incidence of the variable in the recharge of the aquifer, whereas the value 10 expresses the maximum influence on recharge.

The layers of information corresponding to each variable are taken from different sources [33,34,35,36]. The application of the method to our case study is detailed in Appendix B.

3.4. Water Budget Estimation, Trend and Uncertainties

The water budget analysis primarily relied on the modeling investigations conducted by [31,37]. This study aimed to evaluate the water balance and quantify the aquifer storage capacity associated with the Jurassic limestone formations of Jebel Zaghouan.

Refs. [31,37] developed a conceptual deterministic model to convert the precipitation received by the calcareous rock mass into the aggregated discharge fluxes, including springs and galleries. The model’s validity was confirmed through the utilization of meteorological and hydrodynamic data. Ref. [37] accounted for a comprehensive dataset spanning both dry and wet years, with a daily calculation time step.

The model underwent calibration and validation processes, covering the natural dynamics of the system from 1915 to 1927 and validation from 1970 to 1995, inclusive of aquifer exploitation via galleries and wells. The model’s performance was deemed acceptable, with Nash criteria ranging from 54% to 77%.

Figure 7 (produced for the present study) presents a summary of the water budgets for Jebel Zaghouan, encompassing rainfall, spring flow, Real Evapotranspiration (RET), runoff, water budget, and infiltration rate, for both the calibration (Figure 7a) and validation (Figure 7b) periods. These components collectively represent the natural water budget (1915–1927), as outlined in Sagna’s (2000) [37] classical calibration and validation methodology.

3.5. Comprehension of Groundwater Behavior to Drought Conditions

The Standardized Precipitation-Evapotranspiration Index (SPEI) is a drought index that was first introduced by [38]. It is based on the Standardized Precipitation Index (SPI), which was introduced by [39]. SPEI characterizes meteorological drought by taking into account both precipitation and evapotranspiration. The SPEI is particularly suited for arid and semi-arid regions as it integrates temperature alongside precipitation [38]. Unlike SPI, which only considers precipitation, SPEI provides a comprehensive atmospheric water balance. Its reliance on readily available data and simplicity make it a practical tool for drought monitoring in data-scarce regions like ours. Nevertheless, SPEI’s limitations, such as its dependency on evapotranspiration estimation methods, should also be considered. Additionally, SPEI does not directly account for other factors affecting drought conditions, such as soil moisture or vegetation response, which may be captured by other indices like the Palmer Drought Severity Index (PDSI) or the Vegetation Health Index (VHI) [38].

The SPEI is calculated by comparing the difference between precipitation and evapotranspiration with the long-term average for a specific location and time period, using historical data. Typically presented on a scale from −3 to +3, a positive SPEI value indicates that the moisture balance is above average (wet conditions), while a negative value indicates that it is below average (drought).

The Standardized Groundwater Level Index (SGI) is a drought index that characterizes groundwater droughts; it was first introduced by [40]. It is built on the Standardized Precipitation Index (SPI). Unlike meteorological droughts, groundwater droughts are characterized by a decline in groundwater levels, rather than a lack of precipitation. The SGI is calculated using a continuous variable, the groundwater level, which is collected from monitoring wells. Because groundwater level is a continuous variable, there is no need to accumulate it over a specified time period, as it is with precipitation data. The SGI is calculated using a normalization technique that allows the comparison of groundwater levels across different time periods and locations [41,42].

By combining and superposing SPEI and SGI, one can gain a more holistic understanding of drought conditions, considering both surface water and groundwater dynamics [43]. This integrated approach can lead to more informed decision-making in managing water resources during periods of water stress.

3.6. Vulnerability Mapping Using COP

The COP method assesses the vulnerability of groundwater catchment areas by evaluating factors such as topography, soil, land use, lithology, karst features, water table depth, and climatic conditions. The method uses Geographic Information System (GIS) tools, such as QGIS, to compile different combinations of geological and hydrogeological datasets to identify vulnerable areas and create a final vulnerability index. The method focuses on the surface of the groundwater (the water table) by considering the unsaturated zone factor (O), the ability of the underlying soils to protect groundwater, and the impact of diffuse or concentrated infiltration (C) and climatic conditions (P) on the level of protection (Figure A2). The development of the various layers was based on multiple data sources and references [44,45,46]. The processing steps for each layer are detailed in Appendix B.

4. Results

4.1. Water Availability

Due to water shortages in the city of Zaghouan, two additional boreholes were drilled in 2017 and 2018 to supplement the supply. The daily amounts of water extracted from these boreholes were collected with the assistance of SONEDE and CRDA.

Figure 8a shows the temporal evolution of the static level at the piezometer located in the site Ain Ayed (near the dry spring and the new boreholes) as well as the monthly produced amounts (both natural discharge and extracted from boreholes) and Figure 8b shows the temporal evolution of the static level with the monthly rainfall. It is noted that:

• The natural discharge that is produced through the galleries nearly stopped in 2015.
• The static level follows a continuous declining trend since 2002 with an expectation of a slight rise recorded in 2012 following a snow event. Indeed, the recorded groundwater level dropped from −0.7 m (2002) to more than −100 m (2022) (below soil level), the total depth of the piezometer being surpassed making it dysfunctional. The continuous depletion of the groundwater seems to be mainly linked to overexploitation (due to urban expansion) than a state of meteorological drought.

Figure 9 illustrates the variation in monthly calculated SPEI (Zaghouan) and SGI (Ain Ayed) from 2002 to 2020, highlighting two distinct periods. Between 2002 and 2011, several drought events, as indicated by negative SPEI values, correspond to a significant decline in groundwater levels. This period reflects the strong influence of climatic conditions on groundwater depletion. In contrast, from 2011 to 2013, a slight recovery in groundwater levels is observed, indicated by the positive trend in SGI, which correlates with the wetness conditions (SPEI > 2). However, beginning in 2018, severe groundwater drought (SGI < −1.5) is evident despite normal to wet conditions, underscoring a decoupling between climatic conditions and groundwater dynamics, likely driven by non-climatic factors such as overexploitation.

The analysis reveals a consistent decline in groundwater levels, as evidenced by a significant negative trend detected using the Mann-Kendall test (Kendall’s tau = −0.807, p-value < 0.0001, Sen’s slope = −0.163). This non-parametric method identifies monotonic trends in time series data without assuming a specific distribution, comparing all possible pairs of observations to assess whether differences consistently increase or decrease. The results confirm a downward trend, with an acceleration of this decline observed from 2012 onwards.

Prior to 2011, a positive correlation between SPEI and SGI was evident, indicating that increased precipitation led to aquifer recharge, while decreased precipitation resulted in declining aquifer levels (

Table 3. Correlation between SPEI and SGI.

Observation Period	Kendall’s Tau	p-Value
2002–2011	0.157	0.01844
2011–2020	0.066	0.291

). However, starting from 2011, these variables became independent and no longer exhibited correlation (Table 3). The significance of these findings is evaluated based on commonly accepted p-value thresholds, where p-values less than 0.05 indicate statistically significant trends or correlations, providing robust support for the observed changes in aquifer behavior. Despite above-average precipitation levels recorded from May 2018 to February 2020 (SPEI > 1), a continuous decline in groundwater levels was observed without significant recharge, reaching extreme lows (SGI approaching-3).

To assess the impact of climate change, an in-depth statistical analysis of trends in drought characteristics—such as magnitude, duration, intensity, and frequency of events—would be required, as outlined by the WMO (2023) in their guidelines on extreme weather and climate events (WMO-No. 1310) (WMO, 2023: Guidelines on the Definition and Characterization of Extreme Weather and Climate Events (WMO-No. 1310)). In this study, the observed drought trend is primarily hydrological, as evidenced by the continuous decline in groundwater levels, even during wet periods. This pattern aligns with the British Geological Survey’s (https://www2.bgs.ac.uk/groundwater/waterResources/drought_overview.html, accessed on 18 January 2025) definition of groundwater drought, which describes it as the sustained and extensive occurrence of below-average groundwater availability, marked by lower-than-average water levels in aquifers, boreholes, and wells.

4.2. Recharge Estimation

Using APLIS method, the majority of the study area has estimated infiltration rates in the ‘‘moderate’’ category (40–50%) (Figure 10). More than 73% of the study area has an infiltration rate between 40% and 42.5%. This result is compliant with a previous reservoir modelling performed at the test site [31]. Only 6.7% of total area shows an infiltration with a percentage exceeding 45%. Recharge rates issued from the conceptual model [31] for the natural flow period ranged from 30% to 68% (average 42%). For the period corresponding to the functioning of the system via galleries by regulating valves, the infiltration rates varied between 21% and 57% (average 39.4%). The total volume of recharge in millimeters calculated for the total surface of the aquifer is approximately 207 mm, considering the precipitation of mean year (Table 2). This value corresponds to 45% of the median year’s rainfall depth (calculated based on spatial average). Figure 11 shows the spatial distribution of the recharge volume. Results from the new Mediterranean Karst Aquifer Map “MEDKAM” [47] show that recharge is about 162 mm which corresponds to 35% of median annual rainfall depth. MEDKAM gives a recharge value belonging to the range found using with APLIS (21% to 57%) and the conceptual model developed (30% to 68%).

Based on the average recharge rate estimated by the APLIS method, we calculated the safe yield for a year of average, maximum and minimum rainfall over the period of natural fluctuation from 1915 to 1927 (

Table 4. Yield of groundwater.

	Rainfall (mm/Year)	Observed Discharge (Million m3/Year)	Recharge (APLIS) Million m3/Year	Yield 1 APLIS (Million m3/Year)
Average year	467	3.6	3.8	0.2
Wet year	867	6.5	7.0	0.5
Dry year	239	1.9	1.9	0.0

Note: ¹ https://gw-project.org/books/a-glossary-of-hydrogeology/, accessed on 5 January 2025.

).

In average and wet years, the yield is positive (0.2 and 0.5 million m³/year, respectively), indicating excess recharge compared to discharge. This suggests that the aquifer can sustain exploitation without overexploitation. In a dry year, the yield is zero (0.0 million m³/year), meaning recharge equals the observed discharge. This represents a critical situation where there is no surplus to support additional extraction or natural recovery.

Furthermore, not only is the recharge moderate, but the hydrodynamic response is considered slow [48,49]. Ref. [49] assessed the lag time between the precipitations and discharges, using Frank copula. They found a lag of about three to four months. This was in accordance with results of earlier research in the PRIMA-KARMA project (2023) [29], with both the lumped reservoir model KarstMod [50,51] and artificial neuronal networks simulations [52]. Ref. [48] introduced a novel classification system for karst hydrological functioning based on the analysis of 78 spring discharge time series data. The Zaghouan aquifer was classified as C6 according to this typology characterized by low to medium variability of the hydrological response. Determining the hydrological response time or lag time is also important for the comprehension of the hydrodynamic behavior of the karst and can be useful for sustainable management [53].

Despite the uncertainties of measurements of discharges and flow and climatic data processing, APLIS method gave recharge rates of the same range as the values calculated by the model. The results should be checked and validated by other recharge assessment methods in order to contrast the reliability of the results in terms of aquifer management and protection groundwater.

4.3. Water Quality Isotopes and Geochemical Processes

4.3.1. Water Types and Saturation Indices

The Piper diagram (Figure 12) indicates that the groundwater in the study area can be classified into two facies:

• Calcium sulfate water: This is the predominant facies, represented by five groundwater samples from boreholes of the Jebel Zaghouan aquifer (Pristine, Ain Ayed 3 bis, Temple, Ecomusée and Cristaline). The facies results from the presence in the carbonate rocks of some evaporitic minerals (gypsum) easier to be dissolved than carbonates
• Calcium bicarbonate water: This facies is less common, represented by a single sample from the natural spring of Sidi Medien.

Electrical conductivity (EC) of karst springs, easily measured at high temporal resolution and low cost, fluctuates significantly during storms, reflecting internal aquifer behavior and piston flow effects. These fluctuations, influenced by carbonate dissolution and stormwater dilution, make EC a valuable indicator for distinguishing karst systems, identifying runoff composition, and studying aquifer dynamics [54,55].

EC of the first group of boreholes (calcium sulfate water type) ranged between 924 and 1244 µS/cm, with an average conductivity of approximately 1100 µS/cm. In contrast, the conductivity at Sidi Medien is 518 µS/cm, which is less than half of this average (

Table A1. Chemical analysis of water samples (Concentrations in mg/L), Electrical Conductivity (EC) in µS/cm, pH, and Saturation indices.

				(mg/L)								Saturation Indices
Well/Spring	T °C	pH	EC (µS/cm)	Ca	Mg	Na	K	HCO3	Cl	SO4	NO3	Is Calcite	Is Dolomite	Is Gypse	Is Anhydrite
Pristine	22.3	7.22	1244	175	27	40	0.98	216	59	317	5.7	0.28	0.09	−0.65	−0.87
Sidi Medien	22	7.23	518	77	13	21	2.4	192	38	38	19	−0.03	−0.49	−1.82	−2.04
Aïn Ayed 3 bis	21.8	7.99	1092	161	32	50	1.15	188	66	386	3.8	0.94	1.54	−0.61	−0.83
Temple	25	7.6	1178	163	31	70	1.06	193	91	380	2.9	0.57	0.76	−0.62	−0.84
Ecomusée	21.7	7.76	1091	167	32	43	1.13	163	63	409	2.5	0.66	0.97	−0.57	−0.79
Cristalline	22.2	7.45	924	137	28	41	1.07	236	65	227	3.4	0.46	0.58	−0.88	−1.1

). This significant difference suggests that the residence time of water in the karst system at Sidi Medien is considerably shorter compared to the other boreholes [1].

The saturation indices of groundwater samples from the deep aquifer of Jebel Zaghouan, with respect to minerals such as calcite, dolomite, gypsum, anhydrite, and halite, reveal that all samples are undersaturated concerning halite, gypsum, and anhydrite, indicating a tendency to dissolve these minerals (Table A1). Conversely, most groundwater samples are saturated to supersaturated with respect to calcite and dolomite, suggesting a tendency for these minerals to precipitate into a solid phase. This implies that the salt load is not significantly influenced by interactions between water and carbonate minerals. Therefore, the dissolution and precipitation of minerals contribute to the saline charge acquisition in groundwater from Jebel Zaghouan. Notably, the natural spring sample from Sidi Medien is the only one undersaturated with respect to calcite, indicating rapid flow and short residence time conditions within the aquifer [56,57].

4.3.2. Results for Stable Isotopes

Any time series or sampling in an aquifer must be compared to a local or regional isotopic chronicle of monthly or even daily precipitations. However, the response times to rain are much longer in an aquifer, that corresponds to a mixture of rain on an annual or even probably interannual scale. Initially, the focus is on the isotopic value of the theoretical annual recharge, or even interannual if data spans several years. This isotopic value is weighted by the rainfall corresponding to a volume of water likely to recharge the water table. For karst aquifers, where the hydrodynamic response can be extremely fast, interest extends to smaller time periods, such as monthly or even event-driven intervals, depending on the aquifer’s rapid response to recharge.

Although no isotopic measurements of rainfall are available on site, there are two daily isotopic rainfall records in Bizerte (2009–2014) and Tunis (2009–2018) analyzed at LAMA Hydrosciences laboratory that will serve as a reference for comparison with the karst waters of Zaghouan.

For Bizerte the annual weighted oxygen-18 composition of rainfall over six years ranges from −4.12‰ to −6.12‰ V-SMOW with an interannual mean of −5.37 ‰ V-SMOW and for Tunis (ten years) the annual weighted isotopic composition of rainfall ranges from −4.90‰ to −6.72‰ with an interannual mean of −5.64 ‰ V-SMOW. At the monthly level, the isotopic variability, from one month to another, is greater reaching a range of 6 to 8‰ V-SMOW and it can exceed 10‰ in a daily basis. The local meteoric line (MWL) defining the oxygen 18 vs deuterium relationship, calculated for the Bizerte station gives ${\delta }^{2}H=7.02\times \delta {}_{ }{}^{18}O+8.3$. This relationship becomes ${\delta }^{2}H=8\times \delta {}_{ }{}^{18}O+13.9$ when recalculated by setting the theoretical equilibrium slope of 8 (in reference to the global rainfall line GMWL ${\delta }^{2}H=8\times \delta {}_{ }{}^{18}O+10$).

This strong D-excess of 13.9 is found in the daily rains by defining two groups, one close to 10 concerning the rains of Atlantic origin and another higher than 12 and up to slightly more than 20 defining rains of Mediterranean or mixed origin occurring mainly between September and November.

The time series of samples taken at the Temple borehole (Figure 13a) shows, over the investigation period, a certain stability going from −4.96 ‰ to −5.55 ‰ (apart from the samples of 21/05 and 18/06 which seem to have been subject to an evaporation process probably due to poor preservation and storage marked by isotopic enrichment and D-excess below 10). This low temporal isotopic variation shows a significant mixing of waters on a time scale greater than one year.

The isotopic content is consistent with that found in the two long-term stations Bizerte and Tunis. knowing that the latter stations are coastal and of low altitude, and if the conditions can be applied to the area of Zaghouan (massif whose altitude reaches 1200 m), the recharge would be mainly in the low altitudes of the massif, considering the isotopic gradient of altitude (−0.2/100 m for oxygen 18 in the Mediterranean area). To confirm this hypothesis, it would be necessary to make measurements at the base of the massif and in its high part to define the annual contents according to the altitude.

A D-excess ranging between 10.8 and 13.3 demonstrates a predominance of recharge from Atlantic and mixed rainfall. The variability observed also tends to show that the mixture is not totally perfect (Figure 13b). The mixing possibly originates from communications between different reservoirs in the karst more than a direct punctual contribution of rain.

The longer Ain Ayed 3 bis series is more complex to interpret and requires the full year time series to refine the interpretation. The recorded isotopic variability is sufficiently large to show an evolution of water origin over the period, from −5.08‰ to −6.34‰ in oxygen-18 (Figure 13c). A notable fact is that two isotopic groups can be determined, one ranging between −5‰ and −5.8‰ and the other one below −6‰.

This variability also seems to follow a temporal sequence and is also very closely related to the D-excess (R² = 0.989) (Figure 13d). The most depleted values have the highest D-excesses and vice versa (the same similarity is found for the samples issued from the Temple borehole which aligns with the same line). Note also that the least depleted isotopic pole corresponds to a low altitude recharge while the other pole below −6‰ would seem to be related to a higher altitude input. It is thus important to know if there is a direct influence of more or less depleted rainfall events on the water reservoir in the karst.

In the absence of a rainfall isotopic time series sampled over the same area and period, it is difficult to decide. However, the periods where a depleted signal is recorded seem to correspond to the months of heaviest rainfall more likely to infiltrate (colored in green in Figure A3). Another hypothesis would be the mixing of different reservoirs within the karst which could explain this excellent correlation between isotopes and D-excess, and which corresponds in hydrodynamic terms to pressure loading.

4.3.3. Dating

Preliminary results give some indications about the circulation in the karst system but obviously these few results are local, and they cannot be generalized in the whole area.

• Results for Tritium (³H)

Two samples were analyzed (Ecomusée and Temple) and showed very close values (1.7 ± 0.6 TU and 1.6 ± 0.4 TU), that confirms recent circulation flow in the karst which is an obvious element with respect to present recharge in the zone. In the past, with the nuclear test bomb period in atmosphere (1950–1965), ³H has been widely used to estimate recent groundwater residence times. The return of its natural atmospheric content since the end of 1990’s limits the use of ³H to quantify with precision but allows to estimate the presence of young water flow between present time and some centuries depending on the hydrogeological model used (Piston flow or mixing model representing the possible range values). To use this tool, the first step was to reconstruct concentration time series in precipitation during the last 70 years period impacted by anthropogenic input in Tunisia. In North Tunisia this work was done by [58] who have compiled regional ³H data [59] from different Mediterranean stations to build a local rainfall ³H content from 1950 to 2013 including Tunis and Sfax data between 1968 and 2013. To extend this chronicle until 2022 because the lack of data since 2016, it was considered the period between 1999–2015 period where the annual ³H content is in natural steady state (5.22–2.86 TU) to calculate an annual mean of 3.87 TU reported from 2016 to 2022 period.

Two models were used to estimate the dating of the water: the piston flow model and the mixing model [58]. With the Piston model, ³H concentration in groundwater at the time of measurements in 2022 (A₂₀₂₂) is related to the ³H input of precipitation during the preceding year via (1).

$$A_{2022}=P_i{e}^{-kt}$$

(1)

where P_i is the annual mean ³H concentration in precipitation for year i; t is the transit time between year i and 2022 and k is the decay constant of ³H (0.0566/yr).

For samples with values of 1.6 and 1.7 TU in groundwater the estimated age would be 12 or 65 years, the 2 possible period related to anthropic peak of Tritium.

The mixing model used [60] assumes homogeneous tritium deposition over the entire aquifer at any given time, and a constant water reserve each year. The model simulates the aquifer ³H content depending on the renewal rate, which varies according to annual rainfall (2):

$$A_i=\left(1-a_i\right)A_{i-1}{e}^{-k}+a_i{P}_i$$

(2)

with A_i the ³H content in the aquifer during year i, a_i the renewal rate for the year i, P_i the rain ³H content during year i and k the decay constant of ³H.

The ³H content indicates renewal rate between 1.2–1.3% then close to 80 years.

• Results for Carbon-14

Only one sample of carbon-14 was made on the borehole close to Ain Ayed 3 bis (Ecomusée), which shows a participation of much older water (Ecomusée A = 63.4% and 13C = −10.45‰ vs V-PDB). Before calculating the age of water by the classical decay Equation (3), it is necessary to estimate the initial activity (A₀), taking into account the origin of carbon in dissolved bicarbonates, (i) HCO₃⁻ from CO₂ soil gas-biologic pole (A₀ = 100%) and (ii) the mineral pole from rock carbonates. The Gonfiantini mixing model (3) was used [61].

$$A=A_0{e}^{-\lambda t}$$

(3)

With A activity of water sample, A₀ initial activity, $\lambda$ decay constant of 14-carbon, t decay time.

$$A_0={\displaystyle \frac{\left({}_{ }{}^{13}{C}_{\mathrm{T}\mathrm{D}\mathrm{I}\mathrm{C}}-{}_{ }{}^{13}{C}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d}}\right)}{\left({}_{ }{}^{13}{C}_{\mathrm{g}\mathrm{a}\mathrm{z}}-{\epsilon }_{gas-\mathrm{H}\mathrm{C}\mathrm{O}3}\text{-}{}_{ }{}^{13}{C}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d}}\right)}}$$

(4)

with ¹³C_TDIC, ¹³C_solid, ¹³C_gas, ε_gaz-HCO3, respectively carbon-13 of dissolved bicarbonates, carbonate rock = 0‰, CO₂ in the root zone = −20‰ in Mediterranean area, isotope enrichment factor (8.5‰) between bicarbonates and gas depending on annual temperature (19 °C) which gives an age of 2975 years BP. The hydrodynamics of the karst by classical methods does not seem to show a high velocity and the first interpretation of dating with ³H and ¹⁴C confirms it. According to the isotopic data, it appears that the system functions as water reservoirs with highly mixed ancient water and a small amount of newer water that is currently feeding the borehole.

4.4. Vulnerability Mapping and Nitrates

4.4.1. COP Map

The vulnerability map for COP (Conductivity, Overlying layers, and Precipitation) is derived through the multiplication of three individual layers representing the C, O, and P factors, as illustrated in Figure 14. This map provides a comprehensive depiction of the vulnerability of the Jebel Zaghouan karst groundwater system. Upon analysis, it becomes evident that approximately 80% of the area falls within the category of moderate vulnerability. This indicates a noteworthy susceptibility of the karst aquifer to potential contaminant infiltration or over-exploitation, highlighting areas where careful management strategies may be particularly crucial for maintaining water quality and sustainability. This map needs to be improved by better integrating scenario 1 of factor C based on more field visits and remote sensing. These field campaigns should also improve the knowledge of land use and geology. The introduction of the thickness of the area could also be improved. The Sidi Medien Spring used by the local rural population has the highest nitrate concentrations and should be monitored (especially during wet years) by the SONEDE laboratory.

4.4.2. Nitrates

Measured nitrates in different springs and boreholes were mapped for two seasons 2012 and 2022. In 2012 samples were taken from natural springs (Ain el Guelb and Ain Ayed), galleries (Gallerie 44 and Gallerie 47) and two boreholes, whereas in 2022, only the spring of Sidi Medien was still flowing.

One can note that over the time series of 2004–2022, the maximum nitrates concentrations for boreholes and natural springs (except Sidi Medien Spring) were registered in June 2012 (Figure 15). Concentrations ranged between 5.8 mg/L (Ain El Guelb) to 11.8 mg/L (Gallerie 44). In fact, from September 2011 to august 2012, the registered rainfall reached 906.5 mm (537.1 mm from September to December 2011) and snow events were recorded. Inversely, the campaign of 2022 (considered as dry years within 2021, rainfall 224 mm and 248 mm respectively) showed lower values of nitrates, they ranged from 2.5 mg/L (Ecomusée borehole) to 5.7 mg/L for boreholes. These results confirm the impact of rainfall and the unsaturated zone thickness on groundwater quality [62,63] and as in 2022, the piezometric level continued to decline reaching more than 100 m.

The highest value of nitrates (19 mg/L) was registered in 2016 and 2022 in Sidi Medien Spring which is unfortunately not controlled by SONEDE. Unlike the Jebel Zaghouan area, mostly covered with forest and bare lands, the region of Sidi Median is known for agricultural activities (including traditional livestock grazing). Although the nitrate concentration in this spring does not exceed the limits set by Tunisian Standards (NT, 45 mg/L) or the WHO (50 mg/L), it is recommended to monitor and regulate its use for drinking purposes by the rural population [64].

By employing a combination of investigation methods, including the APLIS method for recharge estimation, stable isotope analysis, and vulnerability mapping, this study offers a first and comprehensive understanding of the aquifer dynamics. The results reveal a severe decline in groundwater levels since 2002, with a particularly sharp acceleration noted from 2012 onwards, due to overexploitation linked to urban expansion. Recharge estimation indicates moderate infiltration rates, consistent with previous models, while isotopic analysis shows recharge occurring primarily at lower altitudes with mixed rainfall sources. Tritium and Carbon-14 dating reveal a complex mixture of young and ancient water within the aquifer, enhancing the understanding of its internal circulation. The vulnerability map highlights moderate contamination risks for 80% of the area, underscoring the need for protective measures. Additionally, nitrate levels, influenced by rainfall and agricultural activities, emphasize the importance of continuous monitoring, especially for the Sidi Medien Spring.

5. Conclusions

This investigation conducted a comprehensive assessment of the Jebel Zaghouan aquifer, a critical water source in Northeast Tunisia, addressing groundwater availability and quality. The research also mapped the aquifer’s vulnerability to pollution. Data from nine boreholes and galleries were analyzed to assess historical trends and aquifer sustainability. Recharge rates were evaluated using the APLIS method, and physicochemical analyses from 2004 to 2021 were conducted to understand groundwater quality variations. Stable isotopes, tritium, and Carbon-14 were analyzed to understand the hydrodynamic system and aquifer functioning. Vulnerability mapping was based on the COP method. Emphasizing ongoing monitoring, this investigation aimed to provide decision-makers with insights to improve karstic system management, ensuring long-term sustainability.

The findings reveal significant challenges and opportunities for sustainable management strategies. The analysis of water availability underscores the severe depletion of groundwater levels in the Zaghouan region, primarily due to overexploitation linked to urban expansion. Despite occasional wet periods, the groundwater remains severely depleted, with little to no recharge observed. Despite measurement uncertainties, the APLIS method gave consistent recharge rates, emphasizing the need for further validation to inform effective aquifer management and groundwater protection.

The analysis of groundwater recharge, observed discharge, and safe yield for the period corresponding to the natural functioning of springs highlights the critical influence of rainfall on water resources. In an average year with 467 mm of rainfall, the recharge estimated using the APLIS method is 3.8 million m³/year, while the observed discharge is 3.6 million m³/year, resulting in a modest but positive sustainable yield of 0.2 million m³/year. During a wet year with 867 mm of rainfall, the recharge rises significantly to 7.0 million m³/year, exceeding the observed discharge of 6.5 million m³/year and yielding a surplus of 0.5 million m³/year. Conversely, in a dry year with only 239 mm of rainfall, the recharge drops to 1.9 million m³/year, equaling the observed discharge and leaving no surplus, resulting in a sustainable yield of 0.0 million m³/year.

These findings underscore that groundwater recharge is strongly dependent on rainfall, with a recharge surplus available for sustainable use in wet and average years. However, dry years present a critical scenario where recharge and discharge are balanced, leaving no margin for additional extraction. This variability emphasizes the importance of adaptive groundwater management strategies that consider climatic fluctuations, particularly in dry periods, to prevent overexploitation and ensure the long-term sustainability of water resources.

Isotopic analysis highlights the importance of comparing aquifer time series with local or regional rainfall chronicles. The observed isotopic content aligns with findings from long-term stations, suggesting primarily lower altitude recharge. However, comprehensive data collection and continuous monitoring are essential to refine interpretations and improve understanding of karst hydrodynamics. Dating analysis provides insights into the circulation within the karst system, indicating recent flow based on Tritium analysis and the presence of much older water based on Carbon-14 analysis. These findings support the interpretation that the karst operates with a mixture of ancient and newer water, contributing to borehole supply. Finally, vulnerability mapping highlights significant risk to potential contamination or over-exploitation, emphasizing the need for careful management strategies to ensure water quality and sustainability. Further enhancements and monitoring efforts are recommended to refine the accuracy of vulnerability assessments and inform effective management practices.

For the perspectives, future strategies should focus on implementing integrated monitoring programs, including isotopic analysis, nitrate mapping, and groundwater level monitoring, to gather comprehensive data. It is crucial to enhance the geological understanding of the reservoir through geophysical investigations and geological modeling. Additionally, remote sensing serves as a valuable tool for identifying and monitoring karst features. With ongoing advancements in technology, data analysis, and collaborative research, its effectiveness in studying karst areas is expected to further improve.

Furthermore, enhancing stakeholder engagement with local communities, water management authorities, and research institutions is crucial for effectively designing and implementing management strategies. Adopting adaptive management approaches will allow for timely adjustments based on ongoing data analysis and evolving hydrological conditions. Additionally, addressing conflicts arising from water use and the impact of urbanization, which has led to an unplanned rise in water consumption demand, is essential.