Appraising groundwater potential zones using geospatial and multi-criteria decision analysis (MCDA) techniques in Andasa-Tul watershed, Upper Blue Nile basin, Ethiopia

An integrated approach of remote sensing, geographic information system and multicriteria decision analysis of analytical hierarchical process (AHP) were applied to delineate groundwater potential (GWP) zones in Andasa-Tul watershed, Upper Blue Nile Basin, Ethiopia with an area of 1872 km2. Nine GWP influencing thematic layers comprising lithology, lineament density, geomorphology, slope, soil, drainage density, land use/land cover, rainfall and depth to groundwater level prepared from conventional and satellite data were used. The thematic layers and classes within them were given scale values based on literature and experts’ decision by using Satty’s AHP. The thematic layers have been integrated via their weights/rates using weighted overlay spatial function tool of ArcGIS to provide GWP map. The result shows that GWP map comprises 250.8 km2 (13.4%) very good, 132.9 km2 (7.1%) good, 441.8 km2 (23.6%) moderate, 662.7 km2 (35.4%) poor and 383.8 km2 (20.5%) very poor zones. Validation of the GWP map with existing boreholes and springs yield shows 83.8% agreement indicating good accuracy of the method. The map removal sensitivity analysis result reveals that GWP is more sensitive to lithology (mean variation index, 1.92%) and less sensitive to geomorphology (mean variation index, 0.59%). Similarly, from the single layer sensitivity analysis, lithology and slope are found to be more effective parameters, whereas rainfall and depth to groundwater level are less effective variables. The generated GWP map in the study region can be used as important guideline for locating GWP areas for further exploration, planning and management of the groundwater resources.


