Environ Earth Sci (2017)76:699 DOI 10.1007/s12665-017-7021-y

ORIGINAL ARTICLE

Effect of the spatial variability of transmissivity on the groundwater flow and budget of the Takelsa multilayer aquifer, Tunisia Nesrine Ghouili1,2 • Luı´s Ribeiro3 • Mounira Zammouri1 • Faten Jarraya Horriche2

Received: 19 February 2017 / Accepted: 4 October 2017  Springer-Verlag GmbH Germany 2017

Abstract Groundwater models are vulnerable to uncertainties, especially those that are developed based on incomplete knowledge of the hydraulic parameters. Transmissivity is often the most uncertain parameter in groundwater modeling. This work evaluates the effect of the spatial variability of transmissivity on the outputs of the groundwater flow model of the Takelsa multilayer aquifer. The study was based on a combination of deterministic and stochastic modeling. Groundwater flow modeling and stochastic simulations were made using the PMWIN software (a version of MODFLOW). Initially, the groundwater flow model was developed in the steady state. The criteria used for the calibration process were based on the comparison between the observed and the simulated hydraulic level. Then, a stochastic approach was used and 100 simulations of transmissivity were performed. The simulated transmissivity values were used to recalculate the hydraulic

& Luı´s Ribeiro [email protected] Nesrine Ghouili [email protected] Mounira Zammouri [email protected] Faten Jarraya Horriche [email protected] 1

Faculty of Sciences of Tunis, LR01ES06 Mineral Resources and Environment Laboratory, University of Tunis El Manar, 2092 Tunis, Tunisia

2

Center of Water Research and Technologies, Georessources Laboratory, Technopark Borj Cedria, BP 273, Soliman, Tunisia

3

CERIS, Instituto Superior Te´cnico, University of Lisbon, Av. Rovisco Pais, 1040-001 Lisbon, Portugal

head and the water budget of the Saouaf deep aquifer. Finally, the consequent 100 results obtained by running MODFLOW were analyzed. The results showed variations in the hydraulic head levels, but the flow direction remained the same. The impact of the spatial variability of transmissivity on the water budget shows a range of inflow and outflow rates of Saouaf aquifer, with special relevance to the sea intrusion which may cause a probable deterioration of the quality of the groundwater. Keywords Hydraulic head  Groundwater flow modeling  Stochastic simulation  PMWIN  Seawater intrusion

Introduction Tunisia is influenced by the Mediterranean climate in the North and the Saharan climate in the South. The average rainfall varies from less than 100 mm in the extreme south to more than 1500 mm in the extreme north of the country. Such precipitation regime is naturally associated with considerable spatial variability in the water availability and hence limits the renewable water resources in some regions. Nevertheless, Tunisia has so far satisfied its water needs in the different water consuming sectors. In fact, water resource management—which is based on the rational use of the groundwater resources—is an important tool to control water consumption in the different sectors especially in agriculture—a main water consumer. Evaluating and protecting the available groundwater resources have become a growing necessity (Ben Alaya et al. 2003; Kemper 2004; Raul and Panda 2013; Ebrahimi et al. 2016), particularly because Tunisia, like most of the Mediterranean countries, is influenced by the impacts of

123

699

Page 2 of 18

the changing climate. The predicted impacts include a decrease in precipitation, hence leading to a decline in groundwater recharge and simultaneously to an intense exploitation of water resources. This situation will have negative impacts on water availability and the groundwater management will be more difficult (Nouiri et al. 2015). Mathematical modeling is a useful tool to study and improve groundwater management (Chenini and Ben Mammou 2010). Since the first works of Daniel (1967) on the Djeffara aquifer, mathematical modeling has been used for the quantitative management (Besbes 1976; Ennabli 1980), and mostly to study and evaluate the different problems of Tunisian aquifers such as temporal evolution of groundwater flow (Zammouri et al. 2007) exchange between aquifers, over-pumping, groundwater recharge (Chenini and Ben Mammou 2010), pollution (Kanzari et al. 2012; Zammouri et al. 2013), marine intrusion (Paniconi et al. 2001; Kerrou et al. 2010; Zghibi et al. 2011), and also as decision support element (Haddad et al. 2013; Nouiri et al. 2015). The Northeastern Tunisia (Cap Bon) is one of the regions that contain the most important aquifers, exploited especially for agricultural needs. However, some aquifers including the Takelsa multilayer aquifer remain poorly studied due to the structural complexity of the region. For a better understanding of the functioning of this aquifer, a mathematical groundwater flow model was implemented in this study. The PMWIN (Chaing and Kinzelbach 1998) is one of the worldwide used software for simulating groundwater flow and has been chosen for this study (El Yaouti et al. 2008; Wang et al. 2008; Xu et al. 2011; Zghibi et al. 2011; Park et al. 2014; Haddad et al. 2013; Nouiri et al. 2015; Bouaamlat et al. 2016). This software can provide deterministic predictions while providing also information on prediction uncertainty (stochastic modeling). The construction of groundwater flow models depends on hydraulic parameters such as transmissivity, hydraulic conductivity and storage coefficient. Although these data can be obtained by analyzing pumping test results, their usefulness is limited as they may not reflect the aquifers heterogeneity (Zhang 2001; Pool et al. 2015). The problem of data heterogeneity can be taken into account using a stochastic approach (Li et al. 2009; Bau 2012; Kerrou et al. 2013; Sidiropoulos et al. 2015) as shown by many previous studies (Delhomme 1979; Lin et al. 2000; Liu et al. 2002; Nakaya et al. 2002; Jang and Liu 2004; Zammouri and Ribeiro 2017). The stochastic modeling allows an analysis of the spatial variability of data by using a geostatistical approach (Delhomme 1979). The aim of this research is to study the influence of the transmissivity spatial variability on the groundwater flow model of the Takelsa multilayer aquifer. Transmissivity

123

Environ Earth Sci (2017)76:699

was chosen because it is considered as an important hydrogeological parameter in groundwater flow modeling (Lin et al. 2000), and essential key input is required for the model’s calibration (Ahmed and Marsily 1987). The modeling is supported mainly by groundwater flow directions, and the budget is analyzed with a combination of numerical and stochastic modeling. The model developed in this study is expected to develop our currently very limited knowledge of the studied aquifer system, and hence contribute to a more informed groundwater management of the aquifer region.

