Optimizing Rainwater Harvesting Site and Structure Suitability Using Geospatial, Multi-Influencing Factor, and Analytic Hierarchy Process Approaches

Summary

This paper describes an integrated Geographic Information System (GIS)-based approach using Multi-Influencing Factor (MIF) and Analytic Hierarchy Process (AHP) models to identify suitable rainwater harvesting sites and select appropriate structures in mountainous river basins, demonstrated in the Panjkora basin, northern Pakistan.

Abstract

Rainwater harvesting (RWH) is an essential practice for water conservation, improved water resource management, and water-related hazards mitigation in mountainous areas. Selection of a suitable site and an appropriate structure of the RWH facility are vital for improving water availability and agricultural productivity under all circumstances, particularly due to climate change-related hydrological uncertainty. It is challenging to evaluate and analyze RWH sites across diverse conditions worldwide, particularly in remote, inaccessible mountainous areas where these sites have significant impacts on the environment, society, and economy of the region and downstream. In this study, the Multi-Influencing Factor (MIF), and the Analytic Hierarchy Process (AHP) were applied in Geographical Information System (GIS) using customary and remote sensing (RS) data for selecting a suitable RWH site and appropriate structures in the Panjkora river basin in the Hindu Kush region in northern Pakistan. According to MIF (and AHP) results, the study region has 80.22 (1572.58) km² less suitable area, 1681.99 (1605.69) km² as moderate suitable, 3116.10 (1768.62) km² as suitable, 844.86 (689.15) km² as high suitable and 35.10 (122.61) km²as very high suitable for RWH structures. The resulting maps were validated using Receiver Operating Characteristic and Area Under the Curve (ROC-AUC) tests (MIF score = 0.724 and AHP score = 0.692) to check the accuracy and robustness of the models. This research, presenting results with promising accuracy, will provide new technical insights on the topic for further improvement, suitability and applicability of the models under different hydro-meteorological and physiographic conditions. Overall, both models successfully identified suitable rainwater harvesting sites; however, the MIF model outperformed the AHP model in terms of predictive accuracy and spatial reliability. The proposed GIS-based framework can support sustainable rainwater harvesting planning and water resource management in mountainous watersheds.

Introduction

Water scarcity, including the depletion of both surface and groundwater resources, has become a major challenge in many developing countries¹. Rapid population growth, urbanization, agricultural expansion, complex topography, and changing climatic conditions have intensified pressure on available water resources^2,3. In water-stressed regions, excessive groundwater abstraction to meet domestic and agricultural demands has resulted in declining water tables and reduced long-term water security. Therefore, sustainable water management strategies are urgently required, particularly in regions experiencing rapid population growth, high vulnerability of water resources from climate change, and increasing water demand, such as Pakistan. Among various adaptation strategies, rainwater harvesting (RWH) has emerged as an effective approach for conserving rainfall runoff and supplementing available water resources⁴. RWH systems collect and store rainfall during wet periods for subsequent use, providing a sustainable alternative for addressing water shortages at both local and regional scales^5,6. However, the success of RWH systems largely depends on identifying suitable locations and appropriate harvesting structures. Site selection is a complex process because it is influenced by multiple physiographic, environmental, hydrological, and socio-economic factors that vary spatially across regions^7,8. Traditionally, field surveys have been used to identify potential RWH sites; however, these approaches are often expensive, time-consuming, and impractical for large and inaccessible areas, particularly mountainous regions. Therefore, Geographic Information System (GIS) and Remote Sensing (RS) technologies have increasingly been adopted as efficient alternatives for regional-scale RWH assessment⁹. GIS-based multi-criteria evaluation (MCE) provides an effective framework for integrating diverse thematic layers, including rainfall, runoff, slope, soil characteristics, land use/land cover, drainage density, and socio-economic parameters, to identify suitable RWH zones^10,11. Remote sensing datasets further enhance this approach by providing cost-effective spatial information over large and inaccessible areas, including land cover, geomorphology, soil characteristics, and hydrological features¹². Pakistan is among the countries facing severe water challenges due to increasing water demand, climate variability, and dependence on agriculture-based livelihoods¹³. The country's agricultural economy requires substantial water resources, while declining groundwater levels and irregular precipitation patterns have increased vulnerability to water shortages. Consequently, rainwater conservation and storage have become important strategies for improving water availability, especially in water-deficient and mountainous regions^12,14. Previous studies have demonstrated the effectiveness of GIS and RS-based approaches for RWH site selection by integrating environmental and socio-economic parameters^15,16,17,18. Several researchers have applied various multi-criteria decision-making techniques, including Analytic Hierarchy Process (AHP), Fuzzy AHP, and other weighting approaches, to identify suitable RWH locations under diverse environmental conditions¹⁹. Despite significant progress, the performance of different multi-criteria decision-making approaches may vary depending on regional environmental characteristics, available datasets, and expert-based weighting systems. In particular, limited studies have compared the effectiveness of different approaches under similar watershed conditions. Therefore, evaluating and comparing alternative decision-making techniques is essential for improving the reliability of RWH site suitability assessments. The present study addresses this research gap by comparing the Multi-Influencing Factor (MIF) and Analytic Hierarchy Process (AHP) techniques for identifying suitable RWH sites in the Panjkora Basin, Pakistan. It is hypothesized that incorporating different variables like topography, geological properties, land cover dynamics, drainage-lineament density, soil characteristics and rainfall in GIS using a multi-criteria decision approach can identify the best suitable RWH sites under a mountainous environment. The study provides valuable insights into the applicability of GIS-based decision-support approaches for sustainable water resource planning in mountainous regions.

