**Faten Jarraya Horriche and Sihem Benabdallah \***

Centre of Water Research and Technologies, Geo-resources Laboratory, Techno-park Borj-Cedria, BP 273 Soliman, Tunisia; faten.horriche@topnet.tn

**\*** Correspondence: sihem.benabdallah@certe.rnrt.tn; Tel.: +216-98-404-777

Received: 11 December 2019; Accepted: 20 January 2020; Published: 25 January 2020

**Abstract:** This study was carried out to examine the impact of an artificial recharge site on groundwater level and salinity using treated domestic wastewater for the Korba aquifer (north eastern Tunisia). The site is located in a semi-arid region affected by seawater intrusion, inducing an increase in groundwater salinity. Investigation of the subsurface enabled the identification of suitable areas for aquifer recharge mainly composed of sand formations. Groundwater flow and solute transport models (MODFLOW and MT3DMS) were then setup and calibrated for steady and transient states from 1971 to 2005 and used to assess the impact of the artificial recharge site. Results showed that artificial recharge, with a rate of 1500 m3/day and a salinity of 3.3 g/L, could produce a recovery in groundwater level by up to 2.7 m and a reduction in groundwater salinity by as much as 5.7 g/L over an extended simulation period. Groundwater monitoring for 2007–2014, used for model validation, allowed one to confirm that the effective recharge, reaching the water table, is less than the planned values.

**Keywords:** artificial recharge; groundwater; treated wastewater

#### **1. Introduction**

Water is a critical resource in many Mediterranean countries because of its scarcity and uneven geographical, seasonal, and inter-annual distribution [1,2]. In a country with frequent water stress, several national strategies were established in Tunisia to optimize the management of water resources to meet growing freshwater needs and planning for climate change adaptation [3,4]. The use of reclaimed municipal wastewater was considered one of the main axes in the national water strategy, mainly for agriculture and groundwater artificial recharge [3]. Thus, treated wastewater (TWW) is used to improve groundwater storage and reduce seawater intrusion in coastal aquifers. Previous studies have demonstrated that such alternatives are reasonable when conventional freshwater sources become very limited [5–9]. Artificial recharge using conventional freshwater has already been acknowledged for inducing groundwater level rise and for improving water quality in several aquifers in Tunisia. The Teboulba coastal aquifer registered a rise of groundwater level up to 30 m following the artificial recharge in wells during six years [10]. The artificial recharge in the Zeroud riverbed of Kairouan aquifer led to an increase of water table between 0.2 and 5.25 m for a distance of up to 8 km [11]. A positive impact of the artificial recharge in El Khairat wadi was confirmed by the authors of [12], indicating an increase in water table between 0.4 and 2.6 m.

This study mainly focuses on the Korba aquifer, located in North-Eastern Tunisia. This aquifer registered groundwater decline and water salinization over the past three decades due to high groundwater pumping rates for agriculture uses. This serious threat has motivated many authors to study seawater intrusion by modelling [13–15] or by hydro-geochemical investigations in order to determine the spatial extension of the salinization or to identify the origins and the mechanisms governing its contamination [16–19]. The construction of the Korba–Elmida artificial recharge ponds, using treated municipal wastewater, was among the actions undertaken by the water authority for aquifer storage recovery to overcome groundwater level decline and salinity rise for the Korba aquifer [20].

Numerical modelling and hydro-geochemical investigations were often used by scientists to estimate aquifer storage or to carry out groundwater recovery based on exploratory simulations and scenarios development [8,9,21–24]. They can also help to estimate the potential benefits for constructed structures on hydrologic conditions under a range of management scenarios. Groundwater flow and solute transport models are applied as predictive tools to plan and quantify the impact on the local groundwater, and determine geochemical processes and the resulting recovery efficiency [25,26]. Basic descriptions of various physical and chemical equilibrium for solute transport models are given by the authors of [26], indicating that the development of geochemical transport models or hydrogeochemical models represents still new pursuit, although some mathematical flow models date back to the late 1960s. Thus, solute-transport models are inherently more complex in terms of conceptualization and governing equations, numerical methods, parameter estimation, and boundary conditions, as well as concerns about model complexity [27]. It should be noted that analytical geochemical methods and path way modelling were also used when geochemical data was available to explain groundwater salinity variations and to support groundwater management and preventing salinization [28,29]. Yet, numerical modelling was frequently used for coastal aquifers to simulate seawater intrusion in natural and anthropogenic conditions and to predict also climate change impact [13,26]. In this paper, we evaluate the impact of the constructed artificial recharge ponds on groundwater level and water salinity by using 2D flow and solute transport models. First, we give an overview for the hydrogeological characterization of the artificial recharge site, at a small scale. Secondly, we simulate the predicted impact of the artificial recharge site built on flow and solute transport calibrated models based on several scenarios. Sensitivity tests for the effective artificial recharge are then performed during the artificial recharge period between 2009 and 2014, by using a monitoring network for groundwater levels and salinities.

