A combined methodology for estimating the potential natural aquifer recharge in an arid environment

ABSTRACT

An innovative methodology that combines an indirect physiography-based method for determining the runoff coefficient at a sub-basin scale and a water balance model applied on a daily time scale was developed to calculate the natural groundwater recharge in three watersheds within the Oum Zessar arid area, Tunisia. The effective infiltration was calculated as part of the water surplus by considering the average available water content (AWC) of soil and an average runoff coefficient for each sub-basin. The model indicates that the sub-basins covered mainly by the “artificial” soils of tabias and jessour, characterized by average AWC values greater than 150 mm, did not contribute to natural groundwater recharge over the 10-year period (2003–2012) considered. The estimated volume for the Triassic aquifer amounted to about 4.5 hm³ year⁻¹, which is consistent with previous studies. For the Jurassic and Cretaceous aquifers, the estimated volumes amounted to about 200 dm³ year⁻¹.

KEYWORDS:

• natural aquifer recharge
• hydrogeology
• arid environment
• Oum Zessar
• Tunisia

Editor S. Archfield; Associate Editor F. Tauro

1 Introduction

The areas most vulnerable to environmental changes, particularly hydrological changes, are arid regions (approximately 20% of the world’s land area) that experience severe land and water degradation. These vulnerable ecosystems are subject to great pressure due to rapid population growth, which leads to the uncontrolled use of surface water and groundwater for irrigation (Shen and Chen 2010). The Mediterranean areas of Southern Europe, North-Africa, and the Middle East are affected by dramatic changes in the sustainability, quantity, quality, and management of water resources. An increase in temperature and a decrease in precipitation at the end of the 21st century are forecast by most climate models (García-Ruiz et al. 2011).

Groundwater is essential for meeting fresh water demand worldwide to satisfy urban, agricultural, industrial, and environmental needs, particularly in arid and semiarid regions. Currently, groundwater is the main source of drinking water for at least half the world’s population, which includes a significant part of urban population as well as rural communities (Puri et al. 2008, Gleeson et al. 2012). In developing countries, dependence on groundwater is high owing either to water scarcity or a lack of safe drinking water from surface water supplies (Ghiglieri and Carletti 2010, Owen et al. 2010). In particular, during extended drought events, groundwater also acts as a key strategic resource. Furthermore, groundwater is often inadequately managed in terms of ensuring its long-term sustainability, especially in arid and semi-arid regions. Most of the major aquifers in these regions, particularly dry areas that rely most heavily on groundwater, are experiencing rapid depletion rates due to overexploitation (Famiglietti 2014).

Interactions between climate, geology, morphology, soil condition, and vegetation affect the groundwater recharge process. As potential evapotranspiration generally exceeds rainfall in (semi-)arid areas, the groundwater recharge process depends particularly on the following aspects: highly concentrated rainfall events; accumulation of runoff water in depressions and streams; and rapid percolation through cracks, fissures, or solution channels. The presence of thick alluvial soils normally hinders groundwater recharge, as they allow for high water retention during the wet season. Vegetation subsequently extracts soil water in the next dry season. Conversely, favourable recharge conditions are created by a poor vegetation cover on permeable soil or a fractured porous bedrock near the surface combined with high-intensity rainfall events (De Vries and Simmers 2002). Water budget calculations are important because they provide an accurate estimate of the water resources in an area, which are related to the natural recharge of aquifers.

Groundwater recharge can be estimated by several methods which have been intensively discussed in the literature (Lerner et al. 1990, De Vries and Simmers 2002, Scanlon et al. 2002). Recently, the importance of the groundwater recharge calculation was stated by Smerdon and Drewes (2017), who provide an overview of recent studies exploring the recharge process through different estimation methods and modelling approaches. The different approaches depend on the availability of data and the required level of accuracy. Of these, numerical models are an effective tool that have been developed to compute groundwater recharge using daily balances (Bellot et al. 2001, Samper et al. 2007, Touhami et al. 2013). Spatially distributed modelling is an important instrument for studying the hydrological cycle, and the results of such simulations are particularly relevant to water resources assessment. The following are examples of effective distributed models: Water balance Simulation Model-ETH (WaSiM-ETH) (Klok et al. 2001, Wagner et al. 2009); TOPographic Kinematic APproximation and Integration (TOPKAPI) (Liu and Todini 2002); Gridded Surface/Subsurface Hydrologic Analysis (GSSHA) (Downer and Ogden 2004); and the WetSpass model (Tilahun and Merkel 2009), the latter of which integrates GIS and the water balance method. Spatially semi-distributed models, such as PREecipitation-Runoff-EVApotranspiration Hydrological Response Unit (PREVAH) (Viviroli et al. 2009) and the Soil and Water Assessment Tool (SWAT) (Arnold et al. 1998, Neitsch et al. 2005), use a conceptual process-oriented approach by implementing HRUs as a cost-efficient way of achieving spatially distributed results. In SWAT, HRUs are the smallest spatial units of calculation and are homogeneous in terms of soils, land use, and slope (Gassman et al. 2007, Xu et al. 2009), which allows the simulation of a high level of spatial detail.

The objective of this study was to develop a combined methodology to estimate natural aquifer recharge in three watersheds of the Oum Zessar area, an arid region in south-eastern Tunisia. This study proposes an innovative method that combines an indirect physiography-based method for determining the runoff coefficient first developed by Kennessey (1930) and subsequently modified and tested by several authors in different semi-arid regions in the Mediterranean area (Tardi and Vittorini 1977, Barazzuoli et al. 1986, Spadoni et al. 2010, Ghiglieri et al. 2014, Grillone et al. 2014), and a water balance model applied on a daily time scale as a modification of the model proposed by Allen et al. (1998). This model allows the calculation of effective infiltration I_e through the water balance of the root zone, when a water surplus occurs, by considering crop evapotranspiration under stress conditions. Ultimately, estimates of the regulatory reserves/stock volumes for each aquifer were obtained by overlaying the results of effective infiltration (I_e) over the spatial distribution of the aquifers.

2 Study area

The Oum Zessar study area is located in the province (governorate) of Médenine. It extends from the limit of the Great Erg Oriental in the west to the Mediterranean Sea (Gulf of Gabès) in the east, crossing the Dahar plateau, Beni Khedache Mountains and Jeffara plain (Fig. 1(a)). From ecological, hydrological and socio-economic points of view, this area can be considered representative of arid southeastern Tunisia.

Figure 1. (a) Oum Zessar study area and (b) hydrogeological map.

