Delineation of Suitable Zones for the Application of Managed Aquifer Recharge (MAR) in Coastal Aquifers Using Quantitative Parameters and the Analytical Hierarchy Process

Abstract

Coastal aquifer salinization is usually related to groundwater overexploitation and water table decline. Managed Aquifer Recharge (MAR) can be applied as a measure to reverse and prevent this phenomenon. A detailed literature review was performed to identify the various methods and parameters commonly used to determine suitable sites of MAR application. Based on the review results, a new multi-criteria index (SuSAM) that is compatible to coastal aquifers was developed to delineate suitable zones for MAR application. New parameters were introduced into the index, such as distance from the shore and hydraulic resistance of the vadose zone, while factor weights were determined using the Analytical Hierarchy Process (AHP) and single sensitivity analysis. The applicability of the new index was examined in the coastal aquifer of the Anthemountas basin located in northern Greece. The most suitable areas for MAR application cover 28% of the aquifer’s surface area, while 16% of the area was characterized as non-suitable for MAR application. The new method constitutes the first step of the managed aquifer recharge concept for the delineation of MAR-suitable zones in coastal aquifers.

1. Introduction

Salinization of coastal aquifers has become a global issue in the last decades, influencing socio-economic development, agricultural productivity, and environmental sustainability. Groundwater depletion due to overexploitation is the main cause of coastal aquifer salinization in Greece [1]. Two main salinization processes occur in depleted aquifers: (a) seawater intrusion [1] and (b) salt water upconing [2]. Mapping the vulnerability of coastal aquifers to seawater intrusion [3] and salt water upconing [4] has been proposed as a tool with which to prevent groundwater salinization. More specifically, vulnerability maps depict zones where salinization prevention measures can be applied. Decreasing pumping rates, well reallocation, and crop type changes are commonly recommended to inverse negative groundwater balances. However, in many cases it is not feasible to implement such actions due to low acceptance from farmers, land owners, and other individual users. Therefore, water resources of coastal zones will continue to be influenced, and the salinization phenomenon will spread further [5]. Clearly, more active measures are necessary to prevent further salinization and finally inverse the phenomenon. For instance, increasing the recharge of an aquifer by applying Managed Aquifer Recharge (MAR) procedure and techniques can balance outflows and inflows, thus stabilizing and reinforcing piezometric head. Managed Aquifer Recharge is the purposeful recharge of water to aquifers for subsequent recovery or environmental benefit. It involves methods such as riverbank filtration, stream bed weirs, infiltration ponds, and injection wells in order to increase groundwater storage. However, the following question arises: Where can Managed Aquifer Recharge be successfully applied to counteract the salinization of coastal aquifers? To answer this question, the main principles of MAR should be understood and are therefore briefly presented below.

Due to the existing techniques available, the application of MAR often provides the cheapest form of a new, safe water supply [6]. Additionally, it can be applied by small communities and individuals using locally available materials and manpower. A detailed description of MAR is presented by Page et al. [7], who list the main advantages of the method as:

(a)

MAR can be applied in urban areas,

(b)

Water losses are negligible compared to surface storage,

(c)

MAR requires less surface area than surface storage,

(d)

The reallocation of existing wells is avoided.

If the application of MAR is successful, then brackish aquifers can be restored, the water supply to urban areas can be enhanced, groundwater-dependent ecosystems can be protected, and evaporation losses can be reduced [8].

The quantities of water required for MAR application can be obtained from (a) precipitation collected from rooftops or other demarcated areas; (b) surface water from lakes, rivers, and torrents; (c) treated municipal and industrial wastewaters; (d) other aquifers; (e) stormwater; and (f) mains water. Depending on the water quantities available and the characteristics of the aquifer concerned, the main approaches used to introduce water into the aquifer are (a) spreading methods, (b) recharge shafts, (c) injection wells, (d) induced recharge, and (e) improved land and watershed management. The application of MAR is usually proposed in areas where a decrease in groundwater level occurs; the availability of water is inadequate during the dry season, and groundwater quality is poor.

According to Karanth [9], the feasibility of MAR is driven by the following five factors:

• The availability of sites suitable for MAR application,
• The presence of water sources,
• A favorable hydrogeological environment,
• Optimal hydrodynamic conditions of the aquifer,
• Results of the cost–benefit evaluation.

Comprehension of the favorable hydrogeological environment is essential to ensure the successful application of MAR. More specifically, the criteria of site selection should be first determined. Various site selection factors for MAR application have been examined since the 1970s [10,11]. In this study, a detailed literature review was carried out to identify the optimal parameters and methods used to select suitable sites for MAR application to date. Following the literature review, a site selection index was developed to delineate the suitable zones for MAR application in coastal aquifers. The applicability of this index was then examined in the coastal aquifer of Anthemountas basin (Greece).

Literature Review