#### **2. Materials and Methods**

#### *2.1. Study Area Description*

The Korba aquifer is located in the Cap Bon Peninsula in Tunisia (Figure 1a). It covers an area of 500 km<sup>2</sup> and is bounded by the Mediterranean Sea to the south-east. The average annual rainfall is about 450 mm per year computed for the period 1959–2009 using the national rainfall network over the study area. The study area is characterized by its important agricultural activity, leading to high water demand and a groundwater overdraft [30]. The shallow aquifer is composed of the Plio-Quaternary formation, characterized mainly by sand surmounting the Miocene clay bedrock [30]. Groundwater levels and water quality have been monitored since the 1960s. Early measures showed that water salinity did not exceed 3 g/L [31]. Starting in 1970, the groundwater level declined and water quality degraded as well; an increase of the size of the affected areas was registered mainly in some overexploited localities where the groundwater level declined below the sea water level and the salinity exceeded 5 g/L reaching in some points 10 g/L [31]. This situation was a result of groundwater over abstraction, which generated a reverse hydraulic gradient and seawater intrusion [13–18]. In addition, other authors [18,30] attributed the groundwater salinization to the increase in irrigated agricultural areas, which induced soil leaching and migration of fertilizers to the aquifer.

In order to overcome the problem of groundwater level decline and water quality deterioration in the Korba aquifer, water authorities planned the construction of the Korba–Elmida artificial recharge site, using treated municipal wastewater. The site location was selected after a feasibility study that covered the Cap Bon region, considering several technical criteria related to geologic aspects, depth of water table, uses of groundwater, location of nearby wastewater treatment plant, and other economic constraints [20]. The selected site is located in the north-east of Tunisia, about 300 m north of the Korba wastewater treatment plant and 1.5 km from the coast (Figure 1a). The artificial recharge site, consisting of three recharge ponds with 1500 m2 each, was designed for a recharge rate of 1500 m3/day. In fact, the monitoring of the actual recharge rate, which we realized between 2009 and 2014, ranged between 800 and 1745 m3/day for duration of 0 to 10 h per day [32]. The TWW was transferred from a nearby wastewater treatment plant (WWTP), which received an average inflow of 5000 m3/day with an average salinity 3.3 g/L in the year 2003 [20].

**Figure 1.** (**a**) Location map for the Korba aquifer with the treated wastewater artificial recharge site (TWWARS) and model setup structure and boundary conditions and simplified surface geology, (**b**) scatter diagram calibration of steady state groundwater flow modelling, (**c**) and groundwater level (GWL) map simulated for steady state (1971).

#### *2.2. Hydrogeological Characterization*

The efficiency of the artificial recharge depends on the aquifer hydrogeological characteristics [5,33,34]. The site characterization was based on the field hydrogeological investigations carried out before the project implementation and during artificial recharge period. Soil properties, which control the flow rate of infiltration and the downward percolation, are of special importance to this type of technique. Thus, subsurface formation was identified using well logs and test infiltration results.

A well (W), shown in Figure 2, was drilled up to a depth of 52 m (Figure 2a) and was used to investigate the subsurface formations and to estimate the horizontal transmissivity and the storage coefficient. In order to identify the unsaturated zone, eight complementary drills were performed to check the thickness of sandstone layers. Suitable sites for artificial recharge were identified using four surface infiltration tests (IF), followed by twelve infiltration tests in ditches (PM) (Figure 2b) to the depth of 1.8 m corresponding to the depth of the planned artificial recharge ponds.

**Figure 2.** (**a**) Log for well (W) and (**b**) infiltration velocity map with location for test infiltration in subsurface (IF) and in ditches (PM).

#### *2.3. Modelling Approach*

