Hydrogeochemical modeling for groundwater management in arid and semiarid regions using MODFLOW and MT3DMS: A case study of the Jeffara of Medenine coastal aquifer, South‐Eastern Tunisia

The study of water quality and the quantification of reserves and their variations according to natural and anthropogenic forcing is necessary to establish an adequate management plan for groundwater resources. For this purpose, a modeling approach is a useful tool that allows, after calibration phase and verification of simulation, and under different scenarios of forcing and operational changes, to estimate and control the groundwater quantity and quality. The main objective of this study is to collect all available data in a model that simulates the Jeffara of Medenine coastal aquifer system functioning. To achieve this goal, a conceptual model was constructed based on previous studies and hydrogeological investigations. The regional groundwater numerical flow model for the Jeffara aquifer was developed using MODFLOW working under steady‐state and transient conditions. Groundwater elevations measured from the piezometric wells distributed throughout the study area in 1973 were selected as the target water levels for steady state (head) model calibration. A transient simulation was undertaken for the 42 years from 1973 to 2015. The historical transient model calibration was satisfactory, consistent with the continuous piezometric decline in response to the increase in groundwater abstraction. The developed numerical model was used to study the system's behavior over the next 35 years under various constraints. Two scenarios for potential groundwater extraction for the period 2015–2050 are presented. The predictive simulations show the effect of the increase of the exploitation on the piezometric levels. To study the phenomenon of salinization, which is one of the most severe and widespread groundwater contamination problems, especially in coastal regions, a solute transport model has been constructed by using MT3DMS software coupled with the groundwater flow model. The best calibration results are obtained when the connection with the overlying superficial aquifer is considered suggesting that groundwater contamination originates from this aquifer.

To study the phenomenon of salinization, which is one of the most severe and widespread groundwater contamination problems, especially in coastal regions, a solute transport model has been constructed by using MT3DMS software coupled with the groundwater flow model. The best calibration results are obtained when the connection with the overlying superficial aquifer is considered suggesting that groundwater contamination originates from this aquifer.

Recommendations for water resource managers
• The results of this study show that Groundwater resources of Jeffara of Medenine coastal aquifer in Tunisia are under immense pressure from multiple stresses.
• The water resources manager must consider the impact of economic and demographic development in groundwater management to avoid the intrusion of saline water.
• The results obtained presented some reference information that can serve as a basis for water resources planning.
The current research aims to study the water management of the Jeffara of Medenine aquifer system and more precisely to integrate all available data collected and analysed as part of this study in a model that simulates the aquifer system functioning and serves as a tool for sustainable management of this resource by using professional groundwater software package (Modflow and MT3D code). To achieve this goal, groundwater flow and transport models have been constructed based on previous, hydrogeological studies investigations and data.

