The Charf El Akab aquifer, with an area of 17 km2 and a location 20 km southwest away from the city of Tangier (Figure 1) has been used for quite a long time to provide potable water to the cities of Tangier and Asilah. At the moment, the strategic reserves of this aquifer are used as additional security resources in case of an exceptional drought or unavailability of works currently serving drinking water in the cities of Tangier and Asilah.
At this time, the aquifer has an artificial recharge system that includes 11 injections pits and a large infiltration basin. The recharge works have been carried out in stages since 1958 on favorable sites for this particular operation  .
Between the years 1958 all the way to 1995, the injections were provided from the surface waters of river Mharhar via a supply line  . The implementation of the “9 April 1947 dam”, in 1995, allowed the aquifer to be recharged from the El Hachef treatment plant. The cessation of injections between 2004 and 2011 led to a startling plummet in the piezometric level, which in turn led to the resumption of the artificial recharge in 2011.
The study area is agricultural and forestry for the most part. The habitat is moderately dense and dispersed. In the east of the area, there are three major quarries for exploiting biocalcarenites which are currently at an abeyance. The low human density and the sparse activities present in the area of Charf El Akab, imply the existence of the good water quality that is naturally preserved. Nonetheless, water resources are still at risk and the protection of this aquifer is the main concern.
The bottom groundwater is the only one that could be exploited. Consequently, the protection of this latter against any source of contamination is considered to be the top worry; hence the absolute need to assess and map its vulnerability to pollution.
In this respect, the purpose of this work is to:
- On one hand, to establish the most vulnerable areas by using the DRASTIC method. The reason behind this choice is the ease and the availability of the parameters used by this method.
- On the other hand, to simulate the effect of a possible pollution to determine the perimeter of close protection, using numerical modeling by applying the Arcfem extension under the Arcmap environment.
2. Case Study: Charf El Akab Aquifer
The Charf El Akab aquifer extends along the Atlantic coast in the shape of a closed basin (Figure 1). It forms a system without hydrographic network  . On a marl substratum of Tortonian age, the geological formations are essentially detrital and dated to the Pliocene age.
The main aquifer formations, of Pliocene age, are organized in two levels, separated by an aquitard horizon  :
・ A Lower aquifer formation dominated by the calcarenites, sandstones and sands where the lower water table is located and extends over 15.5 km2;
・ A semi-impermeable middle formation (aquitard) of marls and marly sands;
・ An Upper aquiferous formation rich in sand and fine sandstone, even in calcarenites which contains the upper aquifer and extends over 10.5 km2.
Nowadays, the National Office of Electricity and Drinking Water (ONEE) only utilize the lower water table; this is mainly attributed to the good quality this latter has, on top of the fact that it presents a lesser threat to pollution and a very good hydrodynamic characteristics  .
The total thickness of the aquifer is around 344 m  . In the southern and southeastern part of the aquifer, the upper and middle formations are almost
Figure 1. Location of the study area  , modified.
absent. Subsequently, the lower aquifer behaves like a free layer (section B of Figure 2).
During the period 1953-2013, the average annual cumulative rainfall is 710 mm and the annual natural recharge of the lower aquifer was 1.1 Mm3.
2.1. Piezometric Evolution
The piezometric history of the lower water table is portrayed from the measurements taken at the IRE 30/1 piezometer (Figure 3), which represents the longest and most continuous observation series.
This figure shows fluctuations that reflect the sensitivity of this aquifer to the natural and artificial recharge. As matter of fact, the piezometry demonstrates sometimes a drawdown that can reach up to 7 m (knowing that the reference piezometric level is roughly 9 m for the period (1958-1973). This drawdown recorded in the years 1975, 1995 and 2010 is primarily ascribed to the intensification of the pumping intended for the drinking water supply, occasionally joint with a cutback of the precipitations contributions.
Figure 2. Geological sections of the ENE-WSW and NNW-SSE direction. 1: Upper formation; 2: Middle formation; 3: Lower formation; 4: Tangier pelites substratum (marly unit of Tangier); 5: location of existing boreholes  , modified.
Figure 3. Injection and piezometry relationship in the 30/1 piezometer between 1958 & 2013  .
Alternatively, the artificial recharge contributions exhibit an effect on the piezometry especially during the period 1980-1995 where they helped to uphold a piezometric folding of around 3 m, regardless of an intensification of the pumping.
2.1.2. Piezometric Evolution 2004-2009 and 2009-2013
The piezometric map of the inferior aquifer at the end of the 2004 dry season was established based on the piezometer and borehole data collected from the Lukkos Hydraulic Basin Agency (ABHL). This map illustrates a radial flow groundwater, with a depression cone at the center of the aquifer which represents the capture field of ONEP intended for the Tangier and Asilah’s AEP.
Overall, the piezometric curves are parallel to each other. The hydraulic gradient is quite homogeneous, with an average value of approximately 1.4%. In the South-East zone, the hydraulic gradient shows a larger value of about 3.3%. This boost is attributable to the exploitation by pumping of the aquifer, it’s also noted that the largest drawdown are in the center of the aquifer.
Figure 4 showcases a decrease of the piezometric level of roughly 20 m between 2004 and 2009. This is due to the overexploitation of the groundwater combined with the stop of injections since 2004. The average destocking during this period is about 2 Mm3/year.
Figure 5 displays a slim rise in the piezometric levels during the period of 2009-2013, which amounted to more than 3 m in the aquifer’s central zone. This can be explained by the contribution of the artificial recharge which was resumed in 2011.
Figure 4. Piezometry evolution 2004-2009.
Figure 5. Piezometry evolution 2009-2013.
It’s imperative to mention that the 2009 and 2013 piezometric maps (Figure 4 and Figure 5) maintain the same structure as in 2004. Howbeit, the variation of the hydraulic load is not uniform. In the central zone the decrease of the piezometric level reached 20 m between 2004-2009 (Figure 4), whereas the period 2009-2013 (Figure 5) registered an increase of the piezometric level of roughly 3 m.
3. Environmental Vulnerability
The vulnerability DRASTIC method was developed by the EPA (Environmental Protection Agency) in the United States in 1985 and  , for the purpose of estimating the potential for groundwater pollution  . It allows the assessment of vertical vulnerability based on seven criteria.
The final vulnerability index (Equation (1)) is the weighted sum of the seven parameters according to the subsequent formula (    ):
: Depth to groundwater: Distance to the water table, thickness of the unsaturated zone;
: Recharge rate (net): Refill;
: Aquifer media: Nature of the saturated zone;
: Soil media: Nature of the soil;
: Topography, (slope): Topography, tilt in %;
: Impact of the vadose zone: Nature of the unsaturated zone;
: Hydraulic Conductivity of the Aquifer: Permeability of the Aquifer;
: Rating for the area being evaluated: Rating given to each parameter;
: Importance weight for the parameter: weighting factor assigned to each parameter.
DRASTIC does not provide absolute answers. It distinguishes among the vulnerable and the less vulnerable areas. Each parameter in the model is symbolized with an index, otherwise noted as (r), varying on average from 1 to 10.
Table 1 shows the typical notes assigned to each of the seven parameters  . A weighting factor (w) between 1 and 5 is then applied to the various criteria rating  , which amplifies in accordance with the significance of the parameter in the assessment of vulnerability, with the aim of making their respective importance relative in terms of vulnerability.
The depth map of the water plan is based either on direct measurements of the water level in the boreholes profundity, in proportion to the ground, as well as the piezometers implanted in the area, or indirectly by geophysical exploration in the areas which lack water points capturing the water table  . In our case, it is calculated by subtraction from the piezometric level retrieved during the measurements of the piezometric survey carried out in 2013. It involves a raster interpolation of the topographic levels with a method of reverse distance weighting, with 4 neighboring points.
Table 1. Attribution of notes for DRASTIC model indicators  .
The same concept has been embraced to establish the vulnerability maps for both the aquifer recharge and the thickness of the unsaturated zone. Soil texture at the Charf El Akab level is mostly differentiated by a high amount of fine yellow quartz sand, fine and well graded, finely bioetritic  .
The incline map is determined from the topographic maps of the study area. After reviewing the topographic map of the region, 4 ranges of slopes are accentuated. One of the morphological characteristics of Charf El Akab is the monotony of a flat topography in the center and the west zone, with a tilt that does not exceed 2%. In the South and South-East zone, near the artificial recharge zone, we notice a slight increase in this parameter, which reaches up to 18%.
The nature of the unsaturated zone features two major notations  and  . The interpolation of these notations made it possible to discern three zones, with different vulnerability degrees. The hydrodynamic parameters of the aquifers; the transmissivity, permeability, storage coefficient, and the hydraulic gradient determine the rate of pollutants migration into the aquifer in addition to the residence times of pollutants in the saturated zone  .
The amalgamation of the seven thematic maps presented above allowed us to draw up the intrinsic vulnerability map of the Charf El Akab aquifer. There are three classes with different degrees of vulnerability (Figure 6):
Taking into account the vulnerability assessment criteria by the DRASTIC method (Table 2), the Charf El Akab lower aquifer vulnerability map (Figure 6) uncovers three zones of different degree of vulnerability. The low vulnerability areas are located in the limit sector of the two layers, the zone of medium vulnerability concerns the northern part which coincides with the sectors where the two layers are superimposed (existence of the aquitard) and the dayas, whereas the very vulnerable zone happen to be in the south and south-east (free zone of the lower formation) of the aquifer and in the quarries of biocalcarenites exploitation (currently at a standstill).
4. Mathematical and Numerical Approaches for the Protection Perimeters of Water Intakes
The notion of protection perimeters is extensively discussed in the international literature (  -  ). These perimeters match up to a zoning established near a water capture (source, well, borehole, etc.) intended for human consumption, in order to preserve its quality and prevent its contamination  .
The criteria for demarcating the protection perimeters are of a spatio-temporal nature, specifically: the distance to the well, the drawdown, the transfer time taking into account the rate of degradation of the pollutant, the flow limits, the soil purifying power, etc. The choice of a criterion depends on technical, socio-economic and regulatory considerations  .
The most prominent protection around the world includes the following zoning:
・ Zone I: Immediate protection against physical degradation or direct insertion of pollutants into the capture (10 to 30 m from the catchment site).
Figure 6. Vulnerability map according to the DRASTIC method.
Table 2. Criteria for vulnerability assessment in the DRASTIC method  .
・ Zone II: Close protection against bacterial or viral pollution (defined as a transit time of approximately 50 days to reach the catchment, with a distance of at least 50 m from the structure).
・ Zone III: Remote protection against persistent pollutants (commonly bounded by the entire capture).
The gears used to restrict the protection perimeters are diverse. There are: geological mapping, test pumps, water appraisal, study of the capture environment, tracing  and numerical modeling. It must be noted that, ordinarily, protection is based on purely hydrogeological criteria.
Zones II and III are customarily outlined from analytical or numerical models  , which can be deterministic or stochastic (    ).
For the most part, only the advection process of flow is considered. Numerous authors, however, take into consideration the advective/dispersive processes (      ).
Today, the close protection perimeters are frequently defined from formulas giving the dimensions of the “calling zone” and the representation of the transfer times of the groundwater, in the form of isochronous curves. The calling zone is the part of the piezometric surface influenced by pumping, where the water is intended to converge to capture.
The isochronous curves are the sites of the points surrounding the capture, symbolizing a given transfer time of the groundwater up to the catchment. They can be calculated using relatively simple formulas, taking into consideration only the convective type transfers  .
The utilized formulas are based on the flow laws in porous media and on the principle of mass conservation (   ).
Considering bacteriological pollution, it is estimated that the life span of the pathogenic bacteria in groundwater does not exceed 50 days according to  .
Generally, the isochron 50 is the limit beyond which a particle cannot achieve capture within a time period shorter than a fixed duration, usually of 50 days (   ).
For the sake of calculating the close protection perimeter of Charf EAkab, we will present the application of a numerical approach, based on the finite element method, in order to resolve the diffusivity equation coupled with that of the mass transfer.
It should further be noted that the resolution of the transport equation uses the inverse direction of the Darcy velocities obtained by the resolution of the diffusivity equation. A certain amount of solute is then injected briefly into the pumping wells. The propagation of the pollutant around the latter makes it possible to delimit numerically isochrones including that of the “50 days” which corresponds to the perimeter of close protection (    ).
4.1. Introducing the Arcfem Application (Arc Finite Element Modeling)
The ArcFem extension in the ArcGIS environment, where the central application is ArcMap, comprises several menu bars  . This application is based on the method of finite element as a numerical method for the resolution of the diffusivity and pollutant transport equations.
4.2. Conceptual Model and Hydrodynamic Modeling of the Charf El Akab Aquifer
Commonly, in the Arcgis environment, the data applied to solve the diffuse and transport equations are in the form of thematic maps. The digitization of these cards makes it feasible to classify their contents by themes, therefore, assisting the analysis and the updating of the information. The decomposition of the data into several themes also bestows greater flexibility in the output operations.
According to Holländer (1989) the longitudinal dispersion αL is generally between 12 and 480 m. It ought to be noted that, usually, lateral dispersion embodies more or less 10% of the longitudinal dispersion (     ). In fact, we have given a value of 50 m for the longitudinal dispersivity and a value of 15 m for the transverse dispersivity. In addition, the hydrodynamic and all saturated thickness characteristics of the aquifer are shown in Table 3. In effect, these data are the results given by the hydrodynamic model  . The Charf El Akab aquifer has 11 operational boreholes in 2013, regulated with a total flow of 604 l/s (Table 4).
The domain then undergoes discretization into 3605 triangular mesh elements with 1867 nodes in the horizontal orientation with an average spacing “characteristic length” of 100 m between the nodes (Figure 7). The mesh conception of the study area was generated with refining at the pumping boreholes with a surface ranging from 3000 m2 to 200 m2.
The boundary conditions adopted in the model design are the following:
・ Null flow conditions: corresponding to the limit of the aquifer extension. All limits are at zero flow since the basin acts like a closed reservoir;
・ Dirichlet conditions with imposed potential: these potentials signify the piezometry equivalent to the edges of the aquifer;
・ Newman conditions: the flow rates were assigned to the nodes corresponding to the positions of the ONEP operating boreholes (Figure 7);
・ Flow condition imposed on the areas associated with natural and artificial recharge;
・ The different stages of the modeling were the following:
・ Simulation in steady state of an average situation corresponding to the piezometric levels of the year 2013;
・ Calculation of the darcy speed in each node;
・ Solving the pollutant transport equation and determining the protection perimeter.
Table 3. Input data and synthesis of Charf El Akab properties  .
Table 4. Borehole characteristics (source ONEE).
Figure 7. Mesh domain and application of the initial and the boundary conditions.
4.3. Results and Discussion
4.3.1. Simulation in Steady State
The hydrodynamic modeling, by ArcFem, goes through 3 junctures; these latter are detailed by  .
・ Pre-calculation of the piezometric levels;
・ Calculation of the piezometric levels in each node;
・ Production of iso-value curves of the piezometric levels.
Figure 8 portrays the validation result of the interpolation window parameters. Where it shows the iso-value curves generated with a width-space of 5 m, in addition to the impact of the pumping on the piezometry of the water table.
In the case of a steady state, we conducted the calibration of the model by fine-tuning the natural recharge and permeability. A satisfactory piezometric state is obtained.
Figure 8. Superposition of the calculated hydraulic piezometry and the reference one.
4.3.2. Calculation of the Darcy Velocity in Each Node
The speeds calculated in each node using the “Darcy velocity Calculation” window are employed to resolve the transport equation and to determine the flow direction (Figure 9).
4.3.3. Dimensioning of the Close Protection Perimeter
The close protection perimeter thus calculated is symbolized in yellow in the ArcMap interface with all other layers (Figure 10). Its area equates to 72.76 ha.
The projection of the close protection perimeter on the vulnerability map (Figure 11) shows that it corresponds with the medium vulnerability area and encloses all the operational boreholes. The south of this zone is distinguished by a dense habitation which has the ability to expose the waters of the water table to the pollution through the discharges of the untreated wastewater, at the abandoned wells or septic pits.
Setting up the close protection perimeter will accordingly further protect the water table against pollution.
This study allowed us to recognize three distinct zones based on their vulnerability to pollution:
・ The area of low vulnerability, which characterizes the limits sector of the two layers. This area is feebly vulnerable to potential pollution; it happens to be the artificial recharge zone.
Figure 9. Direction of flows obtained by Darcy’s speeds.
Figure 10. The calculated close protection perimeter.
Figure 11. Overlay of the close protection perimeter and the vulnerability map.
・ The area ofmedium vulnerability: this latter depicts the northern part which coincides with the areas where there is the superposition of the two layers (existence of the aquitard) and the dayas (small lake). It is thus moderately vulnerable to pollution; it focalizes in the capture area.
・ The zone of high vulnerability, this particular zone portrays the south and south-east part (free zone of the lower formation) of the aquifer and in the quarries of biocalcarenites exploitation (presently at a standstill). This exhibits that this part of the water table is especially exposed to pollution. It therefore necessitates proper protective measures.
It’s remarked that at these areas, the thickness of the unsaturated zone is important and plays its role of protection. The aquifer is generally protected, and the groundwater pollution can be caused essentially by untreated wastewater discharges, abandoned wells or septic pits of residential, and waste at abandoned quarries.
The determination of the close protection perimeter, calculated by the ArcFem numerical method in the ArcMap interface, gave us an area equal to 72.76 ha around the operational boreholes. This will strengthen the water table protection against pollution.
The superimposition of the close protection perimeter and the vulnerability map reveals that it coincides with the medium vulnerability zone which encloses all the operational drill holes. The implementation of the close protection perimeter will consequently further defend the water table versus pollution.
In order to guarantee the quality of the water distributed to the population, we recommend protecting the aquifer by simply controlling public discharges (especially in the abandoned quarries), wastewater, agricultural waste, and by implementing protection perimeters around water captures zones. The protection to be accomplished through perimeters is a complementary protection to the more general protection provided by the legislation against spills, discharges, streams, direct or indirect deposits of water or material.
After having identified the vulnerable zones, and outside the immediate protection perimeters that are already put in place by the ONEE, we have proposed the limits of the close protection perimeter on the basis of the hydrogeological conditions of the aquifer system operated. The ArcFem numerical method in the ArcMap interface opens a new numerical modeling area.