Numerous studies exist in the literature, with each proposing indices with different parameters to delineate suitable sites for MAR application. In the 1990s, the delineation of suitable areas was based on geomorphological and geological analysis using remote sensing techniques [12]. Ramsamy and Anbazhagan [13] set priority areas based on geomorphological units and hydrogeological data to determine suitable sites to apply MAR methods such as percolation ponds, pitting, induced recharge, and desiltation of existing tanks. Geophysical methods were also used to determine the optimum hydrogeological sites for placing recharge wells by linking electrical resistivity to aquifer permeability [14]. Saraf and Choudhury [15] used remote sensing and Geographical Information Systems (GIS) to overlay slope, geological, geomorphological, and lineament maps so as to determine a suitable site for recharge basins. The application was performed in a fissured rock aquifer in India, and remote sensing was used to map lineaments and link them with fault zones. Ghayoumian et al. [16] also considered characteristics of the aquifer itself such as thickness, transmissivity, and groundwater quality to localize zones for flood spreading using a Decision Support System (DSS). In the porous aquifer studied, the DSS was used to evaluate the contribution of each parameter present in the optimal location for MAR application. The authors later modified their initial approach using Boolean and fuzzy logic to determine suitable MAR application sites [17]. Hence, based on the selected criteria, unsuitable zones were excluded, and a more precise suitability map was created. Jasrotia et al. [18] used Boolean logic and GIS and introduced, in addition to other parameters, aquifer storativity and specific capacity. Taheri [19] used geophysical methods to determine the lithology of an aquifer to indicate suitable sites for MAR. In this approach, electrical resistivity was linked to the sedimentary formations and their permeability in order to determine the permeable zones. The aforementioned methods were all applied in India and Iran.

After the 2010s, an increasing number of studies was published that introduced new indices and methodological approaches for site suitability in other countries with different management practices and socio-economic cultures. Chenini et al. [20] developed a multi-criteria method in a GIS environment for the selection of optimal sites for MAR application in Tunisia. Saud [21] was the first to introduce an index for site suitability in Saudi Arabia. Chowdhury et al. [22] followed the established parameters to develop an index of site suitability for MAR using the Analytical Hierarchy Process (AHP). In this research, AHP was applied to overcome the subjectivity of weight definition. Site selection for the application of MAR in karst aquifers was also tested in Iran using GIS and fuzzy logic [23]. The use of fuzzy logic contributed to overcome the subjectivity of the parameters classes. Nasiri et al. [24] specified the methodology to delineate suitable zones for flood spreading areas in Iran. The authors used a preference ranking organization method for enrichment evaluations named PROMETHEE, which is a multi-criteria decision analysis method. The SLUGGER-DQL score model was used to assign weights to the parameters influencing site suitability of MAR application in Jordan [25], whereas Mahmoud et al. [26] also included rainfall surplus in their DSS-based model for an application in Saudi Arabia. Genetic algorithms were innovatively introduced into the selection of flood spreading in Iran [27], and Boolean logic, in conjunction with GIS, was selected for the application of Zaidi et al. [28] in Saudi Arabia. In the latter study, a genetic algorithm was included in the analysis of the selected parameters (slope, alluvium thickness, geology, morphology, electrical conductivity, land use, drainage density, aquifer transmissivity, and elevation) and their weightings. Brown et al. [29] combined three indices to define zones of well injection in a karst aquifer in the USA. GIS-based methods for site selection were applied in Sri Lanka [30], Costa Rica [31], and Argentina [32]. Steinel et al. [33] considered existing infrastructure, such as dams, to select sites suitable for the infiltration of captured surface runoff. The optimum surface spreading basin was evaluated using a weighted overlay analysis model in a porous aquifer in the USA [34]. Farhadian et al. [35] used the Nash conflict resolution method to determine suitable sites for MAR application. The Nash method uses an equation to produce optimal levels to resolve conflicts between two or more stakeholders. Ahani Amineh et al. [36] introduced new parameters for the delineation of suitable zones for MAR application, including erosion density and proximity to existing wells. Remote sensing techniques were used to determine the fracture zones of a fissured rock aquifer in India [37], whereas Ghasemi et al. [38] included the hydraulic gradient of the aquifer and distance from rivers to specify suitable zones for MAR application in Iran. Singh et al. [39] included drainage order into their GIS-based application in India, while Christy and Lakshmanan et al. [40] selected sites for percolation ponds according to the permeability values obtained from an electrical resistivity analysis in a coastal porous aquifer. This approach is mainly used to determine local conditions for MAR applications and is site-specific. The parameters and tools of the existing methods used to select sites for the application of MAR according to the available relevant literature have been summarized and are presented in

Table 1. Summarized literature review of available research articles on MAR site selection.

