Integrated Groundwater Flow Modeling for Managing a Complex Alluvial Aquifer Case of Study Mio-Plio-Quaternary Plain of Kairouan (Central Tunisia)

Abstract

In central Tunisia, anthropic activities, such as groundwater abstraction for irrigation, have resulted in excessive groundwater level declines of the Mio-Plio-Quaternary aquifer hosted in Kairouan Plain. Besides, the two dams El Houareb and Sidi Saad’s impoundment since the 1980s has deeply modified the natural process of aquifer recharge. Hence, some studies claim the dam’s instauration of this groundwater depletion; however, some other studies attribute this critical situation to an issue of groundwater management. A multidisciplinary study was carried out to retrace the groundwater flow dynamics for 48 years before and after the dams’ erection and to understand the main factors causing the groundwater depletion. Hence, a conceptual model was developed based on gathering all available data from 114 borehole logs, 10 seismic lines, and 08 petroleum wells. Based on this reconstructed geometry, the groundwater head was simulated using the numerical code Modflow. The model was calibrated in steady-state with reference to the piezometric levels measured in 1969 and in the transient state for the period 1970–2017 and validated for the period 2007–2017. The outputs of the calibrated model show a relevant finding of the decrease of the inflows coming from the rivers’ beds (Zeroud and Marguellil) from 1990 to 2017 by 48%; yet, the pumping rate has increased by 119%. The simulated scenario without dams and maintaining the same withdrawals has shown a groundwater level rise downstream of the plain; yet, in its upstream, the depletion was less intense compared to the current model. However, the case of doing without dams and raising withdrawals from the aquifer has generated a huge decline reaching 22 m near Draa Affane.

1. Introduction

In the past decades, the Kairouan plain basin (Central Tunisia) has faced a serious issue of groundwater resources degradation, a problem that is particularly emblematic in semi-arid and arid regions (e.g., Ethiopia [1], Egypt [2]). Hence, the economic growth of the plain relies entirely on the exploitation of groundwater for both irrigated agriculture and the supply of drinking water.

The large Mio-Plio-Quaternary alluvial aquifer of the Kairouan plain represents a key resource with the highest potential for further agricultural development in the region. For more than 40 years, the basin has been shaped by regional planning strategies including the creation of numerous water and soil conservation works, hillslope lakes, and big and small dams [3]. Before the construction of large dams, this alluvial reservoir used to be recharged by the direct infiltration through the beds of the Zeroud and Marguellil main wadis, especially during flooding events. Undeniably, deep percolation through the stream beds was the main source of groundwater recharge [4,5,6,7].

Since the 1980s, to protect the Kairouan city from hazardous floods and to mobilize groundwater resources essentially for irrigations purposes, two dams have been erected, namely Sidi Saad and El Houareb. These dams have deeply modified the groundwater cycle both in terms of fluxes and processes [8,9] and the floodwater hardly reaches the aquifer [6,10], “erection of large dams and erosion control structures moved and reduced the recharge fluxes [11]. Consequently, the recharge process is limited to El Houareb dam leakage and Ain El Beitha flow infiltration through the karst formation [8,11,12]; and hence, the Kairouan plain aquifer (downstream of the dams) has started experiencing a net groundwater depletion [13,14].

Although the dams’ construction constitutes the only visible factor underlying the water storage degradation, overexploitation and legislation related to groundwater access are pointed out by recent studies as the main factors.

For sustainable management of water resources on a regional scale, this integrated study aims to retrace the dynamic flow paths before and after the impoundment of the two dams as well as to identify the factors that could threaten the availability of groundwater.

In the literature, methods developed to simulate flow dynamics paths of an aquifer and evaluate its flow budget vary as a function of the available inputs data and also of the modeling approach [15]. Numerical modeling could be a helpful tool for the assessment of groundwater resources and the analysis of future management scenarios [2]. The conception of numerical models is based on the available geological information and on the aquifer system conceptualization, taking into account the existence of errors relative to hydraulic parameterization [16].