Protocol

This study involved geospatial, remote-sensing, and field-validation data and did not involve human participants, identifiable personal data, animals, or vertebrate tissues; therefore, institutional human or animal ethics approval was not required.

Study region

The Panjkora river basin is an important physiographic region situated in the eastern Hindu Kush Mountains in the northern Pakistan (Figure 1). Panjkora river is the main river of the basin (113 km long with 5758.27 km² catchment area) and starts as a torrent from ice-capped mountains of Hindu Kush. It joins the River Swat near Chakdara, Dir Lower²⁰. Panjkora River is joined by five significant torrents or streams, including Barawal, Dir, Gawaldai, Jandol, and Kohistan. It extends from 34°39′30′′ to 35°46′1′′ North latitudes and from 71°13′08′′ to 72°22′13′′ East longitude. The region's location and rugged topography significantly influence its climate (mountainous and temperate). The Upper (Kumrat, Thal) region of the basin has a longer winter season and a colder summer season. From November onwards, the temperature drops sharply. However, in Dir Lower (Timergara, Talaash, Maidan, Samarbagh), the temperature is usually above the freezing point from December to February. The warmest months in Timergara are June through August, with average maximum temperatures above 35 °C, whereas June and July are the hottest months in Dir Town (with maximum temperatures of 32.4 °C and 31.5 °C). Monsoon is the source for summer rainfall, whereas the Western Depression brings winter rainfall. The study area is characterized by a high relative humidity throughout the year. Riverine and flash floods²¹ occur (almost) every year, especially in areas up and downstream of Wari. The major crops cultivated in the region include rice, wheat, corn, potato, and onion, while significant fruits grown in the study area are persimmon, orange, apple, walnuts, apricot, plum, loquat, and mulberry.

Figure 1: Study area map of the Panjkora River Basin, northern Pakistan. (A) Location of Khyber Pakhtunkhwa within Pakistan; (B) location of the Panjkora River Basin within Khyber Pakhtunkhwa; and (C) Panjkora River Basin showing the basin boundary, elevation distribution, river network, and major locations within the study area. Please click here to view a larger version of this figure.

Data collection and preparation

