*3.12. Plan Curvature*

Plan curvature explains the geometry of a particular region. It helps understand the way contours intersect the horizontal region and their impact on the slope inclination of a particular zone [70]. It explains the flow of groundwater and helps establish a generalized flow pattern. Plan curvature approximates the inclination of various zones that impacts groundwater recharge through topographical influence [111]. The inclination of the area is marked with a slope that runs from the region of higher inclination towards the lower inclination, thus indicating groundwater flow [100,109]. The region of Islamabad is higher in inclination towards the northern region that goes down towards the southern zones. This is because the northern area is comprised of mountainous regions, and the southern zone consists of high-populous flat regions, establishing a generalized pattern of inclination decrease [111]. The inclination and gentle slopes and the presence of flat surfaces greatly influence groundwater recharge [108]. Therefore, plan curvature has been included as a key factor in this study.

## *3.13. Profile Curvature*

Profile curvatures define the nature of the ground zones under study: linear, concave, and convex. It is defined as the line parallel to the direction of the maximum slope. Patterns might indicate a general linear formation with a defined value approaching zero. A positive value indicates an upward concave profile, while the negative region represents an upward convex profile [70]. The profile curvature helps classify the area into lower or higher waterretention zones depending upon its convexity and concavity. Accordingly, the regions comprising elevated convex profiles within center zones are regarded as less water holding and vice versa [110]. The curvature of the study area is included in this study to assess its effect on the water-retention capability of the zone following related studies [109,110]. Considering the variability of Islamabad's surface in terms of slope and elevation, it is important to consider the influence of profile curvature on groundwater recharge in this region. Therefore, this factor has been used in the current study.

#### *3.14. Distance to Fault*

Faults describe the change in geological composition in a particular zone [106]. These indicate the movement and change a particular rock surface has undergone in a specified period. For example, earthquake-induced faults can indicate rapid geological movement beneath the surface. The parameters of faults can have vast ranges. Distance to faults impacts the flow and occurrence of groundwater [108]. It is important, as it indicates groundwater flow and can highlight the zones contributing to underground-water flow [106,109]. In our study area, Islamabad and nearby regions have more faults that influence the groundwater recharge. Thus, this factor is included in the current study.

#### *3.15. Fault Density*

The magnitude of faults (density) indicates the potential groundwater regions. In a similar study, lineaments such as faults have been reported to impact the groundwater recharge potential zones and are considered key IF [108]. Fault density helps determine the occurrence and movement of groundwater beneath the surface. Many relevant studies have included fault density as a key factor in assessing groundwater recharge potentials [60,65,66]. As previously discussed, Islamabad has higher faults than the rest of the country. Therefore, fault density is included as a key factor in the study.

## **4. Methodology**

The current study follows a four-step approach. In the first step, the relevant thematic layers are identified. First, the thematic layers used for the study were extracted that act as input data for the eventual delineation of recharge zones. These thematic maps present the geographical map of the study region in accordance with the subject matter. The current study utilizes thematic maps for 15 hydrological factors. These include distance to faults, land use, lithology, drainage density, slope, soil, rainfall, plan curvature, fault density, profile curvature, TWI, elevation, aspect (the front-facing direction of a slope), drainage distance, and slope length. These factors were extracted from previous literature [60,66,87,103,113] considering the geological properties of the study area as listed and are discussed in Section 3 of the study.

The thematic maps used in the research were generated at a 1:200,000 scale considering that this would eventually increase accuracy. In addition, the majority of the data sets were available at this scale. The differing scales were later normalized for the sake of uniformity. The digital elevation model (DEM) data are used on a global scale at 30 m × 30 m resolution for topographic analysis. This resolution is highly important, as it contributes to how sharply the objects can be seen in an image. It represents the size of the tiniest feature captured by a satellite sensor or portrayed in a satellite photo. It is commonly expressed as a single number representing the length of one of the sides of a square (grid) [12]. In addition to the normalization of the input data, uniformity is ensured in their format for easy integration of these thematic maps into the GIS platform. For this purpose, the acquired maps are converted into raster form before integration with the GIS.