Introduction
Water is one of the indispensable natural resources for all living things. Without adequate and incessant source of water, life cannot exist in our planet. The survival of 70% animal and 90% plant species is dependent on water (Gilli et al. 2012). This resource is not only used for normal functioning of life but also important for other numerous applications like socioeconomic development (Olea-olea et al. 2020) food security and ecosystem sustainability (Kebede 2013). Recently, different countries of the world have faced water scarcity problem as a result of rapid increase in population, industrialization, urbanization and mismanagement of the resources (Manap et al. 2012;Singh et al. 2017). Despite the availability of water in abundance on earth, it is dominated by the saline water (97.5%). As stated by Jenifer and Jha (2017) fresh water consists only about 2.5% of the water reserve of the world out of which 30% is stored as groundwater (Karamouz et al. 2011). Since two-third of the fresh water is frozen and stored as glacier (Oki and Shinjiro 2006;Elbeih 2015), groundwater consists of 96% of unfrozen and readily available global fresh water (Downing 2014).
Presently, life and any development activities in different countries have become strongly dependent on existence of groundwater resource (Kebede 2013;Erdogan et al. 2019) unlike surface water that can easily be affected by drought (Nofal et al. 2018) and pollution (Diamantini et al. 2018). Groundwater has the advantage of having a wide area of occurrence, buffer climate variability, better natural protection, cheaper development cost, ecological sustainability and demand based application flexibility at field condition (Hancock et al. 2005;Karamouz et al. 2011;Kebede 2013;Erdogan et al. 2019).
To reduce extreme conflict between demand and supply, people are always forced to explore groundwater sources usually using conventional hydrogeological studies to assure water supply reliability. However, climatic variability and complexity of the aquifer systems in Ethiopia make conventional groundwater exploration techniques (for instance geological, hydrogeological, geophysical, etc.) troublesome. Since groundwater is a concealed resource (Jenifer and Jha 2017;Patra et al. 2018), and its occurrence and distribution is controlled by multi-parameters (Fenta et al 2014;Murmu et al. 2019;Nithya et al. 2019), it is imperative to use advanced technology (Oikonomidis et al. 2015;Senanayake et al. 2016) to minimize the risk of missing groundwater potential (GWP) zones in an area of interest. Unlimited need of groundwater resources for various purposes necessitates to use remote sensing (RS), Geographic information system (GIS) and analytical hierarchical process (AHP) of the multi-criteria decision analysis (MCDA) for identification and delineation of GWP zones ahead of detail exploration of the resources using geophysical, hydrogeological and drilling techniques (Fenta et al. 2014).
The developments of RS and GIS technologies coupled with MCDA provide a new scientific inquiry in the area of hydrogeological studies. The RS method by its virtue of generating repetitive, cost effective, synoptic, feasible spatio-temporal and spectral data of extensive area quickly, it becomes an invaluable technique in delivering pertinent information regarding the different variables controlling groundwater occurrence (Aluko and Igwe 2017).
GIS has grown swiftly as a powerful tool for resolving complex problems. It has the ability to integrate, visualize, analyze and manage digital geographic data from the various sources and become eminent working environment for complicated spatio-temporal data. As described by Gupta (2003) and Weng (2010), integration of RS and GIS provides a significant role in groundwater exploration. The purpose of integrating different groundwater occurrence controlling factors in GWP delineation is to attain a complete and clear picture of target area with better confidence and certainty.
MCDA is a technique used to simplify the problems of applying multiple parameters (criteria) by decision makers. It provides various ways in which decision issues can properly be structured and evaluated by decision makers for prioritization of different decision alternatives (Malczewski 2006). MCDA has gotten its worldwide acceptance as effective technique because of its role in dealing with intricate decision problems (Agarwal and Garg 2016). The commonly used MCDA for GWP area studies is Saaty's AHP . Saaty (1980) developed AHP as part of MCDA tool to solve complex decision problems. The problems are hierarchically structured in to goals, criteria, subcriteria and alternatives at different levels. AHP data are derived by a set of pairwise comparisons which in turn helps to obtain the weights of importance of the decision criteria. AHP is used for normalization of assigned weights/rates to thematic layers/classes ) and computation of eigenvectors, maximum eigenvalues and consistency ratios.
Integration of RS, GIS-based MCDA (AHP) combines all geospatial data (thematic layers) and value judgement and provides composite GWP model map. Different authors combined different thematic layers such as lithology, rainfall, drainage density, lineament density, geomorphology, land use/ land cover soil, etc. to produce comprehensive GWP maps. Irrespective of slight variation in the type and weight given to groundwater occurrence governing variables, some authors use integration of RS and GIS technology for GWP zone delineation without MCDA (Fashae et al. 2014;Gumma and Pavelic 2013;Jasrotia et al. 2016;Kumar and Pandey 2016;Magesh et al. 2012;Oikonomidis et al. 2015) and with MCDA (e.g., Abrams et al. 2018;Arabameri et al. 2019;Fenta et al. 2014;Ghosh et al. 2020;Jenifer and Jha 2017;Kaur et al. 2020;Nanda et al. 2017;Patra et al. 2018).
Groundwater exploration in Ethiopia is commonly carried out by conventional geological and hydrogeologic studies followed by geophysical survey irrespective of complexity of the aquifer system. Authors' field experiences of this paper in volcanic terrain have shown that different abandoned (unproductive) wells are commonly available even within the same lithological units. This shows that complexity of the aquifer system in volcanic terrain (Fenta et al. 2020) makes difficult to apply the commonly used conventional exploration methods. Application of RS, GIS-based MCDA for GWP is recent technology in the country and used to minimize problem of conventional groundwater exploration techniques. A few authors have recently used this technology to explore GWP for future sustainable use of the resources (Ahmad et al. 2020a, b;Andualem and Demeke 2019;Berhanu and Hatiye 2020;Fenta et al. 2014;Tolche 2020).
The Andasa-Tul watershed (the study area) is important development corridor for fast growing Bahir Dar City (capital of Amhara Regional State) and different rural villages. Due to population growth pressure, water insecurity is found to be an ever increasing challenge in Ethiopia (Ayenew et al. 2008) particularly in the study region (Bawoke et al. 2019). Bahir Dar, having well pronounced urban expansion, is currently undergoing city structural planning including (North and Northeast boarder of the study region) with different interest and issues of which long term sustainable use of groundwater is one of the most important targets. Different studies show that geospatial based GWP studies have been used in Ethiopia in recent times. Few of the scientific investigations have been commonly conducted in different parts of the country since a decade ago. Ahmad et al. (2019) and Ahmad et al. (2020a) have carried out GWP investigations by coupling GIS with AHP to delineate GWP areas in Wonbera and Beshilo river basins of Ethiopia, respectively. They have prepared and used more or less similar thematic layers (such as soil texture, soil type, soil depth, soil saturated hydraulic conductivity, slope, drainage density, land use, geology, fault, rainfall and stream power index) in the two basins at different times. Their results have shown that the soil texture (sandy loam) and slope (flat) play a vital role for the existence of high GWP areas in both studies. However, the two studies are dominantly dependent on soil thematic layers than others which lead them to conclude soil texture as important thematic layer. These studies also lack validation test to approve the groundwater potential zones. Similarly, Ahmad et al. (2020) has employed related study using fuzzy logic techniques in GIS domain in Jema river basin of Ethiopia using eight thematic layers (Geology, slope, drainage density, distance to fault, curvature, wetness index, soil texture and type and land use). In this study, soil texture (sandy loam soil) with high terrain wetness index value is found to be GWP promising area. Yet, there is no validation test for this study as ground truth evidences to confirm GWP areas.
Integrated approach of RS, GIS coupled with MCDA (AHP) have been used by (Fenta et al. 2014;Andualem and Demeke 2019;Berhanu and Hatiye 2020) in Raya Valley, Guna-Tana landscape and Megech watershed of Ethiopia, respectively. They used seven to nine thematic layers (lithology, lineament density, slope, geomorphology, elevations, drainage density, soil, landuse land cover and elevation) to delineate GWP zones in their respective study areas. They have validated GWP zones with existing groundwater scheme yield data. In addition to the aforementioned thematic layers, depth to groundwater level plays important role for affecting GWP zones. However, all the above authors and any other authors in Ethiopia do not have used depth to groundwater level as a thematic layer for GWP studies using geospatial and multicriterial approaches. Furthermore, there is no any GWP study in the area under investigation using integrated RS, and GIS based multicriterial approaches. Therefore, the main purpose of this study is to delineate GWP zones in Andasa-Tul watershed by employing integrated RS and GIS-based MCDA techniques and validate the result via field observed groundwater sources (spring and borehole yields). The specific objectives are: (1) to prepare different thematic layers that can affect groundwater occurrence, (2) to carry out map removal and single layer sensitivity analysis and identify most sensitive and effective thematic layers for the provision of GWP zones. The study is indispensable for sustainable use of the groundwater resource by applying suitable groundwater management practices. Such study also paves the way for conducting detail geophysical and hydrogeological studies at ease. Thus, the method optimizes locations for future sustainable groundwater development areas.

