Using Geographic Information Systems and Multi-Criteria Decision Analysis to Determine Appropriate Locations for Rainwater Harvesting in Erbil Province, Iraq

Abstract

Water scarcity is a prominent consequence of global climate change, presenting a significant challenge to the livelihoods of wide parts of the world, particularly in arid and semi-arid regions. This study focuses on Erbil Province in Iraq, where the dual effects of climate change and human activity have significantly depleted water resources in the past two decades. To address this challenge, rainwater harvesting (RWH) is explored as a viable solution. The purpose of this study is to make a suitability zone map that divides the study area into several classes based on the features of each area and its ability to collect rainwater. The map will then be used to find the best place to build different RWH structures. Seven different layers are used to make the RWH suitability zone map: rainfall, runoff, land use/cover (LU/LC), soil texture, slope, drainage density, and the Topographic Wetness Index (TWI). Each layer was assigned specific weights through the Analytical Hierarchy Process (AHP), considering its relevance to RWH. Results revealed four suitability classes: very highly suitable 1583.25 km² (10.67%), highly suitable 4968.55 km² (33.49%), moderately suitable 5295.65 km² (35.69%), and lowly suitable 2989.66 km² (20.15%). Notably, the suitability map highlights the northern and central regions as particularly suitable for RWH. Furthermore, the study suggested three suitable locations for constructing medium dams, six for check dams, and twenty-seven for farm ponds, according to the requirements of each type. These findings provide valuable insights for the strategic planning and effective management of water resources in the study area, offering potential solutions to the pressing challenges of water scarcity.

1. Introduction

Water is an essential and indispensable resource that plays a crucial role in all aspects of our lives, particularly agricultural production [1]. Numerous nations around the world are currently dealing with a serious water crisis, particularly those in arid and semi-arid climates like Egypt, Iraq, Iran, and Syria [2]. Climate change, global warming, and population growth have all exacerbated this issue, resulting in severe water scarcity on a global scale [2,3]. The diminishing availability of water can have adverse effects on agricultural land in various regions worldwide, which poses a major obstacle to food production [4]. Given that agriculture is the primary consumer of water, it is likely to face conflicts with other water users in the future due to escalating demand [4,5]. In the future, ensuring an adequate water supply, particularly for agricultural irrigation in arid and semi-arid areas, will be a challenging task [6].

Erbil Province in Iraq has experienced severe droughts over the past two decades, which have significantly impacted its water resources [7]. Water wells in Erbil province are under stress due to both drought conditions and human activities. In some areas, the groundwater level has dropped by over 54%, leading to the complete drying up of some wells, particularly in the southern and central parts. These wells are crucial sources of drinking water and irrigation [8,9]. Insufficient rainfall during the growing season in the past rendered a significant portion of agricultural land unproductive. Limited water resources also challenged villagers engaged in livestock farming, forcing some to migrate elsewhere [10]. Furthermore, the security instability in some parts of the north and south and the mismanagement of water resources contribute to water scarcity in the study area. Therefore, in response to these challenges, it is crucial to explore strategies for reducing water outflow and increasing water retention [11].

Rainwater harvesting (RWH) is a traditional method used to collect, store, and reuse rainwater for various purposes, with an emphasis on supplementary irrigation [12]. It has recently regained attention in many parts of the world as a viable option for water supply. RWH represents a comprehensive approach to supporting agriculture in regions with limited precipitation during the period of crop growth and scarce water resources [13]. For easier and more cost-effective access to water for irrigation purposes, it is important to scientifically identify appropriate locations for RWH [14].

Various methods have been employed to identify potential locations for RWH. Remote sensing (RS) and geographic information systems (GISs) are currently among the most valuable tools for managing ecosystems and natural resources [15]. Additionally, multi-criteria decision analysis (MCDA) plays a significant role in selecting suitable zones for RWH [16]. The integration of MCDA with GIS, which combines spatial data layers, has been widely utilized in the RWH process [17]. Many authors have created or improved MCDA methods in the last few decades. The Fuzzy Analytic Hierarchy Process (FAHP), the Technique for Order Preference by Similarity to Ideal Solutions (TOPSIS), the Analytic Network Process (ANP), the Analytic Hierarchy Process (AHP), and others are some of these. The key distinctions among these methodologies pertain to algorithmic complexity, criteria weighting methods, the approach to representing preference evaluation criteria, handling uncertain data, and the type of data aggregation [18]. Saaty introduced the Analytical Hierarchy Process (AHP) in 1977 [19], which is a notable example of MCDA. AHP is a highly regarded decision-supporting technique for addressing complex problems [15]. It has been recognized as the most suitable decision method for identifying appropriate RWH locations. Numerous studies have employed GIS-based MCDA in various countries to identify suitable locations for RWH. Notable examples include the work of Aziz et al. (2023) in Iraq [10], Modak and Das (2022) in India [15], Jha et al. (2014) in Saudi Arabia [20], and Ezzeldin et al. (2022) in Egypt [21]. These studies have demonstrated that combining RS and GIS techniques with MCDA is the most effective approach in the RWH process.

The aim of this study is to implement a robust methodology for generating a suitability zone map for rainwater harvesting (RWH) in the study area. This approach will be tailored to the existing environmental, economic, and social conditions of the region. The study also seeks to determine optimal locations for constructing different RWH structures, including medium dams, check dams, and farm ponds. The utilization of Remote Sensing (RS), Geographic Information System (GIS), and Multi-Criteria Decision Analysis (MCDA) techniques will be employed to tackle water scarcity issues specifically in Erbil Province, Iraq. In this study, updated data, such as satellite images and climate data, were utilized to prepare thematic layers. This study makes use of the Topographic Wetness Index (TWI) criteria, which Beven and Kirkby developed in 1979 [22] and are crucial for determining the best locations to hold water [23]. TWI criteria can enhance the accuracy of the suitability zone map for RWH. On the other hand, it is important to highlight that this study covers the entirety of Erbil Province, in contrast to previous studies that typically focused on specific parts of the province. The proposed methodology is illustrated through a case study. The research outputs demonstrate the potential of RWH to achieve various objectives, including mitigating the impacts of water scarcity, reviving surface water resources, increasing groundwater levels, and fostering agricultural development. This study provides an initial overview of the potential of RWH in the study area, and the results are intended to assist decision-makers and local officials in future planning and water resource management.

