Site selection for managed aquifer recharge in the city of Kabul, Afghanistan, using a multi-criteria decision analysis and geographic information system

While the success and sustainability of managed aquifer recharge (MAR) strongly depends on many characteristics of the site, it is necessary to integrate the site characteristics and develop suitability maps to indicate the most suitable locations. The objective of this study is to integrate geographic information system (GIS) and multi-criteria decision analysis (MCDA) techniques to identify the most suitable areas for a MAR project in the Kabul city area, Afghanistan. Data for six effective criteria, including slope, drainage density, surface infiltration rate, unsaturated zone thickness, soil type and electrical conductivity, were collected and then a classification map was produced for each criterion in the GIS environment. By applying MCDA techniques, the weights of the effective criteria were obtained. A suitability map was generated from each technique separately based on a combination of all criteria weights and thematic layers. The result of the analytical network process (ANP) method was found to be more precise and reliable compared with that of the analytical hierarchy process (AHP) method. Based on the final suitability map produced from the ANP model, there is 3.7, 15.0, 37.4, 33.1 and 10.3% of the total area that is unsuitable, of low suitability, moderately suitable, suitable and very suitable for MAR application, respectively. As a final result of this work, seven sites have been prioritized based on land use. The integration of multi-criteria decision analysis and GIS is recognized as an effective method for the selection of managed aquifer recharge sites.


Introduction
Kabul is both the capital city and the largest city in Afghanistan, with a population of 4.5 million, and it is the fifth fastest-growing city in the world (Asian Development Bank 2015). In recent years, the city has seen a rapid increase in population, owing primarily to refugee returns and better job opportunities. Rapid population growth and urbanization in Kabul have created huge pressure on the groundwater resources. A US Geological Survey (USGS) report reveals that, during the 10-year monitoring period 2004-2013, water level decline was observed in all but two of 23 monitoring wells in the city of Kabul ranging from a few meters up to 30 m (Mack et al. 2013). Extraction of groundwater at rates greater than natural replenishment results in groundwater level decline. As a result of a lack of water storage and the seasonal variability of river flows, Kabul is among the world's most waterstressed cities (World Bank 2010) as it depends almost entirely on groundwater resources, with the majority of water coming from the shallow aquifer, reaching to 150-200 m(m) in depth. Within urban areas, most households use privatelyowned dug tube wells, while many households, particularly in informal hillside settlement areas, rely on costly tankerdelivered water (World Bank, 2010).
Many studies have been carried out over the past 35 years that have looked into the characteristics and performance of the Kabul aquifer system. This includes the studies from the Danish Committee for Aid to Afghan Refugees, the Japan International Cooperation Agency (JICA), the German development cooperation through KfW, the United States Agency for International Development, the USGS, and the Russian "Passport" study of 1982. A key recommendation of a 5year program funded by JICA (2011) was the need to use "artificial recharge" to strengthen sustainable groundwater yields (JICA 2011). Artificial recharge of groundwater where applicable can be used to help natural replenishment of groundwater and reduce groundwater decline. The application of managed aquifer recharge (MAR) often provides the cheapest form of new and safe water supply (Dillon 2005). A few studies on MAR have been carried out in Afghanistan, of which one modeled MAR using the surplus flow of the Kabul River in the rainy season (Masoom 2018). In addition, an Asian Development Bank (ADB) funded project is currently ongoing to conduct a pilot study on Kabul's MAR. The project will run from 2017 to 2021 and will create pilot recharge systems to test the feasibility of various recharge options such as recharge basins, injection wells, contour banking, and trenches.
Internationally, various studies propose different criteria and various methods to delineate suitable sites for MAR application. From 2010 onwards, an increasing number of studies have been published that have introduced new indices and methodological approaches for site suitability with different management practices and socio-economic cultures. Geographic information system (GIS)-based multi-criteria decision methods for site selection of MAR applications were applied in Tunisia (Chenini et al. 2010), Iran (Sabokbar et al. 2012; Moghaddam et al. 2017), Sri Lanka (Senanayake et al. 2016), Costa Rica (Valverde et al. 2016), and Argentina (Londoño et al. 2016), and in many other locations. Sallwey et al. (2018) performed a review study to increase knowledge of GIS-multi criteria decision analysis (MCDA). They established a database based on 63 studies applying GIS-MCDA in the context of MAR site selection. According to their study, GIS-MCDA studies are increasingly used to locate suitable areas for MAR application. In addition, some studies were conducted to develop a method using GIS-MCDA approach that can be used to locate and determine suitable areas of MAR on coasts, islands, and urban aquifers (Andersson and Engdahl 2017;Kazakis 2018;Aghazadeh et al. 2019).
Most of the previous studies used an integrated GIS-MCDA approach in order to determine the suitable areas for flood spreading and artificial recharge of aquifers. The different methods, such as Boolean logic, fuzzy logic, multiinfluence factor, analytical hierarchy process (AHP) and analytical network process (ANP), were applied and also their applicability was compared together. In order to perform multi-criteria decision analysis, many effective criteria were used in the previous studies. Thematic layers for these effective criteria were prepared and classified in a GIS system and then weighted using MCDA techniques and finally integrated in a GIS environment. The output of the studies was usually a final map of suitable areas for groundwater recharge. Thus, the general steps in all of those studies were data-layer preparation, classification, weighting, and integration into a GIS environment. The results of those studies indicate the efficiency of GIS-based multi-criteria decision making in site selection and its effectiveness in the rapid assessment of large areas in order to determine a suitable site for MAR application. As reviewed, several techniques have been applied for the site selection of groundwater recharge, but there are very few studies that used the ANP approach.
The main objective of this research is to identify suitable sites for the application of MAR in Kabul, Afghanistan, by employing a variety of multi-criteria decision analysis methods and GIS tools. The findings of the study can be used to gain good knowledge and insight into the various multicriteria decision analysis methods used in locating the most suitable sites for MAR application in an urban area.