In the second step, the pre-processing of the thematic layers was performed to ensure uniform projection and resolution. This is followed by the assignment of scores and suitable weightage to each factor. During weightage overlay analysis, the ranking was given for each parameter of each thematic map, and weights were assigned according to the influences (following IF technique) of the feature on the hydrogeological environment of the area coupled with that parameter's contribution toward the groundwater recharge as shown in previous researches [65,108,109].

The IF technique was used to assign scores and get a diverse and error-free outcome. A diversely produced thematic map considers the input from multiple hydrological procedures, thus not relying on a single hydrological process where the outcome can be manipulated and is prone to error. Moreover, due to finer resolution, any errors in the weighted overlay analysis within the ArcGIS were eliminated since such resolutions result in finer interpretation.

The third step involves using ArcGIS to deploy the thematic layers to get the processed images containing the potential zones. In this step, all the scored thematic maps along are integrated by employing the "Spatial Analysis tool" in ArcGIS 9.3, whereby rankings are assigned to all the thematic maps. Then, these weighted thematic maps are overlaid using ArcGIS to highlight the potential recharge zones. In the fourth (last) step, the study area was categorized based on the potential groundwater rechargeability into five different classes: poor, low, medium, high, and best in terms of their capability for the groundwater recharge potential.

Figure 2 shows a flowchart summarizing the methodology used in this study. The associated steps include acquiring the data, converting to raster, preprocessing (confirming projection and resolution coupled with assigning scores and weights), integrating GIS for final output, and categorizing the study area based on groundwater recharge capability. Figure 2 also shows the source of the acquired data. Accordingly, the thematic maps are acquired from Landsat-8 TM Satellite, Aster DEM, and soil and geological maps of Pakistan. The following sections explain the IFs used in this study in detail, their sources, and the procedure for assigning weights to each of these factors. ity. Figure 2 also shows the source of the acquired data. Accordingly, the thematic maps are acquired from Landsat‐8 TM Satellite, Aster DEM, and soil and geological maps of Pakistan. The following sections explain the IFs used in this study in detail, their sources, and the procedure for assigning weights to each of these factors.

*Water* **2022**, *14*, 1824 12 of 30

shown in previous researches [65,108,109].

in finer interpretation.

groundwater recharge potential.

influences (following IF technique) of the feature on the hydrogeological environment of the area coupled with that parameter's contribution toward the groundwater recharge as

The IF technique was used to assign scores and get a diverse and error‐free outcome. A diversely produced thematic map considers the input from multiple hydrological pro‐ cedures, thus not relying on a single hydrological process where the outcome can be ma‐ nipulated and is prone to error. Moreover, due to finer resolution, any errors in the weighted overlay analysis within the ArcGIS were eliminated since such resolutions result

The third step involves using ArcGIS to deploy the thematic layers to get the pro‐ cessed images containing the potential zones. In this step, all the scored thematic maps along are integrated by employing the "Spatial Analysis tool" in ArcGIS 9.3, whereby rankings are assigned to all the thematic maps. Then, these weighted thematic maps are overlaid using ArcGIS to highlight the potential recharge zones. In the fourth (last) step, the study area was categorized based on the potential groundwater rechargeability into five different classes: poor, low, medium, high, and best in terms of their capability for the

Figure 2 shows a flowchart summarizing the methodology used in this study. The associated steps include acquiring the data, converting to raster, preprocessing (confirm‐ ing projection and resolution coupled with assigning scores and weights), integrating GIS for final output, and categorizing the study area based on groundwater recharge capabil‐

**Figure 2.** Flowchart for potential groundwater assessment using integrated remote sensing and GIS techniques.

#### **Figure 2.** Flowchart for potential groundwater assessment using integrated remote sensing and GIS *4.1. Acquisition of Thematic Maps for Contributing Factors*

techniques. Table 3 below enlists the sources for acquiring thematic maps for all the contributing factors. The soil thematic map was generated using the Soil Map of Pakistan [114]. Land use, rainfall, and TWI thematic maps were generated using Landsat 8TM satellite data. Drainage distance, slope, plan curvature, profile curvature, slope length, elevation, drainage density, and aspect thematic maps were generated using ASTER global DEM. Finally, distance to faults, lithology, and fault density thematic maps were generated using data from the geological map of Pakistan on a scale of 1:200,000 [115].


**Table 3.** Acquisition of Thematic Maps for Contributing Factors.