Takelsa multilayer aquifer system Location and climate The Takelsa multilayer aquifer system covers an area of about 300 km2. It is located between latitudes 36400 000 and 36550 000 N and longitudes 10300 000 and 10500 000 E. This aquifer is situated in a region of a typical semiarid Mediterranean climate with very irregular annual rainfall, on average around 518 mm/year. The region is characterized by a rainy winter, with the months of December and January representing the rainiest period with an average value of 81 mm/month and a mean temperature of 12 C. The summer has a mean temperature of 25 C especially in July and August. Maximum values of evapotranspiration are recorded during the summer season (July and August), and the mean annual potential evapotranspiration reaches 882 mm. As a result of its Mediterranean climate, the Takelsa basin is characterized by a very developed hydrographic network especially down in the Mountains. There are Wadi Al Abid in the northeast, Wadi Bezikh in the southwest, and some coastal wadis in the north. For a better management of the surface water resources, two dams have been built at the level of Wadi Al Abid and Wadi Bezikh (Fig. 1). Geology The Takelsa multilayer aquifer extends from the Korbous Mountain to the Western flank of Abdurrahman Mountain, passing through the Takelsa Syncline (Fig. 1). The Takelsa Syncline is a large depression formed mainly by sandstone and clay alternations, locally called by Saouaf formation. The axis of the syncline is oriented NE– SW and changes direction toward the east as it approaches the north and continues into the Mediterranean Sea. The eastern flank is affected by faults (direction N30), indicating an asymmetrical structure of the syncline (Ben Salem 1992). The geological map shows variable dips for

Environ Earth Sci (2017)76:699

Page 3 of 18

699

Fig. 1 a Geographic location, b geological map of the study area, and c cross section ‘‘C1’’ through Takelsa multilayer aquifer

both sides and the central part of the Takelsa Syncline with 30–35 for the western flank, 20–25 for the eastern flank, and horizontal dips for the central part. The Korbous Mountain is a structure that dominates the southeast part of the Gulf of Tunis with its western flank collapsed under the sea and a coastline marked by hot and cold springs. These springs gush out from the sandstone and limestone levels of Fortuna formation (Fig. 1). They are marked by low dips around 20E, but they are sharply accentuated on the eastern part. This structure is affected by several faults (Fig. 1), the main ones being: F1: a great collapse fault, substantially straight with a main direction N20, 70NNW, and F2: a large fault (direction E–W) which separates the mountain into two compartments. There are other syn-sedimentary normal faults with direction N50 and N90 which reflect the distensional tectonic movements during the period of deposits of Fortuna formation. The Eastern flank of this structure is truncated by reverse faults that separate the Korbous Mountain from the Takelsa Syncline. The Abdurrahman Mountain occupies the central part of the Cap Bon (Ben Salem 1992). It corresponds to a slightly asymmetric anticline, with SW–NE direction and with an elliptical shape. It reproduces the same axial disharmony as that of the Takelsa Syncline with additional complications. It is affected by several faults with different directions, which are often marked in the limestone and the sandstone levels, respectively, of Aıˆn Ghrab and Beglia formations of the Miocene.

Hydrogeology The Takelsa multilayer aquifer is mainly formed by sandstone and clay series ranging from the Oligocene to the Quaternary (Fig. 1). Based on the geological ages, four aquifers are distinguished from the top to the bottom: the phreatic aquifer contained in the Quaternary deposits and in the shallow sandstone level of Saouaf formation, two deep aquifers contained in the deep sandstone of the Miocene (Saouef and Beglia), and the Oligocene aquifer (Fortuna). The Oligocene aquifer outcrops in both Korbous and Abdurrahman Mountains. It is formed mainly by sandstone levels with intercalation of marl and white limestone. This aquifer is exploited by about 12 boreholes during the period between 1986 and 2015, but only 7 boreholes were operating in 2015 (DGRE 2015). The existing boreholes are located down Korbous and Abdurrahman Mountain, near the Oligocene outcrops and their depths range between 150 and 395 m. Captured levels which are mainly the sandstone levels are located at a depth of about 200 m with an average thickness of 100 m. In turn, the Miocene aquifer is characterized by sandstone deposits, clay deposits, and alternate levels of clay and sandstone. This reservoir is constituted by two formations locally called by Saouaf and Beglia. The Saouaf formation (Serravalien–Tortonian) is mainly clay formation, with alternating sandstone layers representing the captured aquifers. Almost all boreholes operating in

123

699

Page 4 of 18

Takelsa multilayer aquifer are implanted in this formation. The captured depth varies between 50 and 450 m. The Beglia formation (Serravalien) is mainly formed by sandstone levels. However, this aquifer is less exploited than the Saouaf aquifer. The depth of boreholes capturing this aquifer varies between 68 and 381 m. The sandstone levels were exploited by 208 boreholes in 2015 (DGRE 2015). Finally, the Quaternary reservoir is formed by sand deposits coming from the degradation of sandstone levels of the area. The number of boreholes capturing this reservoir has decreased from 8 in 1986 to 3 boreholes in 2015 (DGRE 2015). It should also be noted that the sandstone levels of the Miocene and the Quaternary reservoir were exploited also by about 900 shallow boreholes (DGRE 2015) mostly from the Miocene levels. According to the Tunisian water resources agency (DGRE), the Takelsa multilayer aquifer has been exploited since the 1980s, and the groundwater is mainly used for irrigation needs with 9.39 9 106 m3 in the year 2015, while 0.14 9 106 and 0.62 9 106 m3/year are, respectively, used for industrial and drinking water supply. The exploitation of the phreatic aquifer by shallow boreholes has increased from 1.61 9 106 m3/year in 1980–12.24 9 106 m3/year in 2015 (DGRE 1980–2015), and the groundwater resources are estimated to be 12 9 106 m3/year in 2015, hence with a registered deficit of about 0.24 9 106 m3/year. For the deep aquifers (Miocene and Oligocene), the exploitation has increased from 1.08 9 106m3/year in 1981 to 9.97 9 106 m3/year in 2015. Since the groundwater resources were estimated to be 3.5 9 106m3/year in 2015, a deficit of 6.47 9 106 m3/year was registered in 2015. The sandstone levels of the Miocene are the most exploited with almost 90% of the total exploitation (Fig. 2). The exploitation of the Miocene levels was nearly stable in the early 80 s, since then it

Fig. 2 Evolution of the exploitation for Takelsa multilayer aquifer

123

Environ Earth Sci (2017)76:699