Study area
The study area (the city of Kabul) is divided into 22 districts and is located nearly 1,800 m above sea level (m asl) within Kabul province in north-eastern Afghanistan. The city of Kabul is located approximately 34°31′ latitude to the north and 69°12′ longitude to the east (Ahmadi and Kajita 2017); the total area is 375 km 2 .
Kabul has a continental climate. Meteorological data have been extremely scarce in recent years; however, the data recorded from 1956 to 1983 by the World Meteorological Organization (WMO) shows a mean minimum temperature of −7.1°C in January and a mean maximum temperature of 32.1°C in July (Broshears et al. 2005). Three main rivers flow through Kabul: the Kabul River, the Logar River, and the Paghman River. Regarding the hydrogeological setting, shallow groundwater represents the main source of water supply in Kabul. The basin is situated at the intersection of three major fault systems of partially translational and extensional character. It comprises three interconnected aquifers, 20-70 m thick, consisting of coarse sandy to gravely detritus originating from the surrounding mountains. The aquifers were deposited by the three rivers flowing through the basin. The coarse aquifer material implies a high permeability. Deeper parts are affected by cementation of pore spaces, resulting in the formation of semidiagenetic conglomerates, causing decreased pore space and interconnectivity and consequently lower well yields. The main groundwater recharge occurs after the snowmelt from direct infiltration from the rivers (Houben et al. 2009). The geographical location of study area is shown in Fig. 1.

Materials and methods
For this study, existing hydrogeological and relevant data of soils, geological/lithological units, structural features, geomorphologic, and climatic conditions of the study area were collected. The overall study concept involved integration of thematic layers (six criteria) of soil type, drainage density, slope, thickness of the unsaturated zone, electrical conductivity and surface infiltration rate using ArcGIS software and application of multi-criteria decision analysis approaches, including Boolean logic, fuzzy logic, weighted overlay, AHP and ANP to find suitable sites for the application of MAR in the study area.
In this study, six criteria were selected by considering the expected scale, accuracy of the work, site conditions, the influence of each criteria, and adequate availability of information. These criteria were selected after reviewing the related literature and checking the availability of data for the study area. According to the literature, these six criteria are commonly used and are the most important criteria for MAR site selection. On the other hand, there were major restrictions regarding the availability of data for other possible criteria in the study area.
ArcGIS software, which supports both vector and raster formats, was used for the calculations. ArcView 10.6 software was employed for visualizing and weighting input maps and for classification and integration of the maps and preparing output maps. Table 1 provides the details on data sets used in this study along with source and period of the data. The different steps regarding the methodology in this study are presented in Fig. 2.

Drainage density
The drainage density is a measure of the length of stream channel per unit area of the drainage watershed (Magesh et al. 2012). Normally, drainage density is computed by dividing the distance of the stream into area coverage, which will generate a drainage density value. Drainage density is an inverse function of permeability (Carlston 1963). A high drainage density value indicates favourability for runoff and therefore indicates a low recharge groundwater potential zone. The drainage density (km/km 2 ) was calculated using Eq. (1): where Li denotes the total length of drainage of all streams in kilometers, and A denotes the area of the study area in square kilometers (Maniar et al. 2019).
The drainage density map of the area was prepared by spatial analysis tools in ArcGIS using a digital elevation model (DEM) of the area and classified into four classes considering suitability for groundwater recharge purposes: 0-0.5 (very suitable), 0.5-1.5 (suitable), 1.5-2 (moderate) and >2 (unsuitable). Generally, the drainage density of the area varies from 0 to 4.6 km/km 2 . It is worth noting that drainage density was included in this study's multi-criteria decision analyses, but this criterion has little impact on the selection of suitable recharge sites because the study area is an urban region, and urbanization changes the