A/A	Authors	Parameters Applied	Methods—Tools	Method of MAR	Aquifer Type	Country
1	Chopra and Sharma [12]	Geomorphological units (e.g., ridge, structural hills, alluvial fans, sand dunes, flood plains, river channels, seasonal rivulets)	Remote sensing	Tapping flood plains	Porous aquifer (inland)	India
2	Ramsamy and Anbazhagan [13]	Drainage density, Aquifer material, Groundwater level, Geology	Priority areas, Remote sensing	Percolation ponds, Pitting, Small dams, Induced recharge, Desiltation of existing tanks	Fissured rock aquifer (mountainous)	India
3	Anbazhagan and Ramsamy [14]	Water level, Electrical resistivity, Thickness of the vadose zone	Geophysical methods	Wells	Fissured rock aquifer (mountainous)	India
4	Saraf and Choudhury [15]	Slope, Geology, Geomorphology, Lineaments	Remote sensing, GIS	Recharge basins or reservoirs	Fissured rock aquifer (mountainous)	India
5	Ghayoumian et al. [16]	Slope, Infiltration rate, Sediment thickness, Transmissivity, Water quality	Decision support system, GIS	Flood spreading	Porous aquifer (inland)	Iran
6	Ghayoumian et al. [17]	Slope, Infiltration rate, Depth to groundwater, Quality of alluvial sediments, Land use	Boolean, Fuzzy logic, Remote sensing, GIS	Not specified	All types (coastal)	Iran
7	Jasrotia et al. [18]	Lithology, Geomorphology, Land use/land cover, Drainage, Hydrologic soil texture, Depth to water table, Transmissivity, Permeability, Storativity, Specific capacity, Infiltration	Boolean logic, Conditional methods, GIS	Not specified	Porous aquifer (inland)	India
8	Taheri [19]	Electrical resistivity	Geophysical methods	Not specified	Porous aquifer (inland)	Iran
9	Chenini et al. [20]	Watershed limit, Drainage, Drainage density, Lithology, Fractured outcrops, Lineament, Permeability, Piezometry	Multi-criteria analysis, GIS	Not specified	All types (inland)	Tunisia
10	Saud [21]	Precipitation, Lithology, Rock fractures, Slope, Drainage, Land cover/use	Remote sensing, GIS	Not specified	All types (inland)	Saudi Arabia
11	Chowdhury et al. [22]	Geomorphology, Geology, Drainage density, Slope, Aquifer transmissivity	Remote sensing, Analytic Hierarchy Process, GIS	Not specified	Fissured rock, Porous aquifer (inland)	India
12	Malekmohammadi et al. [23]	Slope, Geology, Groundwater depth, Potential for runoff, Land use, Groundwater electrical conductivity	Fuzzy logic, GIS	Not specified	Karst aquifer, Porous aquifer (inland)	Iran
13	Nasiri et al. [24]	Slope, Water quality, Geology, Alluvium thickness, Land use, Transmissivity, Geomorphology, Drainage density	PROMETHEE II, Analytic Hierarchy Process, GIS	Flood spreading	Porous aquifer (inland)	Iran
14	Hammouri et al. [25]	Slope, Land use, Geomorphology, Geology, Well density, Water quality, Depth to groundwater, Runoff available	GIS, SLUGGER-DQL score model	Not specified	All types (inland)	Jordan
15	Mahmoud et al. [26]	Rainfall surplus, Slope, Potential runoff coefficient, Land cover/use, Soil texture	GIS, DSS, Analytic Hierarchy Process	Not specified	Not specified	Saudi Arabia
16	Rahimi et al. [27]	Slope, Alluvium thickness, Geology, Morphology, Electrical conductivity, Land use, Drainage density, Aquifer transmissivity, Elevation	GIS, Genetic algorithm, Analytic Hierarchy Process	Flood spreading	Porous (inland)	Iran
17	Zaidi et al. [28]	Slope, Soil texture, Vadose zone thickness, Groundwater quality (TDS), Type of formation, Land use	Boolean Logic, GIS	Not specified	All types (inland)	Saudi Arabia
18	Brown et al. [29]	Density ratio, Effective porosity, Aquifer gradient, Injection time, Storage duration, Dispersivity, Aquifer thickness, Hydraulic conductivity, Water quality	Index, Statistical analysis	Well injection (brackish water)	Karst aquifer (coastal)	USA
19	Senanayake et al. [30]	Rainfall, Lineament, Slope, Drainage, Land use/land cover, Geology, Geomorphology, Soil characteristics	GIS	Not specified	All types (inland)	Sri Lanka
20	Bonilla Valverde et al. [31]	Hydrogeological aptitude, Terrain slope, Top soil texture, Drainage network density	GIS, Boolean logic, Sensitivity analysis	Not specified	All types	Costa Rica
21	Quiroz Londoño et al. [32]	Drainage density, Geomorphologic units, Soil media, Land cover, Slope and aspect	Remote sensing, Fuzzy logic, GIS	Not specified	Fissured rock, Porous aquifer (inland)	Argentina
22	Steinel et al. [33]	Distance to international borders, Distance to wadis, Catchment size, Rainfall, Land cover, Slope, Existing dams, Thickness of aquifer, Depth to water table, Flow gradient, Distance to faults, Groundwater salinity, Groundwater contamination, Distance to roads, Distance to active government wells	Boolean logic	Infiltration of capturedsurface runoff	All types (inland)	Jordan
23	Fournier et al. [34]	Hydraulic conductivity, Existing land use, Composite suitability, Binary mask, Reference source with selected destination	GIS, Weighted overlay analysis model	Surface spreading basin	Porous aquifer (inland)	USA
24	Farhadian et al. [35]	Precipitation, Vegetation, Distance from connected roads, Soil, Distance from rivers, Geology, Slope, Land use	GIS, Analytic Hierarchy Process, Nash conflict resolution method	Not specified	All types (inland)	Iran
25	Ahani Amineh et al. [36]	Source and groundwater compatibility, Source water quality, Storage availability, Groundwater quality (EC), Construction cost, Source water availability, Aquifer characteristics, Demand, Operating cost	GIS, Analytic Hierarchy Process	Surface spreading	Porous aquifer (inland)	Iran
26	Selvarani et al. [37]	Geology, Geomorphology, Slope, Drainage density, Lineament density	Remote sensing, GIS, Analytic Hierarchy Process	Not specified	Fissured rock, Porous aquifer (inland)	India
27	Ghasemi et al. [38]	Hydraulic gradient, Transmissibility, Aquifer thickness, Land use, Minimum area, Distance of supply sites, Distance from highways and freeways, Distance from residential areas, Distance from rivers, Distance from wastewater, Elevation difference	GIS, Fuzzy logic	Not specified	Porous aquifer (inland)	Iran
28	Singh et al. [39]	Slope, Soil, Land use, Drainage order	GIS, Analytic Hierarchy Process	Not specified	Porous aquifer (inland)	India
29	Christy and Lakshmanan et al. [40]	Electrical resistivity	Geophysical methods	Percolation ponds	Coastal porous aquifers	India