These thematic map data were cross-checked using ground surveys for cross-validation. The imagery was visually interpreted to delineate rainfall, land use, and other factors with the help of slandered characteristic image-interpretation elements such as tone, texture, shape, size, pattern, and association using the Landsat 8 satellite data products. These data sets are used for assessing groundwater recharge potential [65,109,111].

#### *4.2. Weightage Assignment via IF Technique*

Weights and rates were assigned to the factors to obtain a final combined recharge potential map. Using the IF technique, the influence of various factors was taken into account, and the level of impact they have on the hydrological aspect of groundwater flow and its occurrence was assessed. A weightage approach was included as used by [65] to assign weightage to the factors that would ultimately define the control they can assert over the groundwater recharge of the study area. The current study follows a similar approach. In assigning weights to the considered factors, five major descriptive levels were plotted for each factor ranging from very high to very low, including some interrelated levels. These weightage values range from 10 to 1 point, i.e., a very high range is assigned a score of 10, and the minimum level is 1 following relevant groundwater studies [113,116]. These weights for each factor were assigned based on their degree of impact on groundwater recharge as extracted from relevant literature [6,10,63,113].

## **5. Results and Discussions**

This section presents the results and discussions in line with the adopted method.

#### *5.1. Spatial Analysis of Considered Key Factors*

Figure 3 represent the resulting thematic maps of the 15 considered factors for the current study area. Figure 3a highlights the wells or water extraction points in the study area. These are primarily located in the residential zones and plain areas of Islamabad. Figure 3b shows the thematic map of rainfall for Islamabad. The resulting map highlights that Islamabad receives ample rainfall. Further, it shows a rhythmic increase in rainfall volume from south to north. The northeast outskirts receive the highest rainfall, consisting of regions from Rawat to Crore Village. Low-rainfall regions are evident in the southwest. Considering the high rainfall in the northeastern regions, there are more chances for more groundwater recharge and high groundwater levels in alluvial plains [64], thus displaying a higher potential for groundwater recharge. Moreover, the map shows that around 44% of the area receives less than 882 mm of rainfall, 16% area receives rainfall between 882–999 mm, 10% area receives rainfall between 999–1116 mm, 9% area receives rainfall between 1116–1233 mm, while 21% of the area receives most rainfall ranging between 1233–1350 mm. This shows that around 40% of Islamabad receives good rainfall. This assessment can help policymakers preserve the natural terrain in the region receiving more rainfall and utilize it for groundwater recharge.

Figure 3c shows Islamabad's thematic layer of plan curvature data. The figure categorizes the regions based on concavity and convexity. The map shows that the northeast region of Islamabad is composed of higher convexity, whereas a systematic decrease in convexity is observed from north to south. This indicates a higher surface and altitude in the north and a gradual decrease towards the south. This heavily contributes to the groundwater flow from north to south, where a gentler slope and plain area can accumulate this water and get recharged. A similar study accounted for alluvial plain and gentle slopes to be more promising for groundwater potential due to large infiltration rates, high porosity, and permeability [116].

Figure 3d shows the thematic layer of soil data for Islamabad, where the region is classified based on soil composition. Soil types impact groundwater flow directly, but they also impact other important phenomena, such as infiltration [117], which ultimately impact groundwater recharge. The soil conditions define permeability, which impacts groundwater infiltration and soil porosity. For example, the calcareous loamy soil is abundant in arid and densely populated areas. Figure 3d shows that the mountainous soil forms the northern edge of Islamabad that receive a decent amount of rainfall. Such soil helps infiltration, enabling the groundwater to flow towards the inner zones. While no definite pattern exists throughout the study area, calcareous soils are mostly reported for various regions.

**Figure 3.** *Cont*.

**Figure 3.** Thematic layers of selected factors for Islamabad's data (part 1), (**a**) well data, (**b**) rainfall data, (**c**) plan curvature, (**d**) soil data, (**e**) distance to fault, (**f**) drainage distance, (**g**) profile curvature, (**h**) TWI, (**i**) slope, (**j**) elevation, (**k**) slope length, (**l**) lithology, (**m**) land use, (**n**) fault density, (**o**) drainage density, (**p**) aspect. **Figure 3.** Thematic layers of selected factors for Islamabad's data (part 1), (**a**) well data, (**b**) rainfall data, (**c**) plan curvature, (**d**) soil data, (**e**) distance to fault, (**f**) drainage distance, (**g**) profile curvature, (**h**) TWI, (**i**) slope, (**j**) elevation, (**k**) slope length, (**l**) lithology, (**m**) land use, (**n**) fault density, (**o**) drainage density, (**p**) aspect.