Placeholder TextSlope
The slope is an important criterion for the classification of potential groundwater recharge zones. A higher degree of slope results in rapid runoff and increased erosion rate with insufficient recharge potential (Magesh et al. 2012), where the alluvial plain, flood plain, or plateau areas are more favorable for groundwater occurrence due to a longer duration of travel time to downstream and provide adequate time for infiltration to increase groundwater recharge. The slope map for the study area was created using the spatial analysis tool in ArcMap and digital DEM data with a cell size of 30-m resolution and pixel depth of 16 bit.
According to the previous studies, area with slope more than 10% is not suitable to be considered for groundwater recharge; therefore, the map was categorized in two classes using a Boolean logic method as: 0-10% (suitable area) and more than 10% (unsuitable area). In other multi-criteria decision methods and techniques applied in this study, the suitable range (0-10%) was further divided in three more classes. In the study region, slope varies from 0 to 31% and therefore the complete slope map is divided into four classes (0-2%, 2-5%, 5-10%, >10%). The major part of the study region falls under the class of (0-2% and 2-5%) but some higher-slope areas can be seen around the hills inside the study area.

Thickness of unsaturated zone
In areas where the water table is close to the surface and groundwater is not being exploited, recharge will cause the water table to rise, and it will become marsh land. In other words, if the water storage capacity of the natural underground reservoir is small, there is lower certainty of successful recharge (Mahdavi et al. 2010); therefore, thickness of unsaturated zone (distance from ground surface to the water table) is an important criterion to find suitable sites for groundwater recharge.
In this study, water-table depth measurements were collected from 100 observation wells. Then an interpolated map of the unsaturated zone thickness was produced using the inverse distance weighting (IDW) method. This map was used to show the variety of unsaturated zone thicknesses in the area of study. According to previous studies, groundwater recharge is unlikely to be successful where the thickness of the unsaturated zone is less than 10 m, but it can be considered a suitable site where the thickness of the unsaturated zone is greater than 10 m. The map was classified into four classes with respect to suitability for groundwater recharge as: very suitable (>30 m), suitable (20-30 m), moderate (10-20 m) and unsuitable (0-10 m).

Soil type/texture
Types of soil play an important role in the infiltration of surface water to recharge the groundwater in the area. The impervious clayey soil causes more runoff, whereas the sand and gravelly soil causes more infiltration (Maniar et al. 2019). A soil map of study area was prepared based on the results of soil investigations conducted in different locations of Kabul metropolitan area. In total, soil data of 110 boreholes were collected from private and governmental organizations and used for preparation of the soil map. Since only the surface soil type was required in this study, data of soil from shallow depths (less than 3 m) was extracted from the reports. The interpolation using the IDW method was performed on soil type data and then a classified soil map was generated which consisted of the four main soil classes found in the area of study-siltyclay, loam, sandy loam, and sand. The ranks of soil are assigned according to their degree of infiltration. The sandy soil has a high degree of infiltration and thus has a higher priority, while the silty-clay soil has the lowest degree of infiltration and thus low priority for groundwater recharge (Rajasekhar et al. 2019).

Surface infiltration rate
The ground surface infiltration rate is vital for groundwater recharge. The main recharge source for groundwater is water infiltration from rainfall and surface water. The areas with a high ground infiltration rate can be recognized as suitable sites for groundwater recharge. Based on the field measurements, no related data for these criteria were available in this study; therefore, the rate of surface infiltration was obtained according to the soil texture-infiltration relationships established by the Food and Agriculture Organization (FAO 1979) in the study area (Ghayoumian et al. 2007; see Table 2 for presentation of the relationship).
A value for the infiltration rate was assigned to each soil data location based on the correlation suggested in Table 2. A

Placeholder text quality of groundwater
The quality of groundwater reveals the amount of chemicals and biological impurities present in groundwater and is a major criterion in specifying water for certain uses (Nasiri et al. 2013). In the present work, electrical conductivity (EC) was employed as a criterion for assessing water quality from a salinity point of view. Since total dissolved solids (TDS) and EC showed the same trend of change, only EC was used as the water quality index in this study. In order to generate the water quality map, the EC data layer constructed from the data from 54 observation wells was utilized. Salinity classification was used to divide the area into four classes on the basis of EC values. For groundwater recharge suitability, classification was done as follows: <1,000 μS/cm (very suitable), 1,000-2,000 μS/cm (suitable), 2,000-3,000 μS/cm (moderate), >3,000 μS/cm (unsuitable).
The thematic map of EC shows that the study area is characterized by suitable to very suitable, due to the presence of low EC values. High EC values are only observed in two small parts of the area: one in the northern Kabul subbasin and the second near the intersection of the Logar River and the Kabul River. Figure 3 illustrates the classified thematic maps for the effective six criteria applied in this study. In addition to the six effective criteria previously discussed, a land-use/landcover map (prepared by the UN-Habitat team) was used in this study, which included 12 classes: agriculture, barren land, buildings under construction, green areas, institutional, residential, road streets, transport, industrial, vacant plots, and water bodies. The amount of infiltration depends on the land-use pattern of the area. It is well understood that vegetation increases infiltration, whereas roads, concrete pavements, and buildings decrease infiltration (Maniar et al. 2009). Landuse/land-cover is a significant indicator of groundwater status and use of groundwater, as well as an important indicator within the reach of areas with groundwater potential. Undertaking artificial recharge projects is feasible in areas with an appropriate density of vegetation coverage because these regions not only recharge water into aquifers but also prevent surface soil erosion. In addition, barren land and vacant plots can be suitable areas for a MAR projects. The landuse/land-cover map was used as a filter layer on the final suitability map to remove areas with land-use restrictions. Only barren land, agriculture and vacant plots were considered suitable land types for MAR application.