.

Site selection constitutes the first step for the application of MAR, based on hydrogeological and morphological parameters, as well as existing infrastructure. Initial research methods of site selection were simple, using just a few parameters such as geomorphological units, geological formations, and groundwater depth. The most commonly used parameters have been slope, water level, and drainage density. As research and technology progressed, GIS-environments allowed the use of a higher number of parameters, including aquifer hydraulics and existing infrastructure. The Analytical Hierarchical Process (AHP) has been the most commonly used technique to define parameter weights. The literature review revealed that site selection indices should be developed according to the specific characteristics of the target aquifer. Additionally, the use of GIS can ensure the use of multi parameters in a wider area, while mathematical processes can overcome the subjectivity of ratings and weight assignments.

In respect to the presented literature review, this research deals with the first element of managed aquifer recharge, i.e., the delineation of MAR-suitable zones within a coastal aquifer. This approach is currently lacking from the existing literature, and hence it could provide a further tool for optimum site selection to apply MAR in coastal aquifers. Thus, a spatial, multi-criteria index was developed by incorporating parameters found in the specific hydrogeological environments of coastal aquifers. The index was then applied to a specific coastal aquifer located in northern Greece; however, it can also be adapted for application in other countries and regions.

2. Methodology

The literature review revealed the most commonly used parameters and tools for the delineation of suitable sites for MAR application. Although many indices have been developed, to the best of our knowledge previous research has not considered coastal aquifer environments. Hence, a multicriteria approach of MAR siting specific to coastal aquifers was developed and applied in a case study. The thematic maps were developed in a GIS environment, while the final index was produced by using overlay techniques.

2.1. Anthemountas Coastal Aquifer

The coastal aquifer of Anthemountas covers an area of 157 km², with a mean topographic slope and elevation of 5% and 65 m, respectively (Figure 1). The water demands of the basin are met with groundwater obtained from the coastal aquifer, while a high population density and intensive agricultural activities have led to overexploitation of the groundwater. Neogene, Pleistocene, and Holocene sediments host the porous aquifer that consists mainly of gravel, sand, and marls. A detailed description of the aquifer can be found in relevant studies [1,41]. The aquifer is found in both confined and unconfined conditions, while negative piezometric head reaches up to 40 m below sea level (b.s.l.) in a variable zone up to 8 km from the coastline. Additionally, the concentration of Cl⁻ reaches 350 mg/L in some areas. Indisputably, the confrontation of groundwater salinization and the progressive recovery of depleted reserves using MAR should be a priority in the coastal aquifer of Anthemountas basin.

2.2. Site Selection Index to Apply MAR

