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

Emna Gargouri-Ellouze, Fairouz Slama, Samiha Kriaa, Ali Benhmid, Jean-Denis Taupin, Rachida Bouhlila

► To cite this version:

Emna Gargouri-Ellouze, Fairouz Slama, Samiha Kriaa, Ali Benhmid, Jean-Denis Taupin, et al.. Comprehensive Assessment of the Jebel Zaghouan Karst Aquifer (Northeastern Tunisia): Availability, Quality, and Vulnerability, in the Context of Overexploitation and Global Change. Water, 2025, 17 (3), pp.407. 10.3390/w17030407. [https://doi.org/10.3390/w17030407](https://doi.org/10.3390/w17030407)

HAL Id: hal-05081134

https://hal.science/hal-05081134v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

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

Emna Gargouri-Ellouze <fnref n="1" />,* ID, Fairouz Slama <fnref n="1" /> ID, Samiha Kriaa <fnref n="1" />,<fnref n="2" />, Ali Benhmid <fnref n="1" />, Jean-Denis Taupin <fnref n="3" />,*, ID
and Rachida Bouhlila <fnref n="1" />

<figure><img src="image_2.png" /></figure>

check for
updates

Academic Editor: Chin H Wu

Received: 9 December 2024

Revised: 21 January 2025

Accepted: 24 January 2025

Published: 1 February 2025

**Citation:** Gargouri-Ellouze, E.; Slama, F.; Kriaa, S.; Benhmid, A.; Taupin, J.-D.; Bouhlila, R. Comprehensive Assessment of the Jebel Zaghouan Karst Aquifer (Northeastern Tunisia): Availability, Quality, and Vulnerability, in the Context of Overexploitation and Global Change. *Water* **2025**, *17*, 407. https://doi.org/10.3390/w17030407

**Copyright:** © 2025 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license
(https://creativecommons.org/
licenses/by/4.0/).

**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.

**Keywords:** karst aquifer; vulnerability mapping; stable isotopes; APLIS method; SPEI; SGI

## 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-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

<footnote n="1">¹ Laboratory of Modelling in Hydraulics and Environment (LMHE), LR99ES19, National Engineering School of TUNIS (ENIT), University of Tunis El Manar, Tunis 1002, Tunisia; fairouz.slama@enit.utm.tn (F.S.); samiha.kriaa@enit.utm.tn (S.K.); ali.benhmid@enit.utm.tn (A.B.); rachida.bouhlila@enit.utm.tn (R.B.)</footnote>

<footnote n="2">² Regional Agricultural Development Commission (CRDA), Zaghouan 1100, Tunisia</footnote>

<footnote n="3">³ HydroSciences, University of Montpellier, CEDEX 05, 34093 Montpellier, France</footnote>

<footnote>* Correspondence: emna.gargouri@enit.utm.tn (E.G.-E.); jean-denis.taupin@umontpellier.fr (J.-D.T.); Tel.: +216-23356050 (E.G.-E.); +33-0467149037 (J.-D.T.)</footnote>
(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:
<table><caption>**Table 1.** Summary of Methods and Applications in Karst Hydrogeology.</caption><thead><tr><td>Method</td><td>Application</td></tr></thead><tbody><tr><td>Geological</td><td>Elucidates the aquifer's geometry and flow paths.</td></tr><tr><td>Geophysical</td><td>Identifies geological structures and fractures.</td></tr><tr><td>Speleological</td><td>Maps conduit networks, providing direct insights.</td></tr><tr><td>Hydrological and Hydraulic</td><td>Establish water balances and determine hydraulic properties.</td></tr><tr><td>Hydrochemical and Isotopic</td><td>Trace water origins, movements, and interactions.</td></tr><tr><td>Artificial Tracers</td><td>Delineate connections and flow velocities [20]</td></tr></tbody></table>

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–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–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, Tunisia) 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.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Location of the Jebel Zaghouan karst aquifer.</figcaption></figure>

### 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 19 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.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Geological context of Jebel Zaghouan.</figcaption></figure>

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.
<figure><img src="image_5.png" /><figcaption>Figure 3. Location map of sampling boreholes and springs.</figcaption></figure>

# 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.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Zaghouan springs production (Million cubic meters per year).</figcaption></figure>

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 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.

<table><caption>**Table 2.** Meteorological stations characteristics.</caption><thead><tr><th rowspan="2">Meteorological Station</th><th rowspan="2">Observation Period of Data Series</th><th colspan="2">Location (WGS 84)</th><th rowspan="2">Mean Annual Precipitation (mm)</th><th rowspan="2">Altitude (m)</th></tr><tr><th>X</th><th>Y</th></tr></thead><tbody><tr><td>Zaghouan PF</td><td>1961 to 2017</td><td>10.129</td><td>36.425</td><td>470</td><td>130</td></tr><tr><td>Mograne CSA SM</td><td>1971 to 2020</td><td>10.090</td><td>36.432</td><td>485</td><td>151</td></tr><tr><td>Zaghouan SM</td><td>1909 to 2016</td><td>10.153</td><td>36.403</td><td>481</td><td>165</td></tr><tr><td>Zaghouan DRE</td><td>1976 to 2017</td><td>10.144</td><td>36.403</td><td>449</td><td>238</td></tr><tr><td>Zaghouan Contrôle</td><td>1915 to 1951</td><td>10.149</td><td>36.395</td><td>497</td><td>241</td></tr><tr><td>Zaghouan Sidi Bou Gabrine</td><td>1991 to 2017</td><td>10.116</td><td>36.375</td><td>492</td><td>677</td></tr><tr><td>Zaghouan Poste Optique</td><td>1912 to 2008</td><td>10.140</td><td>36.441</td><td>516</td><td>945</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Average annual precipitation.</figcaption></figure>

### 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 Migraine 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 Gdhir 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.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Conductivity and Temperature of groundwater samples.</figcaption></figure>

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–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.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) calibration period (1915–1927) water budget, (b) validation period (1970–1995) water budget.</figcaption></figure>