For this study, the data were collected from different sources. Rainfall/ Precipitation data were downloaded from Global Precipitation Measurement (GPM) National Aeronautics and Space Administration (NASA) https://gpm.nasa.gov/missions/GPM from 2014 to 2023. The soil texture data were collected from the Directorate of Soil Survey Khyber Pakhtunkhwa, Pakistan (www.soilconservation.kp.org). The Geological data was obtained from the regional office of Geological Survey of Pakistan (https://gsp.gov.pk/). For collection and calculation of land Scenarios (Land use/Land Cover) Sentinel 2 images were obtained from European Space Agency (ESA) Copernicus Open Access Hub (https://scihub.copernicus.eu/). Sentinel-2B imagery acquired on 10 September 2025 was used for land use/land cover (LULC) mapping. The image was processed and classified using the Maximum Likelihood Classification (MLC) algorithm. A total of 65 training samples were collected across the study area representing seven LULC classes: water bodies, forest, agricultural land, urban areas, bare soil, snow/ice, and rangeland. The prepared training samples were used to perform supervised classification and generate the final LULC map. The classification accuracy was evaluated using an accuracy assessment approach based on validation samples, including overall accuracy and Kappa coefficient. The digital elevation model (DEM) model 12.5 Spatial resolution was acquired from Alaska Satellite facility (ASF) (https://asf.alaska.edu/) on 2/12/2023. The DEM model was further used to generate slope, drainage network, drainage density, and elevation layers. The existing rain water harvesting structures data were collected from relevant departments for cross validation.

All spatial datasets were processed and analyzed using geographic information system (GIS) software (see Table of Materials). The GIS thematic-layer data are provided in Supplementary File 1. All input datasets were projected into a common projected coordinate reference system (CRS) (WGS 1984 UTM Zone 42N) to ensure spatial consistency and accurate area calculation. Raster datasets with different spatial resolutions were resampled and aligned to a common grid using the nearest neighbor resampling method, while maintaining the original spatial characteristics of categorical datasets. The Digital Elevation Model (DEM) with 12.5 m spatial resolution was used as the reference raster for spatial alignment, and all thematic layers were converted into raster format with the same cell size and extent. The study area boundary of the Panjkora River Basin was used as a mask to extract all input layers and maintain a consistent spatial extent for analysis. Missing pixels and areas outside the basin boundary were excluded from the analysis and treated as NoData values. The thematic layers (rainfall, slope, drainage density, lineament density, soil, geology, and land use/land cover) were reclassified into suitability classes using the Jenks Natural Breaks classification method, and corresponding ranks/weights were assigned based on the MIF and AHP approaches. Table 1 shows the data sources.

Table 1: Data sources and characteristics used for rainwater harvesting suitability assessment. Please click here to download this Table.

MIF suitability modeling

Initially, the selection of various parameters is carried out based on literature review¹². To determine suitable locations for RWH, precipitation, lithology, lineament density, drainage density, soil texture, slope, and land use/land cover, were taken into consideration as distinct influencing factors. For this objective, pre-processing of the parameters is carried out to create the parameter's influence scale; the data were then categorized according to their significance to RWH, and the major and minor importance were determined using the multi-influencing factor formula (Equation 1). Table 2 shows the major and minor importance of different factors²² (See the Supplementary File 2)

Table 2: Selected influencing factors and their major and minor influence scores used in the Multi-Influencing Factor (MIF) model. Please click here to download this Table.

The selected Factors were ranked using the relation,

[(X+Y) ÷ ∑(X+Y)] × 100 (1)

where Y denotes the minor effect of factors, and X denotes the major effect. Each factor's major and minor influences are calculated using Equation 1.

The major (X) and minor (Y) influence scores were assigned based on previous studies and the hydrological significance of each factor in controlling runoff generation, infiltration, and rainwater harvesting potential¹². A major influence was assigned to factors having a direct impact on RWH suitability, while minor influence represented indirect relationships among the controlling parameters. The factor weights were calculated using Equation (1) by normalizing the combined major and minor influence scores. The subclass weights were assigned according to their relative contribution to runoff accumulation, infiltration capacity, water retention, and suitability for RWH structures. This approach ensured a transparent and reproducible weighting framework for GIS-based suitability analysis.