Figure 3e shows the thematic layer map of distance to fault for Islamabad. This map categorizes regions with respect to distances to faults. Considering that the faults act as points with more recharge capability, more distance from faults implies less recharge capability and vice versa. In this respect, Figure 3e shows that the major faults are all located on the outskirts of Islamabad. The zones comprising convex geological features and landscapes have nearby faults, whereas the southern regions comprising more land use and less geological convexity comprise low distances to faults. This aligns with several studies that have established patterns with lineaments and groundwater recharge potential [10,118]. Overall, the southern regions with less distance to faults display more recharge potential in the current study area.

Figure 3f shows the thematic layer of drainage distance for Islamabad. It categorizes the study area based on the distance of various zones from the drainage networks. Figure 3f shows that the study area comprises abundant and closely located drainage networks. However, there is no defined pattern for the drainage distances in the study area. Considering that a lesser distance from the drainage pathway displays higher groundwater recharge potential [119], the drainage distance thematic map suggests that the study area has a larger potential for groundwater recharge. Further, there is a well-distributed groundwater flow throughout the region. Figure 3g shows the thematic layer of profile curvature for Islamabad. It highlights the geological characteristics of Islamabad and depicts the concavity and convexity of the region. It is indicated that the outskirts of the northern region are higher in altitude and contribute to the groundwater flow under gravity. A higher profile value indicates a rising elevation, ensuring a systematic flow towards the inner edges with the highest and densest land use in Islamabad. This is in line with a previous study's findings that suggest a higher potential of gentle slopes for groundwater recharge [116].

Figure 3h shows the thematic layer of TWI for Islamabad. The TWI map shows the impact of geology on the hydrological aspects. The outskirts, shaded in deep blue in Figure 3h, show the zones with geological makeup that impact regional hydrology. Following our thematic maps for the land use, well data, and rainfall, the TWI highlights Islamabad's northern outskirts as the areas directly reaching the groundwater. The inner edges with lower index value contribute little to the groundwater flow, while the geological makeup of the outermost skirts contributes greatly to the groundwater flow towards the center, housing the area with the highest and densest land use. A direct relationship between the higher TWI value was also established by another study [120]. Following our findings, a higher TWI value suggests a better groundwater recharge potential in the Islamabad region.

Figure 3i shows the thematic layer of slope data for Islamabad, showing that the northern outskirts of Islamabad have the highest slope. The slope plays an important part in determining the runoff direction of groundwater. The thematic map indicates that 23% of the region has a slope greater than 48 degrees, 38% has a slope ranging from 36 to 48, 16% has a slope ranging from 24 to 36, 9% has a slope ranging from 12 to 24, and 4% of the region has a slope less than 12 degrees. The figure shows that the outskirts of Islamabad in the northern region comprise the highest slopes due to mountains that promote a rapid runoff towards the south. While some water is lost during the runoff, infiltration takes water to the deep soil layers, contributing to recharging the local groundwater table.

Islamabad's outskirts comprise Attock, Wah Cantt, and Taxila in the west; Murree in the northeast; Haripur in the north; Gujar Khan, Rawat, Mandrah, and Kahuta in the southeast; Rawalpindi to the south and southwest; and other Punjab regions in the east. The greater slope in the northern region ensures a flow of water towards the south with the highest settlement and greatest water recharge potential.

Figure 3j shows the thematic layer of elevation data for Islamabad. Islamabad is high on the northern edge due to the mountains that decrease towards the south. The area with residential zones, i.e., the inner edges, and that towards Rawalpindi has higher population density and low elevation. This systematic decrease of elevation contributes directly to the groundwater flow as the water flows under the action of gravity. The higher elevation area also receives greater rainfall, as shown in our rainfall thematic map, ensuring infiltration and surface runoff towards the inner edge. Thus, the area with higher elevation retains rainwater for a lesser time duration and generates more runoff towards the residential areas in Islamabad, in line with published studies [64].