Therefore, in a first step, it is essential to develop a conceptual and groundwater flow model, starting from borehole data, well-pumping tests, and well monitoring, so as to investigate the general groundwater flow behavior and also the main flow budget components as well as to simulate the groundwater head from 1969 to 2017. This needs a straightforward model, which should be able to capture the general groundwater flow behavior within the Kairouan plain aquifer, without introducing complex fine-resolution even though the real system is heterogeneous. The conceptual and numerical models, although simple, have to be robust and hydro-stratigraphically well based. Hence, to satisfy the requirement of simplicity and robustness, a multi-dimensional approach is proposed and illustrated for the Mio-Plio-Quaternary aquifer.

2. Materials and Methods

2.1. Study Area

Kairouan plain (Figure 1) occupies the southern half of the basin of Djebibina-El-Alem-Kairouan-Nsar-Allah. It extends from the town of Kairouan within the North to the Pontian undulations of the Djebels Khordj and BouThadi [5,7]. It is restricted to the west by a series of reliefs with secondary formations (Djebels Baten, Cherichira, El Houareb, Draa Affane, Djebel Cherachine) and to the east by the Sebkhas (salt lakes) El Kelbia, El Hani, and Cherita [5]. Hence, due to its morphology ranging between mountain and plain, this basin, situated in a semi-arid climate where mean annual rainfall from 1978 to 2010 ranges from 400 mm in the north and 200 mm in the south [5,6,17], is characterized by high variability of precipitation in time and space [12,13,17,18]. Recently, [18] have studied the rainfall trends of 07 rainfall stations in Central Tunisia in order to defend or claim the hypotheses of precipitation decrease. This study attests that the regionalization of rainfall decreases in this area is not evident over time and space.

The alluvial plain, foremost, hosting one of the largest reservoirs in central Tunisia, has been continuously considered a very attractive scientific site to study, thanks to its huge collapsing basin, which exceeds quite 700 m in depth [4,13,19,20], and up to 500 m of saturated thickness [7,11]. Bordered on the west by a chain of mountains, it is supplied by two famous intermittent wadis, namely Zeroud and Marguellil, recognized by their seasonal flow [11,14,15]. Yet, not perennial, through the beds of those wadis, a recharge is supposed to result from runoff infiltration, and it takes place during rainy events [12,18,21]. Before the construction of the two dams of Sidi Saad in 1982 and El Houareb in 1989, the two wadis were running through the plain, sometimes feeding the Sabkha Kelbia (a very large endorheic Salt Lake) during extreme floods [17]. Indeed, this Sabkha was considered the natural drain for both surface water and groundwater by [4]. After the devastating floods in 1969, and in the purpose of protecting the Kairouan city, reducing erosions, and mobilizing water resources, [21], little water was actually stored and the leakage in the underneath fractured limestone could not be controlled, inhibiting flood releases and keeping a very localized contribution to recharge [8,12,13]. These hydraulic structures have led to a decline in runoff of around 45% over the period between 1996 and 2011 [22], Since the 1980s, adding to this recharge decline, groundwater abstraction from this thick aquifer has been the major water supply to irrigated agriculture.

Characterized by a complex geological structure and a highly varied lithological cover [23,24] this collapsed basin is filled with detrital sediments of the Plio-Quaternary age [11]. The Mio-Pliocene formations mainly consist of sands, silts, clays, marls, and conglomerate strata which were derived from the erosion of secondary and tertiary age deposits surrounding the basin [11,22,24]. The Plio-quaternary is composed mostly of alternations of silt and sand, further as coarse sand with clays and crusted soil layers [12,17,25]. Those deposit sediments are composed of alternate permeable layers and aquifers (sand, sandstone) and of impermeable layers (clay) [24,26].

2.2. Groundwater Model Development

• Hydrogeological database

