Spatial prioritization of groundwater resources in mountain region of Nepal: a GIS and AHP integrated framework

About one-third of the global supplies of freshwater are found as groundwater which is also a significant component of the hydrologic cycle (Siebert et al. 2010). It has a clear benefit over surface water resources since it is comparatively free of turbidity, pathogenic organisms, biological, chemical, and radioactive contaminants. It is almost abundant everywhere on the land stages and its chemical composition is nearly stable, barring unexpected contamination from factories or mines. Therefore, most of the practitioners and engineers involved in water supply and irrigation are interested in this freshwater reservoir (Barman et al. 2024). Most of the developing countries rely on groundwater because it needs minimal treatment or even no treatment. In Nepal, groundwater is one of the main supplies of potable water. Though it is available in almost every part of the country, it varies in quantity and depth from place to place. The geography of Nepal ranges from a low-lying, 59-m-high terai area to extremely high, 8,848-m-high slopes that are home to a diverse range of lithological units, including high-grade metamorphic rocks and unconsolidated alluvial deposits (Dhital 2015). More than 90% of people living in the highlands rely only on springs for their groundwater needs (Pandit et al. 2024). These resources are utilized for gravity flow water supply systems in the hills and mountains, which serve as the catchment's groundwater storage and are crucial to the Himalayan water budget (Andermann et al. 2012). These springs get recharged mainly with monsoon rainfall (Sharma et al. 2016) and are susceptible to vulnerability (Gurung et al. 2019a, 2019b). Also, these sources are depleting due to human activity and climate change (Tambe et al. 2012). For example, the hill slope hydrology of Nepal's middle mountains has been disrupted by anthropogenic activities such as land use change, watershed degradation, and road network growth (Ghimire et al. 2019). These man-made disturbances have caused spring sources to dry up and disrupt natural flow regimes, especially during the dry season (ICIMOD 2015; Chapagain et al. 2019; Ghimire et al. 2019). Additionally, the 2015 earthquake in Nepal affected the water table in the central Himalayas as numerous spring sources dried up due to the catastrophe (Bricker et al. 2014). As a result, existing springs are under pressure, making it crucial to have a thorough grasp of groundwater storage volumes at various settlement area sites.

In addition to the quantity of groundwater recharged, its quality is equally important. Agricultural and cattle farming activities in the basin may pose significant pollution risks, introducing contaminants such as pesticides, fertilizers, and animal waste (Raju et al. 2015). These pollutants can degrade groundwater quality, impacting environmental health and increasing the costs and complexities of water treatment. Mitigating these issues requires best management practices (BMPs) in agriculture, such as buffer strips, cover crops, and controlled agrochemical use (Jain & Singh 2019). An integrated approach to water resource management, involving collaboration among farmers, policymakers, and water managers, is essential for sustainable and safe groundwater recharge. Effective long-term groundwater management necessitates the proper execution of strategies and policies, which are formulated following a comprehensive evaluation of this critical and fragile resource (Sharma & Shakya 2006).

The essential initial step in groundwater management is to evaluate its available potential. Various techniques exist for assessing groundwater, with traditional investigation methods being particularly time-consuming and costly. These conventional methods require extensive use of geological, hydrogeological, and geophysical survey tools (Sander et al. 1996; Sturchio et al. 2004; Kumar et al. 2005; Israil et al. 2006; Mall et al. 2006; Jha et al. 2010; Singh et al. 2013; Senanayake et al. 2016; Roxy et al. 2017; Asoka et al. 2018; Mukherjee et al. 2018). On the other hand, several methods, including artificial neural networks, decision tree model, Dempster–Shafer model, frequency ratio, fuzzy logic, logistic regression, maximum entropy model, and weights of evidence model have been developed recently to improve the methodological techniques used in groundwater investigation (Kodihal & Akhtar 2024). The AHP is a widely used scientific decision-making technique among the different multi-criteria decision-making (MCDM) techniques because of its structural simplicity, efficacy, efficiency, and consideration of numerous watershed-influencing factors (Franek & Kresta 2014; Rejith et al. 2019; Arunbose et al. 2021; Baghel et al. 2023; Tavana et al. 2023).