2. Materials and Methods

2.1. Study Area

Erbil Province, situated in the northeast of Iraq, covers an area of 14,837 km² and is positioned between 44° and 45° E longitude and 36° and 37° N latitude (Figure 1a). Its internal borders are the governorates of Sulaymaniyah to the east, Kirkuk to the southeast, Salahaddin to the southwest, Ninawa (Mosul) to the west, and Dohuk to the northwest. It shares international borders with Turkey to the north and Iran to the northeast. Erbil Province contains ten administrative districts, as shown in Figure 1b. The population is about 2.25 million, which is mostly Kurdish [24]. The region experiences a Mediterranean weather pattern characterized by arid and semi-arid conditions, with high temperatures during the summer months and cool, damp winters. According to data obtained from the Ministry of Agriculture and Water Resources (KRG) for the period 2014–2023, the average annual rainfall ranges from 250 mm to 1400 mm. Topographically, the northern parts are high and comprise the most famous chains of mountains called Zagros, with the peak of Hasarost being the highest peak in the region with a height of 3607 m above sea level [25]. These heights decrease gradually towards the central portions until they reach the plains in the southern parts, which make up most of the agricultural land in the study area. The soil in the northern areas is shallow to medium chestnut soil that has been created from the original rocks; it has a low potential for agriculture, but it is rich in the natural rangeland, whereas the plain areas consist of dark brown and black soils and are favorable for agriculture due to their great depth, good texture, and high content of organic matter [26]. Wheat and barley are the main crops in the winter season, depending on the rainfall, while many other agricultural crops grow depending on underground water resources during the summer [27].

2.2. Dataset

To prepare the criteria that are fundamental for generating the RWH suitability zone map, the required data were collected from various resources, including the following:

• Remote sensing data, including a Landsat 8 Operational Land Imager (OLI) satellite image dated 14 April 2023, and a Digital Elevation Model (DEM) sourced from the Shuttle Radar Topography Mission (SRTM) with a 30 m resolution, were accessed on 15 June 2023, from https://earthexplorer.usgs.gov [28].
• Rainfall data for 10 years (2014–2023) were obtained from 23 meteorological stations belonging to the Ministry of Agriculture and Water Resources (KRG).
• Soil data were retrieved from the digital soil map of the world (DSMW) developed by the Food and Agriculture Organization (FAO) and the United Nations Educational, Scientific, and Cultural Organization (UNESCO) (FAO-UNESCO 2008) [29], accessed on 10 July 2023.

2.3. Criteria Selection and Preparation

The FAO guidelines, earlier studies, and the availability of data on the study area all played a role in the selection of various criteria used to create the suitability zone map. The criteria that were selected and generated are listed below:

2.3.1. Rainfall

Rainfall is a critical factor in the RWH process. The FAO recommends the adoption of RWH technology (RWHT) in regions with annual rainfall ranging from 100 mm to 1000 mm. In areas receiving less than 100 mm of yearly rainfall, there is no incentive to implement RWH [30]. The rainfall intensity in the study area was assessed using daily rainfall data collected over a 10-year period (2014–2023) from twenty-three meteorological stations within the study area. The geostatistical method known as Inverse Distance Weighted (IDW), commonly used for interpolating rainfall variables, was employed to estimate the spatial distribution of rainfall data. By considering the geographical direction and distance of existing data points, rainfall values for unobserved locations were approximated.

2.3.2. Soil Texture

The texture of the soil plays a significant role in the RWH process [31]. Soils with a higher water-holding capacity are more suitable for RWH [32]. The soil map was obtained from the FAO/UNESCO DSMW (2008) [29], with a scale of 1:5,000,000. However, the Natural Resources Conservation Service (NRCS) of the United States Department of Agriculture (USDA) defined soil characteristics related to water retention and infiltration in 2007 [33]. Based on these characteristics, the different types of soil texture have been put into Hydrological Soil Groups (HSGs).

2.3.3. Land Use and Land Cover (LU/LC)

The LU/LC map was created using Landsat 8 OLI satellite imagery from the year 2023. The classification process involved a supervised image classification approach employing the maximum likelihood algorithm. The image was categorized into six classes: forest, rangeland, barren land, urban or built-up area, agricultural land, and water bodies. To evaluate the accuracy of the classification, the study employed a confusion matrix. Additionally, cross-referencing was conducted using a Google Earth map, random points, and actual ground points.

2.3.4. Runoff Depth

Estimating runoff depth is a crucial step in identifying suitable locations for RWH [6]. Several methods are available for estimating runoff depth, with two commonly used ones being the Rational Method and the Soil Conservation Service Curve Number (SCS-CN). In both methods, it is important to consider the characteristics of the watershed, such as land use, soil type, and antecedent moisture conditions. Additionally, rainfall data—whether observed or estimated—is a key input for runoff estimation. In this study, the SCS-CN approach, originally established by the USDA-SCS in 1972 (now known as the Natural Resources Conservation Service—NRCS) [34,35], was utilized to evaluate potential runoff in the study area. This method calculates direct runoff from rainfall events in a watershed or catchment area for each pixel [36,37]. The equation used to estimate runoff in the study area is as follows:

$$Q=\frac{\left(P-Ia\right)²}{\left(P-Ia\right)+S}\hspace{1em}\hspace{1em}for p>Ia$$

(1)