Gathering all available data starting from geological data (borehole logs, seismic wells) in order to conceptualize the aquifer geometry and define the boundary conditions, to the hydrological data, including groundwater abstraction reports, piezometric signals, climatic data, and pumping tests, a whole complex hydrogeological database has been created. These databases, summarized in Figure 2, provide an essential knowledge basis for characterizing aquifers and assessing their water resources to make valuable decisions [27,28]. However, coming from different sources (

Table 1. Sources of available data for numerical modeling of Mio-Plio-Quaternary aquifer (Kairouan plain).

Data	Sources	Period	Characteristics
Precipitation	INM 3 DGRE	1956–2011 1983–2016	Kairouan station 04 stations
Hydrodynamic parameters.	CRDA 1 DGRE 2	*	43 pumping Test
Recharge	(a) (a) × (b)	1969–1982 1982–2017	*
Piezometer	Thesis [4] CRDA, DGRE reports	1969 1970–2016
Exploitation	* DGRE 2 annual reports Thesis [4]	1970–2016 1969–1970	*

(a) Infiltration Coefficient = 5%, [11]. (b) Runoff loss = 45% after dam’s construction, [22]. * No specific data. ¹ Regional Commission of Agricultural Development of Kairouan. ² General Direction of Water Resource. ³ National Institute of Meteorology.

), this heterogeneous gathered data needed to be harmonized in time and space using GIS tools as they are very powerful for processing, analyzing, and integrating spatial data sets [29,30].

• Model Set up

To deal with a very complex aquifer numerical model, it is imperative to balance what is sufficiently simple and robust to represent reality with a satisfactory degree of reliability [14,31]. Among the most used numerical codes for hydrogeological modeling, the modular finite-difference groundwater flow model MODFLOW, version 2000, [32] developed by USGS (United States Geological Survey) has been used to simulate the groundwater head (GRH) and groundwater flow movement as well as to estimate the water budget [33]. Thus, this code defines the major characteristics of the studied system and simulates ground-water flow in aquifer systems using the finite differences method [32]. The hydrodynamic parameters of the steady state were calibrated based on a [4] piezometric map drawn for January 1969, before the flooding in October and November 1969, which is considered representative of a stable state of the aquifer flows.

The model was designed on a regular grid of 1000 m × 1000 m with a total of 2500 cells. In this study, the numerical model was supposed to be a monolayer aquifer referring to the funding of the 3D conceptual modeling. The top of the layer coincides with the topographic surface; however, the bottom is delineated by the Miocene formation. Obviously, the MODFLOW model calculates the transmissivities and vertical leakage terms from the thickness of the model layers implicitly, given by the top and bottom elevation of the layer [2].

• Aquifer exploitation

Because of the free and uncontrolled access to the groundwater in the Kairouan plain aquifer, even though restriction laws exist, the exact estimation of the plain exploitation seems to be difficult [34]. Indeed, referring to a note from the DGRE (General Direction of Water Resource) annual report from 1969 to 2016, this difficulty comes from the intensification of irrigated agriculture followed by the implementation of deep wells exceeding 100 m in depth especially at the plain upstream (Marguellil). This proliferation of deep wells, coupled with intense pumping flows, has led to overexploitation. Actually, farmers denied this over-usage of the aquifer; however, they have experienced a fall in piezometric levels, “It is hard to ignore because the change has been fast and widespread (a drop of 1–2 m/year)” [14].

The total estimated pumping volume (Figure 3), from shallow and deep wells per year, has been increased spectacularly from 31.3 in 1970 to 111.43 Mm³ in 2016. The second-degree polynomial trend of this increase is about y = 0.0564 x² − 222.79x + 220,208 with a high correlation coefficient (R² = 0.9662). Shallow pumping is defined as the abstraction from shallow wells where its depth was less than 50 m [35].

2.3. Boundary Conditions

The hydrological conditions along the boundaries of the conceptual model determine the mathematical boundary conditions of the numerical model [27].

Indeed, the western geological boundary of the plain (Jbel Cherichira) on the El Houareb threshold (Figure 4), according to [36], was supposed to be connected in-depth with the plain of Kairouan. Indeed, Ref. [37] has evidenced this discharge of the Ain El Beidha aquifer into the Kairouan plain through the hydraulic junction of Plio-Quaternary overlaps. Furthermore, referring to [38], the southern compartment of the aquifer of sabkha El Behira plain converges towards the extreme southwest of the Kairouan plain. Consequently, these two limits along Jbel Cherichira and Jbel Cherachine are considered as a specific head. Hence, in the eastern part of the region, the Kairouan plain is bordered by two topographical depressions corresponding to the two sebkhas: El Kelbia and Sidi El Hani [21]. has displayed the groundwater flow lines shifting to the NNE (North-Northeast) direction towards the topographic depression of the sebkha El Kelbia. This very large endorheic Salt Lake was considered the natural outlet and drain for both surface water and groundwater by [7,35]. On the other hand, Ref. [13] has considered a possible discharge of the plain into sebkha Sidi El Heni and Cherita, thus, a limit of drain has been imposed near each sabkha depression.