increased to 9.34 9 106 m3/year in 2015. This exploitation rise is due to the increase in the number of wells, reaching 208 boreholes in 2015. As a conclusion, the water balance of the Takelsa multilayer aquifer registered a deficit of 6.71 9 106 m3/year in 2015. The interpretation of the results of pumping tests carried out in the available 82 boreholes located in the Miocene aquifer (Saoauf and Beglia formations) provided the transmissivity (T) values. The Oligocene and the Quaternary aquifer were excluded from this study due to the lack of available data. Table 1 shows the basic statistics of transmissivity in Saouaf and Beglia aquifers. By comparing the values of mean and median and the skew coefficient, we concluded that transmissivity is lognormal distributed. The piezometric monitoring network of the Takelsa multilayer aquifer is limited because it is only composed by 9 piezometers in the deep aquifer and by 2 piezometers and 8 shallow boreholes in the phreatic aquifer (Fig. 3). The piezometers and the shallow boreholes used for the monitoring of the phreatic aquifer are located mainly in the Quaternary deposit in the bottom of Korbous Mountain and the Miocene deposits in the center part of the Takelsa syncline. Regarding the deep aquifers, the monitoring network is based by piezometers located mainly in the Miocene deposit in the center of the Takelsa syncline, the eastern flank of Korbous Mountain, and the western flank of Abdurrahman Mountain. The measurements with a biannual frequency started from 1972 with shallow boreholes and from 1997 with piezometers (DGRE 2007). Analyzing Figs. 4 and 5, we can infer the following: For the boreholes implanted in the quaternary deposit we notice that in ‘‘Bor5’’—monitors the aquifer from 1972 onwards—a decrease in the piezometric level till 1998, then a rising tendency (Fig. 4a). This can be related to a decrease in the groundwater exploitation. Regarding

Environ Earth Sci (2017)76:699

Page 5 of 18

699

Table 1 Transmissivity values of the deep Miocene aquifers (m2/s) n

Min

Max

Mean

Median

1st Q

1.6 9 10-2

1.71 9 10-3

6.35 9 10-4

2.83 9 10-4

4.08 9 10-3

1.3 9 10-3

7.10 9 10-4

3.23 9 10-4

1rd Q

SD

Skew coef.

1.6 9 10-3

3.05 9 10-3

3.23

2.03 9 10-3

1.32 9 10-3

1.08

Saouaf aquifer T

62

2.2 9 10-5

Beglia aquifer T

14

6.5 9 10-5

Fig. 3 Piezometric monitoring network of the Takelsa multilayer aquifer

‘‘Bor8’’, located in the boundary limit of Takelsa aquifer near to the water course in a generally unexploited area (Fig. 3), rising tendency is observed, potentially indicating that this area can be a recharge area of the groundwater of the phreatic aquifer (Fig. 4a). The ‘‘Bor 6’’ shows a decrease in the piezometric level that can be the result of a localized exploitation, albeit this conclusion cannot be confirmed due the fact that the monitoring has stopped in this well. For the phreatic aquifer included in the shallow level of Saouaf formation, the monitoring boreholes (Bor2 and 3 and Piezometer 7) indicate a stability of the hydraulic head

level; however, the observed fluctuations can be related to the rainfall variations or to a localized exploitation (Fig. 4b). Regarding the deep aquifer of Saouaf, the piezometers Pz 4 and Pz 6 that monitor these levels and located mainly on the eastern flank of Takelsa syncline reveal a gradual increase from 10 to 20 m. This can be the result of a decompression phenomenon of the aquifer or possible upcoming water from deeper aquifers. However, a general stability of the piezometric level was recorded for the piezometer Pz 5 located on the western flank of the syncline. The observed oscillation is supposed to be the result of the seasonal fluctuations

123

699

Page 6 of 18

Environ Earth Sci (2017)76:699

Fig. 4 Piezometric level for the phreatic aquifer: a Quaternary level, b Saouaf level

(Fig. 5a). Regarding the Pz 8, located almost near to the wadis Bezikh in an area excessively exploited, the piezometric level shows a fluctuation of 20 m. As for the Beglia deep aquifer, the only piezometer (Pz1) monitoring these levels shows a fluctuation between a decrease, an increase, and then stability of the piezometric level; this varied behavior can be the result of seasonal fluctuations (Fig. 5b).

Models and methods Groundwater Flow Model The methodological approach in this work was to firstly create the groundwater flow model and then generate a stochastic simulation based on a geostatistic approach. The PMWIN package (Chaing and Kinzelbach 1998) was chosen for the modeling of the groundwater flow of Takelsa multilayer aquifer because it is well documented and largely tested. It is widely used due to the fact that it is easy to access, less expensive, and capable of simulating large areas. PMWIN use MODFLOW for simulating the flow and MT3D for simulating the transport.

123

The MODFLOW—a modular three-dimensional finite difference groundwater model—was developed by the U.S. Geological Survey in 1989 for the description and the prediction of the groundwater systems behavior (Harbaugh et al. 2000). It resolves by using the finite difference method the groundwater flow equation in the steady and the transient states for multilayer aquifers, which is expressed as follows:     o ohi o ohi Tx;i Ty;i þ þ biþ1 ðhiþ1  hi Þ ox oy ox oy ð1Þ ohi þ qi þ bi1 ðhi1  hi Þ ¼ Si ot where hi, hi?1, hi-1 are hydraulic heads [L] of layers i, i ? 1, and i - 1, respectively; Tx,i and Tx,i is the transmissivity [LT-1] according to x and y of layer i; bi?1 and bi-1 are leakance [T-1] between layers i, i ? 1 and between i and i - 1, respectively; Si is the storage coefficient [L] of layer i; qi is source/sink of water [LT-1] in layer i and t [T] is the time. MODFLOW has been increasingly applied to the description and prediction of the behavior of groundwater systems over the last few years. The ‘‘original’’ version of MODFLOW-88 (McDonald and Harbaugh 1988) or MODFLOW-96 (Harbaugh and McDonald 1996a, b) can

Environ Earth Sci (2017)76:699

Page 7 of 18

699

Fig. 5 Piezometric level for the deep aquifers: a Saouaf aquifer, b Beglia aquifer

simulate the effects of wells, rivers, drains, head-dependent boundaries, recharge, and evapotranspiration. Since the publication of MODFLOW, various codes have been developed by numerous investigators. These codes are called packages, models, or sometimes simply programs. Packages are integrated with MODFLOW, and each package deals with a specific feature of the hydrologic system to be simulated, such as wells, recharge or river. Models or programs can be stand-alone codes or can be integrated with MODFLOW. A stand-alone model or program communicates with MODFLOW through data files. Both codes use MODFLOW as a function for calculating flow fields. The MODFLOW model consists of the main program and various packages (McDonald and Harbaugh 1988). Those used in this study include drainage, evapotranspiration, recharge, and well.

Estimation of recharge The groundwater recharge is one of the important and necessary inputs for the groundwater flow modeling. The WetSpass model (Water and Energy Transfer between Soil, Plants and Atmosphere)—a hydrological physically based model—was used for the quantification of the groundwater recharge. It is integrated into the software ArcView GIS as a raster model (Batelaan and De Smedt 2001). It is one of the widely used models for recharge assessment (AbuSaleem et al. 2010; Al Kuisi and El-Naqa 2013; Park et al. 2014; Tammal et al. 2014; Ghouili et al. 2017). WetSpass allows the calculation of the groundwater recharge, evapotranspiration, and surface runoff by solving first the water balance equation cell by cell and successively, the vegetated area, bare soil, open water, and impervious surfaces. Secondly, the total water balance of a given area is calculated as the sum of the water balance of