The study area
The Andasa-Tul Watershed, which is found in northwest part of Ethiopia (Amhara regional State), is part of Upper Blue Nile Basin. The watershed covers an area of 1872 km 2 and lies in between Adama Mountain (southwest) and the Blue Nile River (northeast). Geographically, it is located between 37°15′26ʺ-37°46′25ʺ east longitude and 11°3′54ʺ-11°36′1 ʺ north latitude (Fig. 1).
The watershed is characterized by a thick succession of Cenozoic volcanic rocks of Oligocene-Miocene flood basalts (the Trap sequence), Quaternary basalt flows and scoria (Asrat 2017). The study region is depicted by predominantly undulating terrain with variety of land forms such as ridges, valleys, mountainous, dissected topography with steep slopes as well as gentle and flat surfaces. The elevation ranges from 1453 m (Blue Nile river gorge) to 3528 m (Adama mountain) above sea level which makes all major rivers such as Andasa, Tul and their tributaries emanated from Adama mountain and drain to Blue Nile river.
Following a significant variation in elevation, the climate of the area varies from warm temperate to cool temperate climatic condition with unimodal rainfall pattern (Alemayehu 2006). The average mean monthly rainfall and temperature of seven meteorological recording stations of (1992-2016) years from Ethiopian Meteorological Agency (EMA) are 1372 mm and 19.4 °C, respectively. The rainy season of Ethiopia including the area under investigation is June to August, whereas the dry season includes December to February. Availability of groundwater sources in the whole watershed is highly variable irrespective of slight variation in rainfall amount. High demand and uncontrolled abstraction of the groundwater cause high stress on surface and groundwater resources, especially in the dry season. The north and northeastern part of the area has been used for irrigation activities by utilizing Andasa and Blue Nile River. The city expansion is well pronounced here and springs and rivers are intensively affected by uncontrolled motor pumps for irrigation purpose.

Datasets and methodology
Acquisition of RS and ancillary data, thematic layer preparation, GIS based MCDA (AHP) and weighted overlay analysis with weighted linear combination (WLC) techniques were used for delineation of GWP zones in the study area. Validity of the model was checked by adding productive boreholes and perennial springs into GWP map. The sensitivity analysis was performed by removing different thematic layers and evaluating the impact of assigned rates and weights on GWP map.

Data sources used and preparation of thematic layers
To demarcate GWP zone in the study area, nine thematic layers (maps) namely lithology, lineament density, geomorphology, slope, soil, drainage density, land use/cover, rainfall and depth to groundwater level were produced using satellite images and other conventional datasets (Table 1). Slope and drainage density thematic layers were produced from Shuttle Radar Topographic Mission (SRTM)-DEM (digital elevation model) at 30 m spatial resolution. The SRTM-DEM (acquired on 11 February 2000) was downloaded from the United States Geological Survey's Earth explorer site (https:// www. earth explo rer. usgs. gov/). The drainage density layer was calculated from drainage network (digitized from toposheet and extracted from DEM) using density calculation in spatial analysis tool of ArcGIS 10.5. The lithological thematic layer of the area was prepared from geological maps of Asrat (2017) and Bawoke et al. (2019) with modification by incorporating field observed data. Twenty-five years (1992-2016) recorded rainfall data of seven meteorological stations (Adet, Bahir Dar, Kimbaba, Merawi, Meshenti, Sekela and Tis Abay) were collected from Ethiopian Meteorological Agency (EMA) of Bahir Dar branch. The rainfall thematic layer was prepared by interpolating point rainfall data of each station using inverse distance weight (IDW) interpolation in ArcGIS. Depth to groundwater level of existing water points (boreholes, hand dug wells) of the study area were collected from secondary sources. All spatially referenced groundwater level data were imported and interpolated using Kriging technique for the preparation of depth to groundwater level thematic layer. The soil and geomorphologic data were extracted from secondary sources of FAO (1998) and ENTRO (2006), respectively. Sentinel 2A images were used to prepare land use/land cover (LULC) thematic layer. The images have 13 spectral bands at 10, 20 and 60 m spatial resolutions. The 10 m resolution bands (band 2, 3, 4 and 8) of the image (acquired on 26 February, 2019 from (https:// www. earth explo rer. usgs. gov/) were used to generate LULC thematic layer. A pixelbased unsupervised image classification was employed to prepare LULC map. SRTM image of the study region was processed by ENVI4.5 to change principal component image (PCI) into Geomatica software readable format for lineaments. Lineaments were extracted from SRTM-DEM image (30 m spatial resolution) using Geomatica software (version 2018) and validated with field collected and measured structural data. The extracted, field measured and validated geological structures were plotted using GeoRose software and displayed in the form of rose diagram. The extracted lineaments were further changed to lineament density thematic layer with spatial analysis tool of line density function in ArcGIS environment.