2.4. Groundwater Water Scenarios Set Up

To detect the impact of the two dams’ impoundment on the groundwater level trends and what the situation would be if the two Sidi Saad and El Houareb dams were not constructed beside the impact of the raise in pumping rates on this decline, two scenarios were defined based on existing data on dams and aquifer withdrawals.

• Scenario 1: The dams were not built and we maintain the same withdrawals values taken from the official data of DGRE (Figure 3).
• Scenario 2: The dams were not built and we raise withdrawals from the aquifer; this raise is equal to the removals volume from dams El Houareb and Sidi Saad based on the reports of DGRE and the study of [7] for irrigation.

3. Results and Discussions

3.1. 3D Geological Conceptual Model Design and Aquifer Geometry Building

Most groundwater problems are addressed with a mathematical model developed from the conceptual model. In developing a conceptual model, relevant hydrogeological data (Figure 2) have been organized, analyzed, and synthesized, with the help of a database tool such as a geographical information system (GIS). Hence, a conceptual model can also be constructed for a generic geologic setting [27].

Owing to its position between two geodynamic domains, namely the central Atlas to the west and the Sahel plain to the East, the complex tectonic structures of the study area [11,25], crossed by numerous faults, make the identification of the lithostratigraphic units difficult. Therefore, it was essential to use seismic data. Thanks to subsurface interpretations given by the correlation between petroleum wells from West to East (direction of Cherachine) to the plain, the cross-section AA’, illustrated in Figure 5, has clearly shown a sequence of sedimentary units consisting mainly of detrital deposits of Oligocene, Miocene Pliocene, and Quaternary ages, overlaying the Upper Cretaceous basement rocks.

Ain Grab and Mahmoud Formations (Mio-Pliocene) are represented by thick sandy and clayey layers, respectively [39]. The thickness of Mio-Plio-Quaternary aquifer (Segui formation) deposits could reach more than 900 m in the KTTIR(KT) synclinal, between Jbel Cherahil and the anticlinal of Nasrallah.

Furthermore, using 10 seismic lines (Figure 1) described by [26], coupled with all available cross-sections and correlation between petroleum wells, the architecture of the productive horizons has been revealed in 3D dimension (Figure 6). These seismic lines provide information on the depth of the Miocene clay considered as the substratum of the Mio-Plio-Quaternary aquifer of Kairouan. In addition, these seismic lines provide information on the aquifers set into a depth exceeding 2000 m, which allows the distinction between the different aquifer layers in the Kairouan plain (Figure 6).

A 3D model has top and bottom as well as side boundaries and conditions [27]. Therefore, building a complex 3D geometry needs a clear and detailed detection of the substratum. Based on the logs of 114 boreholes (Figure 1), 08 petroleum wells, 10 seismic lines, and the 3D conceptual model (Figure 6), the bottom of this collapsed aquifer was the Miocene. The reflectance of the deep oil wells has allowed us to detect and prove the shape of collapsed geometry of the alluvial aquifer reaching more than 900 m in its middle part. Even though the cross-section BB’ has revealed several horizons along the Zeroud wadi, it is still very difficult to separate two aquifers, which may, indeed, support the hypothesis of considering the Kairouan plain as a monolayer with clayey lenses (Figure 7). Confirming these findings, the cross-section correlated between three piezometers (Figure 7) has shown a very large variability between clayey sands and sandy clay, and as a result, the lenticular nature of the aquifer formed by alternating more coarse sands and clays over 500 m thick. Despite the uncertainty about a net separation of productive horizons [34] confirmed by [5,8], especially in the plain upstream, Ref. [4] has considered this aquifer as two productive horizons—deep and shallow—connecting by depth.

3.2. Groundwater Flow Model Calibration and Validation