The AHP technique, in the field of groundwater, was first proposed by Saaty (1980). The integration of AHP with RS and GIS has emerged as an effective approach for accurate and efficient estimation of natural resources (Baghel et al. 2023). This method reduces complex decisions to a series of pairwise comparison matrices which is further reduced to a normalized weight matrix and then it is checked for consistency, therefore reducing the bias in the process of decision-making. The weights assigned for different thematic layers in the AHP are integrated with a GIS and are considered a powerful tool for the efficient assessment of available natural resources. It provides reliable estimation of spatio-temporal information of GPZs within the watershed and is recommended for wider applicability of the technique with some customization (Chandio et al. 2013). The validation methods are slightly different from other methods. Generally, the bore logs and well data have been used by most of the investigators for validation (Berhanu & Hatiye 2020; Serele et al. 2020). However, for mountain region springs, data can be considered for better estimation and validation of the method (Sapkota et al. 2021).

Previous studies suggest that springs reliant on both rainfall and the properties of their recharge zones are drying up at an alarming rate (Adhikari et al. 2021). This is due to growing unpredictability in rainfall patterns (both seasonal and annual) and the impact of human activities on recharge areas. Population growth increases water demand, stressing existing sources like springs, which remain unchanged or have dried up. Unfortunately, inhabitants of the municipality mainly rely on these spring supplies for their domestic, agricultural, and horticultural needs. Therefore, it is imperative to estimate the groundwater potential zones (GPZs) of the study area, as no prior studies have been conducted in this region. This study aims to identify GPZs using AHP techniques integrated with RS and GIS tools, which can serve as a blueprint to mitigate the water shortage problems and effectively manage the existing groundwater. Spring locations will be used to validate the developed model. While this study primarily focuses on groundwater quantity, the quality is equally crucial and has been discussed in further research needs.

The study used primary and secondary data obtained from several sources, including field surveys. The Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) with a resolution of 30 m was downloaded from United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/) and clipped to the study area's boundaries using GIS. This DEM was further processed using GIS to generate thematic layers, including maps of slope, drainage, drainage density, and lineament density. Additionally, LULC maps, soil maps, curvature, TWI, TPI, roughness, and geological maps were sourced from secondary data providers such as the International Centre for Integrated Mountain Development (ICIMOD) and the Soil and Terrain Database (SOTER) and subsequently processed using GIS. Similarly, rainfall data for gauging stations falling in the area were gathered from the Department of Hydrology and Meteorology (DHM). A detailed summary of these secondary data sources and the information they provided is mentioned in Table 1.

The study utilized GIS-assisted multi-criteria decision analysis with analytical hierarchy process (MCDA-AHP) techniques. This method involves the conversion of geographical data (input) into decision outcomes (output). Qualitative information about specific themes and attributes is quantified by creating a pairwise comparison matrix employing Saaty's scale (Saaty 1980). The weight of each influential indicator was allocated by the AHP approach. A weight matrix is used in the AHP method to assign relative weights to the criteria and sub-criteria that are used to evaluate alternatives. This weight matrix takes the form of a square matrix, where each row and column correspond to a criterion or sub-criterion, and the matrix elements reflect the pairwise comparisons made between these criteria. Typically, the construction of the weight matrix relies on expert or stakeholder opinions and judgments. For this study, 11 different parameters were taken into account, drawing on existing research (Arulbalaji et al. 2019; Sapkota et al. 2021; Baghel et al. 2023), data from the Department of Mines and Geology (Nepal), stakeholder viewpoints, and the perspectives of experts and practitioners. The considered parameters are as: landuse/cover (LULC), geology, rainfall, lineament density, drainage density, soil, slope, topographic wetness index (TWI), topographic position index (TPI), roughness, and curvature. The AHP method has been well documented and more details can be found in the available literature (Saaty 1980; Arulbalaji et al. 2019; Baghel et al. 2023). The summary of steps involved in creating a weight matrix in AHP is as follows:

• 1. Identify criteria and sub-criteria: first, identify the relevant criteria and sub-criteria related to the decision. Ensure these aspects cover all vital aspects of the problem, here GPZs mapping.

• 2. Perform pairwise comparisons: compare criteria and sub-criteria pairwise to establish their relative importance. This comparison is done using a scale ranging from 1 to 9, where 1 signifies the same importance, and 9 means extreme importance of one over the other (Table 2) (Saaty 2001).

• 3. Construct the weight matrix: utilizing the pairwise comparison results, a square weight matrix is created as presented in Equation (1). The matrix has the same number of rows and columns as the criteria or sub-criteria. Diagonal elements are set to 1, representing equal importance. Off-diagonal elements are set to the reciprocals of the corresponding pairwise comparison values. For instance, if criteria i and j have a pairwise comparison value of 5, then elements (i,j) and (j, i) in the matrix are set to 1/5 and 5, respectively.

  

  (1)