In this study, a novel index was developed to delineate suitable zones to apply MAR within coastal aquifers. The model provides a multi-criteria analysis of hydrogeological and morphological parameters, as well as existing infrastructure, in a GIS environment. The natural neighbor interpolation method was used to develop the thematic maps. The site (S) suitability (Su) index to apply (A) MAR (M) (SuSAM) comprises the following ten parameters, which can be numerically presented (quantitative parameters): topographic slope (%), shore (distance—m), drainage network (distance—m), depth of groundwater (m), piezometric head (m), vadose zone (log of hydraulic resistance), groundwater quality (electrical conductivity, μS/cm), transmissivity (m²/day), water availability (distance—m), and main roads (distance—m). The parameters were chosen according to their relevance to MAR, as concluded from the literature review above. Additionally, new parameters were added in order to enhance applicability in coastal aquifers. It is worth mentioning that the index was designed for use in coastal aquifer environments, and hence new parameters such as distance from the coast have been included. A thematic map with a pixel size of 25 × 25 m was produced for each parameter. Τhereafter, a rating score was assigned for each factor value on a scale of 0 to 10 (

Table 2. Parameters and rating of the MAR suitability index.

A/A	Parameter	Factor Variable		Rating
		Class	Range
Morphological
1	Slope (%)	Very High High Moderate Low Very low Extremely low	0–2 2–5 5–10 10–15 15–35 >35	10 8 6 4 2 0
2	Shore (distance—m)	Very High High Moderate Low Very low Extremely low	>1000 750–1000 500–750 300–500 100–300 <100	10 8 6 4 2 0
3	Drainage network (distance—m)	Very High High Moderate Low Very low Extremely low	<100 100–300 300–500 500–750 750–1000 >1000	10 8 6 4 2 0
Hydrogeological
4	Depth of groundwater (m)	Very High High Moderate Low Very low Extremely low	>10 10–8 8–6 6–2 2–0 Artesian	10 8 6 4 2 0
5	Piezometric head (m)	Very High High Moderate Low Very low Extremely low	≤ 10 −10–0 0–6 6–10 10–20 >20	10 8 6 4 2 0
6	Vadose zone (log of hydraulic resistance)	Very High High Moderate Low Very low Extremely low	<1 1–2 2–3 3–4 4–5 >5	10 8 6 4 2 0
7	Groundwater quality (electric conductivity—μS/cm)	Very High High Moderate Low Very low Extremely low	<500 500–750 750–1000 1000–1500 1500–2000 >2000	10 8 6 4 2 0
8	Transmissivity (m2/day)	Very High High Moderate Low Very low Extremely low	>100 70–100 30–70 10–30 5–10 <5	10 8 6 4 2 0
Infrastructures
9	Water availability (distance from dams, village/city, waste water treatment facilities—m)	Very High High Moderate Low Very low Extremely low	<500 500–1000 1000–1500 1500–2000 2000–3000 >3000	10 8 6 4 2 0
10	Main roads (distance—m)	Very High High Moderate Low Very low Extremely low	>1000 750–1000 500–750 300–500 100–300 <100	10 8 6 4 2 0

) that covered the following six (6) classes: extremely low, very low, low, moderate, high, and very high. The class ranges were defined based on the literature review, while slight empirical modifications were performed to adapt the index to coastal environments. The final map was obtained using the overlay technique in a GIS-environment and by applying the final SuSAM index (Equation (1)). Figure 2 presents a flow chart of the method followed.

$$\mathrm{SuSAM}=\frac{{{\displaystyle \sum }}_{i=1}^{10}\left(W_i\ast R_i\right)}{{{\displaystyle \sum }}_{i=1}^{10}{W}_i}$$

(1)

in which W and R correspond to parameter weight and rating, respectively.

2.3. Weight Definition and Validation of the Model

The methodological approach applied was based on a multicriteria index incorporating ten (10) parameters each with different influences on the final site selection for MAR application. Defining the weight of each factor is critical to ensure the reliability of the index. The subjectivity involved in weight definition can be overcome by using statistical or structural techniques. In this study, the Analytical Hierarchy Process (AHP) [42,43] was coupled with the single parameter sensitivity analysis [44] to alleviate subjectivity. The AHP approach was performed first and followed by a pairwise comparison test using a 10 × 10 matrix in which diagonal elements are equal to 1. In a pairwise comparison, the higher the parameter value, the higher the influence of that parameter. Hence, the weights of each parameter are produced and applied in the final index. Consistency of the AHP was then checked using the consistency ratio CR (Equation (2)). The CR was calculated as equal to 0.05 (<0.1) and thus verified the consistency of the application.

CR = (CI/RI)

(2)

in which RI is the random index and CI is the consistency index.

Following weight definition using AHP, a single sensitivity analysis was performed to validate the initial weights. The single sensitivity analysis provides effective weights following the application of Equation (3):