| GEOLOGICAL AND HYDROGEOLOGICAL SETTING
The Jeffara coastal plain is located in southeastern Tunisia and northeastern Libya along the Mediterranean coast. It extends from the south of the Skhira region north of Gabes to northwestern Libya at the Tripoli region. In its Tunisian part, this geographical entity corresponds to a vast plain that separates the Dahar plateau from the Mediterranean coast.
The selected area for this study comprises a subzone of the Jeffara called the Jeffara of Medenine. It is bounded by latitudes 33°25′ and 34°03′ North, and longitudes 10°07′ and 11°3 0′ East, it contains the maritime Jeffara (peninsulas of Jorf and Akara-Zarzis), the island of Jerba and the plain of Ben Guerdane (Figure 1).
Due to its geographical location, the study area is characterized by an arid climate. The precipitation irregularities in space and time are mainly due to the influence of the Mediterranean Sea, on the one hand, and the influence of the Sahara, on the other.
The study area is characterized generally by high temperatures indicating a warm Mediterranean-type hot climate (Talbi, 1992;Yahyaoui, 2000).
The evaporation rate in the study area is variable. It reaches a maximum rate of 254 mm/ month during the summer (from April to October). It decreases significantly during the winter season to reach 95 mm in February. However, throughout the year, southern Tunisia generally characterized by high evaporation values (Agoubi, 2012).
Despite the severe climatic conditions, the region of Medenine has a dense ephemeral hydrographic network, articulated around five main wadis (Zigzaou, OumZessar, Zeuss, Sidi Makhlouf, El Morra) which drain the rainwater to the sea in the Gulf of Gabes or the numerous sebkhas and saline depressions on the surface. The surface water resources are minimal because they depend mainly on the seasonal rainfall and the inflow of wadis.
On the other hand, the abundance of closed depressions in the forms of sebkhas in the region may affect the quality of the water infiltrated into the aquifers and, consequently, limit groundwater exploitation (OSS, 2006).
Alternations of continental and marine formations characterize the geology of the Jeffara region due to marine transgressions and regressions that occurred in this sector throughout the geological history of the Northern Sahara (Gabtni, Jallouli, Mickusc, Zouari, & Turki, 2009;Jedoui et al., 2003). The geological ages extend from Paleozoic to Quaternary with a large gap from Paleocene to Eocene.
Structurally, the different tectonic movements have led to the forming of the study area as horsts and grabens and the collapse of the coastal plain of Jeffara, which allows lateral communication between the different aquifer levels that exist in the area.
The aquifer system of the Jeffara of Medenine contains valuable groundwater resources that are widely exploited. It mainly contributes to supply the main cities of the governorate of Medenine with drinking water and contributes significantly to the development of the region's economy.
The main layers constituting this multilayered aquifer system are the sandstone of the Triassic of Sahel el Ababsa, the dolomite-limestone of the Jurassic of Zeuss-Koutine and the sand of the Miocene of Jorf-Jerba-Zarzis.
The previous hydrogeological studies have proved that these three layers are hydraulically connected through the faults in this area. The Trias sandstone aquifer discharges into the Zeuss-Koutine aquifer through the Tajra fault. A part of the groundwater flow joins the overlying Mio-Plio-Quaternary aquifer. The Zeuss-Koutine aquifer discharges into the Jorf-Jerba-Zarzis aquifer (Figure 2; Brini, 2018 andYahyaoui, 2000), and they behave as a single layer aquifer called the deep aquifer of the Jeffara of Medenine, our target in this study, which extends deep beyond the coast and the island of Djerba where we can recognize it far from the shoreline, in offshore oil drillings (OSS, 2006). This connection is confirmed by the configuration of the piezometric map, which indicates a continuity of the piezometric level with a regional flow direction from SW to NE towards the sea, which is the natural outflow of the Jeffara aquifer system. Groundwater salinity analyses show a large range of variation from 700 mg/L in the Triassic aquifer in the upstream to almost 7,500 mg/L in the downstream at the level of Jerba Zarzis Jorf Aquifer (Hamzaoui-Azaza et al., 2013). Groundwater salinization may be due to geogenic factors such as contact with salt deposits, seawater intrusion and upcoming deep natural saline water. It is mostly accentuated and facilitated by anthropogenic activities.

| APPROACH AND METHODOLOGY OF GROUNDWATER FLOW MODELING AND SOLUTE TRANSPORT
The approach followed in this study comprises two main steps; the first step is based on establishing a database integrated into the geographic information system (GIS), and the second step is to develop a numerical model to simulate the aquifer system. The stepwise methodology of groundwater modeling includes the building of the conceptual model, the calibration of the flow model in a steady and transient state, and the calibrated model to simulate various groundwater management plans (Anderson & Woessner, 1992).
A schematic view of the applied methodological approach is represented in Figure 2.

| Conceptual model
The conceptualization of the aquifer system to be modeled influences the reliability of the model and its capacity to reproduce the flow system (Hamzaoui-Azaza et al., 2013;Hassan et al., 2016;Hendrayana, 2017). In this study, the hydrogeological framework for the Jeffara of Medenine aquifer system was constructed from the analysis of the drilled wells data and the lithological logs and the results obtained from hydraulic tests and wells monitoring of groundwater levels. The gathered data assisted with understanding groundwater occurrence and flow processes in the area. They also helped to clarify the vertical exchange of groundwater between the three adjacent layers. The conceptual model considers the two aquifers located in the upstream (Trias and Jurassic) and the Miocene situated in the downstream as a single hydrogeologic unit, forming one continuous aquifer (Figure 3). The Jeffara of Medenine receives an insignificant shallow recharge. In the region of Sahel El-Ababsa, the Triassic aquifer is directly recharged by the infiltration of surface water (Besbes, 2010). A direct recharge estimated at 3% of the total precipitation has been assigned for the calibration of the model in the hole aquifer system made by OSS (2006). In the region of Zeuss-Koutine, where the aquifer gets recharged by the infiltration through the ephemeral stream beds, the recharge rate is estimated at 30% (Fersi, 1979;UNESCO, 1972) of the runoff volume. This recharge is maintained constant with time.
The important groundwater parameters are the horizontal transmissivities and the storage coefficients. The initial transmissivity distribution has been constructed based on the results of field pumping tests. The initial transmissivity values were modified during the model calibration process in a steady state. The storage coefficient values have been obtained from values were initially assigned based on the previous hydrogeological studies (OSS, 2006) and repeatedly adjusted during the transient state model calibration.
Groundwater withdrawal rates recorded by the database of the OSS (2006) and the Tunisian Water Authorities over both national and local scales were quantified and used as good individual abstractions in this model.