• 4. Where M represents a n × n comparison matrix, where n is the total number of parameters selected for comparison, and m_ij denotes the alternatives with respect to a particular criterion; m_ij > 0;1 ≤ i ≤ n;1 ≤ j ≤ n.

• 5. Normalize the weight matrix: normalize the weight matrix by dividing each column by the sum of its elements. This ensures that the weights assigned to each criterion or sub-criterion add up to 1, preserving their relative importance.

• 6. Check consistency: verify the consistency of the weight matrix to confirm the coherence of the pairwise comparisons. If the consistency ratio exceeds 0.1, consider revising the pairwise comparisons to maintain consistency.

• 7. Utilize the weight matrix: once the weight matrix is constructed and normalized, it becomes a valuable tool for calculating overall scores for each alternative based on the criteria and sub-criteria weights. These scores facilitate decision-making in line with the established criteria and sub-criteria.

• 8. Using the Eigenvector technique, the principal Eigenvalue (λ_max) was determined, and the consistency index (CI) was derived from the following equation:

  

  (2)

  where n represents the count of significant factors considered in the research.

Consistency ratio is denoted as CR = CI/RI, where RI is a random index, and its values were given by Saaty's standard as presented in Table 3. Saaty has stated that a CR of 0.10 or lower is considered acceptable to proceed with the analysis. If the consistency value surpasses 0.10, it is necessary to review the judgments, identify inconsistencies, and make corrections. Perfect consistency in pairwise comparisons is denoted by a CR value of 0.

The spring distribution map was used to validate the GPZs map as the availability of spring confirmed the groundwater availability (Sapkota et al. 2021). The spring distribution map was merged with the GPZs map. If areas with higher groundwater potential were associated with a higher number of springs, the produced map was considered to be reliable at this stage. This is because the existence of springs is a clear indication of groundwater availability in those regions. Furthermore, it is essential to verify the output model's accuracy while using the MCDM technique. The receiver operating characteristic (ROC), which calculates the area under the curve (AUC) using inventory information (groundwater level data at numerous sample sites), is a recognized metric for verifying such MCDM models (Pourghasemi et al. 2012; Das 2020). In this study, the number of identified springs was used as inventory data. The maximum level of unbiased accuracy is shown by an AUC value of 1. Higher than 0.7 AUC values are often seen as indicative of appropriate differentiation and, hence, appropriate models (Das et al. 2021). Table 4 shows the AUC value ranges with their test quality meanings.

Eleven thematic maps have been prepared in Arc GIS Pro. Layers were further categorized into different classes and each was given weight as per its relative importance in groundwater potential recharge as discussed earlier.

The distribution and occurrence of groundwater in any terrain are greatly influenced by the geologic setting (Yeh et al. 2016). The different formations of rock layers (lithology) have different water transmissivity and thus hydrological importance of rock is considered while assigning weights to different lithologies. For geology map preparation, geological data at approximately 1:250000 scales prepared by the Department of Mines and Geology of Nepal were digitized and reclassified into 30-m resolution raster data. Details of the geology class and map are presented in Figure 3 and Table 5. The geology map covers four different rock classes such as Slate, Dolomite, and Limestone in one class, encompassing 51.81% of the study area which bears higher water transmissivity and has been assigned higher weights (Arulbalaji et al. 2019; Baghel et al. 2023; Silwal et al. 2023). Similarly, the other two classes of rock, Quartzite and Phyllite cover 14.5 and 3.2% of the study area having moderate to low water transmissivity which is given the second highest weights and the last one is Genesis formation covers 30.52% of the study area, having the lowest water transmissivity rate, which is given very low weight as given in Table 5. The CI and CR values of the geology are found to be 0.03 and 0.033 which is less than 0.1, this confirms that the assigned rank for geology is suitable. The details of the geological parameters are given in Table 5.