where Q is the depth of the runoff (mm), P is the amount of precipitation (mm), S is the maximum amount of water that could be kept after the runoff starts (mm), and Ia is the initial abstraction (mm) that includes all the water that was lost before the runoff started through infiltration, evaporation, and water interception by vegetation. Ia is highly variable but is generally correlated with soil and cover parameters. Through studies of many small agricultural watersheds, researchers approximated Ia as 0.2 [38]. The value of 0.2 for Ia has also been mentioned in previous studies on Erbil province, including Hameed (2013) [26], Babir and Ali (2016) [39], Hameed (2017) [40], and Majeed (2023) [41]. This means that the amount of precipitation (P) is greater than the initial abstraction (Ia) and the amount of water after the runoff starts (S) is suitable for harvesting in the study area. Therefore, in this study, the same value is employed, and the equation can be expressed as:

$$Q=\frac{{(P-0.2S)}^2}{\left(P+0.8S\right)}$$

(2)

The potential maximum retention (S) can be computed by using the estimated CN [42] as follows:

$$S=\frac{25400}{\left(CN\right)}-254$$

(3)

CN is reflecting the surface runoff’s response to a given rainfall, which ranges from 0 to 100. High CN values show that much of the rainfall is transformed into surface runoff and vice versa [32,43]. The CN values were taken from the USDA-NRCS National Engineering Handbook based on the integration of land cover classes and HSGs with respect to hydrological conditions [32,44], reflecting the surface runoff’s response to a given rainfall ranges from 0 to 100. High CN values show that much of the rainfall is transformed into surface runoff and vice versa [32,43].

2.3.5. Drainage Density (DD)

DD is a fundamental geomorphological parameter used to assess the density and complexity of a river or stream network within a specific geographic area. It is a critical morphometric parameter that plays a vital role in the RWH process [10,45]. Areas with higher DD tend to yield greater runoff for harvesting [46]. In this study, the DEM data were employed in ArcGIS 10.8 software to create a flow direction map, which illustrates the movement of water across the landscape, and a flow accumulation map, which highlights areas with a high potential for water flow. By applying a threshold value >1000 to this map, the stream network was identified. DD was subsequently calculated by dividing the total stream length by the size of the study area [47,48,49].

2.3.6. Slope

The slope of the land has a big effect on hydrological aspects, like runoff generation, recharge facilitation, sedimentation modulation, and water flow velocity. In this context, the slope map plays a crucial role in selecting suitable locations for RWH [32,48,50]. The slope ranges of the study area were calculated using SRTM-DEM data projected in the Universal Transverse Mercator (UTM), Zone 38N, and the WGS84 horizontal datum. The resulting slope map was categorized into five slope percentage values following FAO guidelines.

2.3.7. Topographic Wetness Index (TWI)

The TWI measures the balance between water accumulation and drainage on sloped land by looking at the connection between upslope areas and the local slopes [51]. It provides insights into how topography influences water movement, runoff generation, and accumulation [23]. The TWI map was used to evaluate the impact of topography on hydrological processes within the study area. The TWI map was generated using DEM data, which are a representation of the Earth’s surface in a gridded format, and the values were calculated using spatial analyst tools in ArcMap 10.8.2 software. The values were categorized into five groups, with higher TWI values indicating a larger drainage area, implying greater water availability for harvesting, and lower values suggesting the opposite [22,52]. The formula for calculating the TWI is often expressed as:

$$TWI=\mathrm{l}\mathrm{o}\mathrm{g}(\frac{\alpha }{\mathrm{t}\mathrm{a}\mathrm{n}\left(\beta \right)})$$

(4)

where $\alpha$ is the upslope contributing area represents the total area draining through a certain point on the landscape and $\mathrm{t}\mathrm{a}\mathrm{n}\left(\beta \right)$ is the tangent of the slope [23].

2.4. Criteria Prioritization

To generate a suitable zone map for RWH in the study area, seven criteria were identified as thematic layers. Professionals and experts from the area were asked to rate and weight the importance of each criterion using Saaty’s basic scale [53], which goes from 1 to 9, with 9 being the most important and 1 being the least important (

Table 1. The fundamental scale for evaluating the relative importance of criteria in AHP.

The Intensity of Importance Scale	Relative Importance Intensities	Description
1	Equally important	Two activities make an equal contribution to the objective.
3	Moderate important	One activity is slightly preferred over another.
5	Strong important	One activity is greatly preferred over another.
7	Very strong important	One activity is very strongly preferred over another, resulting in its dominance in practice.
9	Extremely Important	The evidence supporting one activity as compared to another is of the highest level of confirmation.
2, 4, 6, 8	Values between two adjacent judgments	Additional subdivision or compromise when required.

Note: Source: [53].

). Experts were selected for the interviews due to their comprehensive understanding of the research area and prior experience in recognizing the significance of RWH. All selected experts have resided in Erbil Province, Iraq, for more than 30 years. Additionally, a literature review was conducted to validate the selection criteria and rankings for accuracy, such as Berhanu and Bisrat 2018 [54], Adham et al. (2022) [32], Alene et al. (2022) [55], Gebremedhn et al. 2023 [56], and Noori et al. 2023 [6].

2.5. Multi-Criteria Decision Analysis (MCDA)

MCDA involves the selection of criteria and decision options [57]. It is a method used to assess the importance of various parameters within a project. The significance of these criteria in determining suitable sites can vary based on their respective weights, and the outcome of the decision-making process heavily depends on these criterion weights. Therefore, ensuring the objectivity of criterion weights is a critical step in the MCDA process [58].