### 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–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><img src="image_10.png" /><figcaption>**Figure 8.** (a) Temporal evolution of the static level at the piezometer, monthly produced amounts and (b) monthly rainfall.</figcaption></figure>

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.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Variation of the monthly calculated SPEI and SGI from 2002 to 2020.</figcaption></figure>

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). 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).

<table><caption>**Table 3.** Correlation between SPEI and SGI.</caption><thead><tr><td>Observation Period</td><td>Kendall's Tau</td><td>p-Value</td></tr></thead><tbody><tr><td>2002–2011</td><td>0.157</td><td>0.01844</td></tr><tr><td>2011–2020</td><td>0.066</td><td>0.291</td></tr></tbody></table>

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%).
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Groundwater recharge distribution by the APLIS.</figcaption></figure>

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).

<table><caption>**Table 4.** Yield of groundwater.</caption><thead><tr><th></th><th>Rainfall (mm/Year)</th><th>Observed Discharge (Million m³/Year)</th><th>Recharge (APLIS) Million m³/Year</th><th>Yield <sup>1</sup> APLIS (Million m³/Year)</th></tr></thead><tbody><tr><td>Average year</td><td>467</td><td>3.6</td><td>3.8</td><td>0.2</td></tr><tr><td>Wet year</td><td>867</td><td>6.5</td><td>7.0</td><td>0.5</td></tr><tr><td>Dry year</td><td>239</td><td>1.9</td><td>1.9</td><td>0.0</td></tr></tbody></table>

Note: <sup>1</sup> https://gw-project.org/books/a-glossary-of-hydrogeology/, accessed on 5 January 2025.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Recharge in mm for a median year of precipitation.</figcaption></figure>

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.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Piper diagram.</figcaption></figure>

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). 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 daily isotopic rainfall records in Bizerte (2009–2014) analyzed at LAMA Hydrosciences laboratory (coll. with ISTEUB-Institut Supérieur des Technologies de l'Environnement de l'Urbanisme et du Bâtiment) that will serve as a reference for comparison with the karst waters of Zaghouan.

In this station 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. 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$ [58]. 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$) [58].

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 [58].

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.
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

**Figure 13.** (a) Stable isotope content for the Temple and Ain Ayed bis 3 boreholes and the global meteoric water line (GMWL); (b) Temporal evolution of D-excess for Ain Ayed 3 bis borehole; (c) Temporal evolution of oxygen-18 for Ain Ayed 3 bis borehole; (d) D-excess vs. oxygen-18 for Ain Ayed 3 bis borehole.

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^2$ = 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 ($^{3}$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), $^{3}$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 $^{3}$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 $^{3}$H data [59] from different Mediterranean stations to build a local rainfall $^{3}$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 $^{3}$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, $^{3}$H concentration in groundwater at the time of measurements in 2022 ($A_{2022}$) is related to the $^{3}$H input of precipitation during the preceding year via (1).

$$A_{2022} = P_i e^{-kt} \quad (1)$$

where $P_i$ is the annual mean $^{3}$H concentration in precipitation for year $i$; $t$ is the transit time between year $i$ and 2022 and $k$ is the decay constant of $^{3}$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 $^{3}$H content depending on the renewal rate, which varies according to annual rainfall (2):

$$A_i = (1 - a_i) A_{i-1} e^{-k} + a_i P_i \quad (2)$$

with $A_i$ the $^{3}$H content in the aquifer during year $i$, $a_i$ the renewal rate for the year $i$, $P_i$ the rain $^{3}$H content during year $i$ and $k$ the decay constant of $^{3}$H.

The $^{3}$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_0$), taking into account the origin of carbon in dissolved bicarbonates, (i) $HCO_3^-$ from $CO_2$ soil gas-biologic pole ($A_0$ = 100%) and (ii) the mineral pole from rock carbonates. The Gonfiantini mixing model (3) was used [61].

$$A = A_0 e^{-\lambda t} \quad (3)$$
With A activity of water sample, A₀ initial activity, λ decay constant of 14-carbon, t decay time.

$$A_0 = \frac{(^{13}\text{C}_{\text{TDIC}} - ^{13}\text{C}_{\text{solid}})}{(^{13}\text{C}_{\text{gas}} - \varepsilon_{\text{gas-HCO}_3} - ^{13}\text{C}_{\text{solid}})} \quad (4)$$

with ¹³Cₜₔ₀, ¹³C₠ₜₔ, ¹³Cₚₜₔ, εₚₜₔ-ₚₜₔ, 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.
<figure><img src="image_19.png" /><figcaption>**Figure 14.** C, O and P factors and COP vulnerability maps of for Jebel Zaghouan.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 15.** Nitrates value maps for the campaigns of 2012 and 2022.</figcaption></figure>

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 Medien 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.

**Author Contributions:** Conceptualization, E.G.-E., F.S., R.B. and J.-D.T.; methodology, E.G.-E., F.S., R.B., J.-D.T., S.K. and A.B.; software, E.G.-E. and F.S.; validation, E.G.-E., F.S. and J.-D.T.; investigation and data curation, E.G.-E., F.S., S.K., A.B. and J.-D.T.; writing—original draft preparation, E.G.-E., F.S. and J.-D.T.; writing—review and editing, E.G.-E., F.S., J.-D.T. and R.B.; supervision, R.B.; project administration, R.B. and F.S.; funding acquisition, R.B. and F.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the General Directorate of Scientific Research in the framework of PRIMA-KARMA project, grant agreement number 01DH19022A.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding authors. The data are not publicly available because they are currently used in theses.