Relative importance based on Saaty’s scale is shown in Table 3.

Table 3: Saaty’s relative importance scale used for Analytic Hierarchy Process (AHP) analysis. Please click here to download this Table.

The thematic levels scores of all the parameters are combined, each subclass score of the MIF parameters is listed in Table 4. Using the reclassification technique, the MIF output is classified into five categories for rainwater harvesting. Finally, maps of the ultimate locations suggested for installing different RWH structures, such as check dams, farm ponds, gully plugs, and other conservation-related structures, are generated and analyzed. Figure 2 shows the methodology framework.

Figure 2: Methodological framework for rainwater harvesting site suitability assessment using GIS-based MIF and AHP models. The framework illustrates the acquisition and processing of field-survey, geology, soil, ALOS PALSAR DEM, ESA, and GPM data to derive thematic layers, including geology, soil, slope, drainage density, lineament density, land use/land cover (LULC), and rainfall. These layers were integrated using the Multi-Influencing Factor (MIF) approach to generate the RWH suitability map, followed by field-based validation to produce the final validated maps Please click here to view a larger version of this figure.

Table 4: Multi-Influencing Factor (MIF)-based ranks and weights of thematic factors and subclasses for rainwater harvesting suitability mapping. Please click here to download this Table.

AHP suitability modeling

Analytic Hierarchy Process (AHP) is an effective technique for handling difficult decision-making situations which also helps the decision-maker set priorities and choose the best option²³. The AHP technique is a systematic framework for organizing and evaluating complicated decisions through the application of mathematics and expert knowledge²⁴. The AHP aids in identifying both the subjective and objective aspects of a decision by simplifying complex judgments through pairwise comparisons and then evaluating the results²⁵. There will inevitably be some disparity because the comparisons are based on subjective or individual viewpoints. By calculating the consistency ratio and removing decision-making bias, the AHP technique provides a useful tool for evaluating the consistency of the decision-maker's judgments, ensuring consistency of perceptions. One of the main benefits of the AHP is the consistency ratio, which quantifies the degree of consistency between paired comparisons of different criteria^26,27,28,29. Geographical data inputs are combined and transformed by the AHP into a decision output. Utilizing Saaty's scale (Table 3), qualitative data on various themes and qualities is transformed into quantitative data by generating a pairwise comparison matrix^30,31. The fundamental process comprises setting the goal, considering and assessing the factors or standards that affect the final decision, and using Saaty's scale to assign a rating to each criteria. To check the consistency of assigned weights, the consistency ratio (CR) as suggested by Saaty²³ was computed using Equations 2 and 3:

CR = CI/RCI (2)

where CI is the consistency index, and RCI is the random consistency index.

The consistency index (CI) is given by the equation:

(3)

where n is the number of criteria and λmax is the major eigenvalue. The consistency index's average value is estimated by the random index.

RWH structure selection

Land cover land use (LULC)

Land use characterizes the use of the land, whereas land cover describes the natural features of the land. Important information about runoff spreading is contained in the LULC³². In vegetation-covered areas, higher absorption and infiltration rates are linked to less runoff, whereas bare land and built-up areas foster high runoff formation^33,34. Sentinel 2b satellite data were used to prepare the land-use/land-cover patterns of the study area. The land use of the Panjkora river basin was classified into seven classes namely; water bodies, forest, crop and agricultural land, urban land, bare land, snow/ice and range land. The suitability weights assigned to the different land-use/land-cover classes were based on their influence on runoff generation, infiltration, and rainwater storage potential. Agricultural land received the highest suitability rating because it generally produces moderate runoff and directly benefits from harvested water for irrigation. Barren land was also assigned a relatively high weight because sparse vegetation and exposed soil surfaces promote greater surface runoff compared with densely vegetated areas. In contrast, forested areas were assigned to lower weights because dense vegetation intercepts rainfall, increases infiltration through extensive root systems, and reduces overland flow. Urban areas and existing water bodies were assigned to lower suitability because they either have limited opportunities to construct additional RWH structures or are already occupied by impervious surfaces or existing water bodies (Figure 3A).