123

699

Page 8 of 18

Environ Earth Sci (2017)76:699

each raster cell. Finally, the water balance components of each area are used to calculate the total water balance. The groundwater recharge is calculated in WetSpass in the last step. It is considered to be the result of the water balance: R ¼ P  S  ET  I

ð2Þ

where R is the groundwater recharge, P the average seasonal precipitation, S the surface runoff, ET the actual evapotranspiration, and I the interception fraction; all with the unit [LT-1]. The model inputs include the land use map, the soil texture map, the digital elevation model, the slope map, the groundwater depth, and the climatic data (precipitation, potential evapotranspiration, temperature and wind speed for winter and summer season). Stochastic simulation Stochastic simulations can be produced using a geostatistical approach. Geostatistics has become an autonomous discipline of the stochastic subsurface hydrology field. It is well beyond the scope of this paper to give an account of the geostatistical theory (more information in Matheron 1970). The praxis in geostatistics encompasses in general two interrelated steps: the structural analysis and kriging or conditional simulation. In the first phase, the spatial structure of the hydrogeological variables is described by a variogram, defined as the expected squared difference between pairs of data values. The function characterizes suitably the continuity of the variable and the cross-correlation between variables. In particular, the range of influence of a sample, i.e., the distance beyond which the variable is uncorrelated, the magnitude of the measurement errors, and the small-scale variability are estimated. Also the behavior of the variogram in specific azimuths highlights the anisotropies if present and evaluates multi-scale effects (nested structures) and trends. For estimation purposes, the experimental functions are fitted with positive definite theoretical models and cross-validated with experimental data. Once the spatial variability is modeled, the kriging algorithm will estimate the unknown value of the parameter in non-sampled locations by a linear combination of the known values (being the distances of the structural type). This phase has been currently related to spatial interpolation operations. The kriging operator has the property of exactitude, meaning that the estimation of the parameter at a location where it has been measured is identical to the observed value. The method always provides a variance of the estimation error which is zero in measured points.

123

Conditional simulation (CS) methodology is an important Monte Carlo type technique because it takes to account for spatial correlation between parameter data and the simulated variables are consistent with the measured ones at well sites. The visualization of several simulations gives an estimate of the uncertainty of the phenomenon. For this reason, each CS is a version of the unknown reality. These realizations are built by stochastic simulation algorithms and are equally probable high-resolution models of the spatial distribution. Each simulation is conditioned in the sense that the resulting realizations honor the hard data values at their locations. Also, the simulated values of the variable have the same statistical distribution and correlation structure as the field measurement. These are good representations of reality, and they measure the spatial uncertainty. Simulation is preferable to kriging in most applications because the latter are designed to minimize a certain local error impact and not designed to reproduce patterns of spatial continuity so as the case of the former method. Furthermore, the alternative images of a particular parameter produced by CS can be processed through a groundwater flow simulator to yield the corresponding measure of uncertainty on the response function, providing answers to such questions as what is the probability of some output variable exceeds a specific threshold value (Ribeiro 1993). The field generator code of PMWIN allows the user to carry out stochastic modeling. For this purpose, the mean value l, standard deviation r, and correlation scales in both I and J directions were used to generate a quantitative description of the transmissivity field. As it is mentioned by Chaing and Kinzelbach (1998), the conditional generation of a single realization proceeds in five steps: •

•

•

• •

The parameter value for each model cell is interpolated from the measurements using the kriging method. The correlation length is determined from the measurements. An unconditional generation is performed using the field generator with the same correlation length (correlation scale). The unconditionally generated values at the measurement locations are used to interpolate values for each model cell by using the kriging method again. The distribution from step 3 is subtracted from the distribution from step 2 yielding kriging residuals. The kriging residuals are added to the distribution from step 1 yielding a realization which has the same correlation length and passes through the measured values at the measurement points.

Environ Earth Sci (2017)76:699

Results and discussion Flow directions The first step in developing a groundwater flow model is the construction of the conceptual model. This step is based on the selection of fundamental assumptions on which the modeling process will be based. The Takelsa multilayer aquifer model is formed of three main layers: the first layer corresponding to the phreatic aquifer—it is formed by the quaternary and the shallow Miocene formations (Saouaf)—and the two other layers correspond to the deep aquifers of Saouaf and Beglia (Fig. 6). The first and second layers are limited by the outcrops of Saouf formation, while the third layer is limited by Beglia outcrops. The Mediterranean Sea borders the aquifers in the north and in the southwest. In the south, the limit is adjusted to the stream of Wadi Bezikh which separates the Takelsa multilayer aquifer from the Grombalia aquifer (Fig. 7). The aquifers’ recharge is provided mainly by infiltration of rainwater through the quaternary deposit and the sandstone outcrops of Beglia and Saouaf. The natural discharge areas are principally the Mediterranean Sea (in the North and the Southwest), corresponding to the down course of wadis draining the superficial aquifers and the evapotranspiration area. During the modeling process, the deep aquifers were considered confined with user-specified transmissivities

Page 9 of 18

699

and the phreatic aquifer as unconfined aquifer. The possibility of vertical exchange between the three layers through semi-confining layers was taken into account. A finite difference grid with 500-m square cells was superposed to the study area. Three types of boundary conditions were used: (1a) Dirichlet boundary condition was imposed in the north and the southwest where a constant head value of 0 m is fixed for the sea level; (2a) Cauchy boundary condition was used to represent the drainage by wadis and; (3a) Neumann boundary condition was applied to represent the groundwater recharge and the pumping by wells (Fig. 7). In turn: •

The groundwater recharge map was obtained by WetSpass and applied to the top layer outcrops (phreatic and the Beglia aquifers).

The land use map was obtained from the local water resources authority (CRDA-Nabeul) database. It was prepared by the union of several maps. The derived land use map reveals five categories of land cover: agricultural land, orchards, shrub, forest, and buildup. However, the agricultural land and the orchards are the dominated land cover of the area (Fig. 8a). Regarding the soil map, it was deducted using four geological maps (1:50 000) of Menzel Bouzelfa; Tazoghrane; La Marsa; and La Goulette. Ten classes of the soil texture were identified with the dominance of sandy clay, sandy loam, and loamy sand (Fig. 8b).

Fig. 6 Conceptual model of Takelsa multilayer aquifer

123

699

Page 10 of 18

Environ Earth Sci (2017)76:699

Fig. 7 Groundwater flow model limits and boundary conditions