Multi-criteria decision analysis methods
Boolean logic is algebra that works with logical relationships rather than numeric relationships. The variables have a logical value of TRUE or FALSE. Boolean logic is an alternative way to find suitable locations for artificial recharge (rather than creating a suitability map) to query the required data. Once all the needed datasets (the thematic layers) have been created, such a query would be used to find all the suitable locations. Probably the simplest and best-known type of GIS model is based on Boolean operations. Robinov (1989) introduced the use of Boolean operations for reasoning with geological maps. In effect, Boolean modeling involves the logical combination of binary maps resulting from the application of conditional operators (Bonham-Carter 1996). Only one or zero values are assigned to each unit area, specifying whether it is satisfactory or unsatisfactory, respectively. The Boolean model consists of AND and OR operators. Based on set theory, the AND operator yields the logical intersection of the two data sets, and the OR operator obtains the logical union of the two data sets. The query data are analyzed based on Boolean logic and entered in the Raster Calculator (from spatial analyst) according to the existing thematic layers as follows: [ "==" means equal, "|" means OR, and "&" means AND).
In Boolean overlay, operation criteria are considered constraints and not factors (Drobne and Lisec 2009). The constraints are assigned a value of either 1 or 0 that are defined as true or false, respectively. Thereafter, they can be combined through either AND or OR operators (Ghayoumisan et al. 2007). To get an output value of 1, the AND operator requires all the input values to be 1. The OR operator only requires one of the input values to be 1 to get an output value of 1 (Bonilla Valverde et al. 2016). In this study, the AND operator was used to exclude unsuitable areas for the application of MAR. Table 3 shows the weight of each effective criterion based on its range, which is divided into two categories: unsuitable (value 0) and suitable (value 1).
Weighted overlay method is one of the most used approaches for overlay analysis to solve multi-criteria problems such as site selection and suitability models. Since the input criteria layers will be in different numbering systems with different ranges, to combine them in a single analysis, each cell for each criterion must be reclassified into a common preference scale such as 1-10, with 10 being the most favourable. An assigned preference on the common scale implies the phenomenon's preference for the criterion. The preference values are on a relative scale. That is, a preference of 10 Thematic maps for six effective criteria: a slope, b surface infiltration rate, c unsaturated zone thickness, d electrical conductivity, e drainage density, f soil type is twice as preferred as a preference of 5. The steps for running the "weighted overlay" tool are as follows: 1. Select an evaluation scale 2. Add input rasters 3. Set scale values 4. Assign weights to input rasters 5. Run the weighted overlay tool A "weighted suitability model" has been developed using GIS techniques for proposing locations suitable for applying groundwater recharge depending on a number of thematic layers and based on the principle of multi-criteria decision analysis. Such models are used for applying a common measurement scale of values to diverse and dissimilar inputs in order to create an integrated analysis. Additionally, the criteria for the analysis may not be equally important. Each individual raster cell is reclassified into units of suitability and multiplied by a weight to assign relative importance to each and finally add them together for the final weight to obtain a suitability value for every location on the map; this can be interpreted by Eq. (2) (Eastman 2001).
where w i = the weight of ith criterion map; x i = score of class of criterion i and S = suitability index for each pixel in the map. In this method, six thematic maps (slope, thickness of unsaturated zone, surface infiltration rate, drainage density, soil type and electrical conductivity) were applied. Each thematic map was classified into four classes with scores of 1, 2, 3 and 4. The score of 1 was assigned for unsuitable classes of the criteria, score 2 for moderate, score 3 for suitable, and finally score 4 for very suitable classes based on effectiveness on groundwater artificial recharge. This method reclassified the range of criteria for the suitable class in Boolean logic into three classes-very suitable, suitable, and moderately suitable. The classified maps were integrated through a weighted overlay tool in the spatial analysis tools of the ArcGIS software. The integrated map uses four classes of suitability for MAR as unsuitable, moderately suitable, suitable and very suitable.
The total weights of each pixel of the final integrated layer were derived from the following Eq. (3): where SL is slope, DD is drainage density, SI is surface infiltration rate, UZ is unsaturated zone thickness, EC is electrical conductivity and ST is soil type. The subscript letter w represents the weight of each criterion, while x represents the score of each class of the individual criteria. Thus, the artificial groundwater recharge index S, a dimensionless quantity that aids in indexing the likely groundwater recharge zones in the area, was estimated using Eq. (3) for each pixel in the final integration layer and regrouped into different four classes with equal class intervals to divide the entire study area into different suitability zones (Chowdhury et al. 2006).
It should be noted that each raster is assigned a percentage of influence according to its importance. The weight is a relative percentage, and the sum of the percentage influence weights must add up to 100. The influence weights in the present study were not assigned for the effective criteria and the influence of all criteria on the site selection for MAR was considered equally. The scores of each class for each criterion are shown in Table 4.
Fuzzy logic method: A fuzzy set is a collection of degrees of membership. Membership of a fuzzy set is expressed on a continuous scale from 1 (full membership) to 0 (no membership at all). Fuzzy functions can be used to separate the maps of a few classes. Then each class is given a membership degree based on their impact on the suitability of the area for groundwater recharge, in the range (0, 1). There are five operators to combine the membership values together: namely AND, OR, algebraic product, algebraic sum and gamma (Mahdavi et al. 2010). Fuzzy logic is generally based on two principles: (1) classifying maps compared together based on importance, (2) weighting classes in every map between 0 and 1 (Ghayoumian 2007). In fuzzy algebraic products, the combined fuzzy membership values tend to be very small with this operator, because of multiplying several numbers with less than 1. The output is always smaller or equal to the smallest contributing membership value and therefore reduced (Bonham-Carter 1996). In this study, the fuzzy algebraic product operator has been used because of its high sensitivity in specifying artificial recharge areas. In fuzzy algebraic products, the combined membership function is defined as Eq. (4) (Tiwari et al. 2017): where μ i is the fuzzy membership function for the ith map, and i = (1, 2, 3 … n) maps are to be combined. Basically, the fuzzy logic system in this study consists of the following steps (Gesim and Okazaki 2018): -Fuzzification: converting the crisp inputs to membership functions which comply with intuitive perception of the system status. -Rules processing: calculating the response from system status inputs according to the predefined rules matrix (control algorithm implementation).
-Inference: evaluating each case for all fuzzy rules.
-Composition: combining information from rules.
-De-fuzzification: converting the result to crisp values.
By using fuzzy logic in GIS, aspects of ambiguity in linguistic variables can be modeled (Benedikt et al. 2002). The fuzzy inference system can be developed in four stages: (1) the definition of linguistic variables, (2) the selection of appropriate membership functions, (3) the definition of rules based on knowledge of the system and (4) the selection of a suitable operator to combine fuzzy sets. The fuzzy inference method was adapted to pixel-based site selection for MAR in this study. The output membership functions in this method are fuzzy values that are de-fuzzified to obtain a crisp output value for each pixel. Raster layers that are prepared in a GIS environment are the input and output of the model. In the output layer, each pixel is assigned a value demonstrating its degree of suitability for MAR (Malekmohammadi et al. 2012). Once the appropriate fuzzy membership value for data criteria is assigned, several reclassified surfaces showing a value from 0 to 1 are generated (ESRI 2014). The next step in applying fuzzy logic is to overlay these surfaces. To complete this step, one of several fuzzy overlay types (AND, OR, Product, Sum and Gamma) must be chosen (ESRI 2014). The maps are then integrated through a fuzzy operation "Product" to model the suitability of areas for MAR application in the study area. The six effective criteria are weighted quantitatively and classified based on the fuzzy logic results in Table 5. Drainage density (km/km 2 ) R a n g e 0 -0.5 0.5-1.5 1. In this study, to prepare a standard fuzzy set as an input, the Gaussian function was used to convert crisp values according to Eq. (5) (Gesim and Okazaki 2018): where, the parameter a is the height of the curve's peak (here in fuzzy logic it is 1), b is the position of the center of the peak and c the standard deviation. After processing input in the inference engine as illustrated in Fig. 4, the result of the process must be converted to a crisp value. By the composition of the crisp values, the membership curves for different classes of thematic layers will be prepared.
Raster layers that are prepared in a GIS environment are the input and output of the model used in this study, and, in the output layer, each pixel is assigned a value demonstrating its degree of suitability for MAR. Thematic raster layers for the criteria were prepared, classified, and weighted. The classification and relative weighting of the thematic layers were based on the literature review, engineering judgments, and expert opinions. Finally, a raster layer is produced as a result. Each pixel in this layer has a value (from 0 to 1) that indicates its suitability for MAR. This output layer is classified into categories such as unsuitable, moderately suitable, suitable and very suitable.
Analytical hierarchy process, first developed by Professor Thomas L. Saaty in the 1970s (Pinto et al. 2017), is a widely  . 4 Structure of the fuzzy logic model (Malekmohammadi et al. 2012) used MCDA technique in the field of water resource engineering (Pinto et al. 2017). Generally, for the AHP calculation, the element of each layer class of individual raster was obtained according to Saaty (1980) which later will be used to calculate the pairwise comparison within the matrix. The pairwise comparison matrix (PCM) was then used for calculating the normalized weight for each thematic layer; thereafter, consistency was calculated to detect the coherence of the judgment matrix. Ultimately, the class weight and consistency ratio (CR) was detected. The steps are as follows: Step 1: Pairwise comparison matrix (PCM) was calculated as Eq. (6): where X nn = indicator of the comparison matrix element.
Step 2: Calculation of normalized weight was determined as Eq. (7): where GM n = geometric mean of nth row of comparison matrix (X), N f = number of factors; GM n was calculated by Eq. (8): Step 3: Consistency ratio (CR) to verify the coherence of the comparison Eq. (9).
The AHP is a concept of measurement by PCM and is based on the judgment of experts to derive priority weights. The PCM is calculated using a weight scale of 1-9, indicating how frequently one map is more important than the others (see Table 6). Based on the PCM, the relative weight matrix and the normalized principal eigenvalue were calculated to determine the percentage of effect of the thematic layers and the classification of the constraints (Rajasekhar et al. 2019).
The integration of the thematic layers depends on the correlation between the layers and their relative significance for the preparation of a suitability map of groundwater recharge. In this study, a 6 × 6 PCM was designed based on the selected criteria, which included slope, surface infiltration rate, thickness of the unsaturated zone, soil type, drainage density, and electrical conductivity. The primary principle for checking the accuracy of comparison is that a consistency ratio (CR), where the value that is not exactly or equivalent to 0.1, shows a satisfactory reciprocal matrix and the ratio above 0.1 demonstrates that the PCM should be changed (Rajasekhar et al. 2019); additionally, the evaluation procedure must be repeated to improve the consistency of the pairwise comparison between the criteria. CR is computed as follows: where CI is consistency index (given from Eq. 10); RI is random consistency index where λ max is principal eigenvalue, n is number of comparisons. At the end, RI can be found from Table 7 based on the number of criteria involved in the analysis.
In this study, AHP analysis and criteria weighting were conducted with the help of Expert Choice (version no. 11) software. Using Expert Choice software, the weightings of all criteria and subcriteria were determined. All layers were analyzed in the ArcGIS environment and converted to raster data. Finally, with the help of ArcGIS software and the Raster Calculator toolbox, thematic layers were overlapped, and a suitability map was obtained and classified.
Analytical network process is one of the most recent MCDA techniques, which was proposed by Saaty (2001), and can be considered as a more recent extension of AHP for decision making, with dependence and feedback that can handle a more complex decision structure (Saaty 2001). AHP is limited to relatively static and unidirectional interactions with little feedback among decision components and alternatives. While AHP has been very popular, ANP is less prominent in the literature (Othman et al. 2011). The basic structure is an influence network of clusters and nodes contained within the clusters. Priorities are established in the same way as in the AHP using pairwise comparisons and judgment. Saaty (2001) suggested the usage of AHP to solve the problem of independence on alternatives or criteria. In addition to AHP, the ANP technique is a general form that allows interdependencies, outer-dependencies, and feedbacks among decision elements in the hierarchical or nonhierarchical structures (Görener 2012). The structural differences between a hierarchy and a network process are pictured in Fig. 5.
The ANP approach comprises four steps (Satty 1996;Chung et al. 2005;Yüksel and Dağdeviren 2007): step 1: model construction and problem structuring, step 2: pairwise comparisons and priority vectors, step 3: super matrix formation, step 4: synthesis of the criteria and the alternative's priorities and selection of the best alternatives.
In this study, the ANP model was established in Super Decision 2.4.0 (2015) software and all the preceding four steps were followed. At the first problem, structure, including clusters and nodes, was created. In the second step, the PCM was arranged exactly according to the matrix from the AHP model because the comparisons of the importance between the criteria are based on the same judgment (see Table 7). Therefore, a matrix of 6 × 6 was produced for all six effective criteria and then the weighting of the criteria was calculated based on the interactions and dependencies between the clusters and nodes of the network model. In the final step, the synthesis of the criteria was done and the weight of each criterion was obtained accordingly (See Table 8).