The MODFLOW code was used to solve the 2D groundwater flow equation in saturated porous media by means of the finite-differences method [35], and the MT3DMS code was used to solve the solute transport equation [36]. Both codes were used within the Processing Modflow package [37]. The equations were resolved within a regular square grid of 25 <sup>×</sup> 10<sup>4</sup> m<sup>2</sup> cell size. The models were calibrated during steady and transient states for the period between 1971 and 1996, by using a manual process based on the trial and error method, referring to measured values in monitoring network. They were thereafter validated for the period 1997–2006. The model area was limited by the boundaries of the Plio-quaternary outcrops (Figure 1a). At the upstream, the Miocene formation outcrop was considered to be part of the modelling area. The bedrock is formed by Miocene clay formation with a thickness ranging between 100 m upstream and 1800 m in the center of the aquifer. The groundwater recharge was mainly assured by rain infiltration with about 7% of the average annual rainfall, as was estimated by the authors of [13,38]. The natural outlet represents the sea, the salt areas (marshes), and the downstream draining areas of the wadis. The coastal limit is presented by a fixed head condition of 0 m, and a fixed water salinity of 38 g/L, which is the average value for the Mediterranean Sea. Near the irrigated areas, salinity of recharged water was increased to consider irrigation return flow as justified by [11,16].

The calibrated models were used to predict the impact of the planned artificial recharge on groundwater level and salinity with a focus on the surrounded region. Several scenarios were used considering hypothesis for natural and artificial recharge and groundwater pumping rates during the period between 2007 and 2050. Model's results were validated during the period 2007–2014 using the actual recharge rates based on monitored values.

### **3. Results and Discussions**

#### *3.1. Characterization of the Artificial Recharge Site*

Hydrogeological setting performed before and after the artificial recharge allowed for the characterization of the subsurface and surface zones. According to the drilled well (W) and the core drills (Figure 2), subsurface formations are mainly sandy with thin sandstone layers. The hydraulic conductivity of the sampled sand from the well, using a grain size analysis, was primarily calculated to 3 <sup>×</sup> 10−<sup>4</sup> m/s. Tests performed for several depths in the core drills showed a range of hydraulic conductivity between 2 <sup>×</sup> 10−<sup>3</sup> m/s for fine sand and 2 <sup>×</sup> 10−<sup>6</sup> m/s for sandstone. Sampling analysis confirmed that unsaturated zone is mainly composed of sand; more than 95% have grain size less than <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>3</sup> m. The pumping test carried out in well W produced an estimate horizontal transmissivity of <sup>4</sup> <sup>×</sup> <sup>10</sup>−<sup>3</sup> m2/s and a storage coefficient varying between 4.5 <sup>×</sup> <sup>10</sup>−<sup>4</sup> and 6 <sup>×</sup> <sup>10</sup><sup>−</sup>4.

The infiltration test results, applied in surface and ditches, gave a vertical infiltration velocity ranging between 1 and 65 m/day depending on the soil type. High values corresponded to sand, whereas low values were linked to sandstone and consolidated sand. The south-east zone indicated infiltration rates higher than 30 m/day, which can be considered as suitable values for artificial recharge, leading to the conclusion that the selected site was well adapted for the construction of the artificial recharge ponds.

#### *3.2. Model Calibration and Validation*

The calibration of the flow model during steady state allowed for the assessment of groundwater recharge ranging between 7% and 11% of annual average rainfall for the Quaternary outcrops. These values are close to those used by previous studies [13–15,38]. The calibration results were satisfactory by comparing the simulated groundwater levels with the observed values (Figure 1b). The coefficient of determination (R2) is about 0.97. The groundwater renewable resources are evaluated as 1.31 m3/s during 1971, as shown in Table 1. Aquifer recharge was provided mainly from recharge at outcrops, and wells' pumping represented the main outlet fluxes. During transient state, the calibrated porosity ranged between 0.05 and 0.35, which are close to the constant value of 0.12 used by [26], and those provided by [14] ranging between 0.04 and 0.25. The general trends of groundwater decline were reproduced by the calibrated model, with a maximum error of +/−1 m (Figure 3a). The results confirmed that the most affected areas were located between Korba, Diar Elhojjej, and Tafeloune agglomerations. The flow model was validated for the period 1997–2006. It reproduced close water level variation for several piezometers as shown in Figure 3a.


**Table 1.** Calculated water balance in steady state (1971).