Figure 3k shows the thematic layer of slope length data for Islamabad that highlights the lengths of slopes in the region. Longer slope lengths are evident on the northern outskirts, while a rhythmic slope length decrease can be observed towards the south. The area with the highest land use comprises regions with lower slope length values. The groundwater flows from the northern sides with the highest slope lengths promoting recharge potentials and infiltration. The gradually decreasing slopes towards the center help with groundwater recharge to meet the requirements of the local population.

Figure 3l shows the thematic layer of lithology data for Islamabad that shows regions with limestone and unconsolidated deposits to be abundant in the area. However, there is no defined pattern, and the data are scattered throughout the region. The concentrated regions are highlighted in red, green, and purple colors in Figure 3l. There is a presence of sandstone in the northeastern region along the dense mountainous regions that continues towards the northwestern region.

Further sandstone and unconsolidated deposits are seen within the areas of highest land use towards the southwest. Past glacial activity has contributed to the unconsolidated deposits in the region due to the weathering of rocks. A previous study also established a pattern between the weathering of rocks towards the increased groundwater recharge potential [121]. An increased recharge was also observed in the area of higher unconsolidated deposits in another study [120]. Accordingly, there is a greater potential for groundwater recharge in the study area.

Figure 3m shows the thematic layer of land-use data for Islamabad, showing areas such as bare land, water bodies, built up, and vegetative regions. Such a map displays the variation of population density and associated water demand throughout the study area [10]. The thematic map for Islamabad indicates that 4% of the region comprises bare land, 36% is built-up region, 51% is vegetative, while 9% of the study area is composed of water bodies. Further, it can be observed that most of the built-up region is around the inner region of Islamabad. This region falls towards the city of Rawalpindi, which has a far greater population density than Islamabad. The runoff from the northern region infiltrates into the groundwater table around these internal regions, where there is a greater need for water.

Figure 3n shows the thematic layer of fault density for Islamabad that highlights geological features induced by the movement of rock bodies. These faults govern groundwater flow following their complex and favorable topography. Accordingly, the fault densities for the area include 43 % area with less than 15 fault density, 6% area ranging from 15 to 41, 35% ranging from 41 to 65, 7% ranging from 65–91, while 9% of the study area has fault density greater than 91. The map indicates that the northeastern edges of Islamabad consist of lower-density faults than the northwestern region, where there are more mountains. The maximum land use is towards the internal regions with no major geological faults. Previous studies have linked fault-dense regions with higher groundwater recharge potential [10]. Thus, there is a higher recharge potential in the northwestern areas of Islamabad.

Figure 3o shows the thematic layer of drainage density for Islamabad that highlights the northeastern regions to have streams or rivers with relatively long lengths. This ensures a deep-water flow towards the inner edges of Islamabad. Thus, the northeastern region contributes majorly towards the groundwater flow in the areas of highest land use. Further, a flow from the northern to the southern edge is seen with the major contribution from the northeastern region. A previous study showed that high-density drainage regions have greater groundwater recharge potential [122]. This is in line with the current study where major water sources contribute to the water recharge. The same has been highlighted by the rainfall thematic map, where the northeastern region receives most of the rainfall and has a high drainage density, thus contributing to the groundwater flow and recharge.

Figure 3p shows the thematic layer of aspect data for Islamabad that categorizes regions based on their compass directions. The aspect map lists out the front-facing direction of regions along with the compass. For example, the major constituting region in the northeast contains southeastern front-facing regions that align with thematic maps of land use and wells in the study region. The dense regions with most residential and commercial zones are in the southeast. The flow from the north region is ensured towards the southeast region. The southeastern compass front directions of the geological regions act as a gentle slope that promotes groundwater recharge in Islamabad [116].

These influencing factors were considered based on a literature review and classified based on their impact on groundwater recharge contribution, i.e., the class at which lesser the groundwater recharge potential would rank lower and vice versa. For example, a higher slope would have lesser groundwater potential, or a lower TWI would mean low water moisture and low groundwater recharge potential; hence, these classes would have lesser weightage.

#### *5.2. Weightage Calculation for Influence Factor (IF) Techniques*