**Acknowledgments:** This Research has been supported by the European Commission through the Partnership for Research and Innovation in the Mediterranean Area (PRIMA) program under Horizon 2020 (KARMA project, grant agreement number 01DH19022A, funded by General Directorate of Scientific Research for the Tunisian team). We extend our gratitude to Regional Agricultural Development Commission (CRDA), Zaghouan, for their indispensable support and collaboration. We would like to acknowledge the National Water Exploitation and Distribution Company (SONEDE), for their invaluable support and contributions to this research, particularly the Zaghouan Water Supply Service (SAZ) for their collaboration and assistance in the field, and the central laboratory of the Ghdir El Golla water complex for the water quality analysis.
**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

## Appendix A

<figure><img src="image_21.png" /><figcaption>**Figure A1.** Example of non-digitized measured hydrograph at Nymphée and Ain Ayed springs from 1929 to 1930.</figcaption></figure>

<table><caption>**Table A1.** Chemical analysis of water samples (Concentrations in mg/L), Electrical Conductivity (EC) in µS/cm, pH, and Saturation indices.</caption><thead><tr><th rowspan="2">Well/Spring</th><th rowspan="2">T °C</th><th rowspan="2">pH</th><th rowspan="2">EC (µS/cm)</th><th colspan="7">(mg/L)</th><th colspan="4">Saturation Indices</th></tr><tr><th>Ca</th><th>Mg</th><th>Na</th><th>K</th><th>HCO<sub>3</sub></th><th>Cl</th><th>SO<sub>4</sub></th><th>NO<sub>3</sub></th><th>Is Calcite</th><th>Is Dolomite</th><th>Is Gypsum</th><th>Is Anhydrite</th></tr></thead><tbody><tr><td>Pristine</td><td>22.3</td><td>7.22</td><td>1244</td><td>175</td><td>27</td><td>40</td><td>0.98</td><td>216</td><td>59</td><td>317</td><td>5.7</td><td>0.28</td><td>0.09</td><td>-0.65</td><td>-0.87</td></tr><tr><td>Sidi Medien</td><td>22</td><td>7.23</td><td>518</td><td>77</td><td>13</td><td>21</td><td>2.4</td><td>192</td><td>38</td><td>38</td><td>19</td><td>-0.03</td><td>-0.49</td><td>-1.82</td><td>-2.04</td></tr><tr><td>Ain Ayed 3 bis</td><td>21.8</td><td>7.99</td><td>1092</td><td>161</td><td>32</td><td>50</td><td>1.15</td><td>188</td><td>66</td><td>386</td><td>3.8</td><td>0.94</td><td>1.54</td><td>-0.61</td><td>-0.83</td></tr><tr><td>Temple</td><td>25</td><td>7.6</td><td>1178</td><td>163</td><td>31</td><td>70</td><td>1.06</td><td>193</td><td>91</td><td>380</td><td>2.9</td><td>0.57</td><td>0.76</td><td>-0.62</td><td>-0.84</td></tr><tr><td>Ecomusée</td><td>21.7</td><td>7.76</td><td>1091</td><td>167</td><td>32</td><td>43</td><td>1.13</td><td>163</td><td>63</td><td>409</td><td>2.5</td><td>0.66</td><td>0.97</td><td>-0.57</td><td>-0.79</td></tr><tr><td>Cristalline</td><td>22.2</td><td>7.45</td><td>924</td><td>137</td><td>28</td><td>41</td><td>1.07</td><td>236</td><td>65</td><td>227</td><td>3.4</td><td>0.46</td><td>0.58</td><td>-0.88</td><td>-1.1</td></tr></tbody></table>
<figure><img src="image_22.png" /><figcaption>Figure A2. Main factors playing a role in the definition of the COP map [44].</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure A3. Monthly rainfall for 2021 and 2022 (orange dot line: threshold of 20 mm, green bars over the threshold).</figcaption></figure>
# Appendix B

## Appendix B.1. APLIS Method

The digital elevation model map (DEM) from Shuttle Radar Topography Mission (SRTM) [33] with 30 m resolution was used for the mapping of altitude and slope. Taking into account the topography, the geological map is grouped into lithological groups with similar hydrogeological characteristics. As part of the project “Map of Water Resources of Tunisia-CRET (Carte des Ressources en Eau de Tunisie)” managed by the Directorate General of Water Resources in 2016, the homogenization of geological maps of the National Office of Mines at the scale 1/500,000, 1/100,000 and 1/500,000 and other petroleum cartographic documents allowed the development of the geological map at the scale 1/200,000.

For the APLIS method, the recharge rates were calculated using QGIS, using this Expression (A1):

$$R = (A + P + 3L + 2I + S)/0.9 \quad (A1)$$

where:

A = altitude, P = slope, L = lithology, I = infiltration landforms, S = soil type.

Since Tunisia is at the same time a Mediterranean and a Saharan country the soils show all the signs of this climatic, morphological and geological diversity. According to the French system of soil classification, the soils of Tunisia are classified as podzols, vertisoils, red Mediterranean soils, calcic-magnesic soils (dominant soils), brown and isohumic soils, saline and hydromorphic soils and also poorly evolved soils. The study area is dominated by complex soils. It includes four soil groups: Cambisols, Leptosols, Arenosols and xerosols, Calcareous regosols and fluvisols [34].

The DEM is introduced into QGIS software, and the altitude map is extracted as output after operation. The height variable is divided into ten series at 300 m intervals, and each level is assigned 1 to 10 points [32]. The altitude of the study area rose from 267 m in the east to 1272 m in the southeast. Based on altitude, only five of ten possible APLIS altitude categories are provided.

The map of slopes in the study area was produced in QGIS software using the DEM, and then was ranked in nine different sequences. The scores assigned to the slope parameters decrease as slope gets steeper which means that increase in slope leads to decrease in groundwater recharge. Approximately 45% of the area has slopes at angles greater than 46%. Only 13% of the study area has slopes less than 21%.

The lithology (L) has been categorized based on the existing local geological, geomorphological, hydrogeological maps and stratigraphic studies. Limestone and dolomite karstified will be considered as the dominant rock formations in the area.