In order to determine the weights of the thematic layers in this study, Saaty’s AHP method from 1977 [19] served as the MCDA method. AHP is a popular method that involves setting up a hierarchy of selection criteria and comparing two items in a matrix pairwise to find out and normalize the weights of each element. The Pairwise Comparison Matrix (PCM) serves as a tool for assessing the relative importance of each criterion in comparison to all others. In this step, the weights of each criterion in the rows were compared to themselves in the columns, where the relative importance is equal and the value “1” is treated as a constant. Subsequently, these weights were compared to all other criteria, where the relative importance varies based on each criterion’s weight. To normalize the PCM, the sum of each column was calculated, and then each cell in the column was divided by its column total to obtain the eigenvectors of the matrix. The average weights for each criterion were found by adding up all the values in the PCM rows and then dividing that number by the number of criteria. In the AHP approach, the consistency ratio (CR) of expert judgments plays a pivotal role in making sound decisions. It is imperative for the decision-maker to closely monitor the consistency of these judgments because inconsistent assessments can potentially lead to erroneous results. Saaty 1984 [53] recommended, through CR calculations, that CR values should ideally be 0.1 or less, indicating a satisfactory reciprocal matrix. Conversely, if the CR is greater than 0.1, it means that the PCM is not consistent and needs to be re-evaluated. In this study, the AHP-CR was computed to assess the consistency of the weights assigned to different layers. This calculation was performed using the following equations:

$$CR=\frac{CI}{RI}$$

(5)

CR is a consistency ratio; CI is a consistency index; it is a factor that measures the consistency of diagonal comparison matrices calculated by Equation (6). RI is a random index; it is a standard value determined by Saaty (1984) [53], and the RI values are variable depending on the number of parameters listed in

Table 2. Random Index values.

Order	1	2	3	4	5	6	7	8	9	10
RI	0	0	0.52	0.89	1.11	1.25	1.35	1.40	1.45	1.49

Note: source: [53].

.

$$CI=\frac{\lambda max-n}{n-1}$$

(6)

CI is the consistency index. $\lambda max$ is the largest eigenvalue of the comparison matrix calculated as the sum of products between each priority value element and the total of columns of the reciprocal matrix. $\lambda max$ was computed in Equation (7). $n$ is the number of criteria.

AX = λ max

(7)

A is the comparison matrix of size for the criteria. X is the Eigenvector of size 1.

2.6. Generating the RWH Suitability Zone Map

To generate the RWH-suitable zone map, a standardization process was applied to the initial criteria layers, which had varying units. This standardization was crucial for enabling weighted overlay analysis. Initially, the criteria layers were converted from vector to raster format using the reclassify tool in ArcMap 10.8. Subsequently, each raster layer was assigned internal values on a scale of 1 to 5, indicating their relative importance. A value of 1 represented the lowest importance, while 5 denoted the highest. These assignments were determined based on expert recommendations, previous studies such as Adham et al. (2016) [59], Tahera et al. (2022) [60], and Surve et al. (2022) [61], and the characteristics of each criterion. The ‘weighted overlay’ tool in ArcMap, commonly used for multi-criteria overlay analysis in site selection and suitability modeling tasks, was then applied. This tool integrated all the raster layers, facilitating the identification of suitable zones within the study area. The resulting map was categorized into four classes: ‘very high suitability’, ‘high suitability’, ‘moderate suitability’, and ‘low suitability’.

2.7. Suitable Sites Selection for Different RWH Structures

Effective RWH relies on prudent land management practices and the construction of appropriate structures [61]. In this study, the resultant suitability zone map was used to determine potential RWH locations within the study area. Factors will have different degrees of impact on finding the best sites for each RWH structure. In this regard, multiple layers were added to the ArcGIS environment based on the characteristics of the RWH structures. As per previous studies, each RWH structure necessitates specific characteristics. For instance, Mahmood et al. (2020) [44] recommended that sites suitable for farm ponds should have a slope of less than 5% and belong to the first or second stream order. The LU/LC should be predominantly agricultural with low infiltration, and runoff depth should be moderate to high. Check dams, following Singh et al. (2009) [62], are best situated in barren lands and riverbed areas with a slope of less than 15%, third-order stream drainage, low infiltration, and moderate runoff. Furthermore, Ibrahim et al. (2019) [63] emphasized that areas with complex topography and extensive drainage networks are ideal for constructing large or medium embankments, relying on a substantial amount of precipitation and runoff as primary prerequisites. Considering these guidelines and the characteristics of each RWH structure, three types of RWH structures were proposed: medium dams, check dams, and farm ponds.

2.8. Validation

The validation process is an important part of making sure that the suitability map made for building different RWH structures is correct and to see how well the methods and techniques used worked. Assessing the RWH suitability zone map is mostly based on a mix of field studies, existing data, and analysis using geospatial technology [15]. The RWH suitability zone map was compared to the study area’s existing functional RWH structures to see how well the method worked. The existing RWH structure points were identified and gathered from the Ministry of Agriculture and Water Resources (KRG), where they represent eleven farm ponds, eleven check dams, and one medium dam. All coordinate points were overlaid on the RWH suitability zone map using the ArcMap platform to determine their respective suitability classes.

3. Results and Discussion

3.1. Thematic Layers

In this study, seven fundamental criteria have been identified that are essential to producing a suitable zone map for RWH. These criteria were meticulously selected based on considerations of data availability and adherence to the guidelines established by the FAO. Notably, these criteria have also been previously utilized in studies with similar objectives conducted in diverse geographical regions. For instance, Ezzeldin et al. (2022) [21] applied them in Egypt, Adham et al. (2022) [32] in Palestine, Wu et al. (2018) [57] in Guatemala, and Yegizaw et al. (2022) [64] in Ethiopia. To create these thematic layers, relevant information was extracted by thoroughly analyzing both remote sensing (RS) data and data collected from field surveys [65]. The criteria employed in this study are outlined below:

3.1.1. Rainfall