W = (Pr × Pw/V) × 100

(3)

in which W is the effective weighting, Pr is the rating value, Pw is the initial weight, and V is the index score.

Sensitivity analysis is usually used to assess the uncertainty in multi-criteria models and to determine the importance of each criterion. The importance of each criterion is quantitatively addressed by the average of the effective weighting. This value can be adopted as a validated weight and assigned to the corresponding parameter. Hence, the validated weights increase the robustness of the multi-criteria model.

The application of AHP overcomes the subjectivity of weight definition of the parameters. Additionally, the validation of weights using sensitivity analysis increases the reliability of the final index. The multi-criteria approach is the most suitable in cases of complex aquifers—similar to the studied one—due to its ability to evaluate large data sets belonging to different parameters. Additionally, the application of sensitivity analysis highlights the less important parameters, which can be excluded in case studies lacking available data. It is worth mentioning that the evaluation of all suggested parameters contributes to a more thorough understanding of the hydrogeological regime.

3. Results and Discussion

In the present study, a novel index (SuSAM) was developed to delineate the suitable zones of an aquifer to apply MAR. The index was customized for the specific hydrogeological conditions of a coastal aquifer. New parameters were included in the concept of site suitability for MAR application, such as distance from the shore and hydraulic resistance of the vadose zone. The distance from the shore was included, because in nearby coastal areas, the groundwater is prone to salinization. The hydraulic resistance of the sediment layers describes the resistance of the vadose zone to vertical water flow. Hence, low hydraulic resistance of the vadose zone favors the application of MAR (e.g., surface spreading). AHP was used to define the weights of each parameter that were then validated using sensitivity analysis.

Table 3. Pair wise comparison of the parameters included in the site selection for MAR application.

Parameter	Topographic Slope	Distance from the Shore	Drainage Network	Groundwater Depth	Piezometric Head	Vadose Zone	Groundwater Quality	Transmissivity	Water Availability	Main Roads	Weights (%)
Topographic slope	1	2	6	4	2	1	2	6	4	8	22
Distance from the shore	0.5	1	4	2	1	0.5	1	4	2	6	12
Drainage network	0.17	0.25	1	0.5	0.25	0.17	0.25	1	0.5	2	3
Groundwater depth	0.25	0.5	2	1	0.5	0.25	0.5	2	1	4	6
Piezometric head	0.5	1	4	2	1	0.5	1	4	2	6	12
Vadose zone	1	2	6	4	2	1	2	6	4	8	22
Groundwater quality	0.5	1	4	2	1	0.5	1	4	2	6	12
Transmissivity	0.17	0.25	1	0.5	0.25	0.17	0.25	1	0.5	2	3
Water availability	0.25	0.5	2	1	0.5	0.25	0.5	2	1	4	6
Main roads	0.12	0.17	0.5	0.25	0.17	0.12	0.17	0.5	0.25	1	2

presents the pairwise comparison of the criteria significance and the parameter weights obtained, while

Table 4. Results of the sensitivity analysis.

Parameter	Effective Weighting (%)
	Min.	Max.	Standard Deviation	Average/Final Weight
Topographic slope	0	43.7	7.9	24
Distance from the shore	0	39.7	5.6	20
Drainage network	0	9.1	1.7	4
Groundwater depth	1.8	19.9	2.8	10
Piezometric head	0	28.4	7.4	7
Vadose zone	0	36.7	6.6	9
Groundwater quality	0	25.5	3.2	12
Transmissivity	0	7.2	1.1	4
Water availability	0	18.5	3.5	7
Main roads	0	9.9	2.1	3

presents the results of sensitivity analysis. The thematic maps produced are shown in Figure 3 and illustrate the spatial distribution of each parameter’s rating score. The geomorphological, hydrogeological, and infrastructural parameters of the SuSAM index are discussed below with focus on their relevance and the spatial distribution of their rating values.

3.1. Geomorphological

3.1.1. Topographic Slope

The parameter of topographic slope is one of the most commonly-used parameters in MAR suitability zones. Flat areas favor the application of MAR (e.g., flood spreading) [24], while steep slopes are a limiting factor for MAR-associated infrastructure. In this study, a digital elevation model with a resolution of 25 m × 25 m was used to produce the slope map (Figure 3). The determination of class range was based on previous studies [24] and related to Demek’s classification [45]. Within the study area, slopes of up to 57% gradient can be found, the highest values being located in the central southern part of the porous aquifer and corresponding to the lowest parameter rating. Convenient zones of shallow slopes suitable for MAR application are located in the lowland part of the basin.

3.1.2. Shore (Distance)