The Dahar is a calcareous plateau with an elevation ranging from 400 to 600 m a.s.l. stretching in a NW–SE direction. It is dissected by small wadis (ephemeral rivers) and ends in the east in a cuesta characterized by a relatively steep escarpment descending abruptly from its crest (Ouessar 2007). The Jeffara is a large plain consisting of crusted Quaternary deposits stretching between the Matmata mountains, the Dahar, the littoral and the Tunisian border with Libya. At its highest point, it has an elevation of 100 m a.s.l.; the plain is crossed by several wadis that drain the mountain zone and, after traversing a long distance on a rather flat slope, flow into the sea or sebkhats (Ouessar 2007). The study area, which is typically agro-pastoral, has experienced a rapid expansion of crop cultivation in flatter areas and marginal rangelands. This phenomenon has had negative effects on native rangelands, reducing the grazing areas at disposal to animals. Desertification is on-going and aggravated by anthropogenic pressures due to changes in socio-economic policies and population growth (Sghaier et al. 2011).

The study area is characterized by an arid Mediterranean climate that is influenced by both the Dahar and the Matmata Mountains (continental arid) and the Mediterranean Sea (Gulf of Gabès) (maritime arid). The average annual rainfall is between 150 and 240 mm per year, and around 20 rain days are recorded each year during the “wet” season, which stretches from November to January (Ouessar 2007, Ouessar et al. 2008). Rainfall is characterized by a high variability, both in time and space. Owing to the high temperatures and low rainfall levels, potential evapotranspiration (ET₀) is very high, reaching a value of 1450 mm in Médenine (Hargreaves method) (Ouessar 2007). Thus, the climatic water balance is almost negative throughout the year (Ouessar 2007, Adham et al. 2016).

The study area mainly consists of Permian to Holocene sedimentary successions that host several aquifers (Ouessar and Yahyaoui 2006, Ouessar 2007) (Fig. 1(b)).

2.1 Jeffara of Médenine aquifer

The shallow aquifer of Jeffara of Médenine covers the entire coastal plain and the valleys in the Dahar domain. It is hosted in Quaternary deposits consisting of continental alluvium. The thickness of this aquifer increases progressively towards the sea, and the total depth of wells exploiting the shallow aquifer varies from 3 to 20 m. The transmissivity varies from 5–9 × 10⁻⁴ to 8–9 × 10⁻³ m² s⁻¹ (Ouessar and Yahyaoui 2006). The recharge of this aquifer is mainly represented by direct infiltration of rainfall, and the infiltration rate from the wadis during floods varies between 30% and 50% of the river flow (Ould Baba Sy 2005, Yahyaoui 2007).

2.2 Mio-Pliocene aquifer

The Mio-Pliocene aquifer is hosted in the Miocene continental deposits, which consist of red clays with intercalations of conglomeratic and sandy levels (Benton et al. 2000). The thickness of these sand deposits increases from 25 to 80 m towards the sea. They are characterized by fine granulometry, and transmissivity ranges from 1 × 10⁻³ to 2 × 10⁻² m² s⁻¹ (Kamel et al. 2013).

2.3 Zeuss–Koutine aquifer

The Zeuss–Koutine aquifer (ZK) in the Jeffara plain includes two main entities separated by the Médenine Fault: the ZK Jurassic aquifer and ZK Senonian aquifer. The ZK Jurassic aquifer consists of a series of dolomite black marls of basal Jurassic. In the ZK Lower Senonian aquifer, the unconfined layer of the lower Senonian and Turonian can be replenished by runoff but also by lateral flow from the ZK Jurassic through the fault (Ouessar 2007). The Zeuss–Koutine aquifer is mainly exploited for drinking water and irrigation. The withdrawals have increased from 206 L s⁻¹ in 1977 to 641 L s⁻¹ in 2014, and this exploitation has produced a continuous piezometric decline of 0.5 to 1 m year⁻¹ on average (Jarray et al. 2017).

2.4 Béni Khédache aquifer

The Béni Khédache aquifer (BK), in the Dahar, includes two carbonate sequences: the condensed series of the upper Triassic to the lower Bathonian and limestones of the upper Jurassic. They are separated by clay and sandstone of the Techout formation. This aquifer can be reached at 200–300 mbgl. It is directly replenished by the infiltration of runoff water from the Dahar plateau or percolation from the eastern cliff (Ouessar 2007).

2.5 Triassic sandstone aquifer

The Triassic sandstone aquifer (Grès de Trias) stretches over a large area between the upper Permian outcrops (Jebel Tébaga) to the north, outcroppings of Dahar to the west, the Jurassic aquifer to the east and, along the Jeffara plain to the south-east up to Tripoli, Libya (Besbes et al. 2005, Ouessar 2007). In the continental Jeffara plain, lower to middle Triassic formations outcrop along the wadis because of the erosional processes affecting the Cenozoic deposits (Arras et al. 2015). Therefore, the piedmont area and the wadis are considered the preferential recharge areas. The transmissivity varies from 0.33 × 10⁻³ to 8.4 × 10⁻³ m² s⁻¹ (Yahyaoui 2007). This aquifer is used mainly for drinking water and irrigation. Since the end of the 1990s, the rate of exploitation has exceeded the rate of replenishment. In 2010, the rate of overuse reached 160%, which resulted in a declining trend in the piezometric level (WADIS-MAR 2016).

In southeast Tunisia, because of water scarcity, various water harvesting techniques have been introduced over many centuries to make the land productive despite its geographical location. Jessour (singular jessr) are hydraulic units consisting of an impluvium, terrace and dyke, generally occupy inter-mountain areas and hill watercourses to intercept runoff and sediment; they are typically arranged in the form of a gradoni, a series of contour ditches used mainly for growing trees. The impluvium or catchment is the area in which runoff water is collected and conveyed, and it is bounded by a natural water divide line. The terrace or cropping zone is the area in which farming is practiced. It is formed progressively by the deposition of sediment, which creates artificial soil (up to 5 m deep) close to the dyke. The dyke is a barrier made of earth and equipped with a central and/or lateral spillway. It is used to block sediment and runoff water (Fleskens et al. 2005, Schiettecatte et al. 2005, Ouessar 2007). The tabia technique is a replica of the jessour system constructed in foothill and piedmont areas, where the slope does not exceed 3% and the soils are relatively deep.

3 Materials and methods

The proposed methodology is represented in the conceptual workflow in Fig. 2, and is described as follows.

Figure 2. Conceptual workflow for the new methodology proposed to estimate the natural aquifer recharge.

This new approach allows researchers to overcome the challenges presented by the critical condition that typically occurs in arid regions, where the magnitude of the recharge rate is small in comparison with other variables, especially evapotranspiration (De Vries and Simmers 2002, Scanlon et al. 2002). In these environments, the water budget should be calculated on a daily time step, as precipitation sometimes greatly exceeds evapotranspiration on a single day, whereas the effects of extreme precipitation events, which are the most responsible for the recharge of aquifers, can be attenuated by averaging over longer time periods. This approach is not a substitute for traditional methodologies relying on expensive field datasets, but it is useful in providing a reasonable, but conservative, natural aquifer recharge appraisal, particularly in ungauged watersheds.