The generated rainfall map reveals that the average annual rainfall falls within the range of 250–1400 mm for the period 2014–2023. Following Modak and Das (2022) [15], the rainfall values were categorized into five groups: very low (250–400 mm), low (400–600 mm), moderate (600–800 mm), high (800–1200 mm), and very high (1200–1400 mm). These categories occupy specific areas within the study region, covering 4201.37 km² (28.32%), 3193.66 km² (21.52%), 2877.27 km² (19.39%), 3810.31 km² (25.68%), and 754.39 km² (5.08%), respectively (

Table 3. Rainfall categories with their ranges and areas.

Rainfall Categories	Average Annual Rainfall (mm)	Area (km2)	Area (%)
Very low	250–400	4201.37	28.32
Low	400–600	3193.66	21.52
Moderate	600–800	2877.27	19.39
High	800–1200	3810.31	25.68
Very High	1200–1400	754.39	5.08
Total		14,837	100

). The rainfall distribution map shows that the lowest average rainfall values are in the southern parts of the study area, gradually increasing towards the north to reach the highest level in the mountainous areas (Figure 2).

3.1.2. Soil Properties

The study region displays diverse soil texture categories, encompassing silty clay, loam, silty loam, sandy loam, sandy clay loam, and clay loam. These groups take up different amounts of the total area, with 2597.69 km² (17.51%), 2594.19 km² (17.48%), 3337.28 km² (22.49%), 1649.88 km² (11.12%), 1438.50 km² (9.70%), and 3219.46 km² (21.70%), respectively (

Table 4. Soil texture classes with their areas.

Class No	Texture	Area (km2)	Area (%)
1	Silty Clay	2597.69	17.51
2	Loam	2594.19	17.48
3	Silty Loam	3337.28	22.49
4	Sandy Loam	1649.88	11.12
5	Sandy Clay Loam	1438.50	9.70
6	Clay Loam	3219.46	21.70
Total		14,837	100

). Notably, silty loam is the most common type of soil texture in the study area. It is mostly found in the central region but can also be found in some parts of the southern region (Figure 3). In addition, the soil texture map shows three Hydrologic Soil Groups (HSGs) based on their water infiltration rates: B, C, and D (Figure 4). These groups occupy distinct regions, covering 4285.84 km² (28.9%), 7398.93 km² (49.9%), and 3152.23 km² (21.2%), respectively (

Table 5. HSG classes with their areas.

Class No	HSG	Area (km2)	Area (%)
1	B	4285.84	28.9
2	C	7398.93	49.9
3	D	3152.23	21.2
Total		14,837	100

).

3.1.3. LU/LC

The supervised image classification process yielded six distinct LU/LC classes. The analysis of the confusion matrix revealed an overall accuracy of 91.6% and a high overall kappa statistic of 90.04, signifying a high level of precision in the classification process. The resulting map illustrates that the predominant land in the study area is agricultural land, covering an area of 4011.22 km² (27.04%), primarily concentrated in the southern regions. Following closely is rangeland, extending over 3880.98 km² (26.16%), with a predominant presence in the northern areas (Figure 5). Interestingly, barren land, occupying 3109.89 km² (20.96%), emerges as the most suitable category for RWH structures. In contrast, water bodies account for only 41.26 km² (0.28%) of the study area. Urban/built-up areas and forests comprise 856.19 km² (5.77%) and 2937.46 km² (19.80%) of the total area, respectively (

Table 6. LU/LC classes and their areas.

LU/LC Class	Area (km2)	Area (%)
Forest	2937.46	19.80
Rangeland	3880.98	26.16
Barren land	3109.89	20.96
Urban/Built-Up	856.19	5.77
Agricultural land	4011.22	27.04
Water Bodies	41.26	0.28
Total	14,837	100

).

3.1.4. Runoff Estimation

The runoff potential map was generated using CN grid values as an empirical parameter. It was possible to estimate the CN values by putting together information from LU/LC, HSG, and precipitation [3]. The results reveal a CN value range of 55 to 100 across the entire study area (Figure 6). The highest CN value observed in the study area is 100, primarily in water bodies. On the other hand, forests in good condition with HSG (B) have the lowest CN value of 55. This is because the forest cover and mostly loam soil texture make it easy for water to pass through [58]. Detailed CN values for each hydrologic soil group and corresponding land use class can be found in

Table 7. CNs according to LU/LC classes and HSGs.

LU/LC Class	HSGs CN
	B	C	D
Forest	55	73	60
Rangeland	63	86	93
Barren land	67	85	88
Urban/Built-Up	62	91	95
Agricultural land	69	83	87
Water Bodies	100	100	100

. The resulting runoff potential map shows four suitability categories: ‘low’ (150–400 mm), ‘moderate’ (400–800 mm), ‘high’ (800–1000 mm), and ‘very high’ (>1000 mm). This classification aligns with earlier studies, such as those by Rajasekhar et al. (2019) [65] and Saha et al. (2021) [66]. In terms of area coverage, the ‘very high’ runoff potential class spans 2752.50 km² (18.6%) of the study area, while the ‘high’ class covers 6781.20 km² (45.7%). The ‘moderate’ and ‘low’ runoff classes occupy 1733.60 km² (11.6%) and 3569.62 km² (24.1%) of the total area, respectively (

Table 8. Runoff potential classes.

Class No	Runoff (mm)	Area (km2)	Area (%)
Very high	>1000 mm	2752.50	18.6
High	800–1000	6781.20	45.7
Moderate	400–800	1733.60	11.6
Low	150–400	3569.62	24.1
Total		14,837	100

). The runoff potential map illustrates that the northern portion, which is characterized by mountains with steep slopes and higher rainfall, is a significant contributor to runoff generation (Figure 7). In contrast, the central study area, mainly composed of flat agricultural and barren lands, predominantly falls under the ‘high runoff potential’ class. The southern region, characterized by flat agricultural lands, is classified under the moderate runoff coefficient’ category. To provide a visual representation, Figure 7 displays the spatial distribution of runoff potential classes across the study area.

3.1.5. Drainage Density