The model was calibrated for steady and transient state conditions as an evolutionary process in which successive adjustments and modifications to the model are based on the results of previous simulations [40]. The transient model was to run until the year 2017 to simulate the dynamic groundwater flow, calibrated, and then validated to predict the aquifer reaction under various scenarios of groundwater management. Validation of the current model was based on the period from 2007 until 2017.

• Steady-state

The calibration of the model in steady state was based on the piezometric map of [4] in January 1969. The flow model was calibrated by adjusting several parameters (permeability, recharge, river stage) within a narrow range of values until the best fit was obtained between the observed heads and simulated heads [2,40]. This parameters adjustment was conducted under steady-state flow conditions based on [2,35] findings.

The model was able to reproduce correctly the observed head in 37 piezometers (Figure 8a) with a high correlation coefficient (R² = 0.93) and acceptable Root Mean Square (RMS = 2.44 m). To ensure a good agreement of the model, the calibration was piloted using automated parameter estimation PEST. As a result, this code has generated dissimilar hydraulic conductivities, compared to initial values. For instance, those initial data are taken from pumping tests interpretation of 43 boreholes and ranged between 4 m/day and 20 m/day. Yet, the calibrated hydraulic conductivities have revealed that the permeability values were unequally distributed over space, suggesting a very heterogeneous sandy clayey medium. It is spanning between 1.5 m/day and 39 m/day, and the outcropping of Daraa affane is the most permeable area in the plain upstream [20]. This variability in the alluvial plain could be a result of its lenticular texture embedded with clayey lenses. Therefore, it is difficult to separate clearly two productive horizons, which is consistent with the findings of [6].

To correctly calibrate the model, reliable quantification of groundwater pumping in 1969 was a challenging key. Indeed, almost all abstraction points were traditional surface wells [4,34,35]. The relative data of these wells were gathered from the study of [4] and DGRE reports (Table 1). The volume of groundwater withdrawals was estimated to be 30.2 Mm³ confirmed by [7]. Moreover, thanks to the available data in [35], the aquifer inputs prior to the dam’s launch have been assessed. Thus, the annual contributions to the aquifer were 6 Mm³ by Zeroud stream beds, of 6.6 Mm³ by the Marguellil, including the underground contribution of 3.6 Mm³ due to a connection with the Aïn Beidha aquifer, at the level of the Marguellil wadi.

• Transient state

The calibration of the current model, in a transient state, was based on the groundwater level of 11 piezometers. These piezometers had been monitored since 1968 [36] and had almost continuous head data until 2016. For 2017, a groundwater head measurement was tackled in November 2017 (14 observation points) in order to ensure reliable data for validation (Figure 8d).

The hydraulic heads and the flow patterns for the alluvial aquifer in 1969, 1990, 2003, and 1994 are shown in Figure 8. Among these years, 1981 was the driest with an annual rainfall of 176 mm, and 1990 was the wettest, with an annual rainfall of 514.5 mm.

Scatter plots of differences between observed and computed groundwater levels show mean residuals (RMS) of 3.11 m, 2.87 m, and 3.22 m, for 1981, 1990, and 2017, respectively (Figure 8b–d). Considering the extent of the study area, these results indicate a reasonable agreement between simulated and observed hydraulic heads, with calibrated hydraulic conductivity values. The coefficients of correlation show acceptable values ranging between 0.81 and 0.87. The general groundwater flow path is from SW to NE, following the general topography (Figure 8), with a mean hydraulic gradient of 2.5 × 10⁻³. In the eastern part of the region, the groundwater flow lines shift to the NNE direction towards the topographic depression of the sebkha Kelbia. This finding was proved by many studies [13,15,21].

After the two-dam erection, the CGWL for 1990 (Figure 8c), compared to CGWL for 1981(Figure 8b), has shown a relevant change, where a decrease of 3 m, especially in the east part of the aquifer, is observed, except near to the El Houareb reservoir, where an increase of 4 m is shown. Indeed, this change in GWL pattern could be justified by the findings of [6,21] that the runoff is stored in the dam reservoir, except for the very exceptional reservoir releases (13 × 10⁶ m³), hence the deep percolation is localized in the vicinity of the El Haouareb reservoir.