The methodology presented here was applied at the watershed scale. The results suggest that it is necessary to enhance the accuracy of the natural aquifer recharge estimates considering the hydrogeological, hydrogeochemical and, isotopic evidence recently obtained as part of the currently unpublished results of the WADIS-MAR project (www.wadismar.eu).

3.1 Watershed delineation

To estimate the water budget, three watersheds were delineated in the study area: Koutine (271.6 km²), Megarine-Arniane (247.4 km²) and Hajar (80.3 km²) (Fig. 3). The Shuttle Radar Topography Mission 1 Arc-Second Global Digital Elevation Model (SRTMGL1 DEM; http://earthexplorer.usgs.gov) released in 2014 was used to define the watersheds, sub-watersheds, stream networks and outlets for each sub-basin. It covers the African continent in 1° × 1° tiles at a 1 arc-second (about 30 m) resolution.

Figure 3. Koutine, Megarine-Arniane and Hajar watersheds.

3.2 Climate data

Monthly rainfall data (1970–2012) recorded by six rain gauge stations (Koutine, Allamet, Ksar Hallouf, Ksar Jedid, Sidi Makhlouf and Béni Khedache) operated by the hydrological service of the Water Resources Directorate of the Ministry of Agriculture (DGRE) were considered (Fig. 1(a)). The only available rainfall and temperature data measured on a daily time step were those from the weather station of Médenine (1977–2012); the data can also be downloaded from the website www.TuTiempo.net. Climatic data from six more weather stations of the National Centers for Environmental Prediction’s Climate Forecast System Reanalysis (CFSR) (1979–2014) (http://globalweather.tamu.edu), located in and around the watersheds, were also considered (Fig. 1(a)). The CFSR climate data are produced using cutting-edge data assimilation techniques (both conventional meteorological gauge observations and satellite irradiances) and highly advanced (coupled) atmospheric, oceanic, and surface-modelling components at a resolution of approx. 38 km. Consequently, the output of CFSR data processing involves various spatial and temporal interpolations (Dile and Srinivasan 2014). The assessment of the reliability of CFSR data compared to measured data was carried out through a statistical analysis to evaluate correlations, trends, and accuracy. For this purpose, the average daily rainfall and temperature data were considered. The average daily data of the CFSR stations were compared to those of the Médenine weather station (the only one with available daily measured data).

3.3 Runoff coefficient, C_k

The mean annual runoff coefficient (C_k) was determined through the procedure described by Ghiglieri et al. (2014) and is given in Fig. 2. For each of the partial runoff components (C_a, C_v and C_p), three series of coefficients are defined according to an interval of values in the annual average aridity index I_a to take the climatic conditions into account. The aridity index I_a is determined according to the De Martonne equation modified by Tardi and Vittorini (1977) (Ghiglieri et al. 2014). Once the I_a is calculated, a series of coefficients are selected to calculate the partial runoff coefficients (C_a, C_v and C_p) (Table 1); then, the C_k is evaluated by their sum.

Table 1. Table of partial runoff coefficients classified by aridity index I_a (modified from Ghiglieri et al. 2014).

Coefficients	Aridity index (Ia)
	Ia < 25	25 < Ia < 40	Ia > 40
Slope, Ca
>35%	0.22	0.26	0.30
10% < S < 35%	0.12	0.16	0.20
3.5% < S < 10%	0.01	0.03	0.05
<3.5%	0.00	0.01	0.03
Land cover, Cv
No vegetation	0.26	0.28	0.30
Grazing land	0.17	0.21	0.25
Cultivation/shrubby	0.07	0.11	0.15
Woods/forests	0.03	0.04	0.05
Permeability, Cp
Very low*	0.21	0.26	0.30
Low*	0.16	0.21	0.25
Medium	0.12	0.16	0.20
Good	0.06	0.08	0.10
High	0.03	0.04	0.05

* class not present in the study area.

Re-produced with permission from IAHS

The slope layer obtained from the SRTMGL1 DEM, where the slope was expressed as a percentage, was classified into four classes (Spadoni et al. 2010, Ghiglieri et al. 2014). The land cover layer was obtained from the land use/land cover map developed by Afrasinei et al. (2017) via a visual interpretation of Landsat 8 satellite images and validated by detailed ground truth data. The original land cover classes were grouped according to the four classes described in Table 2.

Table 2. Classes of land cover, C_v (from Ghiglieri et al. 2014).

Class	Description
No vegetation	Areas lacking vegetation or with sparse vegetation cover, thus having minimum participation to the transpiration process
Grazing land	Areas covered mainly by herbaceous and/or shrub vegetation
Cultivation/shrubby	Units characterized by forest features (even sparse) and by tree and herbaceous cultivations
Woods/forests	Units consisting of real and proper forests, with 75% trees and the rest possibly shrubs and small trees

Re-produced with permission from IAHS

The permeability layer was derived from the official geological map (1:100 000) (Office National des Mines 1996). Each hydrogeological complex was defined by the aggregation of geological formations with permeability conditions similar to those listed in the literature (Swezey 1999, Ould Baba Sy 2005, Ouessar and Yahyaoui 2006, Yahyaoui 2007, Kamel et al. 2013, Chihi et al. 2015) and then classified into three classes according to the runoff coefficient C_k method (Table 3). The spatial distribution of each class for each layer is shown in Fig. 4.

Table 3. Classes of permeability, C_p (modified from Ghiglieri et al. 2014).

Class	K (m s^−1)
Very low*	<10^−9
Low*	>10^−9, <10^−6
Medium	>10^−6, <10^−4
Good	>10^−4, <10^−2
High	>10^−2

* class not present in the study area.

Re-produced with permission from IAHS

Figure 4. Spatial distribution of the classes of each thematic layer (slope, permeability and land cover) for the three watersheds.

The runoff coefficient was calculated using a procedure within the GIS environment (ESRI ArcGIS 10) in which each vector information layer related to the three watersheds (Koutine, Megarine-Arniane and Hajar) is converted to raster format using the same cell size of the SRTMGL1 DEM and then processed by overlay mapping (Fig. 4).

3.4 Multiple regression for aridity index calculation

In this study, we used monthly rainfall data from the six rain gauge stations (Koutine, Allamet, Ksar Hallouf, Ksar Jedid, Sidi Makhlouf and Beni Khedache) and the Médenine weather station, and monthly minimum, maximum, and mean temperature data from the six CFSR weather stations with Médenine data over a 30-year period (1981–2010). To obtain monthly temperature and precipitation maps, a procedure based on the statistical relationship between climate data and the geographic parameters of the study area was adopted and implemented in GIS. Spatial interpolations of monthly climatic variables (average precipitation and average mean temperature) were performed using single and multiple regressions, considering the elevation and distance from the sea obtained from the SRTMGL1 DEM as independent geographic variables (Boi et al. 2011, Canu et al. 2014). The calculation was carried out over an area of 1 600 km² including all the rain gauge and weather stations. The regression equations obtained for precipitation and mean temperature are given in Table 4.

Table 4. Regression equations and statistics. Alt: altitude (m a.s.l.); Dist: distance from the sea (m); Std err.: standard error; Signif. F: significance F.

Regression equations	R^2	Std err.	Signif. F
Precipitation
Pjan = 0.0652Alt + −0.0006Dist + 30.8932	0.92	1.84	0.006
Pfeb = 0.0083Alt + 0.0002Dist + 6.2189	0.90	1.86	0.010
Pmar = 0.0211Alt + 13.6051	0.70	2.74	0.020
Papr = 0.0128Alt + 8.9076	0.85	1.07	0.003
Pmay = 0.0099Alt + 3.9564	0.55	1.78	0.058
Pnov = 0.0029Alt + 0.0001Dist + 16.1887	0.56	2.25	0.197
Mean temperature
Tjan = −0.0069Alt + 12.8988	0.90	0.32	0.001
Tfeb = −0.0049Alt + 14.1049	0.87	0.26	0.002
Tmar = −0.0023Alt + 16.5035	0.88	0.12	0.002
Tmay = −0.0039Alt + 0.0001Dist + 22.4008	0.56	0.36	0.195
Tjun = −0.0027Alt + 0.0001Dist + 25.5459	0.69	0.42	0.095
Tjul = −0.0029Alt + 0.0001Dist + 27.9788	0.51	0.60	0.239
Toct = −0.0035Alt + 24.0322	0.92	0.14	0.001
Tnov = −0.0063Alt + 18.6576	0.90	0.29	0.001
Tdec = −0.0072Alt + 14.3157	0.91	0.31	0.001

These equations were applied to SRTM DEM and the distance from the sea by means of map algebra procedures in ArcGIS to obtain the average monthly precipitation and temperature maps. Other spatial interpolation methods, such as kriging (spherical method) and spline and inverse distance weighting (IDW), were compared using the Geostatistical Analyst extension in ArcGIS; this was done on average monthly data for which the single or multiple regression method was not appropriate (R² < 0.5). Because of the few available stations, their spatial distribution in the study area, and the low autocorrelations of data, the IDW method turned out to be the most suitable for these spatial interpolations. The average monthly and annual precipitation and temperature maps were obtained using map algebra calculations. In addition, monthly precipitation and the monthly mean temperature of the hottest month (August) were used to calculate the aridity index I_a (Tardi and Vittorini 1977) (Table 5).

Table 5. Average annual rainfall and mean temperature and average rainfall and mean temperature in the hottest month (August) for the reference period (1981–2010). P: annual rainfall (mm); p: rainfall in the hottest month (mm); T: annual mean temperature (°C); t: mean temperature in the hottest month (°C).

	P	p (August)	T	t (August)	Ia
Min	148.68	0.08	19.06	29.31	2.51
Max	221.68	6.47	21.17	30.50	4.22
Mean	169.44	2.29	20.57	29.56	3.24

3.5 C_k calculation and annual water budget

The annual average aridity index I_a was calculated for the 30-year period (1981–2010) (Equation 1) by a map algebra procedure considering the same cell size of the SRTMGL1 DEM. As the I_a was less than 25 (Table 5) for each watershed, the analysis considered the series of coefficients related to the class I_a < 25 for each thematic layer defined by the method. The overlay mapping of the three layers by summing the C_a, C_p and C_v values of each cell allowed us to determine the spatial distribution of C_k. The 36 possible combinations of classes and the C_k values were divided into five classes (Ghiglieri et al. 2014). The spatial distribution of the runoff coefficient within each watershed is shown in Fig. 5.

Figure 5. Distribution of runoff coefficient classes for the three watersheds.

Figure 6. Comparison between average monthly rainfall data of the Médenine weather station and average data of the other stations (six rain gauges and Médenine station): (a) histograms (dark blue: all stations data; light blue: Médenine data), and (b) correlation.

Figure 7. (a) Water balance of the root zone and (b) water stress coefficient K_s (modified from Allen et al. 1998).

To calculate the average annual effective infiltration I_e (natural aquifer recharge) over the 30-year period (1981–2010), the standard water balance equation was used:

(1) $I_e=P-\mathrm{E}{\mathrm{T}}_a-R$(1)

where P is the average annual precipitation, ET_a is the annual actual evapotranspiration, and R is the average annual runoff. The annual actual evapotranspiration was calculated using the empirical equation proposed by Turc (1954) and modified by Santoro (1970), which is considered appropriate for arid and semiarid regions in the Mediterranean area (Barroccu and Soddu 2006, Ghiglieri et al. 2006, Giannecchini and Doveri 2016):

(2) $\mathrm{E}{\mathrm{T}}_a=\frac{P}{\sqrt{0.9+\frac{P^2}{\left(586-10T+0.05T^3\right)}}}$(2)

where P is the average annual precipitation (mm) and T is the average annual temperature (°C). By applying Equation (3) (Spadoni et al. 2010, Ghiglieri et al. 2014), it was possible to calculate the average annual runoff R using the runoff coefficient C_k:

(3) $R=\left(P-\mathrm{E}{\mathrm{T}}_a\right)\times C_k$(3)

Compared to the above mentioned studies, where the annual precipitation was greater than the annual actual evapotranspiration, in this study, the annual effective rainfall (P – ET), calculated on an annual basis, was < 0 for most areas, except for a small region in the Dahar in the Koutine and Megarine-Arniane watersheds, thereby preventing an estimate of runoff and effective infiltration from being determined. This is due to the erratic behaviour of rainfall in arid environments, such as southern Tunisia, where high intensity rainfall events are concentrated in a few days. Only after short but intense events does rainfall exceed evapotranspiration, leading to runoff and recharging the aquifers. Such a condition led to the use of daily precipitation and temperature data to estimate the natural recharge of aquifers.

3.6 Water budget on daily basis

Daily precipitation and mean temperature data from the Médenine weather station were considered for the entire area covered by the three watersheds in the daily calculations of the water budget. In fact, because climate data from the Médenine station are the only data available that are measured daily, this station can be considered as representative of the climate of the study area, especially for the Jeffara plain, as shown by the graphs in Fig. 6. The water budget was estimated for each sub-basin of each watershed over a 10-year period (2003–2012) that presented the most complete climate data series.

In this work, a modified water balance model based on the one proposed by Allen et al. (1998) (FAO Paper 56) was applied on a daily time scale. The model considers the effective infiltration I_e (or deep percolation DP) as part of the water surplus in the soil. As in Allen et al. (1998), crop evapotranspiration under stress conditions is considered and depleted water is computed as the result of the water balance of the root zone, as explained in Fig. 7(a). This estimation considered a reference crop with standard characteristics.

Saturation of the soil active layer is reached when the total available water (TAW) content is a maximum. Any additional water gives rise to runoff and deep percolation. If the soil water content decreases below a threshold value, crops begin to experience stress as soil water can no longer be transported quickly enough towards the roots to respond to transpiration demand. Readily available soil water is the fraction of TAW which a crop can extract from the root zone without experiencing water stress:

(4) $\mathrm{R}\mathrm{A}\mathrm{W}=p\times \mathrm{T}\mathrm{A}\mathrm{W}$(4)

where RAW is the readily available soil water in the root zone (mm) and p is the average fraction of total available soil water (TAW; usually ranging between 0–1) that can be depleted from the root zone before moisture stress (reduction in ET) occurs. For a standard crop, a value of 0.50 for p is common (Allen et al. 1998). The effects of soil water stress on crop ET is considered by multiplying the reference evapotranspiration ET₀ by the water stress coefficient K_s. In this study, the ET₀ used to calculate ET_c was estimated using the Hargreaves equation (Hargreaves and Samani 1985). The water content in the root zone can also be expressed by root zone depletion D_r, i.e. the water shortage relative to field capacity (θ_FC). At θ_FC, the root zone depletion is zero (D_r = 0). When soil water is extracted by evapotranspiration, depletion will increase and induce stress if D_r is equal to RAW. After the root zone depletion exceeds RAW (i.e. the water content drops below the threshold θ_t), the root zone depletion is high enough to limit evapotranspiration to less than potential values, and crop evapotranspiration begins to decrease in proportion to the amount of water remaining in the root zone (Fig. 7(b)). For D_r > RAW, K_s is given by

(5) $K_s=\frac{\mathrm{T}\mathrm{A}\mathrm{W}-D_r}{\mathrm{T}\mathrm{A}\mathrm{W}-\mathrm{R}\mathrm{A}\mathrm{W}}=\frac{\mathrm{T}\mathrm{A}\mathrm{W}-D_r}{\left(1-p\right)\mathrm{T}\mathrm{A}\mathrm{W}}$(5)

where K_s is a dimensionless transpiration reduction factor dependent on available soil water (0–1), and the other terms are those described above. In this simplified model, daily soil water balance can be described as follows (Allen et al. 1998):

(6) $D_{r,i}={D_{r,i-}}_1-P_i+R_i+{\mathrm{E}\mathrm{T}}_{c,i}+{\mathrm{D}\mathrm{P}}_i$(6)

where D_r,i is the root zone depletion at the end of day i (mm), D_r,i−1 is the water content in the root zone at the end of the previous day, i – 1 (mm), P_i is the precipitation during day i (mm), R_i is the runoff from the soil surface during day i (mm), ET_c,i is crop evapotranspiration during day i (mm) and DP_i is deep percolation during day i (mm). Thus, it is also true that

(7) $R_i+{\mathrm{D}\mathrm{P}}_i=D_{r,i}-{D_{r,i}}_{-1}+P_i-{\mathrm{E}\mathrm{T}}_{c,i}$(7)

A water surplus (R_i + DP_i) occurs only when ΔD_r,i (D_r,i _D_r,i−1) < P_i – ET_c,i,. The DP_i is the fraction of the water surplus that can effectively infiltrate and recharge the aquifers.

In this study, TAW was considered equal to the Available Water Content (AWC) of soils in the study area, according to Ouessar et al. (2009) (Table 6). For the Koutine watershed, the soil map of the Jeffara region (1:200 000) (Taamallah 2003), modified by Ouessar et al. (2009), was used to account for the “artificial” soils built up behind the water harvesting structures (jessour and tabias) as deposited sediments. For the other two watersheds, the soil map of the Governorate of Médenine (Carte Agricole 1:200 000) (Ministère de l’Agriculture, Environnement et Ressources Hydrauliques 2003), modified to consider soils of jessour and tabias, was used.

Table 6. Available Water Content (AWC) and depth of soils in the study area (modified from Ouessar et al. 2009).

Soil	Description	Depth (cm)	AWC (mm)
AFFL	Outcropping	10	12
CRCG	Calcimagnésiques sur rendzine calcaire (Rendzinas)	20	20
ISOH	Isohumiques brun calcaires tronqués (Calcic Xerosols)	40	39
JESR	Soil on terraces of jessour	200	301
MBEH	Minèraux bruts d’érosion hydrique (Regosols)	20	24
PEAH	Peu évolués d’apport hydrique (Fluviosols)	200	246
PEEH	Peu évolués d’érosion hydrique (Regosols)	20	24
STAB	Soil on terraces of tabias	100	122

For each sub-basin, an average AWC value was calculated as the weighted average of the AWC values of the soil classes. This considered the weights the area of the surfaces covered by each soil type within the sub-basin.

The daily runoff was estimated by multiplying the water surplus by the runoff coefficient C_k. For this purpose, an average C_k value for each sub-basin was considered. Finally, the daily effective infiltration I_e was obtained from the difference between water surplus and runoff.

3.7 Calculation of the average yearly natural recharge of the aquifers

The average annual natural groundwater recharge for the Triassic, Jurassic, and Cretaceous aquifers were estimated over the 10-year reference period (2003–2012). These water volumes represent the yearly renewable water resources as a result of the rainwater infiltration process. They were calculated by multiplying the I_e of each sub-basin by the surface of each aquifer in each sub-basin. For this purpose, the spatial distribution of each aquifer within the three watersheds was defined (Fig. 8) using an official geological map (1:100 000) (Office National des Mines 1996) and 3D hydrogeological model reconstruction data (Carletti 2017). The Triassic aquifer in the plain, in particular, was considered in the natural recharge calculation as the outcropping Mio-Plio-Quaternary deposits have low thickness, lie directly on Triassic sandstones, host local and perched aquifers, and are in hydraulic communication with the underlying Trias aquifer (Yahyaoui 2007).

4 Results

4.1 Climate data assessment

A comparison between CFSR station 336103 and the Médenine station is shown in Fig. 9. This comparison exemplifies the general results related to the assessment of the reliability of CFSR data and that measured by the Médenine weather station. Among the CFSR stations located in the plain, 336103 has the closest elevation (83 m a.s.l.) to the Médenine station (117 m a.s.l.), and both stations have a similar distance from the sea. The results are summarized as follows: (i) the average daily rainfall data from both stations present no correlation and (ii) the mean maximum and minimum daily temperature data present a very high correlation (R² = 0.98–0.99), even though the minimum temperature values of the CFSR station recorded during summer months are underestimated.

4.2 Distribution of c_k in the watersheds

The most represented classes of the runoff coefficient C_k for all the watersheds (Koutine, Megarine-Arniane, and Hajar), are C_k < 0.2 (frequency from 51.6% to 57.6%), 0.2 < C_k < 0.3 (12.5–16.4%), and 0.3 < C_k < 0.4 (17.8–27.1%) (Fig. 5). The coefficients in the first class are mainly concentrated in the Jeffara plain but are also recorded in the Beni Khedache mountain region, which is characterized by cultivated/shrubby areas and mainly associated with jessour where the slope is quite low (< 10%). Those in the second class are primarily located in parts of the plain corresponding to grazing areas (mixed in the Beni Khedache Mountains with cultivated/shrubby areas) with a slope between 10% and 35%. Areas classified in the third runoff class are ubiquitous and characterized mostly by bare soil conditions in the plain and by grazing/pasture lands, with slopes between 10% and 35% in the Beni Khedache mountain area. In addition, this runoff class also matches outcrops of highly permeable rocks only in the mountainous area of the Koutine watershed.

4.3 Statistical analysis of single runoff factors (Ca, Cv and Cp)

A statistical analysis of the results was performed to assess the influence of the individual factors (C_a, C_v and C_p) on the classification of runoff coefficient C_k (Fig. 10(a–c)). Figure 10(a) shows a graphical representation of the distribution of the slope classes for each class of C_k. For all three watersheds, only slope values less than 10% affect the class with the lowest C_k values (C_k < 0.2), and they decrease to 50% for the C_k values between 0.2 and 0.3. The influence of the slope values between 10% and 35% exceeds 70% for the third C_k class (between 0.3 and 0.4) in the Koutine watershed. Slope values greater than 35% account for 50 – 80% of the C_k values which range from 0.4 to 0.5, whereas they influence 100% the values in the fifth C_k class (between 0.5 and 0.6). Figure 10(b) shows that the cultivated/shrubby land-cover class has a decreasing influence with increasing values of C_k. The grazing land class increases from 45% for C_k values between 0.2 and 0.3 to 100% for C_k values between 0.5 and 0.6 in the Koutine watershed. The no vegetation/bare soil class affects mainly the C_k values between 0.3 and 0.5 in the Hajar and Megarine-Arniane watersheds. In the latter watershed, this class reaches 100% for C_k values between 0.5 and 0.6, although this C_k class represents only 0.01% in the C_k spatial distribution. Figure 10(c) shows the influence of the permeability classes. The good permeability class oscillates between 35% and 85% for C_k values up to 0.5 and, only in the case of the Megarine-Arniane watershed, it reaches 100% for C_k values between 0.5 and 0.6. The medium permeability class ranges from 15% to 55% for C_k values lower than 0.5. In the Koutine watershed, it represents 100% of the C_k values between 0.5 and 0.6. The high permeability class, recorded only on the Beni Khedache mountains in the Koutine watershed, generally has a low influence.

4.4 Runoff and effective infiltration

Water budget estimations carried out on a daily basis, over the 10-year period (2003–2012), allowed for the identification of those rainfall events that generated a water surplus, therefore causing runoff and recharge of the aquifers, for each sub-basin of each watershed. A statistical analysis based on these rainfall events and related to the estimation of runoff and effective infiltration for each sub-basin of the three watersheds was carried out. The results included the following: average runoff (R) (mm), average effective infiltration (I_e) (mm), average precipitation (P) (mm), total number of such rainfall events, average number of days for each event, and number of events of 1 day (Supplementary material, Tables S1, S2 and S3). Furthermore, the average AWC of the soils related to each sub-basin was classified into four classes (mm): Class 1: AWC < 50; Class 2: 50 < AWC < 100; Class 3: 100 < AWC < 150 and Class 4: AWC > 150.

The results related to the Koutine watershed are shown in Fig. 11 (and in the Supplementary material, Table S1). Considering the AWC classes of the soils, they are summarized as follows:

• Class 1 (AWC > 150 mm): the sub-basins included in this class are characterized mainly by deep soils (1 or 2 m) corresponding to the “artificial” soils of tabias and jessour. They are located in the plain (K1 and K2) and in the northern part of the watershed (K6, K7, K9 and K10). These sub-basins did not contribute to R and I_e during the reference period;

• Class 2 (100 < AWC < 150 mm): these sub-basins are distributed across the entire watershed. The average R varies from 0.55 to 5.50 mm; the latter value corresponds to sub-basin K21 which is characterized by the highest value of C_k (0.34). The average I_e varies from 1.86 to 17.30 mm. This highest value is related to sub-basin K17, which has almost the same characteristics (AWC, average rainfall, number of rainfall events, etc.) as K21, except for the value of C_k, which is around the half (0.164). Within this category, for those sub-basins having an average AWC between 100 and 120 mm, the rainfall events that produced a water surplus vary between 90 and 100 mm, with an average duration of about 3–4 d; however, only 3 events occurred in 10 years, of which only one event lasted a day. For the other sub-basins included in this AWC class, the average rainfall is as high as 138 mm, with an average duration of 6 d and only 1 event occurred in the reference period;

• Class 3 (50 < AWC < 100 mm): the two sub-basins (K18 and K19) included in this class are mainly concentrated on the Dahar. In this case, the average rainfall that generated R and I_e decreases to 66–68 mm and the number of 1-day events increases to 4–5; and

• Class 4 (AWC < 50 mm): this class includes the sub-basins that contributed the most to R and I_e (K5, K8, K11, K13 and K14). They are concentrated in the central-eastern part of the watershed and are covered by shallow soils (depth between 0.20 and 0.40 m). The average R varies from 2.32 to 5.79 mm, whereas the average I_e varies between 16.00 and 18.58 mm. The total number of rain events that produced a water surplus was between 24 and 42, with average rainfall between 32.96 and 46.68 mm and an average duration of about 2 d. The 1-day rainfall events increased to 18.

Figure 12 (and the Supplementary material, Table S2) shows the results for the sub-basins of the Megarine-Arniane watershed, which are described below:

• Class 1 (AWC > 150 mm): the sub-basins belonging to this class are located in the central-eastern part of the plain and did not contribute to R and I_e during the period 2003–2012;

• Class 2 (100 < AWC < 150 mm): these sub-basins are distributed across the entire area. The average R varies from 0.13 to 4.87 mm and the average I_e varies between 0.45 and 17.70 mm. In general, as AWC values increase, the average rainfall that produced a water surplus increases and the number of events decreases: for AWC values over 120 mm, the average rainfall exceeds 110 mm and only 1 event occurred;

• Class 3 (50 < AWC < 100 mm): these sub-basins are located in the central and northern part of the watershed. The average rainfall that generated R and I_e decreases to 64.23 mm and the 1-day event number increases to 6–7; and

• Class 4 (AWC < 50 mm): only sub-basin M4, which is located in the northern part of the watershed, is included in this class. The average R is 4.79 mm, and the average I_e is 10.34 mm. The average rainfall amount is 31.26 mm, and the total number of rain events that generated a water surplus is 45, with an average duration of about 2 days/event, whereas the 1-day rainfall events are 19.

The results for the Hajar watershed are reported in Fig. 13 (and in the Supplementary material, Table S3), and are summarized as follows:

• Class 1 (AWC > 150 mm): only sub-basin H3 belongs to this class and the same results described above for similar sub-basins of the other watersheds are reported;

• Class 2 (100 < AWC < 150 mm): these sub-basins are located in the central-western part of the watershed. H2 and H4 have almost the same value as C_k, but the latter contributed with a higher value of average I_e (17.23 mm) owing to an unique rainfall event that lasted seven days and was characterized by a cumulative precipitation of about 138 mm;

• Class 3 (50 < AWC < 100 mm): this class includes only sub-basin H1, located in the northern part of the watershed. It contributed with an average R of 5.98 mm and an average I_e of 14.64 mm; 13 rainfall events were recorded, 7 of them being concentrated in 1 day; and

• Class 4 (AWC < 50 mm): No sub-basin belongs to this class.

The annual R (mm) and the annual I_e (mm) were estimated over a 10-year period (2003–2012) for each sub-basin of the three watersheds. In particular, the average R and I_e for each year for each watershed were calculated as a weighted average of the values of each sub-basin considering the areas of each sub-basin as weights. The charts in Fig. 14 show the trends of the average annual I_e for each watershed.

Figure 14(a) shows the trend of the average annual I_e of the Koutine watershed, with values ranging from 65.10 mm in 2003, which was a particularly rainy year (610 mm), to 0.08 mm in 2006 (120 mm). The charts in Fig. 14(b) and (c) show the annual trend of average annual I_e for the Megarine-Arniane and the Hajar watersheds, respectively. The trends are very similar, with the highest value in 2003 (24.35 mm for the Megarine-Arniane watershed and 19.03 mm for the Hajar watershed). For the years 2004–2009, the average annual I_e was practically zero. The higher values for the Koutine watershed, compared to those of the other two watersheds, are due to the watershed having the highest number of sub-basins characterized by AWC < 50 mm, which are those that gave the highest contribution in terms of annual average I_e.

4.5 Average yearly natural recharge of the aquifers

The average annual natural recharge volumes (m³ year⁻¹) for the aquifers are given in Table 7.

Table 7. Average annual natural groundwater recharge (m³ year⁻¹) for each aquifer. AANR: average annual natural recharge (m³ year⁻¹); Area: extent of the aquifer in the watershed (km²).

Watershed	Aquifers
	Triassic		Jurassic		Cretaceous
	AANR	Area	AANR	Area	AANR	Area
Koutine	3 198 250	170.244	119 176	13.971	223 494	78.861
Megarine-Arniane	904 219	186.405	85 399	44.201	12 779	16.784
Hajar	360 196	61.891	16 765	17.393	782	1.000
Total	4 462 665	418.540	221 340	75.565	237 055	96.645

Although the areal extensions of the Triassic aquifer in the Koutine (170.2 km²) and the Megarine-Arniane (186.4 km²) watersheds are comparable (Fig. 8), the average annual groundwater recharge differs consistently (3.2 and 0.9 hm³ year⁻¹, respectively). This is because the five sub-basins (K5, K8, K11, K13 and K14), characterized by AWC < 50 mm and the highest values of I_e in the Koutine watershed, contributed only to the groundwater recharge of this aquifer, whereas the Megarine-Arniane watershed presents only one small sub-basin (M4) with the same characteristics. The contribution to the annual recharge of the other aquifers is 1–2 orders of magnitude lower than that to the Triassic aquifer, and up to 4 orders of magnitude lower in the case of the Cretaceous aquifer in the Hajar watershed. Thus, the total estimate of the average annual groundwater recharge for the Triassic aquifer amounts to about 4.5 hm³ year⁻¹. This value can be considered consistent with the unique reference values of about 4.7 hm³ year⁻¹ (150 L s⁻¹) and 4.9 hm³ year⁻¹ (155 L s⁻¹) for the same aquifer in the study area, as respectively reported by Ouessar (2007) and Yahyaoui (2007).

Figure 8. Spatial distribution of the Triassic, Jurassic and Cretaceous aquifers within the three watersheds.

Figure 9. Comparison between rainfall and temperature data of the CFSR station 336103 and the Médenine station: (a), (b), (c) correlation between average daily rainfall, max temperature and min temperature, respectively; (d), (e) trends of average daily max temperature and min temperature, respectively.

Figure 10. Distribution of (a) C_a, (b) C_v and (c) C_p versus C_k for the three watersheds.

Figure 11. Estimation of average runoff and average effective infiltration, calculated on the basis of the rainfall events generating a water surplus (2003–2012), for each sub-basin of the Koutine watershed.

Figure 12. Estimation of average runoff and average effective infiltration, calculated on the basis of the rainfall events generating a water surplus (2003–2012), for each sub-basin of the Megarine-Arniane watershed.

Figure 13. Estimation of average runoff and average effective infiltration, calculated on the basis of the rainfall events generating a water surplus (2003–2012), for each sub-basin of the Hajar watershed.

Figure 14. Trends of average annual effective infiltration (2003–2012) related to the (a) Koutine, (b) Megarine-Arniane and (c) Hajar (black line: annual rainfall).

5 Discussion

The runoff coefficient method was tested and validated in previous studies (Spadoni et al. 2010, Ghiglieri et al. 2014, Grillone et al. 2014) on semiarid regions in the Mediterranean area, which are typically characterized by a lack of sufficient historical climatic and streamflow datasets that can be used for the implementation, calibration, and validation of hydrological models. By multiplying the runoff coefficient C_k with the effective rainfall (P-ET) in the standard water balance equation on an annual basis, the annual runoff R and, therefore, the annual effective infiltration I_e (natural aquifer recharge) were calculated. This kind of approach is not suitable for ungauged watersheds in arid environments, such as those of the Oum Zessar study area in southern Tunisia, where the water balance is almost negative throughout the year. Indeed, in such climatic conditions, only after very intense rainfall events concentrated in a few days, does the rainfall exceed the evapotranspiration, generating runoff and effective infiltration. To overcome this critical condition, we developed an innovative approach, that combines the runoff coefficient method with a water balance model, which allowed us to calculate the effective infiltration I_e on a daily time scale basis as part of the water surplus in the soil. This calculation was carried out at the sub-basin scale for three watersheds (Koutine, Megarine-Arniane, and Hajar) by considering a mean AWC value and a mean value of the runoff coefficient C_k for each sub-basin. The result, which emerged very clearly, is that sub-basins characterized by average AWC values greater than 150 mm did not contribute to runoff and effective infiltration during the reference period. They are covered mainly by the “artificial” soils of tabias and jessour, the traditional water harvesting structures that are widespread in SE Tunisia. These soils are deep (1 or 2 m) and characterized by high AWC values (122 and 301 mm, respectively); thus, only very intense rainfall events can saturate them and generate a water surplus, especially after dry periods. Conversely, the sub-basins covered by shallow soils (depth between 0.20 and 0.40 m) and characterized by average AWC values less than 50 mm are those that produced the greatest contributions to the effective infiltration. This result emphasizes the delicate issue of the balance between the diffusion of the traditional water harvesting systems that have the primary function of promoting water storage in the soil for agricultural purposes and the need to ensure an adequate natural recharge for aquifers.

Starting from the calculated values of annual effective infiltration I_e and through the reconstruction of the spatial distribution of the Triassic, Jurassic, and Cretaceous aquifers within the three watersheds, the calculation of the average annual natural groundwater recharge for each aquifer was carried out. The results revealed that the Triassic aquifer is characterized by a much larger volume of average annual groundwater recharge than the other two aquifers. This difference does not only depend on its greater areal extension, but also on the fact that all the sub-basins with an average AWC < 50 mm contribute the highest value of average annual I_e to the recharge of this aquifer. This value is consistent with the volumes reported by Ouessar (2007) and Yahyaoui (2007), which represent the only available reference values for this aquifer in the study area. A comparison with the methodology proposed herein was not possible because these authors did not describe the methods they used for the natural recharge estimation in their works. Furthermore, the unavailability of a historical series of piezometric maps of the aquifers for the reference period (2003–2012) prevented us from evaluating the effects of the annual groundwater recharge on the aquifers to assess their overexploitation, as stated in WADIS-MAR (2016).

A shortcoming of the application of the methodology in the Oum Zessar study area was the use of the climate data (daily precipitation and mean temperature data) of the Médenine weather station for the entire area, as they are the only available measured data on daily time steps. In addition, the unavailability of streamflow data measured at the watershed outflow point prevented us from making a comparison between the calculated runoff coefficient C_k values and those based on measured data in order to validate the method. To improve this study and overcome these constraints, it would be appropriate to increase the number of weather stations (at least 5–6) distributed over the entire area. For this purpose, on the basis of our results, it was possible to identify the most suitable locations for new weather stations that should integrate the existing network managed by local institutions, which at the moment, provide, at a maximum, 10-year series of climate data. Concerning the streamflow measured data, it should be necessary to realize at least three stream gauge stations, one for each basin. As a suggestion, the station in the Koutine wadi could be restored and other two could be realized at the edges of the Megarine–Arniane and Hajar watersheds: the former on the wadi on the bridge on the Koutine–Médenine road, and the second one on the bridge in the city of Médenine.

This study assumed that the methodology could be applied to a delimited area corresponding to the three watersheds, even if they cover only a part of the extension of an entire aquifer. In particular, as reported by Besbes et al. (2005), the Triassic aquifer extends to the south-east, beyond the Libyan boundary. The choice was made to refer to a well-known physiographic unit, represented by the watershed, because it has not yet been possible to identify the groundwater divides of each aquifer in the study area. However, among the results achieved by WADIS-MAR, there is some geological, hydrogeochemical and, isotopic evidence which suggests that it is very likely that watersheds and hydrogeological basins do not coincide (Carletti 2017). Three-dimensional hydrogeological model reconstruction highlighted that the tectonic structures in the area significantly affect groundwater flows. A further confirmation comes from the results of the isotopic characterization of groundwater: both the Triassic and Jurassic aquifers are characterized by Tritium concentrations that indicate a mixing between recent and old waters. It could mean a mixing of recently infiltrated water with groundwater characterized by a greater residence time owing to the distance of the recharge areas.

6 Conclusions

The proposed methodology provides a reasonable, but conservative, natural aquifer recharge appraisal in ungauged watersheds, particularly in arid climatic conditions, where the lack of appropriate field datasets prevents the use of hydrological modelling. On the one hand this can be considered its greatest strength, but on the other hand, this methodology shows some weaknesses. The main one is being able to calculate the daily water balance with a resolution that is at most the sub-basin scale, considering average AWC and C_k values for each sub-basin. Indeed, while the C_k is calculated through GIS techniques with the same resolution as the DEM, the runoff and the effective infiltration on a daily time scale are calculated by the water balance model implemented in Excel, and average values per sub-basin are obtained. The implementation of the water balance model in a GIS environment would allow the effective infiltration to be calculated with greater resolution, and therefore, more detailed information on its spatial distribution within the watershed could be obtained. Despite these limits, the application of this methodology in the Oum Zessar study area showed that the land and water management, in particular the wide use of the traditional water harvesting techniques, has had some impacts on the natural aquifer recharge. As a future development of this research, this kind of assessment could be deepened through the implementation of models that are not only hydrological or agronomic but rather Decision Support Systems (e.g. SWAT model). Considering that they include different components (i.e. hydrology, weather, soil erosion, plant growth, land management, and groundwater), they represent a useful tool to simulate several processes at watershed scale on daily time steps. Their implementation needs robust datasets of climatic data and streamflow data, which are essential for the input, calibration, and validation of the models, and are not available at the moment in the study area. Furthermore, the natural aquifer recharge estimation in the study area can be improved considering the hydrogeological basins instead of the watersheds. To achieve this aim, a more detailed geological and hydrogeological study of the area is needed to delineate the hydrogeological basins. In addition, this in-depth analysis will allow us to better understand the groundwater circulation and assess the possible hydraulic connections between aquifers, which imply the presence of a lateral recharge that could significantly integrate the natural aquifer recharge with the effective infiltration.

However, the methodology presented provided natural aquifer recharge estimation which can be used in the future for several hydrogeological applications:

• validating the conceptual model of groundwater flow and the groundwater recharge area delineation;

• supporting the assessment of groundwater vulnerability to the pollution mainly due to human activities and the design of pollution mitigation and aquifer protection measures; and

• implementing managed aquifer recharge (MAR), which is a process of induced groundwater replenishment (Bouwer 2002, Dillon et al. 2009, Fernández Escalante and San Sebastián Sauto 2012, Cosgrove and Rijsberman 2014). Indeed, the estimate of the natural groundwater recharge is very important, in combination with hydrogeological and hydrogeochemical evaluations, to the design of MAR systems and to defining their most suitable location. Furthermore, it allows researchers to evaluate the impact of the implemented MAR systems on the aquifer and thus to assess their efficiency.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary Material

Supplemental data for this article can be accessed here.

Additional information

Funding

This study was carried out within the WADIS-MAR demonstration Project (2011-2016) (www.wadismar.eu) [Grant: ENPI/2011/280-008], funded by the EC under the Regional Programme “Sustainable Water Integrated Management” (SWIM) (www.swim-sm.eu), implemented under the European Neighbourhood and Partnership Instrument (ENPI), and coordinated by Prof. Giorgio Ghiglieri. It formed a part of a PhD study in Earth and Environmental Sciences and Technologies, titled “Trial of protocols and techniques for integrated groundwater management in arid and semi-arid regions to combat drought and desertification” carried out at University of Cagliari. The authors also thank the partial support of the RE-LIVE WASTE Project (www. re-livewaste.interreg-med.eu), funded by the Interreg MED Programme.