The DD map results demonstrate a range of values from 0.21 to 0.94. Following Al-Ghobari and Dewidar 2021 [43] and Setiawan and Nandini 2022 [46], the DD values were put into five groups based on their suitability for RWH: “Very High” (0.62–0.94), “High” (0.47–0.61), “Moderate” (0.34–0.46), “Low” (0.22–0.33), and “Very Low” (0–0.21) (

Table 9. Drainage density classes with their values and areas.

Class	Value	Area (km2)	Area (%)
Very High	0.62–0.94	1391.80	9.38
High	0.47–0.61	2552.12	17.20
Moderate	0.34–0.46	4074.23	27.46
Low	0.22–0.33	4729.72	31.88
Very Low	0–0.21	2089.13	14.08
Total		14,837	100

). The predominant DD category within the study area is ‘low’, covering a substantial area of 4729.72 km² (31.88%). It is widely distributed across all parts of the study area, with a notable concentration in the northern and eastern central regions. Conversely, the ‘very high’ DD category occupies a relatively smaller area, specifically 1391.80 km² (9.38%). It is primarily situated in the central and southern portions of the study area, with a narrow strip in the northern portion (Figure 8). Areas categorized as ‘high’, ‘moderate’, and ‘very low’ DD cover 2552.12 km² (17.20%), 4074.23 km² (27.46%), and 2089.13 km² (14.08%), respectively. Sites with ‘very high’ DD values are particularly suitable for RWH due to their lower infiltration rates and higher surface flow velocities, facilitating more efficient rainwater capture. This observation has been supported by various studies, including those by Balkhair and Rahman (2021) [45], Ahmed et al. (2023) [49], Dragievici et al. (2019) [67], Khudhair et al. (2020) [68], and Alene et al. (2022) [55].

3.1.6. Slope

The slope map shows that slope degrees within the study area vary from 0 to 80%. The slope degrees are classified into five categories based on percentage values: (0–5) nearly level, (5–10) gentle slopes, (10–20) moderate slopes, (20–40) high slopes, and (40–80) very high slopes (

Table 10. Slope classes with their degrees and areas.

Slope Classes	Slope Degree (%)	Area (km2)	Area (%)
Nearly Level	0–5	6952.37	46.85
Gentle	5–10	2342.32	15.78
Moderate	10–20	2759.55	18.59
High	20–40	2555.39	17.25
Very High	40–80	227.37	1.53
Total		14,837	100

). The central and southern parts of the study area predominantly fall within the (0–5%) and (5–10%) slope categories, covering a combined area of 6952.37 km² (46.85%) and 2342.32 (15.78%), respectively. These two categories are generally considered more suitable for RWH, in accordance with findings from Modak and Das (2022) [15] and Adham et al. (2022) [32]. The ‘moderate slope’ category spans 2759.55 km² (18.59%), while the ‘high slope’ and ‘very high slope’ categories encompass 2555.39 km² (17.25%) and 227.37 km² (1.53%), respectively. These steeper slopes are primarily concentrated in the northern portions of the study area (Figure 9).

3.1.7. TWI

The TWI map displays a range of values spanning from 1 to 25, which have been categorized into five distinct classes: ‘Very High’ (20–25), ‘High’ (15–20), ‘Moderate’ (10–15), ‘Low’ (5–10), and ‘Very Low’ (1–5). Among these classes, the ‘Very Low’ TWI class claims the largest area, covering 4266.42 km² (28.75%), while the ‘Very High’ TWI class has the most limited coverage, accounting for just 321.66 km² (2.17%). The ‘Low’ TWI class encompasses 6034.47 km² (40.67%), ‘Moderate’ TWI extends over 2787.57 km² (18.79%), and ‘High’ TWI covers an area of 1426.88 km² (9.62%). It is worth emphasizing that higher TWI values signify a heightened potential for RWH, as they indicate an increased capacity for water accumulation, as elucidated by Berhanu and Bisrat (2018) [54].

Table 11. TWI classes with their values and areas.

TWI Classes	Value	Area (km2)	Area (%)
Excessively high	20–25	321.66	2.17
High	15–20	1426.88	9.62
Moderate	10–15	2787.57	18.79
Low	5–10	6034.47	40.67
Very low	1–5	4266.42	28.75
Total		14,837	100%

and Figure 10 show a comprehensive view of the distribution of TWI values in the study area.

3.2. Criteria Prioritization and the MCDA Process

The prioritization of criteria revealed that the ‘Runoff’ criterion obtained the highest rank with a value of 9, showing its significant influence on the selection of suitable sites for RWH. Following this, ‘Rainfall’ received a value of 7. The criteria ‘LU/LC’, ‘Slope’, and ‘Soil Texture’ were equally important, each assigned a value of 5. ‘DD’ and ‘TWI’ were considered the lowest priority criteria, both having a value of 3 (

Table 12. Criteria used to generate a suitable zone map for RWH and their properties.

Criteria	Criteria Ranking	Unit	Class	Suitability Ranges	Class Value	Weight (%)
Runoff	9	Mm	>1000	Very high	5	38
			800–1000	High	4
			400–800	Moderate	3
			150–400	Low	2
Rainfall	7	mm	250–400	Very High	5	20
			400–600	High	4
			600–800	Moderate	3
			800–1200	Low	2
			1200–1400	Very Low	1
LULC	5	Class	Barren Land	Very High	5	10
			Grassland	High	4
			Cultivated Land	Moderate	3
			Forest	Low	2
			Urban/Built-Up Area	Not Suitable	0
			Snow	Not Suitable	0
			Shadow	Not Suitable	0
			Water Bodies	Not Suitable	0
Slope	5	Degree	0–5	Very high	5	10
			5–10	High	4
			10–20	Moderate	2
			20–40	Not suitable	0
			40–80	Not Suitable	0
Soil Texture	5	Type	Clay Loam	Very high	5	10
			Silty Clay	High	4
			Sandy Clay Loam	Moderate	3
			Silty Loam	Moderate	3
			Sandy Loam	Moderate	3
			Loam	Low	2
Drainage Density	3	Value	Very High	0.62–0.94	5	6
			High	0.47–0.61	4
			Moderate	0.34–0.46	3
			Low	0.22–0.33	2
			Very Low	0–0.21	2
TWI	3	Value	20–25	Very High	5	6
			15–20	High	4
			10–15	Moderate	3
			5–10	Low	2
			1–5	Very Low	1

). To determine the relative importance of each criterion with respect to the others, the PCM, an integral part of the AHP, was utilized. The normalized values obtained from the PCM were used to calculate numerical weights and corresponding percentages for each criterion. These calculations revealed that the ‘Runoff’ criterion carries the highest weight at 38%, while ‘DD’ and ‘TWI’ have the lowest weights at 6%. ‘LU/LC’, ‘Slope’, and ‘Soil Properties’ all share an equal weight of 10%, and ‘Rainfall’ is assigned a weight of 20% (Table 12). Subsequently, after assigning weights, the CR was computed to assess the relative priority of each criterion, following the method employed by Modak and Das (2022) [15]. The results indicated a principal eigenvalue (λ Max) of 7.0135, a CI of 0.0023, and a CR of 0.0017. A CR below 0.1 signifies that the comparison matrix is consistent and the expert judgments are considered acceptable.

