Using Analytical Hierarchy Process and Multi-Inﬂuencing Factors to Map Groundwater Recharge Zones in a Semi-Arid Mediterranean Coastal Aquifer

: Mapping groundwater recharge zones (GWRZs) is essential for planning artiﬁcial recharge programs to mitigate groundwater decline and saltwater intrusion into coastal aquifers. We applied two multi-criteria decision-making approaches, namely the analytical hierarchy process (AHP) and the multi-inﬂuencing factors (MIF), to map GWRZs in the Korba aquifer in northeastern Tunisia. GWRZ results from the AHP indicate that the majority (69%) of the area can be classiﬁed as very good and good for groundwater recharge. The MIF results suggest larger (80.7%) very good and good GWRZs. The GWRZ maps improve groundwater balance calculations by providing estimates of recharge-precipitation ratios to quantify percolation. Lithology, land use / cover and slope were the most sensitive parameters followed by geomorphology, lineament density, rainfall, drainage density and soil type. The AHP approach produced relatively more accurate results than the MIF technique based on correlation of the obtained GWRZs with groundwater well discharge data from 20 wells across the study area. The accuracy of the approaches ultimately depends on the classiﬁcation criteria, mean rating score and weights assigned to the thematic layers. Nonetheless, the GWRZ maps suggest that there is ample opportunity to implement aquifer recharge programs to reduce groundwater stress in the Korba aquifer.