Drainage density

An area's groundwater infiltration and water runoff are described by the drainage density. The subsurface hydrological formation and surface characteristics are both reflected in drainage density. It shows the tightness of channel spacing and the characteristics of the surface material. Runoff decreases with decreased drainage density and vice versa¹². Lower infiltration and lower runoff are generally found in areas with low drainage density, and vice versa. Dense drainage networks are essential for rainwater collection. RWH are better suited to areas with higher drainage densities because they provide a system that allows water to flow and be rapidly conveyed to a point of collection^34,35. The Drainage density of the Panjkora River basin was classified into five classes on the basis of Jenks Natural Breaks classification: 0–9.4907, 9.4907–27.207, 27.207–48.219, 48.219–79.089 and 79.089–161.34 km/km2. Zones with low to moderate drainage densities were assigned a higher weighting value because they are considered as ideal locations for rainwater harvesting (Figure 3B).

Lineament density

Lineaments are linear subsurface features that are typically derived from geological maps and are also visible on satellite images. Lineaments (buried beneath zones of localized or structural weathering) exhibit enhanced porosity and permeability¹². The lineaments were extracted from Landsat 8 image using remote-sensing image-processing software. The line-density tool was employed to generate the lineament raster layer. The lineament density was further classified using Jenks Natural Breaks classification method into five classes: 0.0072-0.406 km/km², 0.406-0.664 km/km², 0.664-0.921 km/km², 0.921-1.33 km/km² and 1.33-2.13 km/km² (Figure 3C).

Soil

Soil texture is an important factor with respect to RWH planning and site selection. The soil's capacity for infiltration is determined by its texture. In general, sandy soils generate low runoff as compared to clayey soil³⁶. The percentages of silt, sand, and clay define the textural class of the soil. Clay soil has poor permeability and can retain the collected water, so areas with medium- and fine-grained soil were often preferred for rainwater collection^8,37. The study region is characterized by five soil textures: Glaciers and snow caps, loamy, non-calcareous clay soil, Loamy shallow non-calcareous soil, Loamy very Shallow soil, and rock outcrops (Figure 3D).

Slope

Infiltration and runoff are significantly impacted by topography⁸. The catchment's variation in slope has a clear impact on how water flows during and after a downpour. Building RWH structures in areas with steep slopes is not cost-effective due to the significant amount of earthwork needed³⁸. For a high RWH potential, a gentle slope is the most suitable location. RWH structures are not durable in areas with steep slopes (slopes greater than 5%)³⁹. Erosion control measures are also considered in areas with steeper slopes⁴⁰ . The slope was calculated in degrees, and the study area was divided into five classes using Jenks Natural Breaks classification: 0°–11.9°, 12°–22.5°, 22.6°–31.8°, 31.9°–42.4° and 42.5°–82° (Figure 3E).

Rainfall

Rainfall is the primary component that produces surface runoff. Rainfall/ Precipitation data, the Global Precipitation Measurement (GPM) data were downloaded from NASA https://gpm.nasa.gov/missions/GPM from 2014 to 2023⁴¹. The GPM rainfall data from 2014-2023 time period and Jenks Natural Breaks classification were used to classify the study area into five classes of rainfall (mm): 49.93–57.014, 57.014–61.773, 61.773–65.262, 65.262–68.646 and 68.646–76.894 (Figure 3F).

Geology