3.3. RWH Suitable Zone

The RWH-suitable zone map delineates four suitability classes: ‘Very High’, ‘High’, ‘Moderate’, and ‘Low’. The results indicate that the ‘Moderate’ suitability class dominates, covering 5295.65 km² (35.69%) of the study area. This class is distributed widely, with a significant presence in the southern part, primarily comprising agricultural land. The ‘High’ suitability class, the second largest, encompasses 4968.55 km² (33.49%), mainly situated in the central area from east to west. These areas include settlement zones and barren lands. The ‘Very High’ suitability class occupies 1583.25 km² (10.67%), primarily in the northern part. In contrast, 2989.66 km² (20.15%) of the total study area is less suitable for RWH, concentrated in the uppermost northern regions covered by forest and the lower southern regions, which mostly comprise agricultural lands (Figure 11 and

Table 13. RWH suitability classes and their areas.

S.n	Suitability Classes	Area (km2)	Area (%)
1	Very High Suitable	1583.25	10.67
2	High Suitable	4968.55	33.49
3	Moderate Suitable	5295.65	35.69
4	Low Suitable	2989.66	20.15
Total		14,837	100

). The analysis highlights that the northern and central areas are particularly ideal for RWH. This RWH-suitable zone map is vital for planning RWH structures and artificial recharge strategies. Harvested rainwater can alleviate pressure on existing water sources by recharging groundwater and serving various purposes, like irrigation and livestock [69]. Additionally, it has the potential to mitigate drought, recharge wells and springs, and reduce groundwater salinity levels [3]. The scientific method used to make this RWH-suitable zone map is also used in earlier research by Jha et al. (2014) [20], Singh et al. (2009) [62], and Tiwari et al. (2018) [70]. Therefore, it serves as a crucial foundation for the development of an efficient and effective water management strategy in the study area.

3.4. Appropriate Sites for RWH Structures

The features of the study area and the building requirements for each type of RWH structure helped find 36 good places for three different types of RWH structures (

Table 14. Appropriate locations for RWH structures.

S.No	Structure Type	Latitude	Longitude
1	Medium dam	36.6373	44.4924
2	Medium dam	36.8123	44.5312
3	Medium dam	36.5872	44.7966
4	Check Dam	36.0178	44.4388
5	Check Dam	35.8459	44.2902
6	Check Dam	36.4203	43.9298
7	Check Dam	36.5783	44.1114
8	Check Dam	36.0895	43.5968
9	Check Dam	36.0797	43.7136
10	Farm Pond	35.6295	43.5572
11	Farm Pond	35.6177	43.6601
12	Farm Pond	35.6884	43.3768
13	Farm Pond	35.8599	43.4097
14	Farm Pond	36.0570	43.5075
15	Farm Pond	36.1631	43.8435
16	Farm Pond	36.1350	43.9112
17	Farm Pond	36.2995	43.7578
18	Farm Pond	36.0185	43.6387
19	Farm Pond	35.8341	44.0945
20	Farm Pond	36.0403	44.3263
21	Farm Pond	36.0454	44.0596
22	Farm Pond	36.0441	44.1104
23	Farm Pond	36.0122	44.1259
24	Farm Pond	36.0038	44.7385
25	Farm Pond	35.9517	44.7632
26	Farm Pond	36.0223	44.6762
27	Farm Pond	36.5956	44.2806
28	Farm Pond	36.5556	44.2965
29	Farm Pond	35.9817	43.9354
30	Farm Pond	35.9711	43.9674
31	Farm Pond	35.9659	43.5463
32	Farm Pond	35.9640	43.6035
33	Farm Pond	35.5290	43.5719
34	Farm Pond	36.2770	43.8913
35	Farm Pond	36.2519	43.7665
36	Farm Pond	36.2472	43.8114

). The results indicate that most parts of the study area are suitable for check dam and farm pond construction. Accordingly, 27 locations for farm ponds and six locations for check dams were proposed, primarily situated in the central and southern portions of the study area from east to west, as shown in Figure 12. All farm ponds are situated on agricultural land with nearly level slopes, adjacent to 1st and 2nd-order streams. Out of the 27 farm ponds, only two sites are in the northern parts of the study area. The potential sites for check dams are located close to agricultural lands and can be useful for supplementing irrigation during the dry season and controlling the speed of flow during stormwater events. Overall, all the check dams are situated in riverbeds with 3rd and 4th stream orders, and the land cover is predominantly barren land. Generally, farm ponds and check dams are both located in highly and moderately suitable zones. The results indicate that the northern parts of the study area are ideal for constructing medium-sized dams due to the high amount of precipitation and runoff potential, in addition to the land cover types and topography. This area exhibits the most complex topography, with a deep drainage network. In this context, three suitable locations have been suggested for constructing medium dams, which fall within the very highly suitable class and are positioned across 4th and 5th-order streams to provide better water access for the dams. These three suggested dams offer a well-distributed solution within the study area and can be utilized for an extended period, providing significant benefits to northern settlements. The outcomes of this research will provide essential guidance for decision-makers responsible for planning reservoir construction within the study area. Equipped with this information, they can make well-informed decisions and ensure the accuracy of their project outcomes. In the end, the research results have effects on many areas, such as better management of water resources, lessening the effects of water shortages, supporting agricultural progress, making sustainable development easier, and improving people’s lives in Erbil Province, Iraq.