In the literature, several applications of site selection for MAR have concentrated on coastal zones [17,29,40]. However, potential salinization of an aquifer due to seawater intrusion near the shore has not been taken into account. Hence, distance from the shore was included as a new parameter in the SuSAM index, similar to the concept of coastal aquifer vulnerability to seawater intrusion [3]. As expected, unsuitable zones for MAR are located close to the shoreline, while more suitable zones occur in the mainland (Figure 3). The distance from the shore was calculated using the multiple ring buffer tool in GIS.

3.1.3. Drainage Network (Distance)

The drainage factor has been widely-used in terms of density corresponding to permeable formations [20,21,22]. In this study, formation permeability was considered within the vadose zone parameter, and the drainage factor used here was the distance (m) from the drainage network. Areas closer to a drainage network have increased water availability, and hence are more suitable for MAR application. However, these zones could be used to establish new infrastructure to collect water, such as small dams. Ahani Amineh et al. [36] also used the distance from the drainage network (rivers) to determine MAR suitability. Figure 3 shows that the drainage network is well-developed in the study area, thus rendering the construction of water-collecting infrastructure feasible within the boundaries of the aquifer.

3.2. Hydrogeological

3.2.1. Depth of Groundwater

Groundwater depth is a critical parameter in the assessment of the MAR-suitable sites [33]. In areas with low groundwater depth, suitability for MAR application decreases due to the possibility of groundwater flooding surface land. Additionally, zones with artesian phenomena are unsuitable for the application of MAR methods. The classes and rating scores used here were determined based on a previous study [18]. In the study site, groundwater depth ranges between 1 and >100 m from the surface. The majority of the study area is characterized by very high suitability to MAR (Figure 3), while the artesian phenomenon does not occur in the Anthemountas basin.

3.2.2. Piezometric Head

Piezometric head has been widely-used in similar studies mainly to describe groundwater flow direction [20]. In this study, the piezometric head parameter was included due to its relation with the salinization process due to seawater intrusion in coastal zones. Negative piezometric head can reverse groundwater flow from the shore towards the mainland rendering the coastal aquifer prone to salinization. Hence, it is suggested that low piezometric head should be favoured in the site selection of MAR in coastal aquifers. In Anthemountas basin, negative piezometric head dominates in the coastal zone due to overexploitation. Hence, higher values corresponded to the coastal zone (Figure 3).

3.2.3. Vadose Zone

The permeability of the vadose zone is a critical parameter for the successful application of MAR. Similar studies usually use soil permeability [32]; however, this is not representative of the entire thickness of the vadose zone. Additionally, upper soil layers are usually treated during the construction of MAR-related infrastructure. Other studies link permeability with geological formations in a quantitative manner [20]. However, porous media is most usually characterized by a high degree of anisotropy. A quantitative approach to consider the permeability of the vadose zone was innovatively introduced into the SuSAM index by using the hydraulic resistance of the sediment layers (Equation (4)). The index is based on the hydraulic conductivity (K) of each of the sedimentary layers and their thickness (d). The relevant data was obtained from previous studies on Anthemountas basin [46], and high hydraulic resistance corresponds to low suitability for MAR application. Figure 3 shows that, based on the vadose zone, suitability is low in the coastal areas, while high suitability is located in the center of the aquifer. The use of this parameter excludes the confined aquifers, while it preconceives the application of methods such as riverbank filtration, stream bed weirs, and infiltration ponds.

$$c={\displaystyle \sum }\frac{di}{Ki}$$

(4)

in which c = hydraulic resistance, di = thickness of the layer, and Ki = hydraulic conductivity of the layer.

3.2.4. Groundwater Quality

MAR application is strongly linked to aquifer groundwater quality. When MAR is applied to prevent and reverse salinization of coastal aquifers, the salinity status of the existing groundwater should be considered. Total dissolved solids (TDS) and electrical conductivity (EC) have been used for this purpose in previous studies [25,33,36]. In the SuSAM index, EC was chosen as the parameter to assess groundwater quality. The use of EC in the proposed index aims to prevent MAR application in salinity-influenced zones and avoid undesirable phenomena. Therefore, MAR will attain higher piezometric head to reverse groundwater flow and hence salinization. The thematic map in Figure 3 shows that the studied aquifer is characterized by moderate to low suitability near the shoreline due to the presence of high EC values. High EC values are also observed in the central part of the basin and are related to the influence of geothermal fluids.

3.2.5. Transmissivity

Hydraulic parameters are also included in attempts to delineate suitable zones to apply MAR. In the literature, site selection has included hydraulic conductivity [34], storativity [18], and transmissivity [27] of the aquifers studied. The thickness of the aquifer has also been included in some studies [33,38]. In this research, transmissivity of the aquifer was incorporated into the SuSAM index, as this parameter includes both the thickness and the hydraulic conductivity of the aquifer. The higher the transmissivity, the higher the suitability for MAR. In Anthemountas, transmissivity ranges from 3 to 430 m²/day. The highest values can be found in the central part of the aquifer, while the lowest values can be found in the south (Figure 3).