The soil texture plays a crucial role in groundwater recharging, and it is one of the important geomorphic components responsible for water transmission (Ibrahim-Bathis & Ahmed 2016; Das 2017). The soil map was obtained from SOTER, cropped using a GIS to the study region, and confirmed through on-site inspection. The proportion of sand, silt, and clay was extracted from the SWAT.mdb file using the recognized soil name. Soil texture was determined using the soil texture calculator triangle and the percentages of sand, silt, and clay. The research area contains four different types of soil such as sandy soil having high water transmission rates (hydraulic conductivity) (>10 cm/h), loam soil having moderate transmission rate (0.2–0.6 cm/h), clay loam soil having low to moderate water transmission rate (0.05–0.2 cm/h), which is dominant soil that covers most of the study area as shown in Figure 3, and clay soil having very low (<0.01 cm/h) water transmission rate (Garca-Gaines & Frankenstein 2015). The above-discussed soil type covers 22.57, 0.12, 55.65, and 21.66% of the total area, respectively. The soil texture, sandy soil has high hydraulic conductivity given a higher weightage value. Similarly, lower weightage is given to clay soil having very low hydraulic conductivity. The assigned weight and rank for soils are found to be satisfactory since CI is 0.05 and CR is 0.053 which is less than 0.10. All the classified soil textures with their coverage area and assigned weights are given in Table 5.

Rainfall is considered one of the major sources of groundwater recharge (Shrestha & Sthapit 2015). Approximately 80% of Nepal's yearly precipitation falls during the summer monsoon, which runs from June to September (Shrestha 2000). For the present study, 30 years (1980–2009) of average annual rainfall data is used. The average annual rainfall ranges from 871.88 to 985.96 mm which are classified into five categories such as very low (871.88–895 mm), low (895–917 mm), moderate (917–939 mm), high (939–961 mm), and very high (961–985.96 mm). The spatial distribution map for the average annual rainfall map is prepared by the inverse distance weighted (IDW) method of interpolation available in Arc GIS as depicted in Figure 3. Groundwater potential depends on the infiltration rate of rainfall and duration of rainfall. Rainfall has a significant impact on groundwater potential availability among various parameters (Kathe et al. 2024). It is high for higher rainfall intensity (Arunbose et al. 2021). Therefore, high weights are assigned for higher rainfall and vice versa as presented in Table 5. The CI and CR of the rainfall parameter are 0.05 and 0.049, respectively. The computed CR value is less than 0.1, thus the assigned rank for the AHP method is valid.

A lineament is a geographical feature that indicates the presence of an underlying geological feature, such as a fault, fracture, joint, or other feature, that leads to higher secondary porosity and permeability (Yeh et al. 2016). These structures are essential for runoff water to enter the subsurface and for the storage and transport of groundwater (Arulbalaji et al. 2019). Maps with hillshade information are useful for identifying possible lineaments and for visualizing topographic features. For this, the lineament information is obtained from four different hillshade maps of different azimuth and altitudes in Arc GIS, and then using the line density tool in Arc GIS, lineament per unit area is obtained which is classified into five categories such as very low (0–0.7 km/km²), low (0.7–1.4 km/km²), moderate (1.4–2.1 km/km²), high (2.1–2.8 km/km²), and very high (2.8–3.6 km/km²). The different class of lineament density covers area of 21.99, 36.41, 26.95, 12.90 and 1.75% of study area, respectively, as shown in Table 5. The higher weight is assigned to high lineament density and lower weights to low lineament density. The CI and CR values are 0.07 and 0.063, respectively. The details of the factors assigned for the AHP method are given in Table 5.

Curvature, which can be concave upward or convex upward profiles, is a quantitative representation of the nature of a surface profile (Nair et al. 2017). Water tends to gather in concave and slow down in convex profiles. The study area's curvature ranges from 35.5 to −34. The data are reclassified and divided into five classes such as −34 to −20, −20 to −6, −6 to 8, 8 to 22, and 22 to 35.5 covering areas of 0.14, 0.39, 99.32, 0.15 and 0.01% of the study area, respectively. The CI and CR values are 0.05 and 0.045, respectively, shown in Figure 4 and Table 5. High curvature values are given high weights, and vice versa.

The degree of variation in elevation between adjacent cells is indicated by the roughness index of a DEM (Riley et al. 1999). The topography's undulations are often expressed by the roughness index. More undulation corresponds with increased roughness and vice versa. A mountainous area with undulating topography is one where weathering and erosion processes gradually change the harsh terrain into a smooth, level surface over time. The roughness of the study area varies from 0.17 to 0.83 which is reclassified into five categories, such as 0.17 to 0.3, 0.3 to 0.43, 0.43 to 0.56, 0.56 to 0.69 and 0.69 to 0.83 which covers 0.46, 11.29, 73.79, 13.97 and 0.49% of the study area, respectively, as shown in Table 5 and Figure 5. In the study, a higher weight is assigned to a low roughness value and vice versa. The CI and CR values are 0.06 and 0.053, respectively.