As a result, an overall decrease is observed from 1969 until 2017. This finding is proved by many studies [13,14]. This decrease is about 19 m. Indeed, this estimation varies according to the study, but almost this decline exceeds 10 m.

• Analysis of GWL fluctuations from 1970 until 2017

To strengthen the analysis of the groundwater flow trends, three piezometers were selected to highlight the groundwater change from 1969 until 2017 (Figure 9). The location of the three piezometers was illustrated in Figure 1.

The findings of head level fluctuations have shown in piezometer (12737/4) (Figure 9A) a net decline in observed head illustrated by a negative slope (−0.4 m/year). From 1970 to 2017, this decline is measured at about 19 m, pointing to an alarming urgent case of water table drop. Towards this decline, farmers of Kairouan plain found a solution: chasing the groundwater level by digging more of their deep wells. On the other hand, in comparison between the observed head and the simulated one in this piezometer, the overall correlation is estimated at about 0.944. However, the mean error was estimated at about +2 m before the dams’ construction; yet, since 1990, the error estimation reaches its maximum (+5 m) in 2017. Even after calibration, this difference between observed and simulated heads is a result of a gap between archives and real exploitation of the aquifer. Indeed, Ref. [35] defend this explication and show that the official pumping rates may be incorrect or that additional individual groundwater withdrawals occurred. As matter of fact, since the proliferation of deep wells, it is very delicate to define even the number of wells.

Generally, dams are expected to raise groundwater levels, especially near the lake [41]; yet, the piezometer 7777/4, located in Draa Affane (near El Houareb dam), observed head level, has manifested such an overall decrease about 16 m from 1987 to 2017 shown by a negative slope (−0.6) with a very high correlation coefficient (0.98). Before 2010, near the El Houareb seepage, this piezometer had displayed a good agreement between observed data and simulated less than +/−2 m. However, within the last 07 years, the difference between simulated and observed data could reach +/− 5 m.

• Water balance in steady and transient state

The hydrogeological budgets computed in steady and transient states are illustrated in the table below (

Table 2. Water balance was simulated by the model in 1969 and 2017 (m³/s).

Year	1969		1990		2017
Parameter	Inflow	Outflow	Inflow	Outflow	Inflow	Outflow
Inflow from Jbel Cherichira and Jbel Cherachine (Figure 4)	0.20		0.17		0.11
Drains		1.42		1.39		0.23
Rivers	1.56		0.8		0.42
Wells (pumping total)	-	0.99		1.59		3.49
Recharge	0.65		0.59		0.23
Groundwater Storage			1.42		2.96
Total	2.41	2.41	2.98	2.98	3.72	3.72

). Total inflow was approximately 2.4 m³ in 1969, shared between recharge from wadis beds in January 1969, inflow from Jbel Cherichira and Cherachine, and recharge from other boundaries. This finding is consistent with the total calculated inflow (2.1 m³/s) by [7]. The total estimated recharge of the Mio-Plio-Quaternary aquifer of Kairouan plain in 1969 was about 2.21 m³/s, where the contribution of the two rivers’ inflow represents 70.5%. This proves undeniably that the stream beds were the main source of groundwater recharge.

The inflow coming from rivers in 1990 has decreased by almost 50% from 1.56 in 1969 to 0.8 m³/s. According to several studies [6,9], this decrease is a result of the two dams’ impoundment where the infiltration from the streams’ beds was blocked. Hence, in 2017, this infiltration from the two rivers is limited to 0.42 m³/s.

• Sensitivity analysis

To ensure the reliability of the model results, sensitivity tests have been carried out. The effect of changes in the sink/source values (i.e., Recharge, Pumping, and Hydrodynamics parameters) on the RMS was assessed through a changing of these values by specific percentages (±50 of the initial values).

The results of the sensitivity analysis are shown in Figure 10. The groundwater system appears to be sensitive to Recharge and Groundwater abstraction more than to hydraulic conductivity. An increase of 50% of recharge has ranged the RMs between 8 m and 12.6 m. However, an increase has ranged the RMS between 12.33 m and 15 m.

3.3. Simulation of Groundwater Scenario Management