3.3. Infrastructures

3.3.1. Water Availability

Water availability is essential for successful MAR application and this parameter corresponds to sources of available water that can be used for MAR. Water sources can be provided by existing dams, wastewater treatment facilities, and villages or cities where buildings can collect rain water. This concept has also been considered as an economic factor for the application of MAR [38]. The wider area of the studied aquifer includes numerous villages and cites, two existing dams, and two wastewater treatment facilities, and hence the buffer zones were determined. The thematic map of Figure 3 shows that a large area of the aquifer is characterized by very high suitability to apply MAR due to high water availability.

3.3.2. Main Roads

Main road networks are a groundwater pollution source, and for this reason have been included in some MAR suitability assessments [33,35]. In this study, distance (m) from main roads was incorporated into the SuSAM index and the buffer zones produced. Lowest suitability for MAR application was obtained when distances from main roads are shorter. On the contrary, the longer the distance from main roads, the greater the suitability for MAR application (Figure 3). However, it is worth mentioning that main roads can be managed, and water be treated to prevent groundwater pollution.

3.4. MAR Suitability Map and Validation of the Index

The final map of MAR suitability in the study area was produced by applying a relative weight to each of the ten parameters used and then overlaying them onto a single map. Based on the AHP results, the parameters of topographic slope and vadose zone were assigned the highest weights, while drainage and transmissivity had the lowest relative weights. Based on the sensitivity analysis, distance from the shore is weighted higher than the vadose zone parameter, while the slope parameter remained the most important in terms of weight value. In addition, the relative weights of drainage network, groundwater depth, transmissivity, and infrastructure increased, and the weight of piezometric head decreased. Decreasing weights correspond to lower parameter influence and vice versa.

The final map produced is presented in Figure 4, and shows the potential suitable sites for MAR application in the coastal aquifer of Anthemountas basin. Three classes (low, moderate, and high) correspond to potential suitability for MAR, while non-suitable sites are located in the southern part of the aquifer. In this area, the steep slopes alone render MAR application unsuitable. Figure 5 shows the distribution (%) of MAR suitable areas within the coastal aquifer. A substantial 16.1% of the porous aquifer is unsuitable for MAR application, while 33% of the area is characterized by moderate suitability. Nevertheless, 28% of the studied aquifer is characterized by high suitability to apply MAR.

The novel SuSAM index developed here deals with the fundamental MAR concept of delineation and pre-screening of potentially suitable sites. Integrated water resource management includes MAR application in order to reinforce the natural recharge of an aquifer [47]. The index is flexible, as criteria can be eliminated in cases in which specific data are not available. Removing parameters with low influence weightings may not significantly influence the final results; however, including additional parameters into the index would ensure a more accurate suitability map. Additionally, in future studies the use of continuous functions could be investigated instead of the discrete values that were adopted in this study.

The next step in this research is a simulation process to test the suitability of the SuSAM index and quantify the results of MAR in the suitable sites [48,49,50]. In coastal zones, multilayer sampling is critically important to define the vertical distribution of salinity [51], while the determination of groundwater residence time [52] could further increase the effectiveness of MAR application. Although the SuSAM index was developed to be compatible with a coastal hydrogeological environment, parameters of agricultural activities and nitrate pollution could be included in the method with the aim of reducing nitrate pollution [53]. An expanded DSS-ΜΑR system could include agricultural planning [54] and groundwater vulnerability maps, which are useful tools for groundwater management [55]. The application of MAR strengthens the concept of integrated water resource management and could also help solve the problem of treated wastewater misuse [56], as well as improve urban water quality [57]. The feasibility and contribution of MAR application is undeniable; however, stakeholders and local populations should begin to accept integrated water resource management as a top priority, including the application of MAR [58].

4. Conclusions

In this study, a multi-criteria index named SuSAM was developed for the selection of suitable sites for MAR application in a coastal aquifer environment. The ten parameters included in the index are all quantitative, and novel parameters, such as distance from the shore and hydraulic resistance of the vadose zone, were introduced to cover the specific conditions of coastal aquifers. Additionally, the method excludes confined aquifers, while it preconceives the application of methods such as riverbank filtration, stream bed weirs, and infiltration ponds. Weight definition and validation were based on the Analytical Hierarchy Process and single sensitivity analysis. Topographic slope and distance from the shore were weighted with the highest values. The method was applied to the coastal aquifer of Anthemountas basin and successfully delineated the suitable areas to apply MAR. The MAR application suitability map showed that 28% of the basin’s surface area can be characterized as very suitable, while 16.1% is non-suitable for MAR application.

The novel index method deals with the first step of MAR application, which is the delineation of suitable sites. The most appropriate MAR method should then be chosen, followed by simulation processes to quantify the results.