The presence of lithologies such as slate, dolomite and sandstone spread in the northwest region to the south as shown in Figure 3 of the geology map, possesses high water transmissivity and has resulted in regions with high to very high GPZs in this area. In the very low and low GPZs regions, rock formations include layers of Gneiss having the lowest water transmissivity. Similarly, the moderate GPZs regions contain Phyllite and Quartzite, which exhibit moderate water transmissivity. Similar research done by various investigators shows that geology has a strong influence on groundwater potential and recharge assessment, for example Andherikhola watershed region of Nepal (Sapkota et al. 2021) and high relief zones of Greece (Oikonomidis et al. 2015). In addition to this Patra et al. (2018) found that the groundwater potential is strongly influenced by geology in the alluvial plains of India.

Similarly, rainfall intensity is also supported by the high water transmissivity nature of lithological units. These regions have high rainfall intensity. Rainfall is most abundant in the northwest part of the watershed and diminishes toward the southeast. Consequently, the areas with higher rainfall exhibit significantly larger GPZs. Arulbalaji et al. (2019) also found that the high rainfall region has high GPZs in the mountains of the southern western ghats of India. Similarly, Baghel et al. (2023) did similar research and found that the high rainfall region of the Mand River catchment of the Mahanadi basin possesses high GPZs.

Additionally, the regions with high GPZs are characterized by low drainage density, which increases in areas with lower GPZs. Low drainage density leads to reduced runoff in streams, allowing most of the rainfall to infiltrate directly into groundwater zones. Similar research done by Diriba et al. (2024) shows that the lowest drainage density results in high infiltration causing high GPZs. However, steep slope will have low infiltration resulting in low GPZs (Diriba et al. 2024). The watershed's terrain slope also plays a role, with very steep slopes covering a small percentage of the area, where fewer springs are found. In contrast, the flat or moderately sloped regions of the watershed have higher GPZs and a greater number of springs as depicted in Figure 6. The results from all the studies mentioned earlier align with our findings. Consequently, it can be concluded that geology, rainfall and drainage density are the most significant factors in groundwater potential zonation in the region. The methodology used has demonstrated its effectiveness in complex, rugged hilly and mountainous areas with hard-rock aquifer systems.

The GPZs map is validated by comparing the results with the number of springs found in the study region for which the GPZs map has been prepared. The GPZs map and the number of springs are overlaid in the GPZs map. The region having high to very high GPZs has a larger number of springs as compared with low GPZs zones as shown in Figure 6. A total number of 56 active spring locations were identified from the field survey. These all springs were plotted on GPZs map. Out of 56 identified springs, 25 fall in very high, 15 in high, 12 in moderate and 4 in low GPZs region with no spring in very low GPZs region.

This study employs AHP and GIS techniques to delineate GPZs in the Doramba rural municipality, a small mountain aquifer in Nepal. It considers 11 thematic layers that directly impact groundwater availability, including drainage density, slope, rainfall, lineament density, LULC, soil type, geology, curvature, roughness, TWI, and TPI. The resulting GPZs map categorizes the study area into five distinct zones of groundwater potential: very low, low, moderate, high, and very high. The north-west and west-south regions predominantly fall within the very high and high zones, covering 39% of the study area, while areas classified as very low represent less than 5%. The majority of the area, 47%, is categorized as having moderate groundwater potential. The accuracy of the GPZs map was verified using spring locations throughout the study area, with an AUC value of 0.707 confirming the model's effectiveness in mapping groundwater potential. This delineated map serves as a preliminary reference for decision-makers to manage and monitor groundwater resources in the future to cope with climate change and other anthropogenic impacts. To assist a user with the practical application of the new AHP_GIS-based GPZs model, a Supplementary User Manual for ‘Spatial Prioritization of Groundwater Resources in Mountain Region of Nepal: A GIS and AHP Integrated Framework’ is prepared and added to the Supplementary Materials.

The study acknowledges the potential for additional groundwater-influencing parameters, such as aquifer properties and geological heterogeneity, recommending their consideration for future research. It also suggests that incorporating the effects of anthropogenic activities and climate change could enhance the development of comprehensive groundwater management plans and policies. Moreover, future research could explore the qualitative assessment of GPZs, mapping zones according to the water quality for drinking, agriculture, and industrial use. Lastly, the application of interpretable machine learning methods for delineating GPZs, although beyond this study's scope, is proposed as a valuable direction for future investigations.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.