The validated model was carried out to simulate the impact of dams’ impoundment and intensive pumping on groundwater level fluctuations. The difference between simulated heads in 2017 and the simulated head in 1969 for the current model and the two scenarios were plotted in the maps below (Figure 11). For the current model, characterized by the influence of the two dams on the recharge pattern and based on the official data of pumping rates, the spatialization of the GWL (Ground Water Level) change between 2017 and 1969 has shown net dewatering ranging between −11 m in the plain upstream and −17 m in the N-East part of the plain.

For Scenario 1, in case of the two dams were not built and the pumping rates remained the same as the official withdrawals data taken from DGRE, the groundwater whiteness rises in the plain downstream with a maximum of 4 m; however, the upstream has shown a decline of GWL about a maximum 8 m, but less than the current situation.

The results of simulating the second scenario without dams and raising the groundwater withdraws with the same percentage that the two dams contributed to irrigation supply show the highest dewatering of the Mio-Plio-Quaternary aquifer ranging between −22 m to −8 m. Indeed, the highest decline in groundwater level is observed near the Draa Affane area, which exceeds 20 m. The lowest decline is observed in the downstream plain.

Even though the impact of the two dams’ erection on the groundwater decline was detected, justified by the raise the GWL level in case of their absence, the impact of the excessive pumping form, the Kairouan plain aquifer is more disastrous on the water drop.

4. Conclusions

The Mio-Plio-Quaternary aquifer of Kairouan plain, within a Mediterranean semi-arid area in Central Tunisia, faces a serious groundwater scarcity. Yet, this collapsed aquifer is a key resource for the regional water supply for irrigation and, to satisfy this supply, an astounding boom in agriculture groundwater use has been perceived during the last 40 years.

Moreover, in order to protect the Kairouan city and to supply water for irrigation, two dams have been erected. Consequently, the impoundment of these dams has deeply modified the process of aquifer recharge. To understand clearly the flow dynamics of the deep percolation and the impact of the two dams on groundwater level trends, a groundwater flow model was developed to simulate different scenarios of water management with and without dams.

The groundwater flow model is based on a 3D geological conceptual model, gathering all available data from 114 borehole logs, 08 petroleum wells, and 10 seismic lines.

After a multi-data analysis, it is still very difficult to separate the two aquifers, which may support the hypothesis of considering the Kairouan plain as a monolayer with clayey lenses against the hypothesis of separated layers. Besides, to simulate the groundwater head from 1969 to 2017, the Modflow code has been used thanks to its capacity to integrate several types of data. Gathering all piezometric data, the model was calibrated in a steady state with reference to the piezometric measurements measured in 1969. The calibration of the transient state was based on the period from 1970 until 2017. The validation was carried out between 2007 and 2017.

The model outputs provide interesting findings, foremost among which is that the main source of Kairouan plain recharge used to be the runoff inflows from the two wadis Marguellil and Zeroud and, since the impoundment of the two dams, this inflow has decreased by 70%. Hence, this reduction may explain the claim of many studies [6,7] about the dams’ impact on groundwater decline. The calibration of the transient state has highlighted the urgent case of the groundwater level (GWL) drop reaching 17 m in a piezometer located in the Zeroud area. The computed scenario management without dams and with the same official withdrawals has shown a groundwater rise downstream of the plain; yet, the decline of GWL persists but is less intense, compared to the case with dams. However, in the case of the absence of dams and with a raise in pumping rates, the GWL decline reaches its maximum of 22 m. Even if the impact of the dams’ erection on GWL decline is detected, the map of GWL changes shows obviously the impact of groundwater exploitation on the GWL of the Mio-Plio-Quaternary aquifer is the key factor to this water table drop.

Even if the outcomes of the model calibration seem to be very decent and quite consistent between observed and simulated heads, there is no guarantee that this is representative of the whole aquifer system due to the heterogeneity of the permeability and recharge process variability over space and time, especially for Plio-Quaternary and Miocene aquifers (e.g., [22]). This strong temporal and spatial variation of the recharge is an inseparable feature of semi-arid regions.

Relying on these findings integrating several data about the Kairouan plain aquifer and its flow dynamics, this study could help managers to predict the faith of groundwater resources within the next decades and to find out innovative strategies in order to balance the water supply and its protection for the next generation.