GIS based analytical hierarchy process (AHP)
AHP developed by Saaty (1980) has become the most widely used MCDA (Fenta et al. 2014;Kousalya et al. 2012) these days. AHP is a well-organized decision making technique which enables groups or individual decision makers to examine their decisions in time of complex situations. It provides a framework for decision makers by simplifying and enhancing natural decision-making processes. The working principle of AHP is based on: (1) breaking down and ordering the problem into hierarchy, (2) providing numerical value on the relative importance of each variable at different levels, (3) making judgements on which variable is the highest priority and influences the output, (4) prioritizing and ranking of preference from the decision alternatives.

Building a hierarchy of decision making elements (decision modelling)
The first step in any AHP needs the formation of hierarchical structuring of different decision making elements at different level of organization (objective, criteria, sub-criteria and alternatives). Modelling of the decision problem in such a way helps the decision makers to increase their understanding of the problem, its context and make alternative approaches to the problem under investigation. The overall goal of the decision is placed at the uppermost levels of the hierarchy, criteria (theme) and the sub criteria (classes) are at the intermediate levels and the decision alternatives at the last level of the hierarchy (Fig. 2).

Establishment of pairwise comparison matrix (PCM)
Following the formation of hierarchical structure, different comparison matrices of all elements at various levels were made (Kousalya et al. 2012). PCM in AHP refers to the system of comparing entities in pairs to evaluate which is preferred and by how much. Comparisons are carried out to determine criteria weighting and also judge the value or score of different classes (sub-themes) within each criterion (thematic layer). The determination in values of nxn PCM (denoted by M) is based on Saaty (1980) who has provided a scale ranging between 1 (equally important) to 9 (extremely important) and also the corresponding reciprocals that is m ij ∈ [1/9, 9] (Table 2, Eq. 1). On this regard, GWP controlling parameters of the area were given values according to their importance for groundwater occurrence. The less desirable entity within the comparison pair is given the inverse (reciprocal) values for less significant compared to the factor which it has been paired (Oikonomidis et al. 2015). Assignment of relative weights to different thematic layers and rates to classes within the same layer were made based on review literature, field experiences and expert's decision (Patra et al. 2018). A PCM for each thematic layer and the rate/rank for class were prepared, and compared with respect to relative importance of scores of element i with element j in m ij matrices (Eq. 1).
Let Eq.
(1) represents nxn square matrix of the form M = m ij where m ij > 0 for every i, j = 1, 2, …, n. A pair wise comparison matrix shows reciprocal provided m ij = 1 m ji (2) for every i, j = 1, …, n (then m ii = 1 for every i = 1, ..., n) (Koczkodaj and Szybowski 2015). Based on Satty's AHP matrix the different criteria (e.g. lithology, lineament, rainfall, etc.) and classes within the criteria were compared with reference to expert's decision, experience and local hydrogeological conditions.

Synthesis of priorities and consistency measurement
Nine matrices for each thematic layer consisting different classes and one matrix for the comparison of all the thematic layers were made and the relative weights/priority vectors were calculated. The main target to use mathematical procedures which perform synthesis involves the computation of eigenvalues and eigenvectors (priority vectors/weights). The relative importance of the various thematic layers and classes within it is dependent on the eigenvector and eigenvalue. Once PCM is generated, eigenvectors, eigenvalues and consistency evaluations (calculation of consistency index (CI) and consistency ratio (CR)) have been computed step by step.
Step 1: The square matrix of M = m ij (Eq. 1) with the element of row i column j was produced and the lower triangular matrix was completed by taking the reciprocal values of the upper diagonal using m ij = 1/m ij .
Step 2: Summation of all column j and normalization of the PCM.
Let Matrix M is a transformed matrix of A= a ij as a result of normalization, it provides normalized PCM of Step 3: Computing the Eigen vector (priority vector/ weights), w = ⌈w i ⌉ by calculating the arithmetic averages of normalized PCM across the row.
Step 4: Calculating the maximum eigenvalue ( max ) from the relation Ax = x , where A is a comparison matrix of nxn, x is eigenvector and max is the eigenvalue. Computation of eigenvalue is carried out by multiplying each column of PCM by the corresponding eigenvector and sum the value of the product across each row to obtain the total sums of the rows. This sum when divided by corresponding eigenvectors (Eq. 3) and averaged the division result gives eigenvalue ( max ) Step 5: Consistency evaluation/measurement. Eigenvalue is the main component for calculation of CI which in turn helps to calculate CR. Checking consistency of judgment is an important undertaking in AHP as human decisions lack perfect consistency in any set of PCM. Consistency check and verification helps to revise the judgment to be consistent with the allowable degree (CR ≤ 10%) but never be negative number as maximum eigenvalue of a matrix will not be less than the value of n, ( max ≥ n ). CR > 0.1 is inconsistent and needs further revision of subjective judgment for consistency. The following steps are followed for evaluation of consistency.
i. Computation of eigenvalue ( max ), step (4) ii. Consistency index calculation: where n is the number of items used for comparison. Following calculation of CI, CR is formulated as where RI refers random consistency index for the same order of matrices. RI is the consistency index of a randomly generated PCM which depends on the number of elements being compared (i.e., size of PCM) and takes on the following values (Table 3).