3.5. Validation

The validation results indicated that all existing RWH structures were classified as ‘successful’ as they conformed to the criteria associated with each structure type. The agreement statements between existing RWH structure sites and the suitability map are presented in

Table 15. Agreement between suitable zone map and the existing RWH structures.

S.No	Latitude	Longitude	Structure Type	Agreement
1	35.87378	43.82843	Check Dam	Agree
2	35.89254	43.84918	Farm Pond	Agree
3	36.14142	44.34366	Check Dam	Agree
4	36.1028	44.21511	Check Dam	Agree
5	36.10235	44.59953	Check Dam	Agree
6	35.90128	44.75567	Check Dam	Agree
7	36.17284	44.38198	Check Dam	Agree
8	35.98385	44.58115	Check Dam	Agree
9	36.2823	44.1516	Farm Pond	Agree
10	36.30324	44.13427	Farm Pond	Agree
11	35.88617	44.76812	Farm Pond	Agree
12	36.10897	44.31211	Farm Pond	Agree
13	36.01644	44.56986	Farm Pond	Agree
14	35.9228	44.87123	Farm Pond	Agree
15	36.16533	44.58236	Check Dam	Agree
16	36.95956	44.34863	Check Dam	Agree
17	36.63323	44.19097	Farm Pond	Agree
18	36.62282	44.19402	Check Dam	Agree
19	36.60145	44.13792	Farm Pond	Agree
20	36.27569	44.28037	Check Dam	Agree
21	36.16454	44.58246	Medium Dam	Agree
22	35.53305	43.41075	Farm Pond	Agree
23	36.52808	43.95224	Farm Pond	Agree

. The findings reveal that 13.04% of existing RWH structures were in areas classified as very highly suitable, 60.87% in highly suitable areas, 21.74% in moderately suitable areas, and only 4.35% in areas classified as lowly suitable. Based on the overlay results, among the 11 check dams, 10 were situated in highly suitable areas, with only one in a moderately suitable area. However, of the 11 farm ponds, five were in highly suitable areas, five in moderately suitable areas, and only one in a lowly suitable area. Regarding the medium dam, the single-identified dam is situated in a very highly suitable area. It is noteworthy that, from the overlaid map, most existing structures fell within the highly suitable zone (Figure 13). The validation results suggest that the produced maps provide a reliable representation of the spatial distribution of suitable areas, aligning with existing rainwater harvesting practices in the study area. Similar validation approaches have been employed in other studies, including Nyirenda et al. (2021) [30], Rajasekhar et al. (2019) [65], Alene et al. (2022) [55], Kumar and Jhariya (2017) [71], and Haile and Suryabhagavan (2019) [72]. In conclusion, these validation results underscore the value of the database and methodology used for generating the RWH suitability zone map in the context of effective water resource management within the study area.

4. Conclusions

Water scarcity presents a significant challenge to communities and ecosystems in arid and semi-arid regions, impacting agricultural productivity and overall sustainability. Over the past two decades, Erbil Province in Iraq has faced issues related to drought waves and water scarcity, particularly in its southern areas. This study underscores the importance of implementing rainwater harvesting (RWH) as a promising and sustainable solution to address the challenges encountered in the study area.

A Geographic Information System (GIS) and Multiple Criteria Decision Analysis (MCDA) were used to create a suitability zone map and look for possible locations for various RWH structures in this study. The study looked at a lot of different factors, such as soil texture with Hydrologic Soil Group (HSG) characteristics, slope, drainage density, topographic wetness index (TWI), and rainfall and runoff potential. The results revealed that only 1583.25 km² (10.67%) of the total area is highly suitable for RWH, mainly in the northern parts of the study area. However, the ‘High’ suitability class covers 4968.55 km² (33.49%), primarily in the central area from east to west. The ‘Moderate’ suitability, dominating the study area at 5295.65 km² (35.69%), is distributed widely, with a significant presence in the central and southern parts. In contrast, 2989.66 km² (20.15%) of the total study area is less suitable for RWH, concentrated in the uppermost northern regions and the lower southern regions. Furthermore, the study indicated that the central and southern parts of the study area are suitable for constructing check dams and farm ponds, while the northern parts are ideal for constructing medium-sized dams. The findings highlight the effectiveness of utilizing remote sensing (RS) data and MCDA in the GIS environment to identify suitable zones and select optimal locations for constructing various RWH structures. The study outcomes propose potential solutions for various issues in the study area, such as water shortages, desertification, and flood risk. The maps that were made can be used to make good decisions and plans for RWH projects. They can help improve the management of water resources and make farming more sustainable in Erbil Province, Iraq, now and in the future. These results may also assist planners in managing rainwater in other regions of the country.

Finally, it’s important to remember that finding the right places for the RWH process based only on environmental and geophysical factors is not enough for it to work. There is a need to intensify research efforts to gain a deeper understanding of these systems’ ability to meet water needs in different contexts, as well as their socio-economic and financial feasibility. Factors such as land ownership, investment and maintenance costs, and labor input should be thoroughly investigated before implementing planned projects.