**Figure 3.** (**a**) Calibration and validation in monitoring wells during transient period (1972–2006). Calibration was established during the period 1972–1996, and validation was performed during the period 1997–2006; (**b**) simulated salinity during transient period (1972–2006); and (**c**) salinity map for 2006.

The calibration of the solute transport modelling was performed using the calibrated flow model for the period 1971–2006. The sea, represented by fixed head and fixed concentration conditions, contributed to groundwater salinization by convection and dispersion. Figure 3b presents simulated salinity results for the transient period. Calibration for salinity was not possible, in part due to the few available measurements but mainly due to the complexity of processes [27]. Simulated salinity for 2006 (Figure 3c) confirms the existence of seawater intrusion along the coastal area between the two water courses (Figure 1), spreading out in the aquifer up to 5 km.

#### *3.3. Artificial Recharge Simulations*

Two simulations (SIM1 and SIM2) were conducted during a forecast period from 2007 to 2050. In SIM1, we maintained the same boundary conditions throughout the simulation period, assuming the continuity of the current conditions of recharge and extraction rates as in 2006 and omitting the additional artificial recharge rate from the treated municipal wastewater. According to this scenario, the model predicted a groundwater decline of less than 2 m near the planned artificial recharge site by the year 2050 and an increase in salinity reaching 12 g/L along the coast.

In SIM2, we considered the same boundary conditions as in SIM1 and assumed an additional recharge flux of TWW of 1500 m3/day with a salinity of 3.3 g/L to be injected in the recharge pools, represented by one cell in the model. The comparison of the simulated groundwater level for 2050 according to SIM1 and SIM2 showed that an increase reaching a maximum of 2.7 m can be achieved (Figure 4a). The influenced area around the site, considering a minimum change of 0.1 m, is around 80 km2. Moreover, the groundwater salinity can be reduced under recharge conditions by a maximum of 5.7 g/L near the site with a recovery area of around 26 km<sup>2</sup> (Figure 4c). It extends over a 6 km distance, including a 17 km<sup>2</sup> of recovery area with a salinity variation of less than 1.5 g/L.

**Figure 4.** Impact of artificial recharge in 2050: groundwater level increase (**a**) SIM2—SIM1 and (**b**) SIM3—SIM1; groundwater salinity decrease: (**c**) SIM1—SIM2 and (**d**) SIM1—SIM3; and (**e**–**g**) groundwater salinity variations in wells P1, P2, and P3.

Results relatives to wells P1, P2 and P3, located, respectively, downstream, upstream, and close to the artificial recharge area as shown in Figure 4c, were used to compare the temporal evolution of groundwater salinity. Figure 4e–g also shows that artificial recharge may reduce salinity, especially near the ponds (P1 and P3). However, its impact was lesser away from the recharging ponds (P2). Salinity decrease reached 4.2, 1.9, and 5 g/L, respectively, for P1, P2, and P3. The salinity decrease occurred during the first ten years of the simulation by an annual decrease of 0.2 g/L (P1) and became insignificant thereafter due to the constant boundary conditions for the aquifer recharge and abstractions. The annual decrease was at a maximum in well P3, located close the ponds.

An additional scenario SIM3 was considered in which only 70% of the artificial recharge rate would reach the water table following the evaporation and losses in the unsaturated zone. This value would be justified for semi-arid regions, and would be close the rates calculated by [38] during the artificial recharge in the plain of Kairouan. As a consequence of this decline, we obtained a lower impact on groundwater level (Figure 4b,d). In fact, the positively influenced area was reduced to 58 km2 and the maximum increase of groundwater levels reached 1.6 m. For the salinity, the influenced area was almost maintained, with a maximum decrease of 4.4 g/L near the recharge ponds.

#### *3.4. Sensitivity Tests*

In this section, we focus on the artificial recharge period between 2007 and 2014. In fact, the artificial recharge was practiced during 2009–2014 with a rate ranging between 0 and 1745 m3/day. Given the uncertainties related to the clogging, the evaporation from the site, and from the unsaturated zone, we investigated whether the measured actual rates contributed effectively to the groundwater recharge by using P3 as an observation well. Two additional simulations (SIM4 and SIM5) were considered. In SIM4, the actual observed rates of the artificial recharge were used, while in SIM5, we considered losses by evaporation and clogging, assuming that only 70% of the measured rates was effective. In both simulations, we maintained an average salinity of 3.3 g/L for the TWW. In fact, the measured salinity of the TWW varied between 3 and 5 g/L. For the starting year of recharge, 2009, the variation of groundwater level for SIM4 is close to the observed values (Figure 5a). As of 2010, SIM5 fits better the observations. These results can be justified by a decrease in the efficiency of the artificial