Each of the land use and the soil texture grids is connected to their attribute tables containing the related parameters. The digital elevation model is frequently a crucial source of information for groundwater recharge assessment. It is required to define the topography of the study area. For the Takelsa multilayer aquifer, the digital elevation model (Fig. 8c) was obtained from the Shuttle Radar Topography Mission dataset of the National Aeronautics and Space Administration (NASA) (Jarvis et al. 2008). The cell resolution was of 100 m for computational efficiency. The slope map was also directly derived from the digital elevation model. The values range from 0 to 44% with a mean value of 6% indicating that most of the study area is relatively flat (Fig. 8d). The initial groundwater depth map was produced using the hydraulic heads map obtained from MODFLOW during the steady state and the digital elevation model. It should be noted that we have used a regional recharge map as input for MODFLOW and it was not a result of WetSpass. The climatic data (precipitation, temperature and wind speed) were obtained from the National Institute of

123

Meteorology (INM) database. The seasonal precipitations were computed from a daily precipitation data measured for the period 1980–1984. Temperature and wind speed were computed using monthly measured values during the same period. For the potential evapotranspiration, it was simulated using the monthly temperate data using the Thornthwaite formula (Thornthwaite 1948). All of the grid data have a cell size 100 m, and they were clipped to give the same boundary as the model (MODFLOW). The simulated groundwater recharge map (Fig. 9) shows values ranging from 0 to 167 mm/year with a mean value of 22 mm/year. The groundwater recharge represents an average of 4% of the precipitation with the highest values observed in areas characterized by sandy loam and loamy sand soil, low slope and the presence of vegetation cover. However, in the forest and shrub area, a low amount of groundwater recharge has been obtained independently of soil types (Ghouili et al. 2017). The infiltration coefficient varying between 0 and 32% is reasonable. Strong values correspond to the recharge rate in wadis generally

Environ Earth Sci (2017)76:699

Page 11 of 18

699

Fig. 8 WetSpass inputs (Ghouili et al. 2017)

characterized by high values that could reach 40% of rainfall in Tunisia. •

•

•

The groundwater exploitation in the steady state was through shallow boreholes capturing the phreatic aquifer (in the center of Takelsa Syncline) and Beglia aquifer (outcrops near the Korbous and Abdurrahman mountains). The phreatic aquifer was exploited by 440 shallow boreholes with a total exploitation of 7.2 9 106 m3/year; however, Beglia aquifer was exploited by 59 shallow boreholes with a total exploitation of 0.2 9 106 m3/year. The transmissivity values introduced for the model were deducted from pumping tests. Average values were used for each aquifer due to the scarcity and the high spatial variability of the measurements. An initial average value equal to 1 9 10-3 m2/s was introduced for the phreatic and Saouaf aquifers and 4 9 10-4 m2/s for Beglia aquifer. The vertical leakage which defines the degree of vertical hydraulic connection between aquifers was deduced by the model’s calibration. We introduced the values of 10-8 s-1 between the phreatic and Saouaf

aquifers and 10-11 s-1 between Saouaf and Beglia aquifers. The reference period for the steady state calibration extends from 1980 to 1984, during which the aquifers were considered to be in an equilibrium state. The trial-and-error method was used in the calibration process; it consists in changing transmissivity until getting hydraulic head values become similar to the values observed at the different measurement points. During the model’s calibration, the transmissivity of the phreatic and Saouaf aquifers was reduced to 5 9 10-4 and 4 9 10-4 m2/s, respectively. The calibration results indicate a good calibration with a correlation coefficient of 0.9 (Fig. 10b). For the phreatic and Saouaf aquifer, the observed hydraulic head varies between 0 and 85 m, and high values are observed in the eastern part. The main flow direction is from the east coming from Abdurrahman Mountain (Fig. 10a, c). For the Beglia aquifer (Fig. 10d), the hydraulic head values are between 0 and 95 m. The flow direction is from Abdurrahman and Korbous Mountain. In general and for the three aquifers, the flow accumulates in the center part of Takelsa synclinal, and then it

123

699

Page 12 of 18

Environ Earth Sci (2017)76:699

Fig. 9 Simulated groundwater recharge rates (Ghouili et al. 2017)

is divided into a part joining the discharges area (the Mediterranean Sea) in the north and the southwest, while the other part flowing to the southeast toward Grombalia aquifer. For Beglia aquifer the drainage by wadis in the north is clearly shown by the calculated piezometric map. Water budget The total water budget is estimated to be 15.3 9 106m3/ year (Table 2). The main outputs of water come from the pumping (48%) with annual rates of 7.4 9 106 m3/year, the evapotranspiration (26%) with 4.1 9 106 m3/year, the drainage by the wadis (18%) with 2.7 9 106 m3/year, and the leakage to the Mediterranean Sea (7%) with 1.1 9 106 m3/year. Another output of water is through the vertical exchange between the different aquifers. The main input of water for the whole system is direct recharge.

123

The total water budget for the phreatic aquifer is estimated to be 17.7 9 106 m3/year, with 73% of water that joins the aquifer comings from the recharge and 28% coming from Saouaf aquifer, while 41% of water is extracted in wells, 30% goes to the Saouaf aquifer, and 23, 5, and 2% are lost, respectively, by evapotranspiration, drainage by wadis and the leakage to the Mediterranean Sea. The total water budget for the Saouaf aquifer in the steady state is about 5.6 9 106 m3/year with the main input of water being from the exchange with the phreatic (95%) and Beglia aquifer (5%). The output of water occurs through the exchange with two other aquifers and the leakage to the Mediterranean Sea (6%). Finally, the main input of water for the Beglia aquifer is from recharge (86%) and exchange with the Saouaf aquifer (13%), while the main output is, respectively, by the drainage by wadis (65%), the leakage to the Mediterranean Sea (18%), with Saouaf aquifer (9%), and pumping (7%).

Environ Earth Sci (2017)76:699

Page 13 of 18

699

Fig. 10 Simulated hydraulic head for: a the phreatic aquifer, c Saouaf aquifer and d Beglia aquifer. b Scatter plot of observed and calculated hydraulic head for the phreatic aquifer. The black

points show the location of piezometric observation boreholes used in the model calibration process

The total water balance estimated by the model during the steady state is lower than the renewable resources of the system evaluated by the Tunisian authorities, estimated to be 20 9 106 m3/year. This difference can be due to the global approach used for the calculation. After the construction of the model and the validation of the results, a geostatistical approach was used to examine the coherence of the transmissivity values used for the calibration of the model in the steady state. All the boreholes exploiting the Takelsa multilayer aquifer are localized in the central part of Takelsa syncline and capture sandstone and clay–sandstone levels of the Saouaf formation. These levels

are exploited by about 208 boreholes (DGRE 2015), but only 62 include transmissivity data. Impact of T spatial variability on groundwater flow and budget To determine the spatial variability of transmissivity, a variographic analysis was firstly performed with the 62 values. The lognormal distribution of transmissivity obliged to work with log transmissivity. The obtained experimental omnidirectional variogram displayed a transition phenomenon with null nugget effect. A spherical model was fitted with a range = 800 m and sill = 0.418