The region is one of the most productive agricultural areas in Tunisia, contributing about 15% of the country's total agricultural production [9]. Major crops grown in the study area are vegetables (potatoes, tomatoes, strawberries, etc.), vineyard plantations and citrus [8,11]. The region is characterized by a semi-arid climate (average annual rainfall: 420 mm) with Mediterranean influence. Average monthly high temperature occurs in July (37 °C) and the lowest average monthly low is 6.1 °C in January. There are ephemeral streams known as wadis, which originate from the Sidi Abderrahmen anticline and receive seasonal runoff. The most important wadis in the region (i.e., Chiba, M'laabi and Lebna) are obstructed by upstream dams. Consequently, groundwater has become the primary water resource in downstream areas.
Three main geological formations constitute the Korba aquifer system (Figure 2a) [2,5,6,10,40]: 1. The Tyrrhenian Quaternary made of arenitic limestones overlying conglomeratic units with a total thickness ranging 10-50 m, which forms a 1.2 km wide strip parallel to the Mediterranean coastline; 2. The Pliocene age corresponds to marine sediments deposited in the Dakhla syncline in the Northern part of the Korba City. The subsurface lithological details collected from some private wells show that vertical variations consist of yellow sands with alternating clay and sandstone layers; 3. Upper Miocene (the Somâa sands) laying on Medium Miocene (sandstones and marls of the Saouaf Formation), which constitutes the base of the aquifer. This formation contains thick fine sand layers of continental origin, including conglomeratic layers.
The main compressive structures that bend in the northeast-southwest direction ( Figure 2a) were generated during the Atlasique phase of the Upper Miocene [5,6,10,40]. These structures include: (1) Djebel Sidi Abderrahmene anticlinal, an elongated antiform parallel to the atlasic folds along the northeast-southwest direction; (2) Dakha synclinal, which contains the Korba coastal aquifer; (3) Menzel Temime anticlinal, which stretches along the Northeast-Southwest direction, beginning in the north of Nabel due to seismicity and ending at Korba. This anticlinal is influenced by inverse The region is one of the most productive agricultural areas in Tunisia, contributing about 15% of the country's total agricultural production [9]. Major crops grown in the study area are vegetables (potatoes, tomatoes, strawberries, etc.), vineyard plantations and citrus [8,11]. The region is characterized by a semi-arid climate (average annual rainfall: 420 mm) with Mediterranean influence. Average monthly high temperature occurs in July (37 • C) and the lowest average monthly low is 6.1 • C in January. There are ephemeral streams known as wadis, which originate from the Sidi Abderrahmen anticline and receive seasonal runoff. The most important wadis in the region (i.e., Chiba, M'laabi and Lebna) are obstructed by upstream dams. Consequently, groundwater has become the primary water resource in downstream areas.
The Tyrrhenian Quaternary made of arenitic limestones overlying conglomeratic units with a total thickness ranging 10-50 m, which forms a 1.2 km wide strip parallel to the Mediterranean coastline; 2.
The Pliocene age corresponds to marine sediments deposited in the Dakhla syncline in the Northern part of the Korba City. The subsurface lithological details collected from some private wells show that vertical variations consist of yellow sands with alternating clay and sandstone layers; 3.
Upper Miocene (the Somâa sands) laying on Medium Miocene (sandstones and marls of the Saouaf Formation), which constitutes the base of the aquifer. This formation contains thick fine sand layers of continental origin, including conglomeratic layers.  [2].
A distensible regime dominated the last phase of the Upper Miocene and a majority of Pliocene causing subsidence of the entire South-West region (Nabeul and Hammamet). This subsidence continued during the Upper Pliocene and Quaternary when the sinking trenches were formed. Normal faults developed during Tortonian synsedimentary in the Somâa Formation along the northwest-southeast and west/northwest-east/southeast directions. These lineaments dip 40-60° in the southeast direction.
The Korba unconfined aquifer system, assigned to the Plio-Quaternary unit, occupies the coastal plain but might also be locally semi-confined due to less permeable lenses of large lateral extent (Figure 2b). The saturated thickness of the aquifer varies from 3-10 m in the east toward the coastal area to a maximum of 50 m in the west near Somâa village with a Miocene marl Formation. The western part of the aquifer may contain sandstone lenses with different thicknesses and lateral variations, creating a third Miocene reservoir with considerable water resources [2,10,11]. The unconfined aquifer system receives infiltration from rainfall, that is, the main natural recharge source. The recharge is estimated to be more than 10% of the 420 mm/year average annual rainfall but also from lateral inflow from the surrounding relief, dunes and wadis during the rainy season [7,10,11]. Other important sources of recharge are irrigation return flow and lateral recharge by leakage (or through faults) from the underlying Miocene sandstone, which are not yet evaluated [2].
Since 1960s, groundwater exploitation in the Cap-Bon area has increased from 270 wells pumping 4 mm 3 in 1962 [10] to more than 13,000 wells pumping 68 mm 3 in 2014 [11]. There is a congestion of wells in the coastal part of the plain, which draws landward the freshwater-saltwater interface [8].  The main compressive structures that bend in the northeast-southwest direction ( Figure 2a) were generated during the Atlasique phase of the Upper Miocene [5,6,10,40]. These structures include: (1) Djebel Sidi Abderrahmene anticlinal, an elongated antiform parallel to the atlasic folds along the northeast-southwest direction; (2) Dakha synclinal, which contains the Korba coastal aquifer; (3) Menzel Temime anticlinal, which stretches along the Northeast-Southwest direction, beginning in the north of Nabel due to seismicity and ending at Korba. This anticlinal is influenced by inverse faults; and (4) Diare el Hajjej anticlinal and Taffeloun synclinal.
A distensible regime dominated the last phase of the Upper Miocene and a majority of Pliocene causing subsidence of the entire South-West region (Nabeul and Hammamet). This subsidence continued during the Upper Pliocene and Quaternary when the sinking trenches were formed. Normal faults developed during Tortonian synsedimentary in the Somâa Formation along the northwest-southeast and west/northwest-east/southeast directions. These lineaments dip 40-60 • in the southeast direction.
The Korba unconfined aquifer system, assigned to the Plio-Quaternary unit, occupies the coastal plain but might also be locally semi-confined due to less permeable lenses of large lateral extent ( Figure 2b). The saturated thickness of the aquifer varies from 3-10 m in the east toward the coastal area to a maximum of 50 m in the west near Somâa village with a Miocene marl Formation. The western part of the aquifer may contain sandstone lenses with different thicknesses and lateral variations, creating a third Miocene reservoir with considerable water resources [2,10,11]. The unconfined aquifer system receives infiltration from rainfall, that is, the main natural recharge source. The recharge is estimated to be more than 10% of the 420 mm/year average annual Water 2020, 12, 2525 5 of 27 rainfall but also from lateral inflow from the surrounding relief, dunes and wadis during the rainy season [7,10,11]. Other important sources of recharge are irrigation return flow and lateral recharge by leakage (or through faults) from the underlying Miocene sandstone, which are not yet evaluated [2].
Since 1960s, groundwater exploitation in the Cap-Bon area has increased from 270 wells pumping 4 mm 3 in 1962 [10] to more than 13,000 wells pumping 68 mm 3 in 2014 [11]. There is a congestion of wells in the coastal part of the plain, which draws landward the freshwater-saltwater interface [8]. The spatial distribution of the piezometric heads measured in 1962, 1977, 2014 and 2018 are presented in Figure 3. The piezometric maps of 1962 and 1977 (Figures 2b and 3a) show a smooth potentiometric surface located above sea level before the major perturbations of surface water regime. The growth of groundwater abstraction during the 1970s caused the hydraulic heads in the Korba aquifer to approach zero piezometric level. The groundwater table decline formed a concentric depression in the region of Diar El-Hojjaj, mainly due to intensive exploitation coupled with the reduction of recharge by damming the wadis (i.e., Chiba (1963) and M'laabi (1964)).
Water 2020, 12, x FOR PEER REVIEW 5 of 28 smooth potentiometric surface located above sea level before the major perturbations of surface water regime. The growth of groundwater abstraction during the 1970s caused the hydraulic heads in the Korba aquifer to approach zero piezometric level. The groundwater table decline formed a concentric depression in the region of Diar El-Hojjaj, mainly due to intensive exploitation coupled with the reduction of recharge by damming the wadis (i.e., Chiba (1963) and M'laabi (1964)). The piezometric map of 2014 ( Figure 3c) shows a groundwater cone of depression reaching 12 m below sea level observed 3 Km from the Mediterranean shoreline. In addition, a localized piezometric depression of up to 5 m below sea level is observed in the western part of the study area [8]. The latter reveals (Figure 3d) a general water table decline (at least 10 m) during the 56-year period, creating local cones of depression that increase in size from year to year. The greatest The piezometric map of 2014 ( Figure 3c) shows a groundwater cone of depression reaching 12 m below sea level observed 3 Km from the Mediterranean shoreline. In addition, a localized piezometric depression of up to 5 m below sea level is observed in the western part of the study area [8]. The latter reveals (Figure 3d) a general water table decline (at least 10 m) during the 56-year period, creating local cones of depression that increase in size from year to year. The greatest piezometric depression (up to 12 m below sea level) is observed in the western part, several kilometers from the Mediterranean shoreline in a region that is prone to vertical saltwater up-coning. Nearer to the shoreline, the piezometric level contours become increasingly negative, reversing the hydraulic gradients toward the central part of the aquifer, which accelerates saltwater intrusion from the Mediterranean Sea.

Material and Methods
We identified suitable sites for groundwater recharge through a knowledge-based analysis of eight factors, namely rainfall, lineament density, slope, drainage density, land use/land cover, lithology, geomorphology and soil type layers [41][42][43][44]. These factors were examined independently for groundwater recharge zoning to estimate recharge rates. Figure 4 illustrates the methodology used to delineate suitable sites for groundwater recharge.
Water 2020, 12, x FOR PEER REVIEW 6 of 28 piezometric depression (up to 12 m below sea level) is observed in the western part, several kilometers from the Mediterranean shoreline in a region that is prone to vertical saltwater up-coning. Nearer to the shoreline, the piezometric level contours become increasingly negative, reversing the hydraulic gradients toward the central part of the aquifer, which accelerates saltwater intrusion from the Mediterranean Sea.

Material and Methods
We identified suitable sites for groundwater recharge through a knowledge-based analysis of eight factors, namely rainfall, lineament density, slope, drainage density, land use/land cover, lithology, geomorphology and soil type layers [41][42][43][44]. These factors were examined independently for groundwater recharge zoning to estimate recharge rates. Figure 4 illustrates the methodology used to delineate suitable sites for groundwater recharge.

Preparation of Input Database
A range of field data and remotely sensed data were collected from various government agencies to map suitable artificial recharge areas. Table 1 summarizes the data types, sources and application. Primary data include digital elevation model (DEM) available from the Shuttle Radar Topography Mission (SRTM) and Sentinel-2A (S2A) multispectral imagery. Rainfall maps were generated by linear kriging interpolation method on a Surfer platform using a long-term average point rainfall for six pluviometry stations in the Cap-Bon over 37 years (1980-2017) [45].
The slope factor was calculated using the 30 m DEM based on the maximum rate of change in value from each cell to neighboring cells [46]. The data related to drainage density were generated indirectly from the slope data. Drainage density (Dd) is defined as the closeness of spacing of stream

Preparation of Input Database
A range of field data and remotely sensed data were collected from various government agencies to map suitable artificial recharge areas. Table 1 summarizes the data types, sources and application. Primary data include digital elevation model (DEM) available from the Shuttle Radar Topography Mission (SRTM) and Sentinel-2A (S2A) multispectral imagery. Rainfall maps were generated by linear kriging interpolation method on a Surfer platform using a long-term average point rainfall for six pluviometry stations in the Cap-Bon over 37 years (1980-2017) [45]. The slope factor was calculated using the 30 m DEM based on the maximum rate of change in value from each cell to neighboring cells [46]. The data related to drainage density were generated indirectly from the slope data. Drainage density (D d ) is defined as the closeness of spacing of stream channels [47] estimated as the measure of the total length of the stream segment of all orders per unit area using Equation (1) [48]. The drainage frequency value of the area ranges from 0 to >4 and the map indicates high moderate, medium and low drainage frequency.
where, ΣD i is the total length of all streams in stream order i (km) and A is the watershed area (km 2 ). The lineament density of the Korba region was prepared from a combination of Sentinel-2A (S2A) multispectral imagery, DEM and black and white photographs using computer-aided supervised extraction. To achieve optimum interpretation with reference to frequency, the "sliding window" method was applied and the number of lineaments in each unit area or pixel was counted [49,50]. The thematic layer for lineament density (L d ) can be defined as the total length of all recorded lineaments divided by the area of the catchment under consideration as follows [51]: where, L i is the length of the ith lineament, ΣL i is the total length of all lineaments in km and A is the area of the grid (km 2 ). Lineament raster was created using false color composite in spatial analysis tool in ArcGIS 10.1 followed by additional editing of the boundary divide line and road features to ensure the quality of the extracted lineaments. The frequency value ranges from 0 to over 2 km/km 2 , providing a basis for evaluating the lineament frequency and their distribution. A lineament frequency contour map was prepared along with a synoptic classification indicating very low, low, moderate, high and very high frequency value.
The spatial soil groups and geomorphological map were derived from the soil map of La Carte Agricole of Nabeul governorate at 1:50,000 for the year 2005, which is an official spatial data source. The vector layer was converted to a grid and then reclassified according to recharge weighted rating [51,52]. The land-use/land cover (LU/LC) raster (30 m resolution) was obtained by processing satellite images from the Landsat 8 Sensor on a scale of 1/500,000 by ENVI software. The lithological map for the Korba area was obtained using a subsurface geology map and several well drilling profiles of the Cap-Bon region [53]. The Korba area is mainly covered by six types of lithological units that were ranked based on available porosity and permeability values of aquifer rocks. Sedimentary features, such as sand and sandstone, generally demonstrate very high porosity (20-31%) and permeability (500-7000 mD) and they are capable of accumulating and storing large volumes of groundwater [54][55][56].

Analytic Hierarchy Process (AHP)
The AHP [57,58] facilitates decision making by organizing perceptions and judgments into a multi-level hierarchical structure that accounts for the forces that control a decision through assessing multiple factors [59]. The AHP is used as a complex multi-criteria decision-solving method to determine the weights assigned to different thematic layers and their respective features. This technique breaks down the multi-criteria decision problem into a hierarchy based on a pair-wise comparison of the importance of different criteria and sub-criteria within the judgment matrix [60][61][62]. The hierarchy allows analysis to consider each of several properties separately [54] for the selection of potential recharge zones from competing sets of parameters. The method was implemented in four steps to delineate potential recharge zones in the Korba study area: (1) selecting factors that influence groundwater recharge zones (2) developing pairwise comparison matrix, (3) estimating relative weights and (4) assessing matrix consistency.

Selection of Factors Influencing Groundwater Recharge Zones
In the first step of the AHP, each factor that influences recharge was given a score between 1 and 9, depending on its significance compared to other factors in pairwise comparisons [58]. For this, a standard Saaty's 1-9 scale was used ( Table 2) to describe the relative influence of parameters, where score 1 denotes equal influence of parameters and score 9 denotes extreme influence of a parameter on groundwater recharge compared to the other parameters [58]. ranked based on available porosity and permeability values of aquifer rocks. Sedimentary features, such as sand and sandstone, generally demonstrate very high porosity (20-31%) and permeability (500-7000 mD) and they are capable of accumulating and storing large volumes of groundwater [54][55][56].

Analytic Hierarchy Process (AHP)
The AHP [57,58] facilitates decision making by organizing perceptions and judgments into a multi-level hierarchical structure that accounts for the forces that control a decision through assessing multiple factors [59]. The AHP is used as a complex multi-criteria decision-solving method to determine the weights assigned to different thematic layers and their respective features. This technique breaks down the multi-criteria decision problem into a hierarchy based on a pair-wise comparison of the importance of different criteria and sub-criteria within the judgment matrix [60][61][62]. The hierarchy allows analysis to consider each of several properties separately [54] for the selection of potential recharge zones from competing sets of parameters. The method was implemented in four steps to delineate potential recharge zones in the Korba study area: (1) selecting factors that influence groundwater recharge zones (2) developing pairwise comparison matrix, (3) estimating relative weights and (4) assessing matrix consistency.

Selection of Factors Influencing Groundwater Recharge Zones
In the first step of the AHP, each factor that influences recharge was given a score between 1 and 9, depending on its significance compared to other factors in pairwise comparisons [58]. For this, a standard Saaty's 1-9 scale was used ( Table 2) to describe the relative influence of parameters, where score 1 denotes equal influence of parameters and score 9 denotes extreme influence of a parameter on groundwater recharge compared to the other parameters [58].

Pairwise Comparison Matrix
The AHP method integrates and transforms spatial data (input) into decision (output), where qualitative information of individual thematic layers and features are converted into quantitative scores based on Saaty's scale. Then, a pairwise comparison matrix (PCM) [19,26,29] is constructed (Equation (3)) using Saaty's scores obtained in the previous step. In the PCM, the matrix column is constructed based on a descending order of parameter influence on recharge. The first element is assigned a score of 1 when compared to itself (see Table 3). Other elements of the rows are filled using the actual Saaty's scores when a more influential parameter is compared with a less influential parameter or the reciprocal of the Saaty's scores score when a less influential parameter is compared to a more influential parameter.
where A is a pairwise comparison matrix where element Xnn denotes relative significance of a parameter for recharge compared to another parameter. Table 3 provides the PCM for the parameters examined in this study. Lithology was selected as the first parameter of the matrix because has a higher influence on recharge potential compared to

More important
Water 2020, 12, x FOR PEER REVIEW 8 of 28 ranked based on available porosity and permeability values of aquifer rocks. Sedimentary features, such as sand and sandstone, generally demonstrate very high porosity (20-31%) and permeability (500-7000 mD) and they are capable of accumulating and storing large volumes of groundwater [54][55][56].

Analytic Hierarchy Process (AHP)
The AHP [57,58] facilitates decision making by organizing perceptions and judgments into a multi-level hierarchical structure that accounts for the forces that control a decision through assessing multiple factors [59]. The AHP is used as a complex multi-criteria decision-solving method to determine the weights assigned to different thematic layers and their respective features. This technique breaks down the multi-criteria decision problem into a hierarchy based on a pair-wise comparison of the importance of different criteria and sub-criteria within the judgment matrix [60][61][62]. The hierarchy allows analysis to consider each of several properties separately [54] for the selection of potential recharge zones from competing sets of parameters. The method was implemented in four steps to delineate potential recharge zones in the Korba study area: (1) selecting factors that influence groundwater recharge zones (2) developing pairwise comparison matrix, (3) estimating relative weights and (4) assessing matrix consistency.

Selection of Factors Influencing Groundwater Recharge Zones
In the first step of the AHP, each factor that influences recharge was given a score between 1 and 9, depending on its significance compared to other factors in pairwise comparisons [58]. For this, a standard Saaty's 1-9 scale was used ( Table 2) to describe the relative influence of parameters, where score 1 denotes equal influence of parameters and score 9 denotes extreme influence of a parameter on groundwater recharge compared to the other parameters [58]. The AHP method integrates and transforms spatial data (input) into decision (output), where qualitative information of individual thematic layers and features are converted into quantitative scores based on Saaty's scale. Then, a pairwise comparison matrix (PCM) [19,26,29] is constructed (Equation (3)) using Saaty's scores obtained in the previous step. In the PCM, the matrix column is constructed based on a descending order of parameter influence on recharge. The first element is assigned a score of 1 when compared to itself (see Table 3). Other elements of the rows are filled using the actual Saaty's scores when a more influential parameter is compared with a less influential parameter or the reciprocal of the Saaty's scores score when a less influential parameter is compared to a more influential parameter.
where A is a pairwise comparison matrix where element Xnn denotes relative significance of a parameter for recharge compared to another parameter. Table 3 provides the PCM for the parameters examined in this study. Lithology was selected as the first parameter of the matrix because has a higher influence on recharge potential compared to

Pairwise Comparison Matrix
The AHP method integrates and transforms spatial data (input) into decision (output), where qualitative information of individual thematic layers and features are converted into quantitative scores based on Saaty's scale. Then, a pairwise comparison matrix (PCM) [19,26,29] is constructed (Equation (3)) using Saaty's scores obtained in the previous step. In the PCM, the matrix column is constructed based on a descending order of parameter influence on recharge. The first element is assigned a score of 1 when compared to itself (see Table 3). Other elements of the rows are filled using the actual Saaty's scores when a more influential parameter is compared with a less influential parameter or the reciprocal of the Saaty's scores score when a less influential parameter is compared to a more influential parameter.
where A is a pairwise comparison matrix where element X nn denotes relative significance of a parameter for recharge compared to another parameter. Table 3 provides the PCM for the parameters examined in this study. Lithology was selected as the first parameter of the matrix because has a higher influence on recharge potential compared to the other factors. Thus, lithology was assigned the value 8. Land use/land cover was selected as the second most important parameter influencing recharge followed by slope, geomorphology, lineaments, rainfall, drainage and finally soil parameter in a descending order of influence. Each parameter in the selected set was assigned a Saaty's score based on its influence on recharge relative to lithology.

Estimating Relative Weights
Weights were assigned to the variables based on 'expert' opinion to estimate the relative importance of variables compared to other variables and to quantify the relative influence of each variable on recharge [26,60,61]. The layers were assigned weights derived by normalizing the pair comparison matrix (NPCM) [62,63]. The NPCM elements were computed by dividing thematic element values by their corresponding total column values from the PCM (Equation (4)) (see Table 4): where X ij is normalized pair-wise matrix value at i th row and j th column, C ij is the value assigned to each criteria at i th row and j th column and L ij is the total values in each column of the pair-wise matrix. Subsequently, a standard weight was calculated for variable i by dividing each normalized pairwise matrix elements by criterion number (N) (Equation (5)) [64,65].
where W i is standard weight. Then, Eigen vector and eigenvalue calculations help determine the percentage of effect of the thematic layers and the classification of the constraints ( Table 5). The eigenvector was calculated by dividing column elements by the column sum in Table 4. The principal Eigen vector was obtained by averaging across the rows to quantify relative weights of each parameter [57,59,66]. A consistency vector was obtained by multiplying two different matrix values from selected thematic layers (Equation (6)), namely, pair-wise comparison matrix value and normalized pair-wise matrix value [3,4,67].
where, λ is the consistency vector. The sum of eigenvalues called principal eigenvalue (λ max ) is a measure of matrix deviation from consistency [66]. According to Saaty [57], for a pairwise comparison matrix to be consistent it must have a principal eigenvalue (λ max ) greater than or equal to the number of the parameters considered (n). The principal eigenvalue of 8.62 was obtained for the 8 × 8 matrix (see Table 5), which was used for the calculation of consistency index.

Assessing Matrix Consistency
Consistency is checked based on Consistency Index (CI) and Consistency Ratio (CR) [48,50,53,55,56]: where, CI: consistency index, λ max : highest matrix Eigenvalue, n: number of variables (thematic layers), CR: consistency ratio and RI: Random Index value based on the number of variables.
A perfectly consistent decision maker should always yield CI = 0 but small values of inconsistency may be tolerated if the CI < 0.1 [57,68]. We obtained an acceptable CI value of 0.08. Also, if the CR is greater than 0.1, the pairwise comparison judgments must be re-evaluated. With a matrix of eight variables, the RI is 1.41 (see Table 6). The applied weighting yielded a CR of 0.05, which shows that the weights (Table 5) assigned to GIS thematic layers of parameters are consistent [58][59][60].

Multi-Influencing Factors (MIF)
MIF is a widely used MCDM technique for environmental management [12,34,36]. The method is useful for mapping groundwater recharge sites by assigning an appropriate weight to various factors and feature classes based on the influence on groundwater flow and storage ( Figure 5).
The ranks and weights are assigned to each factor and different classes based on their relative contribution to groundwater recharge. We specified major and minor influential relationships among the factors ( Figure 5). Each relationship was weighted according to its influence on recharge. A large weight indicates that the factor has a large influence on groundwater recharge. The representative weight of a factor influencing recharge potential is the sum of all weights from each factor. Major effects were assigned a weight of 1.0 whereas minor influential relationships were given a weight of 0.5 [55][56][57]. Furthermore, the classes with no effect on groundwater recharge were assigned a weight of zero. The weights for each factor were summed to obtain the recharge potential factor weight [34,35]. The estimated weight for each influencing factor was obtained as a percentage using Equation (9): (Major e f f ect + Minor e f f ect) (Major e f f ect + Minor e f f ect) × 100.
The results are shown in Table 7.
Water 2020, 12, x FOR PEER REVIEW 11 of 28

Multi-Influencing Factors (MIF)
MIF is a widely used MCDM technique for environmental management [12,34,36]. The method is useful for mapping groundwater recharge sites by assigning an appropriate weight to various factors and feature classes based on the influence on groundwater flow and storage ( Figure 5). The ranks and weights are assigned to each factor and different classes based on their relative contribution to groundwater recharge. We specified major and minor influential relationships among the factors ( Figure 5). Each relationship was weighted according to its influence on recharge. A large weight indicates that the factor has a large influence on groundwater recharge. The representative weight of a factor influencing recharge potential is the sum of all weights from each factor. Major effects were assigned a weight of 1.0 whereas minor influential relationships were given a weight of 0.5 [55][56][57]. Furthermore, the classes with no effect on groundwater recharge were assigned a weight of zero. The weights for each factor were summed to obtain the recharge potential factor weight [34,35]. The estimated weight for each influencing factor was obtained as a percentage using Equation The results are shown in Table 7. Total ∑18.5 ∑100 Figure 6 illustrates the GIS analysis procedure to map the spatial distribution of groundwater recharge using the weighted overlay tool of Weighted Overlay analysis tool of ArcGIS 10.1 based on the eight thematic layers and their corresponding percentage influence on recharge [69,70]. Values in the input raster layers were reclassified into a common evaluation scale of 1 (very good), 2 (good), 3   Figure 6 illustrates the GIS analysis procedure to map the spatial distribution of groundwater recharge using the weighted overlay tool of Weighted Overlay analysis tool of ArcGIS 10.1 based on the eight thematic layers and their corresponding percentage influence on recharge [69,70]. Values in the input raster layers were reclassified into a common evaluation scale of 1 (very good), 2 (good), 3 (moderate), 4 (poor) and 5 (very poor). This was done by multiplying the cell values of each factor class by the factor weight and summing the resulting cell values to produce a map of potential recharge zones [71,72], as summarized in Equation (10): where GWRZ is the groundwater recharge zone, W i is the weight of each thematic layer, R i is the rating of each class of each thematic layer, LI is the lithology of the aquifer, LU/LC is the land-use/land-cover, SLO is the slope, GM is the geomorphology, LD is the lineament, RF is the rainfall, DD is the drainage density and SL is the soil cover. The subscripts c and w refer, respectively, to the factor class of a thematic layer and its percent influence on recharge [46,56]. It means that each layer class rank (each individual class has a rank) was multiplied by the weight from Saaty's assumption of each layer (each layer has a unique weight) to obtain the position of layers in the overlay analysis to map GWRZ.
LD LD , + RF RF + DD DD + SL SL , where GWRZ is the groundwater recharge zone, Wi is the weight of each thematic layer, Ri is the rating of each class of each thematic layer, LI is the lithology of the aquifer, LU/LC is the landuse/land-cover, SLO is the slope, GM is the geomorphology, LD is the lineament, RF is the rainfall, DD is the drainage density and SL is the soil cover. The subscripts c and w refer, respectively, to the factor class of a thematic layer and its percent influence on recharge [46,56]. It means that each layer class rank (each individual class has a rank) was multiplied by the weight from Saaty's assumption of each layer (each layer has a unique weight) to obtain the position of layers in the overlay analysis to map GWRZ.

Figure 6.
A schematic of thematic layer overlay analysis using geographic information system (GIS) to demarcate groundwater recharge potential zones.

Sensitivity Analysis
The GWRZ index is sensitive to the factor weight and the rating value for each class of each layer [73,74]. Sensitivity analysis reveals how much each thematic layer and its weight influence the output map. A map removal sensitivity analysis (MRSA) and a single parameter sensitivity analysis (SPSA) were conducted to examine the effect of each of the eight factors' weights on the final GWRZ.

Map Removal Sensitivity Analysis (MRSA)
MRSA was conducted to determine the most or the least influential thematic layer for groundwater recharge zones [30]. The sensitivity index, S, was computed using Equation (11): where GWRZ and GWRZ' are the output of groundwater recharge map index of all the thematic layers and when one of the thematic layers is removed, respectively. N is the number of the full thematic layers used to compute the GWRZ and n is the number of thematic layers used to compute GWRZ'.

Figure 6.
A schematic of thematic layer overlay analysis using geographic information system (GIS) to demarcate groundwater recharge potential zones.

Sensitivity Analysis
The GWRZ index is sensitive to the factor weight and the rating value for each class of each layer [73,74]. Sensitivity analysis reveals how much each thematic layer and its weight influence the output map. A map removal sensitivity analysis (MRSA) and a single parameter sensitivity analysis (SPSA) were conducted to examine the effect of each of the eight factors' weights on the final GWRZ.

Map Removal Sensitivity Analysis (MRSA)
MRSA was conducted to determine the most or the least influential thematic layer for groundwater recharge zones [30]. The sensitivity index, S, was computed using Equation (11): where GWRZ and GWRZ are the output of groundwater recharge map index of all the thematic layers and when one of the thematic layers is removed, respectively. N is the number of the full thematic layers used to compute the GWRZ and n is the number of thematic layers used to compute GWRZ'.

Single Parameter Sensitivity Analysis (SPSA)
The single-parameter sensitivity analysis was performed to evaluate the effect of each factor in the groundwater recharge potential index. This sensitivity test was conducted because the numerical scores assigned to the parameters, which are the basis of parameter weightings, are essentially arbitrary [11]. The initial arbitrary weights are replaced with effective weights [22,75] and the sensitivity index was calculated according to the Equation (12): where; W i is the weight for each layer, R i is the rating value for each class of each layer, GWRZ is the groundwater recharge zone.

Verification/Validation of Groundwater Recharge Zone
An important step for any model is validation to assess the performance by establishing the relationship between groundwater recharge level class and GWRZ map [76][77][78][79]. Several methods can be used to validate groundwater recharge maps, such as receiver performance analysis, curve area, groundwater yield estimation of wells during field visits and a comparative study between the water level profile and groundwater recharge zones, among others [78,[80][81][82]. In this study, two methods were used to validate the obtained groundwater recharge zones, namely: Receiver Operating Characteristics (ROC) curve and groundwater level monitoring and assessment at 25 shallow wells. The groundwater yield data were collected during several field investigations or obtained from local management agencies [83].
The ROC curve was applied as a mathematical technique to determine the map accuracy [78,80,81]. Using SPSS-statistical software, ROC curves were obtained by considering cumulative percentage of area (on the x axis) and the cumulative percentage of number of wells (on the y axis) [84]. The area under the curve (AUC) was calculated to quantify the ability to predict the occurrence or non-occurrence of pre-defined "events" [83][84][85][86]. The resultant AUC values range from 0.5 (i.e., a random prediction) to 1 (i.e., an excellent prediction) [81]. The quantitative-qualitative relationship between the AUC value and prediction accuracy can be excellent (0.9-1), very good (0.8-0.9), good (0.7-0.8), average (0.6-0.7) and poor (0.5-0.6) [85]. To illustrate the extent to which the occurrence or non-occurrence of pre-defined events was correctly predicted, the AUC method was verified through frequency percentage for AHP and MIF techniques. A comparative study was carried out between the water level profile and groundwater recharge zones through two water level profiles for two sections in the study area AA' and BB' (see Figure 8), for AHP and MIF methods, respectively. Table 8 reports the final weights of the factors governing GWRZs, which are described below. Lithology (LI): LI governs the porosity and permeability of aquifer rocks [18], which in turn influences the occurrence and distribution of groundwater recharge through physio-mechanical properties that control the ability of the aquifer materials to transmit water and the rate at which groundwater flows [27]. The lithology of the study area consists of quaternary alluvium (12% of the surface area; 52.56 km 2 ), sandy clay (22%; 96.36 km 2 ), yellow sand with clay (6%; 26.28 km 2 ), tyrrhenian sandstones (3%; 13.14 km 2 ), sandy sandstone (3%; 13.14 km 2 ) and silty sand (54%; 236.52 km 2 ) with ratings varying from 1 to 6 ( Table 8 and Figure 7a). The quaternary alluvial sediments are identified along the wadis across the study area. The lithology is predominantly silty sand of Dakhla syncline and sandy clay with sandstone intercalations of Pliocene Formation. Laterally, the facies changes southward to more differentiated yellow sand with clay, sand and sandstone layers. The higher ratings (6 and 5), apply to 60% of the area and are assigned to the sediments of sandstone, gravel and sand located at the "Sable of the Somâa Formation" and in the central part of the study area around Diarr El Hojjaj and Menzel Horr villages. These delineated zones are characterized by high permeability index, which implies that rainwater can easily infiltrate [11], providing high hydraulic conductivities (e.g., ranging from 10 −6 to 10 −3 m/s) while hydraulic conductivities at the Somâa Formation are smaller [10]. The Tyrrhenian sandstones form the most permeable unit, followed by the central part of the Pliocene (alluvial deposits), the Pliocene, the Somâa sand and finally, the early Miocene Formations [2].   Land use/Land cover (LU/LC): LU/LC has a major influence on the occurrence and development of groundwater in a terrain [26,87]. Agricultural lands typically allow more infiltration due to pore spaces in the soil, which trap and hold the water in the roots, providing a pathway for water to percolate into the surface by loosening up the rock and soil. By contrast, built-up and barren lands reduce infiltration due to loss of permeable surface and increased runoff potential. Therefore, the areas with agriculture and water bodies are considered as good sites for groundwater recharge, while settlements and barren lands have poor groundwater recharge potential [73]. A major portion of the study area is under agriculture (69.42%; 304.06 km 2 ) followed by barren land (16.58%; 72.62 km 2 ), rural/urban settlements (8.65%; 37.89 km 2 ) and river (5.35%; 24.43 km 2 ) (Figure 7b). Slope (SLO): SLO is a characteristic of local and regional relief, which is an important factor because it influences the water retention, intensity of infiltration [49,52], aquifer recharge and groundwater movement. The higher the ground slope, the smaller the infiltration rate will be due to larger runoff potential [12,27,58]. The study area has five slope classes (Figure 7c), namely very low (0-5%), low (5-10%), moderate (10-15%), high (15-20%) and very high (>20%). The topography layer displayed a gentle slope (0-5%) over most of the study area (87%; 381.06 km 2 ) which has been assigned a weight score of 6 with Very Good potential for artificial recharge. The slope increases from east to west due to the presence of the Abderrahman Mountain range. Slopes (5-10%) are considered Good for groundwater storage. The third category (i.e., Poor) includes two slope ranges of 10-15% and 15-20%. Gently sloping Very Poor recharge areas cover 8.76 km 2 (2%) of the study area, while the area under steep slope is negligible. A higher priority was given to a nearly level to gentle slope categories followed by moderate and steep slopes (Table 8).

Evaluation of Factors Governing Groundwater Recharge Zones
Geomorphology (GM): GM of an area reflects various landforms, which includes their description, species and physical processes that help assess groundwater recharge potential and evaluate possible groundwater areas [20,23,25,50]. Geomorphic features of the study area were categorized into five units namely: shallow flood plain and beach, shallow buried pediplain, moderately buried pedipalin, deep buried pediplain and sedimentary high ground (Figure 7d).
Moderately buried pedipalin covered 9.01% (39.5 km 2 ) of the study area. It comprises gently undulating plains covered with sand layers of continental origin including conglomeratic layers and clay lenses that are moderately favorable for groundwater recharge. Deep buried pediplain and sedimentary high ground cover 4.32% (18.92 km 2 ) and 7.59% (33.24 km 2 ) of the study area, respectively and are characterized by highly sloping topography, smaller amounts of infiltration and high surface runoff characteristic of a poor recharge zone. Lineament (LD): LD density networks such as fractures, joints and faults increase porosity and therefore play an important role in groundwater movement and high groundwater recharge potential [26,27,33]. The LD density, extracted from satellite images and remote sensing, is an indicator of the degree of fracturing of the rock, meaning that groundwater recharge potential is high near lineaments. While high lineaments frequency indicates very high recharge potential due to the presence of recharge pathways, low frequency does not necessarily translate into very low recharge potential. In this study, the lineament density was subdivided and ranked into five classes: Lineament with densities ranging from 0.0-5.0 km/km 2 dominated the study area (57.22%; 250.62 km 2 ). Lineament density ranging from 0.5-1 km/km 2 and higher than 1 km/km 2 covers an area of 45.10 km 2 (32.09%) and 33.78 km 2 (10.69%), respectively. The rose diagram shows the dominance of NE-SW and NW-SE directions of lineaments (Figure 7e). These conjugate lineament directions are characteristic of Tunisian Atlas domain and they increase reservoir permeability through creating interconnected lineaments/fractures. Rainfall (RN): RN is typically the primary source of groundwater recharge in which the water infiltrates into subsurface through fractures and soils. It governs the amount of runoff that would be available to capture in recharge basins to increase infiltration. The mean annual rainfall ranges from 400 to 450 mm in the lowland areas along the coastline, while it ranges from 450 and 500 mm in the mountainous areas. An increase of annual rainfall with altitude is detected. The resulting rainfall map was classified into two major classes (Table 8 and Figure 7f): 450-500 mm/yr and 400-450 mm/yr. Generally, 75% of the 420 mm/yr average annual rainfall occurs during the wet season (from September to March), while summers are usually dry. From the map of rainfall and piezometric levels (see Figure 3), higher-altitude areas have greater recharge potential than lower altitudes. In this study, the rainfall factor was assigned a weight of 0.8 and 0.5 for MIF and AHP, respectively, in the total groundwater recharge potential.
Drainage density (DD): DD is another factor that affects the movement of water and groundwater recharge [44,48,55]. It is a measure of how a watershed is drained by stream channels. DD is influenced by numerous factors, including resistance to erosion, infiltration capacity, vegetation cover, surface roughness and runoff intensity index and climatic conditions [16,18,19]. Areas with high drainage density have less recharge rate, whereas low drainage density areas have a high recharge rate and can directly influence the groundwater recharge. Five drainage density categories were identified in the study area (Figure 7g), namely 'very good' (3-4 km/km 2 or more), 'good' (2-3 km/km 2 ), 'moderate' (1-2 km/km 2 ) and 'poor' (0.1-1 km/km 2 ). High drainage was recorded in the high relief near to Somâa village in the southwestern part. A higher priority was given to a low drainage density category followed by medium and high drainage density categories.
Soil (SL): soil texture controls the percolation and infiltration of surface water into the aquifers, influencing groundwater recharge through properties such as porosity, structure, adhesion and consistency [26,44]. Based on soil texture data, the study area has three main soil types, namely gravel and sand, loam-sandy loam and clay-clay loam (Figure 7h). Sandy loam is the dominant soil textural class in the study area (69.72%; 305.37 km 2 ). Gravel and sandy soils cover a major portion of the northern and central parts (27.18%; 119.05 km 2 ). Nearly 3% of the clay loam soil partially covers north eastern and southwestern portions of the study area. The ranks of soils were assigned according to their degree of infiltration [21,22,46]. The gravel and sandy soil has a high degree of infiltration and, therefore, has a higher influence, while the clay and clay loam soil has the lowest degree of infiltration and, therefore, is given the lowest influence (Table 8).

Analytical Hierarchy Process
The identification of GWRZ was accomplished based on the rates and weights of the eight thematic layers according to Equation (13): Results of the GWRZ mapping are shown in Figures 8a and 9. The Korba area has 70 km 2 of very good groundwater recharge zones (about 16% of the area). The GWRZ map indicates the high groundwater recharge potential zone situated mostly near the shallow flood pediplain zone. This is mainly attributed to the availability of the high quaternary alluvial sediment in plain terrain along with the rivers course where there are high to very high lineament density (1->2 km/km 2 ) with gravel and sand soil type. The GWRZ designated as good is mostly found in the north and the central part of the Korba area, which covers 53% (123 km 2 ). The good potential zone is distributed across irrigated agricultural lands (e.g., vineyards, citrus, strawberries, potatoes, tomatoes, etc.) and grassland on gentle slopes (0-5%) and flat buried pediplains with sandy loam soil type. The moderate zone covers 11.5% (50.37 km 2 ) of the area, scattered along the north-west direction where the lineament density is high (>2 km/km 2 ). The poor 15.8% (69.2 km 2 ) and very poor 2.9% (0.37 km 2 ) groundwater recharge areas are found especially in the south-eastern portion of the study area. This may be due to the fact that the area exhibits high sedimentary ground of mountainous geomorphology, relatively less water bearing lithological Formations (Somâa sand with clay) and steep topography (>10% slope category).

Multi-Influencing Factors Techniques (MIF)
The GWRZ index using the MIF technique was calculated using Equation (14): Very high GWRZ (134.47 km 2 or 30.7%) is identified especially in the center of the study area, along the Mediterranean coastline and along the Chiba and Lebna river channels (Figures 8b and 9). The geomorphology of this region consists of shallow pediplain with shallow quaternary alluvium deposit, which allows water to infiltrate. However, the river channels increase the drainage density, routing the surface runoff downstream, which reduces groundwater recharge. Good GWRZ is about 219 km 2 (50%) covering the north, central and southwestern parts. The lithological Formations with the highest groundwater recharge potential are the alluvial deposits and silt-sand. Furthermore, the good GWRZ is attributed to flat buried pediplains with sandy loam soil type and the relatively high amounts of rainfall on low slopes. Moderate GWRZ (42.05 km 2 or 9.6%) is mainly concentrated in the western part and scattered patches throughout the study area. The low potential zones (poor and very poor) are identified in the north-western and southwestern parts, respectively, due to higher slope and unfavorable geology (Miocene age deposition) and geomorphological condition (high ground) where hard rock terrains are present. The areas covered by poor and very poor GWRZ are about 35.04 Km 2 (8.07%) and 6.92 km 2 (1.58%), respectively.

Map Removal Sensitivity Analysis (MRSA)
In Map Removal Sensitivity Analysis (MRSA) for AHP and MIF methods, the thematic layers were removed from the mapping process one at a time, using the remaining layers to map GWRZ. The computed sensitivity index (S) values are reported in Table 9, indicating that the most sensitive layers for GWRZ map are LI, with the mean variations index (S MVI ) of 1.62% and 1.54%, LU/LC, with S MVI of 1.24% and 1.26% and SLO with S MVI of 1.04% and 1.01%, for AHP and MIF, respectively. This is essentially due to the relatively high theoretical weight assigned to these layers (Table 8). In addition to the relatively high weight of LI (33% for AHP and 22% for MIF), this layer has high rating values in almost every sub-area and it is thus considered the most influential thematic layer for mapping the GWRZ in the Korba area. Likewise, the removal of LU/LC (AHP weight: 23%; MIF weight: 19%) and SLO (AHP weight: 16%; MIF weight: 13%) has a great impact on the GWRZ map.
The fourth and fifth layers are GM (S MVI for AHP: 0.74%; S MVI for MIF: 0.71%) and LD (S MVI for AHP: 0.64%; S MVI for MIF: 0.69%). The least influential factors for the output map determination were the RN thematic map (S MVI = 0.48% and 0.29% for AHP and MIF, respectively) and DD thematic map (S MVI = 0.31% and 0.29% for AHP and MIF, respectively).
The lowest MRSA values, in AHP and MIF, were 0.27% and 0.21%, respectively, for the SL layer. Thus, the sensitivity order of the thematic layers for the AHP method is LI > LU/LC > SLO > GM > LD > RN > DD > SL, which is consistent with the magnitudes of the theoretical weights in the AHP method. However, the sensitivity order of variables is slightly different in MIF technique (LI > LU/LC > SLO > GM > LD > RN = DD > SL). The variation of the sensitivity index is governed by the collective impact of the large data range, the mean rating score and especially the weight assigned to the thematic layer. Overall, the MRSA indicates that almost all the factors are necessary for GWRZ mapping in the Korba aquifer (Table 9).

Single Parameter Sensitivity Analysis (SPSA)
Single Parameter Sensitivity Analysis (SPSA) analysis confirms that most of the study area has good to very good groundwater recharge potential and all the parameters have a significant impact on recharge. The analysis was performed to better understand the contribution of each factor to the groundwater recharge map [88]. In this section, the "effective" weight of each factor was compared with the theoretical weight assigned to it by the AHP and MIF methods (Table 10) in order to determine the "most effective impact" factors. As the factor's effective weight increases, the higher its sensitivity compared to other factors with a lower effective weight. For AHP, groundwater recharge tends to be most sensitive to LU/LC, LI, SLO and GM layers (Table 10) because their mean effective weights, 86%, 66%, 44% and 39%, respectively, are higher than their corresponding theoretical weights, 23%, 33%, 16% and 11% (see Table 7). On the other hand, LD, RF, DD and SL are the least effective, with a mean effective weight of 25%, 25%, 10% and 8% respectively, which are still higher than their corresponding theoretical weights 7%, 5%, 3% and 2%. For MIF, too, the results (Table 10), indicate that LU/LC is the most effective factor, with a mean effective weight of 71%, compared to its theoretical weight (19%). The GM (46%), LI (44%), RN (40%), LD (39%) and SLO (36%) also show a higher effective weight compared to their theoretical weights, 13%, 22%, 8%, 11% and 13%, respectively. The rest of the Factors such as DD (28%) and SL (24%) exhibit less effective weights than others.

Verification/Validation of Groundwater Recharge Zones
The hydrogeological data in combination with the hydrological data provide new information to advance understanding of key hydrological processes controlling groundwater flow and recharge in the shallow Korba aquifer. Lithology is the most important component in determining groundwater recharge potential in the study area due to the nature of the geological formations and how they influence infiltration rates in the surficial materials. Pumping rates are expected to be correlated with parameters such as transmissivity (high transmissivity allows high discharge) or depth of the groundwater table (shallow depth facilitates groundwater abstraction). Ten years of with 10 years of measured well discharge data from 25 monitoring wells were used to cross-validate the GWRZs. Locations of the monitoring wells were overlain on the map of groundwater recharge potential to evaluate the results from the AHP and MIF techniques (Figure 8a,b). The maximum and minimum depths of the monitoring wells are 28 m and 5 m, respectively and the average depth is 9 m. The evaluation approach is deemed adequate given the generally shallow depths of groundwater wells.
The accuracy of the AHP method to map groundwater recharge potential in the Korba area was about 73% whereas the MIF method's accuracy was 67%, indicating good agreement between the two approaches. In relation to topography, the good and very good groundwater recharge zones coincide with shallow and moderate water level depth locations, which can be explained by the low thickness of the unsaturated zone (ranging from 3.5 to 10 m) composed of thick fine sand layers of continental origin. Areas with small depth to groundwater table are located where there is high infiltration. The good recharge zone is mostly found in wet regions located in the central part of the study area with intensive irrigated cultivations on highly permeable soil where the number of wells is proportional to the number of farms. These nearly flat alluvial areas with large hydraulic conductivities (ranging from 10 −6 to 10 −3 m/s) and flood pediplain have good groundwater recharge potential. The poor and very poor recharge zones are found where water level is deep to very deep ( Figure 8). These areas exhibit diverse geomorphological variations (numerous structural hills) where the western most areas are characterized by very steep escarpment and low soil porosity, resulting in low recharge potential in the upland plain.
Eight very high discharge wells (20-30 m 3 /h) were located in the very good recharge zones obtained from both AHP and MIF methods. These wellfields are located primarily in the sandstone and quaternary alluvium in the central part of the study area, which has the highest transmissivity zones of about 10 −2 m 2 /s. Three of five high discharge (10-20 m 3 /h) wells fell within good recharge zones and two were in moderate zones. Further, three medium discharge (5-10 m 3 /h) wells were located in the moderate recharge zone except one that fell within the good zone according to the MIF generated map. For AHP method, five of six low discharge (1-5 m 3 /h) wells exist in the poor zone and one well in the moderate zone. For MIF technique, two of the six low discharge wells occur in the moderate recharge zone. All three wells with very low discharge were found in the very poor zone except one well in moderate zone. Finally, the three very low discharge wells (<1 m 3 /h) fell within very poor recharge zones except one that was located in the moderate zone based on the MIF method.
The groundwater recharge maps were also validated based on AUC. The area under the ROC curve (AUC) is an indicator of model quality. Generally, a good model has an AUC value of 0.7-0.9, while an excellent model has values over 0.9. According to Figure 10, the AHP method outperformed MIF (AHP AUC: 75.6%; MIF AUC: 70.4%), although both methods were reasonably accurate for mapping groundwater recharge in the Korba area. Overall, the validation results increase confidence in the applied methodology as a useful framework for rapid assessment of groundwater recharge to inform siting artificial recharge structures and other groundwater management efforts.

Groundwater Recharge Estimation
Quantifying groundwater recharge rate, which is difficult to measure directly [66,70,72,75], improves understanding of watershed scale hydrologic processes. For the Korba case study, groundwater recharge has been estimated to be up to 10% of the 420 mm/year mean average rainfall [7,10]. Ennabli [10] provides an excellent mathematical description of the problems related to artificial recharge value from precipitation (ranging between 15% and 25%). Furthermore, a 3D numerical model allowed Kerrou et al. [2] to estimate the areal recharge in the range of 8% and 30% of the annual rainfall (33.6-126 mm/yr), which is enough to counteract the intrusion of about 7 mm 3 /yr of seawater into the aquifer. However, groundwater overdraft reduced the submarine groundwater discharge of 16 mm 3 /yr, causing saltwater front to advance inland. Based on the works of Zghibi et al. [8,11], the precipitation recharge contributes 43.5% of the inflow to the aquifer system, depending especially on the frequency and the timing and intensity of the rain events, soil type and geological conditions. within very poor recharge zones except one that was located in the moderate zone based on the MIF method.
The groundwater recharge maps were also validated based on AUC. The area under the ROC curve (AUC) is an indicator of model quality. Generally, a good model has an AUC value of 0.7-0.9, while an excellent model has values over 0.9. According to Figure 10, the AHP method outperformed MIF (AHP AUC: 75.6%; MIF AUC: 70.4%), although both methods were reasonably accurate for mapping groundwater recharge in the Korba area. Overall, the validation results increase confidence in the applied methodology as a useful framework for rapid assessment of groundwater recharge to inform siting artificial recharge structures and other groundwater management efforts.

Groundwater Recharge Estimation
Quantifying groundwater recharge rate, which is difficult to measure directly [66,70,72,75], improves understanding of watershed scale hydrologic processes. For the Korba case study, groundwater recharge has been estimated to be up to 10% of the 420 mm/year mean average rainfall [7,10]. Ennabli [10] provides an excellent mathematical description of the problems related to artificial recharge value from precipitation (ranging between 15% and 25%). Furthermore, a 3D numerical model allowed Kerrou et al. [2] to estimate the areal recharge in the range of 8% and 30% of the annual rainfall (33.6-126 mm/yr), which is enough to counteract the intrusion of about 7 mm 3 /yr of seawater into the aquifer. However, groundwater overdraft reduced the submarine groundwater discharge of 16 mm 3 /yr, causing saltwater front to advance inland. Based on the works of Zghibi et al. [8,11], the The GWRZ maps were used to derive recharge-precipitation ratios for the two calculation methods (Table 11). A simple estimation of water volume (W) in the subsurface media was performed for the AHP and MIF methods by calculating the volume of annual precipitation (P) × recharge ratio ×% of area, where P = 175.56 106 m 3 /yr (420 mm/yr) and S (surface) = 438 km 2 :  Thus, based on the results of the AHP method 27.73% (116.5 mm/yr) of the annual rainfall recharges the Korba aquifer whereas the MIF method indicates that recharge is about 32.60% (137.8 mm/yr) of the annual rainfall. These estimated recharge values can be further verified using tracer techniques and water-table fluctuations method to inform groundwater conservation planning.

Conclusions
Two MCDM methods, namely analytical hierarchy process (AHP) and multi-influencing factor (MIF) were applied using a geospatial analysis framework to delineate groundwater recharge zones (GWRZ) in the Korba unconfined aquifer in northeastern Tunisia. The factors that were considered include lithology, land use/land cover, slope, geomorphology, lineament density, rainfall, drainage density and, soil type. The AHP method classified the GWRZs as very good (16%), good (53%), moderate (11.5%), poor (15.8%) and very poor (2.9%). The classification based on the MIF method suggests that 30.7% of the study area has very good recharge potential while 50% is considered good. The remainder of the Korba area has moderate (9.6%), poor (8.1%) and very poor potential (1.6%). The GWRZ maps were used to estimate recharge-precipitation ratios to quantify percolation (i.e., 27.7% (116.5 mm/yr) based on AHP and 32.6% (137.8 mm/yr) based on MIF). The groundwater recharge maps were verified using field yield data and water level depth from 25 shallow wells. While the accuracy of the GWRZs generated by both MCDM methods was satisfactory, the validation phase indicates that the AHP outperformed the MIF in the Korba region (AHP accuracy: 75.6% compared to MIF accuracy: 70.4%). The results are valuable for groundwater management and planning artificial recharge to mitigate groundwater tale decline and salt water intrusion in this coastal aquifer. The presented case study illustrates how groundwater recharge potential can be assessed using AHP and MIF to support sustainable groundwater management.