After obtaining the individual thematic maps for each of the contributing factors, these factors were integrated to obtain a potential holistic map that highlights the recharge potential of Islamabad. Accordingly, weights and rates were assigned to the 15 key factors. For incorporating the mutual influence of the factors, rate values were assigned to them. Two points were given for every major effect, while one point was given to the corresponding factor for each minor effect. The cumulative weightage of both major and minor effects was considered for calculating the relative rate, as shown in Table 4. Table 4 shows that factors such as lithology influence six of its fellow factors majorly. It has a noticeable impact on the lineament, drainage, land/use, slope, and soil types. Thus, it has been assigned a value of 2 six times (2 × 6 factors).


**Table 4.** Relative rates and scores for each potential factor.

Similarly, other factors have also been assigned their respective rate values using the same approach. Overall, the major effect (*A*) and minor effect (*B*) are summed for all factors, and their cumulative sums are calculated for each factor to get the proposed relative rates. The cumulative proposed relative rates sum up to 160. Using this value, the normalized relative rates are calculated, where the proposed relative rate of each factor is divided by the cumulative proposed related rates and multiplied by 100 using Equation (1). The values are rounded off to the nearest integer.

*Normalized relative rates* (*Y*) = *Proposed relative rates* (*<sup>A</sup>* <sup>+</sup> *<sup>B</sup>*) *Commulative Proposed relative rates* (Σ(*<sup>A</sup>* <sup>+</sup> *<sup>B</sup>*)) <sup>×</sup> 100 (1)

> After the assignment of rate values, the next step is to assign weights. In this process, five major descriptive levels are plotted for each factor ranging from very high to very low, including some interrelated levels as shown in Table 5. Factors contributing majorly, such as rainfall, can be seen as very dominant in relevant studies [108,109] and in abundance in the southeastern regions of the study area and thus were assigned higher weights. In contrast, factors such as profile curvature were assigned a lower weightage, as the area followed a rhythmic curvature, and the influence of curvature was not dominant in terms of groundwater flow, as evident from Figure 3 (previously shown).

> With a weightage of 8.1%, rainfall is a dominant factor in the southeastern parts of the study area. Plan curvature data indicate a slight shift in curvature as seen from the thematic map and thus were assigned a weightage of 3.1%. The higher curvature would result in a greater flow of water beneath the surface [66]. Soil is the primary factor that controls seepage and the associated groundwater recharge [117]. Thereby, it was assigned the highest weightage (8.1%). Likewise, faults being the primary indicator of geographical movement (earthquakes or tectonic) over the years indicate a weaker and vulnerable zone suspectable to the greater flow of groundwater channels beneath the surface. It adds greatly to the groundwater recharge and was hence assigned a weightage of 6.8%. Drainage distance, profile curvature, and TWI were assigned weights of 6.2%, 3.1%, and 5%, respectively.

> The data obtained from thematic maps do not indicate an abrupt or dominant effect of these considered geographical features (key factors) over the study area, thus acquiring a lower weightage in our study area. The slope indicating the natural flow of water towards the lower altitude area was assigned a weightage of 8.1%. Elevation and slope length were assigned the weightage of 6.8% and 3.7%, indicating the flow towards lower-elevated areas and the flow speed. Accordingly, the lower the speed, the greater the infiltration and vice versa [122]. Lithology has been assigned a weightage of 10%. It indicates the rock characteristics that dictate the water flow beneath the surface in channels and streams. Land use is another primary factor that was assigned 10% weightage. It has been utilized by several related studies [60,66]. Finally, fault and drainage densities and aspects indicated the magnitude of faults, drainage networks, and front-facing direction of slopes signifying the flow of groundwater beneath the surface and were assigned weights of 6.2%, 9.3%, and 5%, respectively, in this study.

> After the assignment of rates and weights, the % influencing score was calculated using Equation (2). The % influencing score is defined as the percentage of factor effect on recharge potential (%) and is shown in Table 5 for each factor, where *X* is the normalized weight from 1 to 10, and *Y* is the rate from 1 to 10.

$$\% \text{ of fluorancing score} = \frac{\text{Total Weightage } \Sigma (X \times Y)}{\text{Grad Total Weight} \left(GTW\right)} \times 100 \tag{2}$$



**Table 5.** Weight evaluations of factors influencing potential recharge capacity.


*Water* **2022**, *14*, 1824


*Water* **2022**, *14*, 1824

#### *5.3. Final Combined Recharge Potential Map 5.3. Final Combined Recharge Potential Map*