The physical composition of a watershed and the amount of soil it produces are greatly influenced by the geology of the area. Geological features control the flow of water into subsurface aquifers⁴⁰. Sedimentary and metamorphic rocks are the two major types of rocks found in the current study area. Lithology was broadly divided into Lower Paleozoic rocks, Carboniferous sedimentary rocks, Cretaceous sedimentary rocks, Mesozoic intrusive and metamorphic rocks, Triassic rocks, undivided Paleozoic rocks, undivided Paleozoic rocks and undivided Precambrian rocks, and undivided Silurian rocks. The availability and storage capacity are significantly influenced by the type of lithology; certain rocks have the ability to percolate surface water and replenish the aquifer⁴¹. On the other hand, some rocks allow water to pass through and help recharge subsurface water. Lithology strongly controls runoff generation through its effects on permeability, porosity, and infiltration capacity. In the Panjkora Basin, compact metamorphic rocks generally exhibit lower primary porosity and permeability than unconsolidated or highly porous sedimentary deposits. Consequently, rainfall is less likely to infiltrate and more likely to generate surface runoff, making these formations more suitable for surface rainwater-harvesting structures. In contrast, sedimentary formations containing coarse-grained or sandy materials generally permit greater infiltration, thereby reducing surface runoff available for storage. Therefore, higher suitability weights were assigned to metamorphic rocks, whereas relatively lower weights were assigned to sedimentary formations. Figure 3G shows the geological map of Panjkora river basin. All data are available in Supplementary Files 1 and 3.

Figure 3: Spatial distribution of the thematic factors used for rainwater harvesting site assessment in the Panjkora River Basin. (A) Land use/land cover, (B) drainage density, (C) lineament density, (D) soil texture, (E) slope, (F) rainfall, and (G) geology. Different colors represent the respective classes of each thematic factor. Please click here to view a larger version of this figure.

Results

Thematic map reclassification for MIF model

Land cover land use and RWH

The land use of the Panjkora river basin was classified into seven classes, namely, Water bodies, Forest Land, Crop and Agriculture land, Urban land, Bare Land, Snow and ice and range land Regions where RWH is vital to life have been assigned a high weight based on literature. Agricultural and barren land were assigned higher suitability because these land-cove...

Discussion

This study compared two widely used GIS-based multi-criteria decision-making approaches, namely the Multi-Influencing Factor (MIF) and Analytic Hierarchy Process (AHP), for identifying suitable rainwater harvesting (RWH) sites in the mountainous Panjkora River Basin of northern Pakistan. Although both models successfully delineated potential harvesting zones, noticeable differences were observed in their spatial predictions and validation performance. The ROC-AUC analysis demonstrated that the MIF model (AUC = 0.724) per...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
ArcGIS Desktop	Esri	Version 10.3.1	Used for GIS-based spatial data preprocessing, thematic-layer generation, reclassification, weighted-overlay analysis, and suitability mapping.
Digital elevation model (DEM), 12.5 m	Alaska Satellite Facility (ASF) Distributed Active Archive Center	N/A	12.5 m spatial-resolution DEM used as the reference raster and to derive slope, drainage network, drainage density, and elevation layers; acquisition date reported as 2/12/2023 in the manuscript and should be made unambiguous.
Geological data	Geological Survey of Pakistan	N/A	Geological/lithological data used to generate the geology thematic layer.
Global Precipitation Measurement (GPM) precipitation data	National Aeronautics and Space Administration (NASA)	2014-2023 dataset	Precipitation data used to generate the rainfall thematic layer. The exact GPM product/version and temporal statistic should be specified by the authors.
Global Positioning System (GPS) field-validation points	Not specified	N/A	Field-validation locations used for cross-verification and ROC-AUC assessment. GPS receiver/app manufacturer and model were not specified in the manuscript.
Landsat 8 imagery	Source not specified	Product/scene ID not specified	Used for lineament extraction. Authors should provide the source repository, acquisition date, product/scene ID, spatial resolution, and preprocessing details.
Sentinel-2B satellite imagery	European Space Agency (ESA), Copernicus	Acquired 10 September 2025	Used for land-use/land-cover mapping with supervised Maximum Likelihood Classification. Authors should add the exact product/scene identifier and processing level.
Soil texture data	Directorate of Soil Survey, Khyber Pakhtunkhwa, Pakistan	N/A	Used to generate the soil-texture thematic layer. Dataset/map edition or identifier was not specified in the manuscript.