recharge pools due to clogging phenomena. In fact, we noticed that infiltration velocity decreased during this period. Before the artificial recharge operation, the infiltration velocity was more than 30 m/d. In June 2010, we measured an infiltration velocity of about 1.3 m/d. It was the result of the accumulation of suspended particles in subsurface soil, causing a decrease of the saturated hydraulic conductivity. This decrease of infiltration velocity is consistent with those reported by the authors of [39] during infiltration of TWW in a 3 m diameter column. The decrease of infiltration rates from 3.3 to 0.08 cm/h was due to the development of a clogging in the uppermost layer. The assumption of considering an effective rate of 70% of the recharge rate is thus justified.

For the period 2011–2014, we registered a trend for a general decrease of measured groundwater levels. Sensitivity tests were run in order to reproduce this decrease by reducing the infiltration rate linearly over the simulation period, starting with 70% in 2009. The best simulation, SIM 6 presented in Figure 5a, corresponds to an infiltration rate decreasing linearly with a negative slope of 0.12.

Regarding salinity, the observed variations do not match the simulated values of SIM 4, 5, and 6 (Figure 5b). However, they are much closer to the results given by SIM 2, which corresponds to the projected artificial recharge with a rate of 1500 m3/day. Given that this volume is hard to achieve, we elaborated further sensitivity tests on the salinity of TWW by reducing the average salinity input value. The best simulation, SIM 7, corresponds to the same artificial recharge rates as for SIM 6 with a TWW salinity of 1 g/L. Simulated salinity results are much closer to the observed values. Given that the monitored salinity for the TWW varied between 3 and 5 g/L, we can assume that some of the salinity was retained by the unsaturated zone.

**Figure 5. Well P3** (**a**) Observed and simulated variations of groundwater level; (**b**) Observed and simulated variations of groundwater salinity.

#### **4. Conclusions**

The Korba–Elmida artificial recharge pond using TWW was a pilot site to investigate the recharge impact on groundwater level and salinity. The hydrogeological characterization allowed for the calculation of vertical infiltration velocity, which reached more than 30 m/day near the artificial recharge site. Sand formation was identified as the most suitable for infiltration. The planned recharge rate of 1500 m3/day, and applied to an area of 3000 m<sup>2</sup> (two ponds), is equivalent to 0.5 m/day, which is very low compared to infiltration capacity of sand. Thus, there is still the potential to increase the rate of the artificial recharge to enhance the impact on groundwater level and salinity.

The groundwater flow and transport models allowed for the prediction of the artificial recharge impact on groundwater. The results of the simulations showed that, in the long-term, the artificial recharge with a rate of 1500 m3/day would induce a maximum recovery of groundwater level of up to 2.7 m and a maximum relative decrease in salinity of 5.7 g/L. The influenced area encompassed up to 26 km2 around the site, extending about 10 km belt along the coast. Hence, it was found that the implementation of an artificial recharge site would reduce groundwater decline, improve water salinity, and reduce seawater intrusion given that the water abstraction stays at the actual level. The results of the simulations may have been affected by predictive rainfall, artificial recharge, and abstraction, considered constant for both simulations. The sensitivity tests based on modifying the recharge rate and the reclaimed water salinity allowed for the calculation of different variabilities for groundwater level and water salinity.

The sensitivity of the models to the effective artificial recharge was performed during the period 2007–2014 using the affective artificial recharge rates and the observed groundwater level and salinity. We concluded that effective artificial recharge is less than the planned values. This can be justified by the evaporation and losses in the unsaturated zone and by the clogging phenomena in the ponds. Further, the impact of the artificial recharge is insignificant in terms of salinity.

**Author Contributions:** Conceptualization, F.J.H.; methodology, F.J.H. and S.B.; software, F.J.H.; validation, F.J.H. and S.B.; formal analysis, F.J.H.; investigation, F.J.H.; resources, F.J.H.; data curation, F.J.H.; writing—original draft preparation, F.J.H. and S.B.; writing—review and editing, F.J.H. and S.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding and the APC was funded by the Tunisian Scientific association "Association Eau et Dévelopement".

**Acknowledgments:** The authors would like to thank the staff of the water authority DGRE and its local representatives for their support and interest in this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


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