weight from 1 to 10, and *Y* is the rate from 1 to 10.

% ൌ

*Water* **2022**, *14*, 1824 23 of 30

After considering rate assessment, different layers of recharge potential were superimposed in the ArcGIS tool. As a result of the integration of the 15 contributing factors, the final combined potential map was generated, which highlights the overall recharge potential of Islamabad, as shown in Figure 4. The resulting map generated with the help of influencing factors' relative rates categorizes the region into five descriptive levels based on the rechargeability. These descriptive levels include "best", "high", "medium", "low", and "poor", each with a distinctive color. After considering rate assessment, different layers of recharge potential were super‐ imposed in the ArcGIS tool. As a result of the integration of the 15 contributing factors, the final combined potential map was generated, which highlights the overall recharge potential of Islamabad, as shown in Figure 4. The resulting map generated with the help of influencing factors' relative rates categorizes the region into five descriptive levels based on the rechargeability. These descriptive levels include "best", "high", "medium", "low", and "poor", each with a distinctive color.

After the assignment of rates and weights, the % influencing score was calculated using Equation (2). The % influencing score is defined as the percentage of factor effect on recharge potential (%) and is shown in Table 5 for each factor, where *X* is the normalized

ℎ ሺ ൈ ሻ

ℎ ሺሻ <sup>100</sup> (2)

**Figure 4.** Potential groundwater recharge zones in the study area. **Figure 4.** Potential groundwater recharge zones in the study area.

From the output thematic map (Figure 4), it is evident that the eastern region of the study area is the most suitable for groundwater recharge. Accordingly, it is highlighted to be the "best" region. This region received the highest rainfall as per the previously pre‐ sented maps. This is in line with previous studies that argued that the higher the rainfall, the greater the groundwater recharge and vice versa [116,123]. Moreover, it can be ob‐ served from Figure 4 that the groundwater recharge potential decreases as we head to‐ wards the western side of Islamabad. A decreasing pattern for groundwater recharge is seen as we move from east to west in the study area. Most of the mountainous region is located towards the northeast of Islamabad, receiving the highest rainfall and having higher slopes, inducing rapid runoff. Towards the center and to the west, the slope length decreases, thus indicating a higher recharge potential, as gentle slopes were attributed to higher recharge potential [122]. From the output thematic map (Figure 4), it is evident that the eastern region of the study area is the most suitable for groundwater recharge. Accordingly, it is highlighted to be the "best" region. This region received the highest rainfall as per the previously presented maps. This is in line with previous studies that argued that the higher the rainfall, the greater the groundwater recharge and vice versa [116,123]. Moreover, it can be observed from Figure 4 that the groundwater recharge potential decreases as we head towards the western side of Islamabad. A decreasing pattern for groundwater recharge is seen as we move from east to west in the study area. Most of the mountainous region is located towards the northeast of Islamabad, receiving the highest rainfall and having higher slopes, inducing rapid runoff. Towards the center and to the west, the slope length decreases, thus indicating a higher recharge potential, as gentle slopes were attributed to higher recharge potential [122].

Table 6 presents the data of each category shown in graphical form in Figure 4 and gives the exact portions of the study area having best to worst recharge capability. It shows that the area labeled under the "best" comprises 136.8 km<sup>2</sup> , covering 15% of the study area. Similarly, an area of 191.52 km<sup>2</sup> falls under our map's "high" classification, covering 21% of the study area. Another 35% of the region collectively serves as a competent region (preferred) for groundwater recharge. The moderate zone covers 218.88 km<sup>2</sup> of area, covering 24% of the study area. In contrast, the potentially poor and low zones make up 13% and 27% of the area, i.e., 118.56 km<sup>2</sup> and 246.24 km<sup>2</sup> , respectively.


**Table 6.** Classification of potential recharge areas.