The Infiltration Landform (I) item is set by assigning priority infiltration forms to areas observed during field work. Two rankings 5 and 1 were assigned to Preferential Infiltration variable. Points 1 and 5 are awarded to areas with minimum and maximum Preferential Infiltration, respectively. This assumption is based on observations and investigations on the field. In Mhimi's investigation (2020) [35], findings indicate a relatively lower occurrence of dolines in Jebel Zaghouan, the southwestern sector of the study. The primary concentration of dolines is observed in the area defined by the effluents of the hydrographic network of the Zriba plain, particularly within the southern regions of Jebel el Guebli and Jebel el Leri.

The Soil (S) parameter values were assigned based on a site soils map using soil taxonomy classifications [36]. Four types of soils (S) were classified in the study area: (1) Cambisols, (2) Calcareous regosols and fluvisols (3) Arenosols and xerosols, (4) Leptosols.
## Appendix B.2. COP Method

### Appendix B.2.1. The Overlying Layers: O Factor

Groundwater is protected by the characteristics and thickness of the unsaturated zone [42]. Four major layers make up the overlying: topsoil, non-karstic rocks, subsoil, and unsaturated karstic rocks [43]. But only the lithological layers ($O_L$) and soils ($O_S$) are considered to be of hydrogeological significance. Figure A3 demonstrates how the soil subfactor considers texture, particle size, and thickness. The lithology subfactor considers the thickness of each layer (m), the restricted circumstances, and the geological lithology at the surface ($l_y \times c_n$). In Figure A3, the value for each category is displayed. To produce the O factor for Jebel Zaghouan the pedology and geology maps from the national information system "carte agricole" (PNCA 2004, Programme National Carte Agricole. Société Tunisienne d'Ingénierie), along with nine available borehole and one piezometer completion reports, provided information about the thickness of the unsaturated zone.

Os factor is composed of two layers related to texture and thickness. In Jebel Zaghouan, three texture classes (Loam, silt and clay) and three soil depth classes (very shallow, shallow and deep) were identified. Using the pedology map and the borehole completion reports, 6 lithology types were identified in the study area (Table A2).

<table><caption>**Table A2.** Values assigned to the lithologies identified in the studied catchment area (where "ly" refers to the type of lithology and "cn" refers to the degree of confinement according to the COP method).</caption><thead><tr><th>Lithology</th><th>Nomenclature in the Pedology Map</th><th>Value ly</th><th>Value cn</th></tr></thead><tbody><tr><td>Marls</td><td>102</td><td>1000</td><td>1.5</td></tr><tr><td>Non fissure limestones rocks</td><td>105</td><td>1000</td><td>1.5</td></tr><tr><td>Sands and gravels</td><td>108</td><td>10</td><td>1.5</td></tr><tr><td>limestones encrustation</td><td>119</td><td>500</td><td>1.5</td></tr><tr><td>Gravels</td><td>120</td><td>30</td><td>1.5</td></tr><tr><td>Gravels</td><td>122</td><td>30</td><td>1.5</td></tr><tr><td>Fissured carbonated rocks</td><td>102,105</td><td>3</td><td>1.5</td></tr></tbody></table>

The $O_L$ map is obtained by applying the following Formula (A2):

$$O_L = [l_y \times m] \times c_n \qquad (A2)$$

Finally, the O map is obtained by summing the $O_S$ and the $O_L$ maps.

### Appendix B.2.2. The Concentration of Flow: C Factor

The C factor is used to describe the groundwater's reaction to concentrated infiltration. The level of protection offered by the unsaturated zone during recharge situations is determined differently. The following two scenarios are possible:

* Scenario 1 involves a karstic catchment with limited protection against rapid infiltration. Concentrated surface runoff near the base of slopes or quick infiltration through swallow holes are the main sources of recharge. The Jebel Zaghouan case study was analyzed to gather information about potential sinkholes or swallow holes using inquiries with the water authority and a local speleology expert from the AREZ association (https://arez.tn/, accessed on 6 February 2023). However, it was found that the caves and holes, when placed within the hydrological network did not seem clearly connected to rapid groundwater infiltration. Tracer tests could not be performed due to the drought of groundwater, so scenario 1 was not considered. Further field
and remote sensing investigations are required to identify these infiltration regions and mechanisms.

* Scenario 2 is different from the focused recharge through swallow holes and occurs in the zone of rest for autogenic recharge. The C factor is comprised of three variables: surface features (*sf*), slope (*s*), and vegetation (*v*). In this context, the *sv* factor is considered in contrast to scenario 1.

## Appendix B.2.3. Slope and Vegetation

The slope is extracted with QGIS from DEM in percentage, and reclassified into 4 categories (≤8%, 8% < *s* < 31%; 31% < *s* < 76%; *s* > 76%), which were assigned weights accordingly. Based on the land use maps from "carte Agricole", land use was divided into two main types, mainly "no vegetation" (including bare soils and rocks, grassland, or low sparse vegetation) and "vegetation" [44].

Considering Scenario 2, the slope and vegetation parameter *sv*, is applied in the opposite sense of scenario (1) since runoff is the major factor preventing water from seeping to the water table instead of flowing through swallow holes as the slope is higher and vegetation is lacking.

## Appendix B.2.4. Surface Features

The surface layer *sf* parameter indicates whether karst features are present or not. For our study case the *sf* value was assigned to that corresponding to the fissured carbonates in a permeable surface layer (=0.75). The final C score is obtained by multiplying *sf* and *sv*.

## Appendix B.2.5. The Precipitation Factor: P

The *P* Factor represents the climatic conditions in the catchment area. It is calculated by adding two sub-factors, *P*<sub>Q</sub> and *P*<sub>I</sub>, which determine the amount and frequency of annual precipitation respectively. *P*<sub>Q</sub> reflects the quantity of precipitation, with a range of 0.2 to 0.4, while *P*<sub>I</sub> assesses the intensity of precipitation, which is the ratio of the amount of precipitation and the number of rainy days. This sub-factor has a range of 0.2 to 0.6. The *P* factor considers that the higher the precipitation and its intensity, the more susceptible the studied area is to recharge and potential vulnerability. The influence of precipitation volume and annual recharge on groundwater vulnerability is described by the *P*<sub>Q</sub> subfactor. It corresponds to a historical series of wet years' mean annual precipitation [42]. *P*<sub>Q</sub> map was performed based on the precipitation data of wet years.

Ordered quantile normalization (ORQ) enabled wet years to be identified for each weather station [45,46] and SPI approach. ORQ is a conversion method to transform a vector with any distribution into one that fits to a normal (Gaussian) distribution. Positive transformed data indicates years with precipitation greater than the median, while negative values indicate years with precipitation less than the median. Using the SPI approach, a year is classified as wet if the value of the transformed data is greater than 1. Then, applying the Standard Gaussian quantile function (the inverse of the Standard Gaussian cumulative distribution) derived the threshold defining a wet year. Therefore, the wet years' mean annual precipitation ranges, for wet years from 710 to 857 mm, for the stations located in the study area (Table A3).
<table><caption>**Table A3.** Wet years' mean annual precipitation and altitude.</caption><thead><tr><th>Meteorological Station</th><th>Wet Years' Mean Annual<br>Precipitation (mm)</th><th>Altitude (m)</th></tr></thead><tbody><tr><td>Zaghouan PF</td><td>783</td><td>130</td></tr><tr><td>Mograne Csa SM</td><td>710</td><td>151</td></tr><tr><td>Zaghouan SM</td><td>780</td><td>165</td></tr><tr><td>Zaghouan DRE</td><td>790</td><td>238</td></tr><tr><td>Zaghouan Contrôle</td><td>785</td><td>241</td></tr><tr><td>Zaghouan Sidi Bou Gabrine</td><td>857</td><td>677</td></tr><tr><td>Zaghouan Poste Optique</td><td>808</td><td>945</td></tr></tbody></table>

The precipitation/altitude gradient was calculated using (A3).

$$P(Z) = P(130) + (Z - 130) \times \frac{9}{100} \quad (A3)$$

where $P(Z)$ is the annual precipitation at an altitude $Z$ on the catchment in meters (Figure A4); $P(130)$ is the reference precipitation, 9 mm is the precipitation gradient per 100 m elevation. A raster layer of precipitation distribution was created from the trimmed DEM using the altitude variation, as per (A3).

<figure><img src="image_24.png" /><figcaption>**Figure A4.** Precipitation spatial distribution.</figcaption></figure>

$P_O$ values were given to classes based on the average of the wet years' mean annual precipitation of rain gauges located at the base of the mountain.

The intensity of precipitation $P_I$ parameter describes the temporal distribution of precipitation within a time interval. For the number of rainy days, the daily precipitation measurements at high-altitude stations are irregular; therefore, our calculations were based solely on the data from other stations, giving greater importance to the two Zaghouan and Mograne CSA SM stations, since they are monitored by the National Institute of
Meteorology. Similarly, to precipitation, a gradient for the number of rainy days (NRD) as a function of altitude (A4), Figure A5 was established.

$$NRD(Z) = NRD(130) + (Z - 130) \times \frac{16.4}{100} \qquad (A4)$$

<figure><img src="image_25.png" /><figcaption>Figure A5. Number of rainy days according to altitude.</figcaption></figure>

## References

1. Ford, D.; Williams, P. *Karst Hydrogeology and Geomorphology*; John Wiley & Sons: Hoboken, NJ, USA, 2007.

2. Kalhor, K.; Ghasemizadeh, R.; Rajic, L.; Alshawabkeh, A. Assessment of groundwater quality and remediation in karst aquifers: A review. *Groundw. Sustain. Dev.* **2019**, *8*, 104–121. [CrossRef]

3. Water Resources Mission Area. Karst Aquifers. Retrieved from USGS, 2021. Available online: https://www.usgs.gov/mission-areas/water-resources/science/karst-aquifers (accessed on 22 January 2023).

4. Bakalowicz, M. Karst and karst groundwater resources in the Mediterranean. *Environ. Earth Sci.* **2015**, *74*, 5–14. [CrossRef]

5. Nerantzaki, S.D.; Nikolaidis, N.P. The response of three Mediterranean karst springs to drought and the impact of climate change. *J. Hydrol.* **2020**, *591*, 125296. [CrossRef]

6. Lai, G.G.; Padedda, B.M.; Ector, L.; Wetzel, C.E.; Lugliè, A.; Cantonati, M. Mediterranean karst springs: Diatom biodiversity hotspots under the pressure of hydrological fluctuation and nutrient enrichment. *Plant Biosyst.* **2019**, *154*, 673–684. [CrossRef]

7. Siegel, L.; Goldscheider, N.; Petitta, M.; Xanke, J.; Andreo, B.; Bakalowicz, M.; Barberá, J.A.; Bouhlila, R.; Burg, A.; Doummar, J.; et al. Distribution, threats and protection of selected karst groundwater-dependent ecosystems in the Mediterranean region. *Hydrogeol. J.* **2023**, *31*, 2231–2249. [CrossRef]

8. Sivelle, V.; Jourde, H.; Bittner, D.; Mazzilli, N.; Tramblay, Y. Assessment of the relative impacts of climate changes and anthropogenic forcing on spring discharge of a Mediterranean karst system. *J. Hydrol.* **2021**, *598*, 126396. [CrossRef]

9. Fathi, S.; Sjåstad Hagen, J.; Haidari, A.H. Synthesizing existing frameworks to identify the potential for managed aquifer recharge in a karstic and semi-arid region using GIS multi-criteria decision analysis. *Groundw. Sustain. Dev.* **2020**, *11*, 100390. [CrossRef]

10. Stevanović, Z.; Marinović, V.; Krstajić, J. CC-PESTO: A novel GIS-based method for assessing the vulnerability of karst groundwater resources to the effects of climate change. *Hydrogeol. J.* **2021**, *29*, 159–178. [CrossRef]

11. Bagherzadeh, S.; Kalantari, N.; Fadaei, N.A.; Derakhshan, Z.; Conti, G.O.; Ferrante, M.; Malekahmadi, R. Groundwater vulnerability assessment in karstic aquifers using COP method. *Environ. Sci. Pollut. Res.* **2018**, *25*, 18960–18979. [CrossRef] [PubMed]
12. Hssaisoune, M.; Bouchaou, L.; Sifeddine, A.; Bouimetarhan, I.; Chehbouni, A. Moroccan groundwater resources and evolution with global climate changes. *Geosciences* **2020**, *10*, 81. [CrossRef]

13. Ries, F.; Lange, J.; Schmidt, S.; Puhlmann, H.; Sauter, M. Recharge estimation and soil moisture dynamics in a Mediterranean, semi-arid karst region. *Hydrol. Earth Syst. Sci.* **2015**, *19*, 1439–1456. [CrossRef]

14. Valdes-Abellan, J.; Pardo, M.A.; Jodar-Abellan, A.; Pla, C.; Fernandez-Mejuto, M. Climate change impact on karstic aquifer hydrodynamics in southern Europe semi-arid region using the KAGIS model. *Sci. Total Environ.* **2020**, *723*, 138110. [CrossRef] [PubMed]

15. Parisi, A.; Monno, V.; Fidelibus, M.D. Cascading vulnerability scenarios in the management of groundwater depletion and salinization in semi-arid areas. *Int. J. Disaster Risk Reduct.* **2018**, *30*, 292–305. [CrossRef]

16. Richts, A.; Vrba, J. Groundwater resources and hydroclimatic extremes: Mapping global groundwater vulnerability to floods and droughts. *Environ. Earth Sci.* **2016**, *75*, 926. [CrossRef]

17. Šariri, S.; Valić, D.; Kralj, T.; Cvetković, Ž.; Mijošek, T.; Redžović, Z.; Karamatić, I.; Marijić, V.F. Long-term and seasonal trends of water parameters in the karst riverine catchment and general literature overview based on CiteSpace. *Environ. Sci. Pollut. Res.* **2024**, *31*, 3887–3901. [CrossRef]

18. Gammoudi, S.; Chkir, N.; Boughattas, N.; Hamdi, M.; Arraouadi, S.; Zouari, K. Assessment of urban groundwater vulnerability in arid areas: Case of Sidi Bouzid aquifer (Central Tunisia). *J. Afr. Earth Sci.* **2020**, *168*, 103849. [CrossRef]

19. Liu, W.; Jiang, L.; Liu, B.; Liu, R.; Xiao, Z. Monitoring the evolution process of karst desertification and quantifying its drivers in the karst area of Southwest China. *Environ. Sci. Pollut. Res.* **2023**, *30*, 123259–123273. [CrossRef]

20. Dar, F.A.; Jeelani, G.; Perrin, J.; Ahmed, S. Groundwater recharge in semi-arid karst context using chloride and stable water isotopes. *Groundw. Sustain. Dev.* **2021**, *14*, 100634. [CrossRef]

21. Goldscheider, N. Overview of methods applied in karst hydrogeology. In *Karst Aquifers—Characterization and Engineering*, Stevanović, Z., Ed.; Springer International Publishing: Cham, Switzerland, 2015; pp. 127–145. [CrossRef]

22. Andreo, B.; Carrasco, F.; Durán, J.J.; Jiménez, P.; LaMoreaux, J.W.L. *Hydrogeological and Environmental Investigations in Karst Systems*; Environmental Earth Sciences; Springer: Berlin, Heidelberg, 2015. [CrossRef]

23. Redhaounia, B.; Ilondo, B.O.; Gabtni, H.; Khomsi, S.; Bédir, M. Electrical resistivity tomography (ERT) applied to karst carbonate aquifers: Case study from Amdoun, Northwestern Tunisia. *Pure Appl. Geophys.* **2016**, *173*, 1289–1303. [CrossRef]

24. Ayadi, Y.; Mokadem, N.; Besser, H.; Khelifi, F.; Harabi, S.; Hamad, A.; Boyce, A.; Laouar, R.; Hamed, Y. Hydrochemistry and stable isotopes (δ18O and δ2H) tools applied to the study of karst aquifers in Southern Mediterranean Basin (Teboursouk Area, NW Tunisia). *J. Afr. Earth Sci.* **2018**, *137*, 208–217. [CrossRef]

25. Hamed, Y.; Hadji, R.; Ahmadi, R.; Ayadi, Y.; Shuhab, K.; Pulido-Bosch, A. Hydrogeological investigation of karst aquifers using an integrated geomorphological, geochemical, GIS, and remote sensing techniques (Southern Mediterranean Basin—Tunisia). *Environ. Dev. Sustain.* **2023**, *26*, 6943–6975. [CrossRef]

26. Ascott, M.J.; Macdonald, D.M.J.; Black, E.; Verhoef, A.; Nakohoun, P.; Tirogo, J.; Sandwidi, W.J.P.; Bliefernicht, J.; Sorensen, J.P.R.; Bossa, A.Y. In situ observations and lumped parameter model reconstructions reveal intra-annual to multidecadal variability in groundwater levels in Sub-Saharan Africa. *Water Resour. Res.* **2020**, *56*, e2020WR028056. [CrossRef]

27. MacDonald, A.M.; Lark, R.M.; Taylor, R.G.; Abiye, T.; Fallas, H.C.; Favreau, G.; Goni, I.B.; Kebede, S.; Scanlon, B.; Sorensen, J.P.R.; et al. Mapping groundwater recharge in Africa from ground observations and implications for water security. *Environ. Res. Lett.* **2021**, *16*, 034012. [CrossRef]

28. Tarpanelli, A.; Paris, A.; Sichangi, A.W.; O'loughlin, F.; Papa, F. Water resources in Africa: The role of Earth observation data and hydrodynamic modeling to derive river discharge. *Surv. Geophys.* **2023**, *44*, 97–122. [CrossRef]

29. The KARMA Project. The KARMA project. Retrieved from Karst Aquifer Resources Availability and Quality in the Mediterranean Area, 2023. Available online: http://karma-project.org/ (accessed on 10 March 2023).

30. Jemmali, N.; Rddad, L.; Souissi, F.; Carranza, E. The ore genesis of the Jebel Mecella and Sidi Taya F Ba (Zn Pb) Mississippi Valley-type deposits, Fluorite Zaghouan Province, NE Tunisia, in relation to Alpine orogeny: Constraints from geological, sulfur, and lead isotope studies. *C. R. Geosci.* **2019**, *351*, 312–320. [CrossRef]

31. Djebbi, M.; Besbes, M.; Sagna, J.; Rekaya, M. Les sources karstiques de Zaghouan. Recherche d'un opérateur pluie-débit. *Sci. Tech. Environ.* **2001**, *13*, 125–128.

32. Andreo, B.; Vías, J.; Durán, J.J.; Jiménez, P.; López-Geta, J.A.; Carrasco, F. Methodology for groundwater recharge assessment in carbonate aquifers: Application to pilot sites in Southern Spain. *Hydrogeol. J.* **2008**, *16*, 911–925. [CrossRef]

33. NASA Shuttle Radar Topography Mission (SRTM). Shuttle Radar Topography Mission (SRTM) Global. Distributed by OpenTo-pography, 2013. Available online: https://doi.org/10.5069/G9445JDF (accessed on 9 January 2023).

34. Mtimet, A. Soils of Tunisia. In *Soil Resources of Southern and Eastern Mediterranean Countries*; Zdruli, P., Steduto, P., Lacirignola, C., Montanarella, L., Eds.; CIHEAM: Bari, Italy, 2001; pp. 243–262, (Options Méditerranéennes: Série B. Etudes et Recherches; n. 34).
35. Mhimdi, A. Apport de L'imagerie Satellitaire Pour la Détection des Systèmes Karstiques (Les Dolines) Dans la Région de Zaghouan. Master Thesis, Faculté des Sciences de Tunis, University Tunis El Manar, 2020. Available online: http://www.biruni.tn/v-en/library-catalog.php?ei=28 (accessed on 8 January 2023).

36. Baillie, I.C. Soil Survey Staff 1999, Soil Taxonomy. Soil Use Manag. **2001**, *17*, 57–60. [CrossRef]

37. Sagna, J. Study and Modeling of Zaghouan Karst Sources. Master Thesis, National Engineering School of Tunis, University Tunis El Manar, Tunis, Tunisia, 2000.

38. Vicente-Serrano, S.; Beguería, S.; López-Moreno, J. A multiscalar drought index sensitive to global warming: The Standardized Precipitation Evapotranspiration Index (SPEI). *J. Clim.* **2010**, *23*, 1696–1718. [CrossRef]

39. McKee, T.B.; Doesken, N.J.; Kleist, J. The relationship of drought frequency and duration to time scales. In Proceedings of the 8th Conference on Applied Climatology, Anaheim, CA, USA, 17–22 January 1993; pp. 179–184.

40. Bloomfield, J.P.; Marchant, B.P. Analysis of groundwater drought building on the standardised precipitation index approach. *Hydrol. Earth Syst. Sci.* **2013**, *17*, 4769–4787. [CrossRef]

41. Peterson, R.A. Finding optimal normalizing transformations via bestNormalize. *R J.* **2021**, *13*, 310–329. [CrossRef]

42. Peterson, R.A.; Cavanaugh, J.E. Ordered quantile normalization: A semiparametric transformation built for the cross-validation era. *J. Appl. Stat.* **2020**, *47*, 13–15, 2312–2327. [CrossRef] [PubMed]

43. D'Oria, M.; Balacco, G.; Todaro, V.; Alfio, M.R.; Tanda, M.G. Assessing the impact of climate change on a coastal karst aquifer in a semi-arid area. *Groundw. Sustain. Dev.* **2024**, *25*, 101131. [CrossRef]

44. Vías, J.M.; Andreo, B.; Perles, M.J.; Carrasco, F.; Vadillo, I.; Jiménez, P. Proposed method for groundwater vulnerability mapping in carbonate (karstic) aquifers: The COP method. *Hydrogeol. J.* **2006**, *14*, 912–925. [CrossRef]

45. Daly, D.; Dassargues, A.; Drew, D.; Dunne, S.; Goldscheider, N.; Neale, S.; Zwahlen, F. Main concepts of the "European approach" to karst-groundwater-vulnerability assessment and mapping. *Hydrogeol. J.* **2002**, *10*, 340–345. [CrossRef]

46. Othman, J. GIS Applications in Karst Hydrogeology: APLIS for Recharge Estimation and COP for Vulnerability Assessment: Application on the Qachqouch Spring Catchment in Lebanon. Master Thesis, Lebanese University-Faculty of Sciences, Beirut, Lebanon, 2021.

47. Xanke, J.; Goldscheider, N.; Bakalowicz, M.; Barberá, J.A.; Broda, S.; Chen, Z.; Ghanmi, M.; Günther, A.; Hartmann, A.; Jourde, H.; et al. Mediterranean Karst Aquifer Map (MEDKAM), 1:5,000,000. Berlin, Karlsruhe, Paris, 2022. Available online: https://doi.org/10.25928/MEDKAM.1 (accessed on 15 January 2023).

48. Cinkus, G.; Mazzilli, N.; Jourde, H. Identification of relevant indicators for the assessment of karst systems hydrological functioning: Proposal of a new classification. *J. Hydrol.* **2021**, *603*, 127006. [CrossRef]

49. Zhou, C.; Nooijen, R.; Kolechkina, A.; Gargouri-Ellouze, E.; Slama, F.; van de Giesen, N. The uncertainty associated with the use of copulas in multivariate analysis. *Hydrol. Sci. J.* **2023**, *68*, 2169–2188. [CrossRef]

50. Mazzilli, N.; Guinot, V.; Jourde, H.; Lecoq, N.; Labat, D.; Arfib, B.; Baudement, C.; Danquigny, C.; Dal Soglio, L.; Bertin, D. KarstMod: A modelling platform for rainfall-discharge analysis and modelling dedicated to karst systems. *Environ. Model. Softw.* **2019**, *122*, 103927. [CrossRef]

51. Slama, F.; Gargouri-Ellouze, E.; Faydi, T.; Cinkus, G.; Jourde, H.; Bouhlila, R. Rainfall-Discharge Modelling of the Djebel Zaghouan Aquifer Using KarstMod. Eurokarst2022, Malaga, Spain, 2022. Available online: https://www.researchgate.net/publication/376140976_Rainfall-discharge_modelling_of_the_Djebel_Zaghouan_aquifer_using_KarstMod (accessed on 5 March 2024).

52. Gargouri-Ellouze, E.; Slama, F.; Ouedraogo, T.A.; Bouhlila, R. Statistical Evaluation and ANN Modelling of the Discharge Spring Response for Djebel Zaghouan Karstic Aquifer. Eurokarst2022, Malaga, Spain, 2022. Available online: https://www.eurokarst.org/wp-content/uploads/2022/06/Eurokarst-program-2022_Updated-Thu24.pdf (accessed on 5 March 2024).

53. Wang, Z.; Wu, R.; Huang, K.; Qiu, Y.; Li, Z.; Lv, Y.; Wan, J. Structure identification of a karst groundwater system based on high-resolution rainfall-hydrological response characteristics. *Environ. Sci. Pollut. Res.* **2022**, *29*, 26922–26935. [CrossRef] [PubMed]

54. Ravbar, N.; Engelhardt, I.; Goldscheider, N. Anomalous behavior of specific electrical conductivity at a karst spring induced by variable catchment boundaries: The case of the Podstenjšek spring, Slovenia. *Hydrol. Process.* **2011**, *25*, 2130–2140. [CrossRef]

55. Chang, Y.; Hartmann, A.; Liu, L.; Jiang, G.; Wu, J. Identifying more realistic model structures by electrical conductivity observations of the karst spring. *Water Resour. Res.* **2021**, *57*, e2020WR028587. [CrossRef]

56. Karimi, H.; Raeisi, E.; Bakalowicz, M. Characterising the main karst aquifers of the Alvand basin, northwest of Zagros, Iran, by a hydrogeochemical approach. *Hydrogeol. J.* **2005**, *13*, 787–799. [CrossRef]

57. Wu, P.; Tang, C.; Zhu, L.; Liu, C.; Cha, X.; Tao, X. Hydrogeochemical characteristics of surface water and groundwater in the karst basin, southwest China. *Hydrol. Process.* **2009**, *23*, 1012–1022. [CrossRef]

58. Ben Ammar, S.; Taupin, J.D.; Ben Alaya, M.; Zouari, K.; Nicolas, P.; Khouatmia, M. Using geochemical and isotopic tracers to characterize groundwater dynamics and salinity sources in the Wadi Guenniche coastal plain in northern Tunisia. *J. Arid Environ.* **2020**, *178*, 104150. [CrossRef]
59. IAEA. GNIP Maps and Animations, International Atomic Energy Agency, Vienna, 2001. Available online: https://nucleus.iaea.org/wiser (accessed on 28 January 2025).

60. Leduc, C.; Taupin, J.D.; Le Gal la Salle, C. Estimation de la recharge de la nappe phréatique du Continental Terminal (Niamey, Niger) à partir des teneurs en tritium. C. R. Acad. Sci. Paris **1996**, *323*, 599-605.

61. Salem, O.; Visser, J.H.; Dray, M.; Gonfiantini, R. *Groundwater Flow Patterns in the Western Libyan Arab Jamahiriya Evaluated from Isotopic Data*; International Atomic Energy Agency (IAEA): Vienna, Austria, 1980.

62. Wang, Z.-J.; Yue, F.-J.; Lu, J.; Wang, Y.-C.; Qin, C.-Q.; Ding, H.; Xue, L.-L.; Li, S.-L. New insight into the response and transport of nitrate in karst groundwater to rainfall events. *Sci. Total Environ.* **2022**, *818*, 151727. [CrossRef] [PubMed]

63. Wang, Z.-J.; Yue, F.-J.; Wang, Y.-C.; Qin, C.-Q.; Xue, L.-L.; Li, S.-L. The effect of heavy rainfall events on nitrogen patterns in agricultural surface and underground streams and the implications for karst water quality protection. *Agric. Water Manag.* **2022**, *266*, 107600. [CrossRef]

64. Troudi, N.; Tzoraki, O.; Hamzaoui-Azaza, F.; Melki, F.; Zammouri, M. Evaluation of Groundwater Quality Using Nitrate Pollution Index and the Potential Health Risk Method in Guenniche Basin of Northern Bizerte (Tunisia, North Africa). In *Recent Advances in Environmental Science from the Euro-Mediterranean and Surrounding Regions*, 4th ed.; Ksibi, M., Sousa, A., Hentati, O., Chenchouni, H., Velho, J.L., Negm, A., Rodrigo-Comino, J., Hadji, R., Chakraborty, S., Ghorbal, A., Eds.; EMCEI 2022, Advances in Science, Technology & Innovation; Springer: Cham, Switzerland, 2024. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.