Results and discussion
In this study, the city of Kabul, with an area of about 375 km 2 , has been investigated for MAR site selection. Based on the literature review and available data for this area, effective criteria have been selected. To determine the most suitable areas for MAR in the study area, electrical conductivity of groundwater, unsaturated zone thickness, slope, drainage density, soil type, and surface infiltration rate have been considered. For each criterion, the raster thematic layer with a pixel size of 30 × 30 m was prepared in a GIS environment; additionally, the land use factor was used to create the constraint layer on the suitability map. Among the different land uses, agriculture, barren land, and vacant plots are practicable for MAR. In order to select the suitable areas for MAR in the area, different methods and techniques including Boolean logic, weighted overlay, fuzzy overlay, AHP and ANP were applied and weighting the effective criteria was done accordingly.
Boolean logic analysis was the simplest and best-known method in comparison with other methods that were used in this study. The result is a Boolean TRUE or FALSE map for the locations that meet or do not meet the given criteria. In  order to prepare the map, each criterion was divided into two classes-unsuitable and suitable, based on literature review, engineering judgment and range of data. Then, value 0 was assigned for unsuitable class and value 1 was assigned for suitable class of the criteria (see Table 3). The result of the Boolean model for identifying suitable locations for MAR is shown in Fig. 6. Areas in red are not suitable and areas in green are suitable locations. The difference between querying the data based on Boolean logic and creating a suitability map using the weighted model is that, in the Boolean TRUE or FALSE map, there are no areas of medium suitability. If the analysis needs more flexibility, a suitability map should be created, where every location (cell) has a suitability value. As shown in the map, suitable areas are distributed across the study area where the criteria have a suitable range with a true (1) value. In this method, all of the criteria have equal importance on site selection, therefore the effectiveness of each criterion is the same as others for the analysis which is not an accurate enough approach. This map, on the other hand, only divides the area into two categories-unsuitable and suitable-presenting an overall approach to the site's suitability.
In fuzzy logic analysis, fuzzy membership values were given to different features of effective criteria according to their importance in groundwater recharge. All the maps were integrated by overlay techniques using GIS to delineate suitable areas for MAR. The surface infiltration rate of more than 45 mm/h is assigned with fuzzy membership value of 0.95 and has very high priority. Similarly, the fuzzy values of 0.75 and 0.35 indicate comparatively suitable and moderately suitable, respectively, and for the unsuitable class the fuzzy value of 0.01 is assigned. Likewise, the slope criteria ranging from 0 to 2% received a fuzzy membership value of 0.7, indicating greater recharge potential, whereas higher slopes of 2-5, 5-10, and more than 10% received fuzzy membership values of 0.5, 0.3, and 0.01, respectively. The final weights of all effective criteria are shown in Table 5. Thematic maps for all effective criteria were illustrated based on fuzzy membership in Fig. 7.  The suitability map generated using the fuzzy overlay method is as shown in Fig. 8b. The area with very suitable potential is located in the west, similar to the suitability map of the weighted overlay method. Generally, the map prepared from the fuzzy overlay method has a significant similarity to the generated map from the weighted overlay method and there is an excellent correlation between these two maps.
In the weighted overlay method, each thematic layer was prepared and classified based on scores of 1-4: 1 for unsuitable, 2 for moderate, 3 for suitable and 4 for very suitable class (see Table 4). Then, all thematic layers were combined using the weighted overlay tool in the ArcGIS software with respect to their degree of influence which was assigned equally for all criteria in this study. The final suitability map from the weighted overlay analysis is illustrated in Fig. 8a. This map consists of four classes that range from unsuitable to very suitable. This method has a good correlation with the suitability of the Boolean method, but with a more detailed classification. Results show that only a limited area in the west part is categorized as a very suitable area, and suitable and moderate classes can be found in different parts of the study area, mostly in the central and northern parts.
In the AHP method, the geospatial analysis of the thematic layer was performed using weighted overlay analysis and the resultant map was classified based on the reclassifying method. Each layer has a subcriterion weight. After determination of weights for the criteria, a PCM was developed for weighting the subcriteria. Subcriteria are classified into four classes as defined in the weighted overlay method. Each class was compared to others based on its importance to groundwater recharge and weighted according to Saaty's scale and then the normalized weight of each class (subcriteria) was obtained accordingly.
Because there were six criteria in this study, the random consistency index (RI) was obtained as 1.24, as shown in Table 7. Consequently, the consistency ratio for the pairwise comparisons matrix is calculated as 0.05; therefore, a PCM is acceptable and evaluations are sufficiently consistent in comparing the importance of the effective criteria. Table 9 shows the weightage of all criteria and subcriteria based on the AHP technique. According to this table, the weightage of surface infiltration rate, soil type, slope, thickness of the unsaturated zone, electrical conductivity and drainage density were 34, 26.5, 19.9, 12.1, 4.4, and 3.1%, respectively. For the selection of suitable sites in the study area, the weighting of effective criteria reveals that surface infiltration rate is the most important and drainage density is the least important.
After overlapping thematic layers of all criteria, the final suitability map was divided into five classes-unsuitable, low suitability, moderate suitability, suitable, and very suitable. The final suitability map produced from the AHP model is as shown in Fig. 8c. According to this suitability map, 5, 20.6, 42.1, and 24.3% of the total area are unsuitable, low suitability, moderately suitable and suitable, and 7.9% are very suitable.
As mentioned in section 'Materials and methods', in the ANP model, the same PCM was applied as designed in the AHP model. Table 9 shows the weightages of the effective criteria calculated using the ANP method. According to the ANP method results, the most important criterion for selection for MAR sites is soil type with a weight of 28.24%, followed by surface infiltration rate (25.39%), slope (23.28%), drainage density (10.8%), thickness of the unsaturated zone (9.02%), and electrical conductivity.
The suitability map generated based on the results of the ANP model is presented in Fig. 8d. According to this map, 3.7% of the total area is unsuitable, 15% is of low suitability, 37.4% is moderately suitable, 33.1% is suitable and 10.3% is very suitable for MAR.
As seen from Table 9, different normalized weights for the criteria have been achieved for the AHP and ANP methods. Water can better infiltrate the ground surface if there is an even/gentle-slope surface in comparison with a steep and rough surface (Nouayti et al. 2019). Furthermore, infiltration is faster in gravel and sand soil than in silty clay soil, and infiltration is slower in areas with dense drainage than in areas with lower drainage density. Hence, slope, soil type and drainage density affect the rate of surface infiltration and it can be concluded that the surface infiltration rate has a relationship with these three criteria. This interdependency between criteria was applied to the network structure in the ANP method. As a result, the normal weight of the effective criteria obtained by the AHP and ANP methods differ. The highest weight in the AHP method is assigned to the surface infiltration rate but the highest weight in the ANP method belongs to the soil type. On the other hand, drainage density has the lowest weight or lowest degree of importance in the AHP method, but the lowest weight (lowest degree of importance) in the ANP method was assigned to electrical conductivity. This occurred because the AHP method does not allow one to account for a relationship and the interdependence of criteria. Figure 9 shows the differences between the weightage of the AHP and ANP methods in a bar chart, whereby it can be concluded that when criteria have interdependency, the ANP method gives a more accurate and reliable result in comparison to AHP for site selection purposes (Singha and Pasupuleti 2020). It should be mentioned that, because of the difference between the weightage of criteria from the AHP and ANP methods, the suitability maps obtained from these two methods have different percentages for the areas of suitability classes (Fig. 10). As shown in Fig. 10, the unsuitable class covers 5% of the total area in the AHP map but in the ANP map, the unsuitable class covers 3.7% of the total area. On the other hand, 7.9% of the total area is recognized as very suitable area in the AHP map but 10.3% of the total area is very suitable in the ANP suitability map.
Land-use plays an important role in selecting suitable sites for artificial recharge and evaluating the feasibility of whether any structure could be constructed on the selected site. In this study, barren land, vacant plots and agriculture areas have been selected as appropriate land-use types for MAR application (Fig. 11a). Planning for MAR sites may face less land acquisition issues for these types of land.
Therefore, the area's land-use map was overlapping by the ANP method's suitability map as a filter for deleting areas with a land-use restriction. For this purpose, very suitable areas obtained from the ANP method were clipped from the generated suitability map. And, barren land, vacant plots, and agricultural areas were clipped from the land-use map. To find the final suitable areas for MAR application, both clipped layers were overlaid together in a GIS environment and then seven prioritized MAR suitable sites, named A, B, C, D, E, F and J, were assigned considering appropriate land use (Fig. 11b).
The seven prioritized sites for MAR identified in this study were compared with the proposed MAR sites from the current pilot project (Fig. 12). From this comparison, it can be found that three of the prioritized sites (B, D and E) are matched with the MAR sites of the current pilot project (sites 3, 4 and 6). In addition, prioritized sites for MAR were assessed through site visits, which confirmed their suitability and the results of this study; therefore, MCDA techniques integrated with GIS are an effective method for the selection of MAR sites.