123

699

Page 14 of 18

Environ Earth Sci (2017)76:699

Table 2 Water budget of Takelsa multilayer aquifer in the steady state (9 106 m3/year) Input

0.8 0.7

Output

Phreatic aquifer

0.6

12.8

Pumping

7.2

Drainage by wadis

0.8

Down exchange

4.9

Losses by evaporation Leakage into Mediterranean sea Total

5.3 4.1

17.7

0.5

Variogram

Recharge

0.3

0.3 17.7

0.2

Saouef aquifer Upper exchange

5.3

4.9

Down exchange

0.3

0.4

5.6

5.6

Leakage to Mediterranean sea Total

0.1

0.3

0 0

2.5

Pumping

0.2

Drainage by wadis

1.9

Upper exchange

0.4

0.3

2.9

2.9

Leakage into Mediterranean sea Total

0.5

The multilayer aquifer Recharge Pumping

15.3 7.4

Drainage by wadis

2.7

Leakage into Mediterranean sea

1.1

Losses by evaporation Total

4.1 15.3

15.3

(Fig. 11). Then the ‘‘Field Generator’’ module was used to produce 100 likely transmissivity fields on a grid of 500 9 500 m using the following parameters: mean = - 3.34, standard deviation = 0.72, and correlation length = 800 m in both directions. Finally, MODFLOW was run with the 100 transmissivity fields and 100 hydraulic head maps, and water budgets were produced for the Saouaf deep aquifer. Figure 12 shows the results of five simulations selected among the 100 generated simulations. The comparison of the results of the conditional simulations revealed that the histograms and experimental variograms were similar. Moreover, the average and the standard deviation of every simulation were comparable to the average and to the standard deviation of the original 62 transmissivity values. The obtained hydraulic head varied between a minimum of - 1.04 m and a maximum of 80 m with a 49.5 m mean. The hydraulic heads maps showed that the spatial variation of the transmissivity caused some changes in the hydraulic

123

1000

2000

3000

4000

5000

6000

7000

Lag Distance

Beglia aquifer Recharge

0.4

Fig. 11 Experimental variogram of LOG10 (T)

heads of the Saouaf deep aquifer, but the general flow direction remained the same. The changes manifest mainly in the eastern part of the aquifer, with a fluctuation in the hydraulic head also in the southwest in the discharge area (limit with the sea). A lowering of the groundwater level has been noticed in the simulations 25 and 50, while in simulations 1 and 100 a flow of water from the sea was manifested (Fig. 12). More importantly, however, is to analyze the effect of the spatial variability of transmissivity on the water budget of Saouaf aquifer, especially in the quantity of water exchanged between this aquifer and the two others (phreatic and Beglia aquifer), and the risk of seawater intrusion from the Mediterranean Sea caused by leakage. In Fig. 13, five simulations among the 100 ones generated to represent the exchanged flow between the different aquifers are showed. For the phreatic aquifer, water percolates downward over all the area except in the northeast. In some simulations (e.g., simulation 25 and 50) the quantity of water exchanged between the two aquifers decreases. For Beglia aquifer, the water percolates upward to Saouaf aquifer principally in the south. Tables 3 and 4 represent a descriptive statistics of the volume of water exchange between the Saouaf aquifer with the Mediterranean Sea and between the Saouaf aquifer and the phreatic and the Beglia aquifers. We can infer the following: (a)

The quantity of water out flowing to the sea reaches a maximum of 0.3 9 106 m3/year (Table 3), while the inflow rate is in most cases null, but sometimes reached a small value estimated to 0.03 9 106 m3/ year. Although it is a small rate, it is a warning sign

Environ Earth Sci (2017)76:699

Page 15 of 18

699

Fig. 12 Simulated LOG10 (T) and hydraulic head maps, LOG10 (T) histograms, and experimental LOG10 (T) variograms

Fig. 13 Exchanged flow map between phreatic and Saouaf aquifer (a) and between Saouaf and Beglia aquifer (b) for different simulations: negative values represent inflow rate and positives values outflow rate (flow rate in 9 106 m3/year)

123

699

Page 16 of 18

Environ Earth Sci (2017)76:699

Table 3 Descriptive statistic for the exchange between Saouaf aquifer and the sea (9 106m3/year) Statistic

In

Out

Min

0.00

0.001

0.291 - 0.028

(b)

Out-in

Max

0.03

0.291

Median

0.0001

0.041

0.039

1st Q

0.00

0.017

0.085

3rd Q

0.002

0.085

0.017

Table 4 Descriptive statistic for the exchanged flow between Saouaf aquifer (2), phreatic aquifer (1) and Beglia aquifer (3) (9 106 m3/ year) Statistic

In(1-2)

Out(2-1)

In(3-2)

Out(2-3)

Min

5.81

5.69

0.213

0.298

Max

7.88

7.76

0.291

0.435

Median 1st Q

6.95 6.68

6.77 6.48

0.254 0.239

0.370 0.356

3rd Q

7.26

7.08

0.265

0.388

for the necessity of pumping control in this region. With continuous over-pumping, the aquifer will be contaminated by sea water intrusion.

(c)

(d)

The quantity of water exchanged between the Saouaf aquifer and the phreatic aquifer is more important than the volume exchanged with the Beglia aquifer: almost 75% (& 7.5 9 106 m3/year) of water that joins Saouaf aquifer comes from the phreatic aquifer and 75% (& 0.27 9 106 m3/year) of water is wasted toward the Beglia aquifer. The quantity of water entering to Saouaf aquifer from phreatic aquifer reaches a maximum of 7.88 9 106 m3/year (Table 4), with a median of 6.95 while the outflow rate reaches a maximum of 7.76 9 106 m3/year with a median of 6.77 9 106 m3/year. The quantity of water entering to Saouaf aquifer from Beglia reaches a maximum of 0.291 9 106 m3/ year (Table 4), with a median of 0.254 9 106 m3/ year while the outflow rate reaches a maximum of 0.435 9 106 m3/year with a median of 0.377 9 106 m3/year.

Probability distribution functions can now be calculated using the results of the net quantities (OUT-IN) of the water budget. Figure 14 displays the probability distribution functions of volumes of water interchanged between Saouaf aquifer and the phreatic aquifer, Beglia aquifer, and the

Fig. 14 Probability distribution functions of net interchange rates between Saouef aquifer and: a phreatic aquifer, b Beglia aquifer, and c the Mediterranean Sea

123

Environ Earth Sci (2017)76:699