The results show that around more than half (51%) of the total area of Islamabad does not have sufficient recharge capability, and the city is dependent on only 35% of the total area to fulfill the city's demand for groundwater for daily life usage. This can be taken into consideration by local authorities when planning to meet the local water requirements and groundwater recharge. The city planners and policymakers should take mitigation steps and devise strategies to preserve most of this 35% of the land to avoid any further damage to the already fragile water condition of the city. The information devised from this final groundwater potential zones map can help resolve the long due water shortage issues in various sectors of Islamabad and nearby areas through efficient management and preservation of groundwater resources in the area. Compared to the previous studies [36–42], this study addresses the research gap of applying this methodology in a non-coastal region and modifies it by using thematic maps of larger spatial scale and the DEM data of smaller resolution to refine the accuracy of the process. All the previous published research used the one-time dataset and map the output. However, these do not depict the true representation of the groundwater recharge. This is because the considered datasets may change temporally, needing more datasets to overcome this limitation. Hence, this study used the annual mean for all datasets, which change with respect to season or time. Secondly, previously published research used limitedly influencing datasets that might not present the actual situation of the study. In the current research, all the contributing factors were analyzed and used to consider the entire situation. Accordingly, the model gives reliable actual output. Moreover, the study also considers more contributing factors than the previous studies to further enhance the accuracy of the output. The research presents a holistic approach that gives comparatively improved results and can be applied to other regions as and when required.

## **6. Conclusions**

Considering the constant increase in groundwater demand in Islamabad with increasing population growth, the decreasing groundwater level has become a matter of concern for the local authorities. This study attempts to develop a groundwater potential recharge zone map of the study area of Islamabad, Pakistan, to help the policymakers devise efficient policies for mitigating this problem.

The methodology involves the integration of RS and GIS to develop a map that highlights the groundwater recharge potential in the study area. In our scenario, 15 key factors were selected based on their contribution to the recharge. These include soil, land use/land cover, drainage distance, slope, rainfall, plan curvature, distance to faults, profile curvature, TWI, elevation, slope length, lithology, fault density, drainage density, and aspect. Thematic maps were generated and overlayed using GIS. A holistic map was devised at the end, comprising input from 15 of the influencing factors and their weights to produce a weighted map. The resulting map categorizes the region into five different descriptive levels, namely poor, low, medium, high, and best, based on the groundwater recharge potential. The results showed that 13% of the area falls in the poor-recharge-potential category, 27% area has a low potential, 24% has medium potential, 21% has high potential, and 15% has the best chance of recharging the groundwater table. Overall, around 35% of the study area is suitable for groundwater recharge, and more than half is unsuitable for such purposes.

This study provides a holistic model with more accurate results than the previous studies by introducing a comparatively greater number of factors and employing the thematic maps of larger spatial scale and DEM data of a smaller resolution. The current study paves the way for future infrastructure development by the concerned authorities to meet the water demand of Islamabad and preserve the precious natural terrain with high recharge potential.

The study is limited in terms of the factors considered. Further, it is restricted to a single region in a developing country for testing purposes. Moreover, considering that this study was limited in terms of the unavailability of geophysical data for the case study area, future researchers can conduct further research by including the geophysical and field data from multiple regions. This can help in carrying out the subsurface groundwater modeling as well as 3D modeling of the targeted study area. Further, similar studies can be conducted for larger nearby regions and developed countries to help move toward global sustainability goals and tackle climate change effects. The effects of vegetation on recharge can also be investigated in the future.

**Author Contributions:** Conceptualization, A.M., B.A. and F.U.; methodology, A.M., B.A., N.K. and F.U.; software, A.M., B.A., N.K. and F.U; validation, A.M., B.A., N.K., F.U., H.A., E.E.H. and A.H.A.; formal analysis, A.M., B.A. and N.K.; investigation, A.M., B.A., N.K. and F.U.; resources, A.M., B.A., N.K., H.A., E.E.H. and A.H.A.; data curation, A.M., B.A., N.K. and F.U.; writing—original draft preparation, A.M., N.K., B.A. and F.U; writing—review and editing, F.U., H.A., E.E.H., A.A.A. and A.H.A.; visualization, A.M., B.A., N.K., F.U. and A.A.A.; supervision, A.M., B.A., F.U. and A.H.A.; project administration, A.M., B.A., F.U. and A.H.A.; funding acquisition, A.H.A. and H.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Taif University Researchers Supporting Project number TURSP 2020/252, Taif University, Taif, Saudi Arabia.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data can be shared upon reasonable request.

**Acknowledgments:** The authors appreciate Taif University Researchers Supporting Project number TURSP 2020/252, Taif University, Taif, Saudi Arabia for supporting this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