Combining thematic layers using weighted linear combination (WLC) techniques
The computed weights to each thematic layer were multiplied by the rates of classes for preparation of GWP map. All thematic layers, which were prepared with the same coordinate system and areal extent, were added and overlaid with WLC techniques (Malczewski 1999). The output potential map was reclassified into different zones.
The aggregating equation that has been used for computation of GWP using WLC method is displayed in Eq. 7. Different authors (e.g. Fashae et al. 2014;Fenta et al. 2014;Ghosh et al. 2020;Kaur et al. 2020) have used the equation for calculating groundwater potential index where GWP is groundwater potential, W i is weight for each thematic layer and R i , rates for the classes within a thematic layer derived from AHP. Figure 3 shows schematic diagram displaying overall methodological approaches for generation of GWP zones.

Sensitivity analysis
Sensitivity analysis measures variations and effectiveness of inputs parameters, and provides indispensable evidence on the effects of weights/ rates assigned to each thematic/subthematic layer on the output map. In this study, sensitivity analysis was performed to evaluate the effects of each thematic layer on the GWP map of an analytical model. Map removal (Lodwick et al. 1990) and single parameter (Napolitano and Fabbri 1996) analyses were employed for GWP mapping in this study. The map removal sensitivity analysis was carried out by removing one or more thematic layers (maps) and producing a new GWP each time for the remaining layers by overlaying each other. This technique has focused in verifying the relevance or impact of the different parameters (thematic layers) in assessing the GWP condition and helped to point out the most/least prominent thematic layer(s) in the computation of GWP map. The map removal sensitivity analysis has been worked out in terms of index of sensitivity using the following equation: where S is sensitivity analysis index as a result of removal of one map; GWP is the original GWP index obtained by computing all the nine thematic layers; GWP′ is the GWP index found by excluding one thematic layer at a time; (N = 9) and (n = 8) are the respective number of thematic layers used to calculate GWP and GWP′.
The single parameter sensitivity analysis was implemented to scrutinize the impact of each of the thematic layers on the GWP value. In this technique, effective or real weight of each parameter was compared with theoretical weight (assigned weight) allotted to the same layer in GWP map. The effective weight of each thematic layer (parameter) is expressed by the following equation: Moderately Experience and judgment slightly-to-moderately favor one activity over another 5 Strongly Experience and judgment strongly or essentially favor one activity over another 7 Very strongly An activity is strongly favored over another and its dominance is showed in practice 9 Extremely The evidence of favoring one activity over another is of the highest degree possible of an affirmation 2, 4, 6, 8 Intermediate values between two judgments Used to represent compromises between the preferences in weights 1, 3, 5, 7 and 9 Reciprocal When activity i relative to j is assigned, one of the above numbers, the activity j compared to i is assigned its reciprocal where W is the effective weight of each thematic layer, X r and, X w are the rates and weights of the thematic layers, respectively, GWP is the groundwater potential index as stated in (Eqs. 7, 8).

Results and discussion
The thematic layers Nine thematic layers have been produced for investigating groundwater prospective areas in Andasa-Tul watershed. The thematic layers include lithology (LT), geomorphology (GM), and lineament density (LD), and slope (SL), soil (SO), drainage density (DD), land use/land cover (LULC), rainfall (RF) and depth to groundwater level (DW). The PCM was constructed based on order of importance (LT > LD > GM > SL > SO > DD > LULC > RF > DW) ( Table 2) of each thematic layer to GWP and the corresponding weights were calculated and displayed as shown in Table 4. The weights of each thematic layer influencing the spatial distribution and extent of GWP zone derived from PCM vary from 0.35 (highest) to 0.02 (lowest) corresponding to lithology and depth to groundwater, respectively. The weights for the remaining thematic layers are found in between these values (Table 4).The consistency ration in Table 4 is less than 0.1 (10%) indicating that the judgment in PCM of the thematic layers is acceptable (Saaty 1980).