Mediterranean Sea, respectively. With these functions we can evaluate, given a specific volume of water, the probability of this value to be surpassed. For example: (a) from Fig. 14a we conclude that 80% is the probability that the amount of water gaining by Saouaf from the phreatic aquifer is higher than 0.126 9 106 m3/year and; (b) from Fig. 14b we conclude that 80% is the probability that the amount of water gaining by Saouaf from the Beglia aquifer is higher than 0.105 9 106 m3/year. While analyzing Fig. 14c, we can conclude about the risk of seawater intrusion. The probability distribution functions show 20% probability that the volume of water surpasses 0.11 9 106 m3/year. These results suggest a need to reduce the pumping for this aquifer and consequently, leading to major groundwater management problems and hindering any kind of sustainable development in the region.

Conclusions In this work, stochastic conditional simulations were conducted to generate multiple replicates of aquifer properties and groundwater flow patterns. A ‘‘classical’’ deterministic model of groundwater flow (finite difference method) was used to simulate groundwater flow for a given spatial distribution of transmissivity. The combination of deterministic and stochastic hydrodynamic modeling allows a better characterization of the spatial heterogeneity associated with the uncertainty of groundwater flow parameters given the limited available dataset. As main result of this study, we show that the uncertainty of the results in terms of groundwater flow can be relatively important and that the probabilistic approach is a good means to estimate this uncertainty. A more precise management of the Saouaf aquifer resources is recommended. Also, it is necessary to optimize the exploitation of the aquifers in the region without exhausting their reserves. The monitoring network should become denser and cover all the area of the aquifer and the monitoring and the quality control of the coastal sections is also suggested. A hydrochemical investigation will be examined in a follow-up publication by coupling hydrochemical, isotopic, and multivariate statistical analysis to improve the groundwater quality of the aquifer and for the evaluation of its suitability for drinking and agriculture purposes. Acknowledgements The authors wish to thank the anonymous reviewers of this article for their useful comments and suggestions. Gratitude is also extended to the Resources Water Direction of Tunis (DGRE) and the Regional Direction of Agriculture and Water Resources of Nabeul (CRDA Nabeul).

Page 17 of 18

699

References Abu-Saleem A, Al-Zu’bi Y, Rimawi O, Al-Zu’bi J, Alouran N (2010) Estimation of water balance components in the Hasa Basin with GIS based WetSpass model. J Agron 9:119–125 Ahmed S, Marsily DG (1987) Comparison of geostatistical methods for estimating transmissivity using data on transmissivity and specific capacity. Water Resour Res 23:1717–1737. doi:10.1029/ WR023i009p01717 Al Kuisi M, El-Naqa A (2013) GIS based Spatial Groundwater Recharge estimation in the Jafr basin, Jordan –Application of WetSpass models for arid regions. Revista Mexicana de Ciencias Geolo´gicas 30:96–109 Batelaan O, De Smedt F (2001) WetSpass: a flexible, GIS based, distributed recharge methodology for regional groundwater modelling. In: Symposium proceedings ‘‘Impact of human activity on groundwater dynamics’’, the Sixth IAHS Scientific Assembly, Maastricht, The Netherlands, IAHS Publ.269, pp 11–17 Bau AD (2012) Planning of groundwater supply systems subject to uncertainty using stochastic flow reduced models and multiobjective evolutionary optimization. Water Resour Manage 26:2513–2536. doi:10.1007/s11269-012-0030-4 Ben Alaya A, Souissi A, Tarhouni J, Ncib K (2003) Optimization of Nebhana reservoir water allocation by stochastic dynamic programming. Water Resour Manage 17:259–272 Ben Salem H (1992) Contribution a` la connaissance de la ge´ologie du Cap Bon: stratigraphie, tectonique et se´dimentologie. The`se 3e`me cycle. Faculte´ des sciences de Tunis. Universite´ de Tunis II. Tunisia Besbes M (1976) E´tude hydroge´ologique du bassin de Moulare`sRedeyef sur mode`les mathe´matiques. Paris: Centre d’informatique ge´ologique-E´cole nationale des mines de Paris 6. France Bouaamlat I, Larabi A, Faouzi M (2016) Hydrogeological investigation of an oasis-system aquifer in arid southeastern Morocco by development of a groundwater flow model. Hydrogeol J. doi:10. 1007/s10040-016-1409-8 Chaing WH, Kinzelbach W (1998) Processing modflow: a simulation system for modeling groundwater flow and pollution, p 334 Chenini I, Ben Mammou A (2010) Groundwater recharge study in arid region: an approach using GIS techniques and numerical modeling. Comput Geosci 36:801–817 Daniel MJ (1967) Etudes des re´percussions d’une exploitation de dure´e sur un syste`me aquife`re saharien. Chronique d’hydroge´ologie 11. BRGM Delhomme JP (1979) Spatial variability and uncertainty in groundwater flow parameters: a geostatistical approach. Water Resour Res 15:269–280 DGRE (1980–2015) Situation de l’exploitation des nappes phre´atiques entre, 1980 et 2015. Direction Ge´ne´rale des Ressources en Eau. Tunisia DGRE (1981–2015) Annuaire d’exploitation des nappes profondes de 1981a` 2015. Direction Ge´ne´rale des Ressources en Eau. Tunisia DGRE (2007) Etude d’optimisation des re´seaux de suivi des ressources en eau, Re´seau de surveillance pie´zome´trique et de Qualite´ des eaux souterraines, re´gion Nord-Est.Direction Ge´ne´rale des Ressources en Eau. Tunisia Ebrahimi H, Ghazavi R, Karimi H (2016) Estimation of groundwater recharge from rainfall and irrigation in an arid environment using inverse modeling approach and RS. Water Resour Manage. doi:10.1007/s11269-016-1261-6 El Yaouti F, El Mandour A, Khattach D, Kaufmann O (2008) Modelling groundwater flow and advective contaminant transport in the Bou-Areg unconfined aquifer (NE Morocco). J Hydro Environ Res 2:192–209. doi:10.1016/j.jher.2008.08.003

123

699

Page 18 of 18