Conclusions and recommendations
This study was carried out to identify the most suitable areas for managed aquifer recharge (MAR) application based on various multi-criteria decision analysis (MCDA) techniques. The overall study concept involved the integration of six thematic layers of soil type, drainage density, slope, thickness of the unsaturated zone, electrical conductivity, and surface infiltration rate. The multi-criteria decision analysis methods used were Boolean logic, weighted overlay, fuzzy logic, AHP and ANP.
The Boolean method's output suggests two classes for the area's suitability for MAR: unsuitable and suitable. The detailed classification of the suitable areas was obtained through fuzzy overlay and weighted overlay methods. The suitability map of these two methods has a significant similarity in identifying the most suitable area by classification through four According to the fuzzy logic and AHP method, surface infiltration rate is the most important criterion and drainage density is the least effective criterion for the selection of MAR sites in the study area. On the other hand, based on the ANP method, soil type is the most effective criterion and electrical conductivity is the least effective criterion for MAR site selection.
Generally, a very suitable area for MAR is located in the western part of the study area in the upper Kabul subbasin close to the upstream of Paghman River. Predominant land uses are barren, vacant and agriculture land in the eastern part of the study region.
Suitable areas for application of MAR are located in the southwest, western, northern, and central parts of the study area but the larger portion of suitable areas is in the western and central parts. Based on the surface infiltration map and soil type map, these areas mainly contain alluvium and colluvium deposits with high permeability characteristics and are located adjacent to the Paghman and Kabul rivers.
The results of the present study can serve as guidelines for planners, decision makers and hydrogeologists for planning future artificial recharge projects in the study area in order to ensure dependable water supply and sustainable groundwater utilization on a long-term basis. The databases and the results obtained can also be useful in developing a conceptual model of the study area for future modeling studies. In conclusion, the integration of MCDA and GIS is recognized as an effective method for modeling the uncertainties and reducing the errors that are associated with data classification in GIS for selection of MAR sites.  Fig. 10 Percentage of suitability classes based on AHP and ANP methods Fig. 11 a Land-use/land-cover map, b prioritized sites for managed aquifer recharge (MAR) Future work can concentrate more on the suitable sites proposed in this study, as well as include and evaluate more effective criteria for the environment, economic and social factors, and water-source availability and accessibility, which have not been covered in this study.