Lithology
Lithology plays a fundamental role in governing the spatial distribution and occurrence of groundwater. The porosity, size of pore space and the ease at which the pore spaces are interconnected controls storage and permeability of geologic medium that in turn affect the availability of groundwater in the area of interest. The main lithologic units found in study area consists Quaternary basalt (13.85%), Quaternary deposit (2.16%), Middle basalt (18.36%), Lower basalt (46.37%) and Upper basalt (19.25%). These lithologic units have been given weights (rates) on the basis of hydraulic properties (hydraulic conductivity, transmissivity, storativity and yields observed from pumping test, lithologic log (well completion reports) of the area. The (9) W = X r X w GWP , assigned order of importance for groundwater occurrence is QB > QD > MB > LB > UB, respectively. The spatial distribution of lithologic unit (Fig. 4a) and the assigned weights are depicted in Table 5. Despite its small areal coverage (13.85%) in the study area, Quaternary basalt is found to be the most important groundwater reservoir (aquifer). The existence of high potential aquifer in Quaternary basalt has also been witnessed by Kebede (2013) as he stated discharge of Araki spring in excess of 120 l/s in adjacent Infranz catchment of Tana sub-basin and Jiga spring discharging (400 l/s) of other area of the same lithologic unit in the Upper Blue Nile basin.

Lineament density
Identifying and characterizing the different types of geological structures is an important part of groundwater potential studies. The occurrence of groundwater is affected by orientation of different geologic structures in the area under investigation. Digital elevation model (DEM), satellite image interpretation, field observation and measurement of different structural features help in identifying the geological structures (faults, fractures and dykes) that exist in the study area. These structures, which have been traced as linear features (lineaments), play a significant role for groundwater occurrence and movements as they serve as a conduit or barrier for transmission and storage of groundwater. These lineaments show different orientations as depicted in Fig. 5. However, the most dominant orientations as shown in the rose diagram are NW-SE, NNE-SSW and NE-SW and this is possibly related to intra-plateau uplift, rifting and subsidence associated with Lake Tana formation (Chorowicz et al. 1998) in the adjacent basin. The increase in lineament density favors highest groundwater potential zone (Hammouri et al. 2012;Fashae et al. 2014;Bhuvaneswaran et al. 2015). The water movement between surface and subsurface is controlled by dykes and faults (Ahmed and Sajjad 2018). The lineament density of the area has been classified into five classes bounded between 0.0 km/km 2 (lowest) and 2.1 km/km 2 (highest) where the highest lineament density is intended to provide best groundwater potential zone (Table 6, Fig. 4b). Figure 4b reveals that more lineament density is found in the center close to north boundary of the area.

Geomorphology
The present geomorphologic set up of the study area is the combined result of complex interaction between geologic, tectonic and paleo-climatic variables (Xue et al. 2019). Geomorphology is affected by prevailing geologic structure, underlying geologic formation and surficial processes (Rajaveni et al. 2017) and plays an essential role on the occurrence, distribution and movement of groundwater (Fenta et al. 2014;Patra et al. 2018). Geomorphological mapping entails the assessment and characterization of the different geomorphic units. The main geomorphologic units existed in the study area based on ENTRO (2006) include plateau (14.20%), bottomlands (3.15%), minor scraps (18.08%), minor valleys (12.33%), foot slopes (29.57%), and complex landforms (10.70%), mountains (11.96%), respectively. Their order of importance and spatial distribution has been shown in Table 7 and Fig. 4c Slope Slope is one of the key factors controlling the rate of groundwater recharge and affects GWP of an area. Flat and gentle slopes provide long resident time and favor more infiltrations compared to steep slopes which facilitates run off or less infiltration (Fashae et al. 2014;Ibrahim-Bathis 2016;Jasrotia et al. 2016;Rajaveni et al. 2017). The inverse relationship between variation in slope and corresponding GWP is emphasized by different authors (Oikonomidis et al. 2015;Patra et al. 2018;Kaur et al. 2020). The slope thematic layer of the study area was produced from 30 m spatial resolution of SRTM-DEM and classified into six classes. The classification is based on Berhanu et al. (2013) classification scheme which ranges from flat (0-3) to extremely steep (≥ 50) percent (Table 8). The result shows 12.14% flat, 28.52% gentle, 23.29% moderate, 18.67% steep, 11.05% very steep, 6.33% extremely steep (Table 8 and Fig. 4d).

Soil
Soil plays important role for the spatial variation in quantity and quality of groundwater of an area. The variation is attributed to the type and nature of the soil covering the upper layer of the aquifer. The recharge (infiltration) rate of soil is highly affected by its porosity, size of pore space, adhesion, texture and permeability (Jasrotia et al. 2016;Patra et al. 2018). The soil map of the study area was prepared from the soil database of major soil classification of Ethiopia based on FAO (1998).The assignment of values to each soil group and calculations of the weight/rate was based on information from WRB (2014) and Driessen et al. (2001).The soil consists of eight major soil groups which are Eutric Leptosols (30.01%), chromic luvisols (7.81%), lithic leptosols (0.22%), Haplic Nitisols (5.08%), haplic luvisols (0.43%), Haplic Alisols (27.42) and Eutric Vertisols (28.73%). It has been shown from Table 9 and Fig. 6a that

Drainage density
The study area is dominated with dendritic and parallel drainage patterns reflecting existence of highly weathered geologic formation (regolith and soil), alignment of parallel ridges and geologic structures. The drainage density, which is a measure of length of drainage line to area of the basin, has an inverse relation with permeability and is indirect indicator of GWP (Magesh et al. 2012;Agarwal and Garg 2016;Rajaveni et al. 2017;Patra et al. 2018). The area shows more drainage density in the eastern part and southern end (near the boundary), but less in the northern and western portions (Fig. 6b). The drainage density varies between 0 km/km 2 and 1.5 km/km 2 and it is classified into five classes (0-0.3, 0.3-0.6, 0.6-0.9, 0.9-1.2, 1.2-1.5) km/km 2 with areal percentage of 16.83, 33.37, 35.39, 12.52, 1.89, respectively. Calculated rates are presented in Table 10 and the spatial distribution of drainage density is displayed in Fig. 6b.

Land use/land cover (LULC)
The land use/land cover plays an essential role in availability of groundwater (Berhanu and Hatiye 2020;Jasrotia et al. 2016;Kaur et al. 2020) as it affects the different hydrologic components (infiltration, surface runoff, interception and evapotranspiration). Forest and vegetated lands reduce runoff and enhance infiltration, whereas roads, parking lots, built-up and barren lands hinder infiltration and facilitate runoff (Kaur and Rishi 2018). The types and spatial extent of LULC identified in the area (Fig. 6c) include forest (16.8%), shrubs (31.9%), grassland (5.0%), and cultivated land (42.2%), built up (4.1%) and their order of importance (rate) is presented in Table 11. Results show that the highest and lowest rates are given for forest land and built-up, respectively.

Rainfall
Rainfall is essential hydrologic variable that controls the availability of potential groundwater (Nanda et al. 2017;Murmu et al. 2019;Kaur et al. 2020). The sustainability of groundwater is dependent on the amount of recharge   (Kotchoni et al. 2019). The mean annual rainfall which varies from 1271 to 1441 mm per year is classified into five classes (Table 12, Fig. 6d). Since rainfall has direct relationship with GWP, the highest weight (rate) is given to highest rainfall value and vice versa (Table 12). Figure 6d has revealed that rainfall amount generally decreases from west to east direction suggesting variation in elevation and other climatic variables.

Depth to groundwater level
The groundwater level is a dynamic surface in which its spatio-temporal variation is highly related with seasonal change (unconfined aquifers) and change in pressure (confined aquifer). The depth to water table for (unconfined) is very close to ground surface (shallow) during heavy rainfall season, but relatively deep in the dry season. About 200 groundwater level data were interpolated by kriging in ArcGIS and the output thematic layer was generated. The depth to groundwater level varies between 2.15 and 24.5 m below the ground surface. The thematic layer has been classified into five classes with lowest depth to groundwater level refers to high GWP and vice versa. Similar studies that used depth to groundwater level as parameter for GWP delineation were conducted by Patra et al. (2018) and Kaur et al. (2020). The spatial distribution of depth to groundwater raster thematic layer and its weight/rate are displayed in Fig. 7 and Table 13, respectively.

Groundwater potential map
Integrating different groundwater controlling variables for discovering available potential groundwater resource is imperative to optimize the choice of location for future groundwater development and management. All generated thematic layers were combined based on their weights and rates using WLC of GIS spatial analysis tool. The GWP map was prepared by computations of all the layers using the following equation: Fig. 6 Thematic layers: a soil, b drainage density, c land use/land cover, d rainfall GWP = 0.35 × LT + 0.22 × LD + 0.14GM + 0.10 × SL + 0.07 × SO + 0.05 × DD + 0.03 × LULC + 0.02 × RF + 0.02 × DW, (10) where GWP = groundwater potential, LT = lithology, LD = lineament density, GM = geomorphology, SL = slope, SO = soil, DD = drainage density, LULC = land use/land cover, RF = rainfall, DW = depth to groundwater level.
The GWP model output values which vary from 0.04 to 0.46 have been classified into five distinct zones using natural breaks classification mechanism in ArcGIS environment. The classified GWP values were then labeled into 'very poor' (20.5%), 'poor' (35.4%), 'moderate' (23.6%), 'good' (7.1%) and 'very good' (13.4%). Figure 8 reveals that 'very good' GWP zone is dominantly localized in the northern and western boarder (Sebatamit, Yigomahuletu, Yinesa, Yisala, Tatek and Berta Geberie areas) signifying the existence of Quaternary basalt with the lowest slope class (0-3%). The 'very poor' zone on the other hand is dominantly found in the southern part of the area (Gobe, Ayibar, Ambesit, Debal, Abyotefrie etc.) characterized by existence of upper basalt in mountainous area with little lineament density and moderate to extremely steep slopes. The 'good' GWP zone is found dispersed in the central and western parts (Zemene Hiwot, Foreswega, southern part of Yigoma Huletu, Yemeket) over different lithological units with relatively high to very high lineament density and flat to steep slopes. The 'moderate' and 'poor' GWP zones, consisting nearly 60% of the area, are found intermixed and occupying northeast to southwest aligned zone.

Validation of groundwater potential map
For the purpose of comparing yields of water point data (borehole and spring) versus predefined GWP zones, different categories or classification mechanisms together with consideration of site specific situations have been adopted for this work. Irrespective of lack of general standard classification system, groundwater yield data can be categorized into different zones with respect to site specific conditions. A study by Tuinhof et al. (2011) in Sub-Saharan Africa highlighted aquifer productivity classification into six groups as very high yield (> 20 l/s), high yields (5-20 l/s), moderate yields (2-5 l/s), low moderate yields (0.5-2 l/s), low yields (0.1-0.5 l/s) and very low yields (< 0.1 l/s). Similarly, Enideg (2012) classified yield of volcanic aquifer of central Ethiopia as high productive (5-20 l/s), moderately productive (2-5 l/s) and low productive (0.05-0.5 l/s). Groundwater study by SOGREA (2013) in Tana and Beles Sub-basins, on the other hand, classified aquifer yields into four classes as low (0-3 l/s), moderate (3-6 l/s), high (6-20 l/s) and very high (> 20 l/s) GWP zones. Comparing all aquifer productivity (groundwater yield) classification schemes, SOGREA's classification is selected for this study considering proximity of the area and geological settings to the present study. As a result, the low category of SOGREA is assigned to represent 'poor' and 'very poor' zones of this paper for validation purpose. Similarly, the very high, high and moderate classes refer to very good, good and moderate zones, respectively. The GWP zones were validated with (N = 37) groundwater points consisting 31 borehole and 6 spring yield data collected from field and water well drilling organizations and superimposed on the GWP map. It is found that yield in the 'very good zone' varies between 10 and 75 l/s with mean value of 34.43 l/s. The good and moderate potential zones show spatial variability in the range of 6-30 and 1.5-6 l/s with mean yield of 14.6 and 3.26l/s, respectively. Similarly, very poor to poor potential zones possess 0.04-4 l/s which give rise to mean yield of 0.88 l/s. Figure 8 shows that 31 out of 37 water points coincide with their respective potential zones making 83.8% accuracy (agreement) signifying the reliability of the method. The validation of groundwater potential zones using water point yield data has also been employed by other researchers (Fenta et al. 2014;Kumar et al. 2020;Murmu et al. 2019;Panahi et al. 2017).
The validation map portrays the concentration of very high yield water points (> 20 l/s) in the 'very good' zones and this is attributed to existence of Quaternary basaltic aquifer in the lowest slope values (Fig. 8). In contrary, very low and low groundwater potential zones (0-3 l/s) are characterized by existence of moderate to extremely steep slopes, high drainage density, and high depth to water level, mountains and complex landforms. Table 14 portrays the statistical result of map removal sensitivity analysis. The result summarizes the contribution of each thematic layer in affecting the model output (GWP map). Despite existence of differences in mean variation index (MVI), the removal of each layer has an impact on the output map suggesting all the thematic layers used in this AHP technique play their own specific role for GWP demarcation. Highest variation sensitivity index with (MVI = 1.92) is found by removing lithology layer and this is attributed to the assignment of relatively highest theoretical weight (35%). It is also found that GWP map is highly sensitive to slope (MVI = 1.50%) indicating that the slope is one of the important factors in affecting the distribution of GWP zones. Intermediate variation indexes have been obtained by removing rainfall (1.04%), depth to groundwater level (0.99%), soil (0.95%), LULC (0.85%) and lineament density (0.81%). Drainage density and geomorphology are comparatively least sensitive parameters with MVI of 0.62 and 0.59%, respectively.    (Table 15). Very closer or similar values between mean effective weights and empirical weights are observed in rainfall and depth to water level, respectively (Table 15). These two parameters have been found to be relatively less effective for GWP mapping. The single layer sensitivity analysis result shows that the lithology (effective weight = 26.38%) and slope (effective weight = 21.36%) are the most effective parameters in GWP mapping. This is in agreement with the map removal sensitivity analysis as the two thematic layers are most influential there. However, similar study by Fenta et al. (2014) in the same country but different study region showed slope as less effective thematic layer suggesting dependency of effective weights on local conditions and study regions. Such dependency and sensitivity variations of GWP with respect to different cases and study regions has also been highlighted by Panahi et al. (2017)

Conclusion
Identification and demarcation of GWP zone is imperative for future sustainable use of the resources. The spatial variation of groundwater in Andasa-Tul watershed which has a total area of 1872 km 2 was studied using nine GWP controlling thematic layers prepared from conventional data sources and satellite images. These thematic layers include lithology, lineament density, geomorphology, slope, soil, drainage density, land use/land cover, rainfall and depth to groundwater level. GIS-based AHP was employed to  calculate weights and rates of the thematic layers. Thematic layers were combined using WLC of GIS spatial analysis tool. The aggregated results of all thematic layers produce different GWP zones designated as 'very good', 'good', 'moderate', 'poor' and 'very poor' potential zones. Very good GWP zone has been found in the northern and western parts close to the boundary of the area especially in the Quaternary basalt. The poor and very poor groundwater potential zones are localized in the southern, central and northeastern mountains and complex land forms. Validations of GWP zones with groundwater yield (l/s) data (Table 16) provide efficient result with 83.8% agreement.
Single layer sensitivity analysis result discloses the slight disparities between the empirical and effective weights. The effective weights of lithology (26.38%), lineament density (14.34%) and geomorphology (10.37%) show smaller values compared with empirical weights of the same thematic layers of (35%, 22%, and 14%), respectively. However, in all remaining thematic layers effective weights are generally greater than the empirical weights. Very close similarities between effective and empirical weights of LULC, rainfall and depth to water level have been observed. It is generally found that lithology and slope are more effective, whereas rainfall and depth to groundwater are less effective thematic layers for GWP mapping.
The result of this study could become an important input for further detail groundwater explorations (geophysical and hydrogeological) as it predefines potential areas. It also helps to manage and protect groundwater sources for future sustainable use of the resources. In general, the study of GWP zones by the use of integrative RS, GIS based MCDA (AHP) is an essential time and energy saving technology which has to be used ahead of detail geophysical exploration and water well drillings. Since the method is generic in nature and can be implemented in a wide variety of environments, it can be used for the same or similar purposes in other areas/ regions.