| Reference state and software choice
The year 1973 was chosen as the reference year for the model calibration under steady-state conditions. This choice was justified because this date corresponds to the beginning of piezometric measurements in the study area and, accordingly, the availability of an observed piezometric map. On the other hand, no significant groundwater abstraction in the region has taken place before 1973.
The governing equation of groundwater flow (Anderson & Woessner, 1992;Bear & Verruijt, 1987) is: where h is hydraulic heads (L) of layers i,; T x, and T y are components of transmissivity (LT −1 ) parallel to the x-and y-axis within layer i, respectively; S is the storage coefficient (1) of layer i; q i is source/sink of water (LT −1 ) in layer i; x, y are Cartesian coordinates (L); and t (T) is the time.

| Boundary conditions
The description of the boundary conditions is part of the definition of the structure of the model and, therefore, strongly conditions its validity.
Mathematically, the boundary conditions correspond to the geometry limits and the values of variables or their derivatives at the limits (Konikow, 1996).
In physical terms, there are three types of boundary conditions (Holzbecher & Sorek, 2005;Konikow, 1996): 1. Conditions with imposed potentials: this condition characterizes the groundwater supply or outlet areas. 2. condition with imposed flows is the case of direct infiltration of rainwater, pumping, or natural supply from neighboring aquifers. 3. Imposed drain condition: it is shown schematically by a fictitious mesh with imposed potential representing the connection between groundwater table and the drain. The data used to establish this model brings together all the (piezometric) measurements known since 1973. The boundary conditions envisaged in the steady state are (Figure 4): 4. Conditions of potentials imposed on the eastern and western limit of the system: the value of potentials is deduced from the piezometric map for the year 1973. 5. The northern and southeastern boundaries are drawn by following the current lines to agree with zero flow limits.

| Groundwater recharge and discharge
In the arid and semiarid zone, groundwater recharge is mainly done by direct recharge, which varies between 2% and 5% of rainfall and by infiltration of runoff into the wadis beds estimated at 30% of the runoff volume (Besbes, 2010).

F I G U R E 4 Model area gridding and boundaries
For the Sahel Ababsa region, 5% of the annual rainfall was assigned as the surface recharge. At the Zeuss-Koutine region, where recharge is by infiltration from wadi beds, the recharge rate estimated to be 30% of the runoff volume (Ben Salah, 2016).
3.5 | Results and discussion of flow modeling in steady state 3.5.1 | Flow model calibration under steady-state conditions The primary purpose of this calibration is to reproduce the groundwater flow conditions in the Jeffara Medenine and provide confidence limits to the model parameters, such as transmissivity and recharge. Model calibration is the process whereby model parameter values are adjusted and refined to provide the best match between measured and simulated values of hydraulic heads and flows.
The calibrated model uses the values of the hydraulic parameters, wells, and boundary conditions depending on the initial field conditions. Thus, we tried different zonation of transmissivity to obtain a distribution of hydraulic head close as much as possible to the measured head values in 1973. Figure 5a represents the matching between the model simulated and the field measured hydraulic head contours. It shows that the density of the more significant number of measurements is located in the eastern part of the aquifer system (i.e., Triassic sandstone sand Zeuss-Koutine aquifers), which allowed obtaining in this zone a better degree of certainty of the measured piezometric levels.
The contour maps of simulated groundwater levels versus observed groundwater levels indicate relatively good agreement between the measured and calculated piezometric contours lines by the steady state ( Figure 5a).