Ennabli M (1980) E´tude hydroge´ologique des aquife`res du Nord-Est de la Tunisie pour une gestion inte´gre´e des ressources en eau. The`se de Doctorat es-sciences naturelles. University of Nice. France Ghouili N, Jarraya Horriche F, Zammouri M, Benabdallah S, Farhat B (2017) Coupling WetSpass and MODFLOW for groundwater recharge assessment: case study of the Takelsa multilayer aquifer (Northeastern Tunisia). Geosci J 21:791–805. doi:10. 1007/s12303-016-0070-5 Haddad R, Nouiri I, Alshihabi O, Mabmann J, Huber M, Laghouane A, Yahiaoui H, Tarhouni J (2013) A decision support system to manage the groundwater of the Zeuss Koutine aquifer using the WEAP-MODFLOW framework. Water Resour Manage 27:1981–2000. doi:10.1007/s11269-013-0266-7 Harbaugh AW, McDonald MG (1996a) User’s documentation for MODFLOW-96 an update to the U.S. geological survey modular finite-difference ground-water flow model: U.S. Geological Survey Open-File Report 96-485, p 56 Harbaugh AW, McDonald M.G (1996b) Programmer’s documentation for MODFLOW-96, an update to the U.S. Geological Survey modular finite-difference ground-water flow model. U.S. Geological Survey Open-File Report 96-486, p 220 Harbaugh AW, Banta ER, Hill MC, McDonald MG (2000) Modflow2000 the U.S. geological survey modular ground-water model: user guide to modularization concepts and the groundwater flow Process. US Geological Survey. Reston p 121 Jang CS, Liu CW (2004) Geostatistical analysis and conditional simulation for estimating the spatial variability of hydraulic conductivity in the Choushui River alluvial fan, Taiwan. Hydrol Process 18:1333–1350 Jarvis A, Reuter HI, Nelson A, Guevara E (2008) Hole-filled SRTM for the globe Version 4, available from the CGIAR-CSI SRTM 90 m Database (http://srtm.csi.cgiar.org/) Kanzari S, Hachicha M, Bouhlila R, Battle-Sales J (2012) Characterization and modeling of water movement and salts transfer in a semi-arid region of Tunisia (BouHajla, Kairouan)—Salinization risk of soils and aquifers. Comput Electron Agric 86:34–42 Kemper EK (2004) Preface groundwater-from development to management. Hydrogeol J 12:3–5. doi:10.1007/s10040-0030305-1 Kerrou J, Renard P, Tarhouni J (2010) Status of the Korba groundwater resources (Tunisia): observations and three-dimensional modelling of seawater intrusion. Hydrogeol J 18:1173–1190. doi:10.1007/s10040-010-0573-5 Kerrou J, Renard P, Cornaton F, Perrochet P (2013) Stochastic forecasts of seawater intrusion towards sustainable groundwater management: application to the Korba aquifer (Tunisia). Hydrogeol J 21:425–440. doi:10.1007/s10040-012-0911-x Li W, Lu Z, Zhang D (2009) Stochastic analysis of unsaturated flow with probabilistic collocation method. Water Ressour Res 45:W08425. doi:10.1029/2008WR007530 Lin YP, Lee CC, Tan YC (2000) Geostatistical approach for identification of transmissivity structure at Dulliu area in Taiwan. Environ Geol 40:111–120 Liu CW, Jang CS, Chen SC (2002) Three-dimensional spatial variability of hydraulic conductivity in the Choushui river alluvial fan Taiwan. Environ Geol 43:48–56 Matheron G (1970) La the´orie des variables re´gionalise´es et sans applications. Les cahiers du Centre de Morphologie Mathe´matique de Fontainebleau McDonald MG, Harbaugh AW (1988) A modular three-dimensional finite-difference ground-water flow model: techniques of water-

123

Environ Earth Sci (2017)76:699 resources investigations of the United States geological survey. Book 6. Chapter A1. p 586 Nakaya S, Yohmei T, Koike A, Hirayama T, Yoden T, Nishigaki M (2002) Determination of anisotropy of spatial correlation structure in a three-dimensional permeability field accompanied by shallow faults. Water Resour Res 38:35(1)–35(14) Nouiri I, Yitayew M, Maßmann J, Tarhouni J (2015) Multi-objective optimization tool for integrated groundwater management. Water Resour Manage. doi:10.1007/s11269-015-1122-8 Paniconi C, Khlaifi I, Lecca G, Giacomelli A, Tarhouni J (2001) Modeling and analysis of seawater intrusion in the coastal aquifer of eastern Cap-Bon, Tunisia. Transp Porous Media 43:3–28. doi:10.1023/A:1010600921912 Park C, Seo J, Lee J, Ha K, Koo MH (2014) A distributed water balance approach to groundwater recharge estimation for Jeju Volcanic Island, Korea. Geosci J 18:193–207. doi:10.1007/ s12303-013-0063-6 Pool M, Carrera J, Alcolea A, Bocanegra ME (2015) A comparison of deterministic and stochastic approaches for regional scale inverse modeling on the Mar del Plata aquifer. J Hydrogeol 531:214–229. doi:10.1016/j.jhydrol.2015.09.064 Raul KS, Panda NS (2013) Simulation-optimization modeling for conjunctive use management under hydrological uncertainty. Water Resour Manage 27:1323–1350. doi:10.1007/s11269-012-0240-9 Ribeiro L (1993) Evaluating the uncertainty of leakage rates of Tejo and Sado groundwater system. Geostatistics Tro´ia 92:695–707. doi:10.1007/978-94-011-1739-5_55 Sidiropoulos P, Mylopoulos N, Loukas A (2015) Stochastic Simulation and management of an Over-Exploited aquifer using an integred modeling system. Water Resour Manage 29:929–943. doi:10.1007/s11269-014-0852-3 Tammal M, Kili M, El Gasmi El H, Mridekh A, El Mansouri B (2014) Modeling multi-aquifer system of Tadla basin and plateau of phosphates. Int J Innov Sci Res 6:172–180 Thornthwaite CW (1948) An approach toward a rational classification of climate. Geogr Rev 38:55–94 Wang S, Shao J, Song X, Zhang Y, Huo Z, Zhou X (2008) Application of MODFLOW and geographic information system to groundwater flow simulation in North China Plain, China. Environ Geol 55:1449–1462. doi:10.1007/s00254-007-1095-x Xu X, Huang G, Qu Z, Pereira SL (2011) Using MODFLOW and GIS to assess changes in groundwater dynamics in response to water saving measures in irrigation districts of the Upper Yellow River Basin. Water Resour Manag 25:2035–2059. doi:10.1007/ s11269-011-9793-2 Zammouri M, Ribeiro L (2017) Analyzing the effect of transmissivity uncertainty on the reliability of a model of the northwestern Sahara aquifer system. J Afr Earth Sc. doi:10.1016/j.jafrearsci. 2017.02.034 Zammouri M, Siegfried T, El-Fahem T, Kriaˆa S, Kinzelbach W (2007) Salinization of groundwater in the Nefzawa oases region, Tunisia: results of a regional-scale hydrogeologic approach. Hydrogeol J 15:1357–1375. doi:10.1007/s10040-007-0185-x Zammouri M, Jarraya-Horriche F, Odo BO, Benabdallah S (2013) Assessment of the effect of a planned marina on groundwater quality in Enfida plain (Tunisia). Arab J Geosci 7:1187–1203. doi:10.1007/s12517-012-0814-0 Zghibi A, Zouhri L, Tarhouni J (2011) Groundwater modeling and marine intrusion in the semi-arid systems (Cap-Bon, Tunisia). Hydrol Process 25:1822–1836. doi:10.1002/hyp.7948 Zhang D (2001) Stochastic methods for flow in porous media: coping with uncertainties. Academic Press, San Diego, p 368