| Water balance
The results of the water balance of the entire model domain are shown in Table 1. Budget terms are expressed in m 3 /s, and the water-budget components are boundary condition, recharge, and pumping.
We can deduce the importance of the lateral contribution and the weak contribution of the recharge due to severe climatic conditions.
With a very negligible exploitation rate (considered to be zero in this year), the system discharges a significant input towards the sea (0.49 m 3 /s).
The steady-state model calibration allows the evaluation of the various aquifer system components. The calculated water balance for each hydrogeological unit shows that the three layers are interconnected and confirm the conclusions of the previous geochemistry study (Hamzaoui-Azaza et al., 2013). This vertical exchange of water explains the local sources of mineralization in each aquifer and the relationship of these three units. The annual groundwater pumping required to meet the water demands in the region is deduced from the DGRE (General direction of water resources) and OSS (DGRE, 1973(DGRE, -2013OSS, 2006). The recharge is generally low and considered to be constant over the simulation period. For the transient state conditions, three parameters determine the aquifer's behavior, groundwater recharge, pumping and storage coefficient. The recharge was constant with time since no annual measurements of runoff are available. The pumping rates were assigned to the model in each well for a water extraction period of 33 years . The calibration, under transient conditions, concerned only with the storage coefficient values. Few data about this parameter are available and its zonation was assigned during the calibration process. The total modeled pumping duration was divided into 33 stress periods.

|
Overall, the computed and observed heads display a satisfactory matching from 1973 to 2015 ( Figure 6). The transient model results show that the most significant drawdowns are recorded in the East of the study area (Triassic sandstone aquifer), where there is the most significant number of wells intensively exploited to meet the growing needs of the region in drinking water and agriculture.
The highest drawdowns (10-20 m) are recorded eastern part of the study area at the Triassic sandstone aquifer that is the most exploited to meet the growing needs for drinking and agricultural purposes.

| Water balance
The water budgets of the Jeffara aquifer obtained from the groundwater flow model simulation results are presented in Table 2; The total water budget over the entire aquifer reveals a great balance between inflows and outflows of water from 1973 to 2015.
From the comparison between the flow components calculated in 2015 with that of the steady-state conditions, the subsurface outflow towards the sea has diminished because of

| Management simulations
Predictions comprise those model simulations that provide the outputs to address the questions defined in the modeling objectives. The time over which predictions are made with a calibrated model should also be related to and limited by the historical record's length. Konikow (1996) stated that a reasonable guideline is to predict only a time comparable to the matched period. After successful calibration, the model can be used to predict the aquifer's evolution under different constraints and at different times. These results can help resource managers develop their strategy for the sustainable exploitation of groundwater in the region (Dassargues, 1997;Hendrayana, 2017).
Due to the constraints of our study region (continually increasing population, water deficit) and to meet growing industrial, agricultural, tourist and domestic needs, the current calibrated transient model was also used to predict the effects of future groundwater withdrawal on hydraulic heads for two different scenarios.

| Management scenario no 1
The first scenario simulates 35 years (2016-2050) of pumping, considering that the conditions that were existed in the year 2015 kept unchanged. However, the effect of current withdrawals becomes more visible even if they remain constant in the future. The drawdown is significant everywhere, but the Triassic sandstone aquifer is the most affected since it records the highest values. In the Jerba region, the drawdown increases by 7 m compared to the year of 2015 while in the Triassic sandstone, the drawdown reaches 30 m.

| Management scenario no 2
The second scenario predicts the aquifer's behavior when running the model up to 2050 when the discharge rates assigned in the first scenario are doubled. The run of this predictive simulations showed a significant decrease in the piezometric levels because of the additional withdrawals. The drawdown will reach 15 m in the coastal area and exceed 35 m in the Triassic sandstone layer.

| SOLUTE TRANSPORT MODELING
Demographic growth, severe climatic conditions, economic and social development in the Medenine region requires an increase in the exploitation of resources to meet all these needs. A drawdown in the piezometric levels was recorded, accompanied by a degradation of the water quality. For all these reasons, it was very necessary to study the impact of the overexploitation of this resource on groundwater quality by applying the solutes transport model, which simulates the Spatio-temporal distribution of the water salinity over the entire Jeffara of Mdenine aquifer system. The groundwater flow models were used as a prerequisite for a solute transport model.
Reliable prediction of contaminant movement can be made if hydrodynamic dispersion and chemical reactions that affect the dissolved chemicals in groundwater can be accurately represented in a systematic model.
For this study, the three-dimensional multispecies transport model MT3DMS package has been used (Chen et al., 2013;Zheng & Wang, 1999).
The MT3DMS code has a modular structure similar to MODFLOW, which allows for the simulation of advection, dispersion/diffusion, source/sink mixing, and chemical reactions (Maliva & Missimer, 2012).
The first step is to obtain realistic initial salinity distribution; the contaminant transport model is developed based on the steady-state hydrodynamic model. It assigns representative values to each of the water bodies, fluxes as well as boundary conditions.

| Solute transport process
The transfer of solutes into groundwater is governed by the convection-dispersion equation, which translates the principle of conservation of solute mass at the local scale, or any point and at any time in the field of study (Bouhlila, 1999). This equation is written:

| Construction of the steady-state solute transport model
In this stage of the research, only an overall concentration, identified as the total dissolved solids (TDS), will be retained.
The conceptualization of the solute transport model is identical to that of the groundwater flow models. It requires to define the initial distribution of the model parameters (effective porosity and dispersivity) and the boundary conditions. The first step is to reconstitute the observed initial state and then calibrate the transport model over a historical period.

| Boundary conditions
The upstream boundaries of the model domain, which represents the subsurface groundwater inflow boundary, are assigned as fixed concentrations based on the observed initial concentration distribution in the study area.
At the boundary with the sea, the concentrations were defined according to the groundwater analyses' data.
Transport conditions are closely related to the direction of flow. The piezometric level adopted is that resulting from the calibrated steady-state flow model.
In our simulation, we imposed an initial concentration fixed at t = 0; C = C0 = 1.5 g/L over the entire area of the water table. This concentration represents a value lower than all the concentrations measured in the aquifer.

| Effective porosity
This parameter is decisive for calculating the sufficient velocities at any point and specifying the convection and dispersion flows. It is also particularly crucial for the distribution of the steadystate and transient levels because it controls all the processes of mixing water and the resulting concentrations. Data about this parameter are rare. The initial porosity distribution is based on available data and previous modeling efforts (OSS, 2006).

| Dispersivity
The dispersivity values are not available in the field. Consequently, the values assigned for this model correspond to that used in the OSS model (OSS, 2006).

| Salinity of the recharge rate
The recharge was not of great influence on the aquifer's salinity, so an average salinity value of 1 g/L has been estimated for the Triassic sandstone and Zeuss-Koutine dolomite-limestone layers and a zero value for the Miocene sandy aquifer.

| Results and discussion of the solute transport model calibration
The solute transport model calibration consists in reproducing the observed history of groundwater contamination in the study area. Much of the flow model output was used as input to the solute transport model. The first stage is to reproduce the initial state of TDS observed in 1973. Figure 7 shows the salinity distribution in the year 1973 extrapolated according to the availability of the observed data. This will serve as a reference for the calibration of the solute transport model to calculate a more consistent initial concentration by flushing the aquifer in the state over 100,000 years.
Due to the lack of recorded data over some parts of the study area, especially in the southern zone, the solute transport model's calibration was restricted to reproduce the distribution of the historical TDS values.
Several simulations have been carried out. The values of porosity and dispersivity have been modified to obtain a salinity distribution, which is as close as possible to that of the reference year 1973.
But the construction of the model developed for the steady flow regime seems insufficient to explain the salinity of the waters of the Jeffara of Medenine aquifer system.
Moreover, neither the recharge nor the exchange between the three layers allows reaching salinity values of the order of 6 and 7 g/L in the region of Jerba-Zarzis.
Thus, to properly calibrate the transport model, a new initial flow model must be designed to include potential sources of solutes that can help to recover salinity in the area.
F I G U R E 7 Spatial distribution map of salinity for the year 1973 According to OSS (2006), these potential sources are: 1. The convective exchanges between the Miocene aquifer and the superficial aquifer were not included in the first conceptual model. 2. The contribution of the Triassic fossil layer to the Miocene aquifer through th Medenine fault. 3. Marine intrusion is a source of degradation of water quality.
Hypothesis 1: Exchange with the lower Evaporite Triassic (fossil) aquifer The initial state of flow has been changed in the hydrodynamic model to consider the exchange with the Lower Triassic aquifer through the Medenine fault.
A new condition (General Head Boundary) has been added to the south of the aquifer. It assumes that the high salinity in the Jerba-Zarzis region comes from its exchange with the Lower Triassic aquifer. The new General Head Boundary condition is shown in Figure 8.
After several calibration attempts by changing the transmissivity and the recharge rate, we succeed in reconstructing the piezometric reference map. The initial piezometry obtained will be used in the transport model.
The transport model was constructed using the same initial conditions of the first model and adding the new general head boundary condition. A concentration of 10 g/L was assigned to this limit.
To obtain salinity values close to those of reference, the sensitivity of the model was studied. The porosities, the longitudinal dispersivities and their relationship with the horizontal and F I G U R E 8 New General Head Boundary condition vertical dispersivities were modifiedill, all the results obtained are close to the calculated salinity map shown in Figure 9.
This indicates that the salinity of groundwater of the Miocene of Jorf-Jerba-Zarzis aquifer is not entirely attributable to the waters of the Lower Triassic.
Hypothesis 2: Exchange with the superficial aquifer The initial flow state in the hydrodynamic model has been changed by adding a new condition (General Head Boundary), representing the convective exchange of the system with the overlying superficial aquifer. The extension of this connection is shown in Figure 10.
After several simulation runs in which the transmissivity and recharge have been modified, the model has reproduced the reference piezometric map that can be used in the transport model.
By adding the new condition in the transport model with a concentration value of 7 g/L in the Jerba-Zarzis region and after several calibration attempts in which the dispersivity and the porosity have been modified, satisfactory calculated TDS values close to those measured in the observation wells were obtained (Figure 11).
The calculated TDS map obtained by the new condition is shown in Figure 12. However, the calculated TDS map does not match very carefully for the whole modeled area, but it can be accepted with its little discrepancies due to the lack of data in some parts of the study area.
F I G U R E 9 Salinity map, in g/L, calculated by adding the lower Triassic aquifer The porosity values used for the calibration of this model ranging from 0.006 to 0.04 and the distribution of longitudinal dispersivity varies between 250 m and 1,500 m. The correlation coefficient between the observed and calculated TDS equal to −0.54 while the variance of the error calibration equal to 5 g/L indicate an acceptable calibration quality.
F I G U R E 1 1 Correlation between measured and calculated concentration values F I G U R E 1 0 Exchange with the superficial aquifer

| CONCLUSION
The management of aquifers includes both qualitative and quantitative aspects. The best way to optimize the exploitation of groundwater resources is to take maximum advantage of all available data relating to an aquifer and combine it with the appropriate physical laws to build a model.
Thus, the main goal this study was to develop a hydrogeochemical model of the Jeffara of Medenine aquifer system to simulate its functioning following natural and /or anthropogenic forcing. MODFLOW code is used to establish this model.
The calibration of the mathematical model under steady-state conditions (1973) used to refine the spatial distribution of transmissivity, to restore the groundwater level at each point and to evaluate the water balance of the aquifer system, despite the lack of data and information concerning the geometry and the hydrodynamics parameters in some places.
The reliability of the model was assessed by comparing the simulated and observed values of the hydraulic head measured in 1973. A satisfactory concordance between measured and calculated levels has been registered, especially in areas where there is sufficient information about the hydrodynamic properties (transmissivity and piezometric).
The simulated groundwater flow model evaluates the water budget of the three modeled aquifers and quantifies groundwater flow exchange between them, which confirms the chemical relationship between these aquifers, as mentioned by many previous works.
The transient model has refined the spatial distribution of the storage coefficient of the aquifer. The water balance resulting from this calibration of 2015 has highlighted a gradual fall F I G U R E 1 2 Salinity map, in g/L, by adding the superficial aquifer in piezometry in response to increasing exploitation. The most substantial drawdowns are recorded in the Triassic sandstone aquifer of Sahel el Ababsa in the East of the study area.
The various scenarios for the exploitation of groundwater in the Jeffara of Medenine have shown that the aquifer system is in a water stress situation. The model results showed that the Jeffara of Medenine is subjected to overexploitation conditions, especially in the East of the study area, where most boreholes are concentrated to meet the needs of the region, especially in drinking water. These areas at risk of water shortage require rapid action for reasoned management. All these results will help the managers and planners set up a priority pumping program that considers the aquifer's resources.
Monitoring of the quality must accompany the assessment of the quantity of water resources in Jeffara of Medenine; this is why a model of solute transport has been investigated using MT3DMS software coupled with Modflow.
The results obtained seem to be insufficient to explain the orders of salinity magnitudes in the Miocene aquifer.
For this reason, new "General Head Boundary" conditions that reflect the Miocene aquifer's communication with the Lower Triassic aquifer and the superficial aquifer have been added to the flow model.
The salinity of the Miocene Jorf-Jerba-Zarzis aquifer was not entirely attributable to the waters of the Lower Triassic. The best result obtained is the one that considers the connection with the overlying superficial aquifer. The salt transport model gives salinity values close to those measured.
The results of the hydrodynamic model and that of the solute transport show that an immediate integrated management plan remains necessary to limit the deterioration of water quality and the drawdown of the piezometry.