Evaluation of Integrating SWAT Model into a Multi-Criteria Decision Analysis towards Reliable Rainwater Harvesting Systems

Abstract

Rainwater harvesting (RWH) has been recognized as one of the most reliable and efficient methods for water supply, especially in arid and semi-arid regions (ASARs) facing freshwater scarcity. Nevertheless, due to the inherent uncertainty of input data and subjectivity involved in the selection of influential parameters, the identification of RWH potential areas is a challenging procedure. In this study, two approaches for locating potential RWH sites were implemented. In the first approach, a frequently-used method of the multi-criteria decision analysis and geographic information system (MCDA-GIS) was utilized, while, in the second approach, a novel strategy of integrating the soil and water assessment tool (SWAT) model as a hydrology model into an MCDA-GIS method was proposed to evaluate its performance in locating potential RWH sites. The Mashhad Plain Basin (MPB) was selected as a case study area. The developed potential RWH maps of the two approaches indicated similar patterns for potential RWH areas; in addition, the correlation coefficient (CC) between the two obtained maps were relatively high (i.e., CC = 0.914) revealing that integration of SWAT as a comprehensive hydrologic model does not necessarily result in very different outputs from the conventional method of MCDA-GIS for RWH evaluation. The overlap of developed maps of the two approaches indicated that 3394 km² of the study area, mainly located in the northern parts, was identified as high-potential RWH areas. The performed sensitivity analysis indicated that rainfall and slope criteria, with weights of 0.329 and 0.243, respectively, had the greatest sensitivity on the model in the first approach while in the second approach, the criterion of runoff coefficient (with weights of 0.358) had the highest impact. Based on results from the identification of the potential locations for conventional RWH techniques, pond and pan techniques are the most proper options, covering high-potential areas of RWH more effectively than other techniques over MPB.

1. Introduction

Freshwater scarcity has become a pivotal issue in sustainable development [1,2,3], especially in arid and semi-arid regions (ASARs) where communities are encountering water scarcity problems, not only in agricultural and industrial sectors, but also for satisfying domestic water demands [4,5,6,7]. Most countries in North Africa and the Middle East are facing water scarcity and its related consequences [8]. Moreover, various factors, such as accelerated climate change, population growth, intensive agricultural activities, and industrialization, as well as continuous pollution of water resources, will exacerbate freshwater scarcity and shortages in the future [9,10]. Therefore, more studies are needed in order to focus more specifically on water conservation planning and management, and the sustainable use of potential water resources [5,8,11].

Rainwater harvesting (RWH) can be considered as one of the most cost-effective and environmentally-friendly water conservation methods, especially in ASARs, which simultaneously address water scarcity issues and alleviates groundwater over-extraction [2,12,13,14]. From a historical perspective, archeological evidence found in Southwest Asia revealed that RWH has been practiced since the Neolithic Age (around 10,000 BC to 4500 BC) [11,15]. Moreover, examples of RWH structures were discovered in Tunisia, Greece, China, and historical Palestine, tracing back thousands of years [2,16]. In general, all schemes for collecting and storing rainwater to supply safe and inexpensive water for different purposes, such as domestic, industrial, and agricultural purposes are identified as RWHs [17,18,19].

In developing countries, where domestic water utilization is low, RWH can meet a large amount of water demand [20,21]. Furthermore, harvested water can be managed and diverted for other purposes, such as groundwater recharge, downstream flood mitigation, soil moisture improvement, irrigation, and livestock purposes [1,22,23]. However, identifying RWH potential areas/zones could be challenging, since several criteria need to be considered, including hydrology, climatic condition, topography, and soil parameters, to maximize water availability, particularly in ASARs [24,25].

Many methods have been applied to incorporate contributing RWH factors to locate potential RWH areas. In this regard, the integrated application of multi-criteria decision analysis (MCDA) and geographic information systems (GIS) has been frequently used for a wide range of objectives in water resource research, including potential RWH assessment [26,27,28]. In Saudi Arabia, geospatial and MCDA methods were implemented for identifying potential areas for some RWH structures [29]. In Egypt, potential RWH and storage areas were determined using remote sensing and decision support systems [30]. In South Africa, in-field RWH and ex-field RWH suitability maps were developed by combining ecological, physical, and socio-economic factors [20]. There are other studies that have been done in which assessment of potential RWHs was conducted through the incorporation of biophysical factors (e.g., slope, soil type, land use) [31,32].

As a multi-criteria decision-making (MCDM) system, the analytical hierarchy process (AHP) was introduced by Saaty [33,34]. It is a favorable decision-supporting technique for solving multiplex problems [35]. Particularly, it has been recognized as the most applicable decision method for the identification of potential RWH areas, along with the GIS platform [36]. In the AHP technique, based on the knowledge of experts and mathematics, decisions are organized and analyzed in a structured procedure. The primary idea behind AHP is to use hierarchical symbols to represent the components of any issue and show links between them [35]. Many studies have made extensive use of AHP in order to identify possible RWH sites [37].

The majority of related studies indicate that the utilization of practical approaches, as well as the selection of proper/suitable criteria, are the two main factors in the successful implementation of geospatial and MCDA methods in mapping potential RWH areas [29,38]. In this study, an assessment of potential RWH through the incorporation of biophysical parameters into MCDA, within a GIS environment, as a frequently used approach, is conducted. Furthermore, the novel approach of integration of SWAT as a hydrological model with MCDA and GIS is adopted to evaluate potential RWHs. The integrated approach of hydrological modeling with MCDA and GIS toward the identification of potential RWH areas in ASARs has already been addressed in numerous studies [25,39,40,41]. The primary reason for utilizing a hydrological model is to obtain, either the potential runoff generation depth and coefficient map, or to incorporate the obtained layer into MCDA to enhance the RWH site selection criteria. Additionally, studies revealed that the combination of SWAT with GIS and MCDA reduces the likelihood of inherent biases in GIS-based MCDA.

The runoff coefficient is considered as one of the most influential parameters in RWH assessment, accounting for major contributing factors, such as rainfall, topography, land use, and soil texture and structure [29,42]. Meanwhile, in most studies, a runoff coefficient map is generated by utilizing the widely-used method known as the soil conservation service-curve number (SCS-CN) method [43,44,45]. The SCS-CN method calculates runoff depth based on land cover, hydrologic soil group, and antecedent soil moisture content. While the SCS-CN method is simple, easily applicable, and conceptually stable, it is constrained by some factors, such as the basin area (less than 8 km²), low infiltration capacity, and high rainfall depth [44,46].

On the other hand, in ungagged basins, the SWAT model [47] can be considered as a suitable option for hydrological simulations, particularly in developing countries lacking rainfall and runoff records [48,49,50]. Data availability for most parts of the world, along with other advantageous features, including ease of representation, use of spatially available data, and capability of result presentation through GIS, increase the potential of integrating SWAT with the MCDA-GIS approach for the identification of potential RWH areas [51,52]. Moreover, the SWAT model, which requires a higher number of inputs, is more comprehensive compared to the SCS-CN, and its application can lead to a more accurate runoff coefficient map. In addition, the SWAT model includes, not only all factors used in the SCS-CN method, but also the independent interaction of influential factors, including infiltration and evaporation in rainfall–runoff modeling, which enables capturing more variability rather than relying on a single CN as a lumped parameter [48].

This paper evaluates the impact of incorporating a hydrological model on improving the identification of suitable RWH locations. To enable the evaluation, two approaches, one including hydrological model and the other one excluding the hydrological model, are examined. The integrated MCDA-GIS approach was used as the first approach. The second approach, which incorporates hydrological modeling, was the combination of the SWAT model with MCDA and GIS. MPB is located in northeast Iran and struggles with water supply, especially in recent years. A parameters sensitivity analysis was performed to reduce the subjectivity of parameters selected using expert judgment [48]. In addition, it identified the most sensitive criteria for locating potential RWH zones in MPB. In the end, the most applied RWH structures in ASARs were evaluated for MPB to find the best practices for high-potential areas for RWH.

The inclusive objectives of this research are to: (1) develop and compare the MCDA-GIS approach with the MCDA-GIS integrated with SWAT approach to identify potential RWH zones in the MPB; (2) implement sensitivity analyses on the two approaches to identify the most contributing criteria for potential RWH zone selection; (3) identify the potential location for three conventional RWH techniques (i.e., ponds and pans, terracing, and percolation tank) over MPB. The combined approach of the SWAT model, MCDA, and GIS was evaluated for the first time in ASARs. This paper is categorized as follows: in Section 2, we describe the study basin and the data and methods used in this study. Section 3 presents the results of the evaluation. Section 4 describes conclusions and recommendations for future studies.

2. Materials and Methods

2.1. Study Framework

In this study, a GIS-based MCDA framework has been developed and implemented at a large basin scale for the identification of potential RWH areas. Figure 1 represents the conceptual methodological framework of this study. In the present study, two approaches were developed and compared. In the first approach, five biophysical factors: rainfall, slope, land use, soil type, and soil depth, were incorporated into an MCDA for the evaluation of potential RWH zones. While in the second approach, four factors: slope, rainfall, and soil type, along with the runoff coefficient map derived from the SWAT model were used. The selected criteria in both approaches were weighted using AHP and then normalized and reclassified using GIS techniques. Afterward, different sub-criteria were scored and rasterized. Then, weighted linear combination (WLC) was used in the GIS environment to develop the final potential RWH maps.

2.2. Study Area and Data Collection

The methodology was applied to MPB, located in the northeast of Iran, with an area of 9762 km², between 35°59′ N to 37°03′ N latitude and 60°06′ to 58°22′ E longitude (Figure 2). The climate of the study area is semi-arid to arid, and the average monthly temperature ranges from 11.6 to 26.7 °C [48]. A review of 19 rain gauge stations during 30 years (1979–2010) indicated an annual average rainfall of 303 mm, which varies from 166 to 486 mm [53]. In general, a V-shaped monthly rainfall pattern is observed. January, February, and March are the wettest months, while the rainfall trend decreases until July, the driest month of the year, and increases again until December. Climate change has increased the temperature trend, while no considerable changes in rainfall patterns are observed [54]. The annual average evapotranspiration ranges from 236 to 310 mm [55]. In addition, the mean elevation of the basin is 1487 m above sea level (ASL), which varies from 856 to 3247 m ASL; and the average land slope in the study area is 16.2%. From the geological point of view, the main constituents of the study area are ultramafic, granitic, and metamorphic rocks. In addition, the northern parts of the study area are sedimentary zones, mainly comprised of limestone, conglomerate, dolomite, and gypsiferous marls, and the higher parts of the MPB include discontinuous loess deposits [56].

2.3. Data Collection

The available data were obtained from various governmental organizations and other direct sources, such as field surveys and RS (i.e., satellite imagery) data; data collection focused on the main influential factors to identify potential RWH areas, including maps of rainfall, soil, and land use. Furthermore, data processing and an assortment of supporting techniques were conducted to perform the data analysis. All the generated/available data (GIS layers) were geo-referenced to the MPB geographical coordinate system (WGS_1984_UTM_ Zone_40N).

2.4. Data Processing and Analysis

2.4.1. Criteria Map Generation

Slope map generation was performed using the ASTER DEM surface and the 3D Analyst tool in the ArcGIS environment. The hydrological soil map was reclassified into three soil groups, according to the soil texture and characteristics (infiltration and soil grading). The land use data were merely descriptive and were not readily available to be incorporated directly into the model; the geological field mapping along with the satellite images were used for better land-use data preparation.

The SWAT model was used due to model compatibility for regions lacking proper available data. The SWAT model incorporates various contributing factors in the runoff process and lessens the inherent biases involved in the MCDA-GIS approach [48].

ASTER DEM was utilized in the SWAT model to delineate basin, sub-basins, and hydrological response units (HRUs). Then, weather data for the period of 1981 to 2010 (i.e., CFSR (climate forecast system reanalysis), daily rainfall, minimum and maximum temperature collected from the GSMaP (global satellite mapping of precipitation)) was imported to build the basin model and enable hydrological simulation. Calibration and validation processes were neglected for two main reasons; the first reason is the inadequacy of observed runoff data to examine the model simulations and the comparative and/or relative survey of the results because in the comparative environment, calibration and validation processes are obviously not necessary. Therefore, the non-calibrated model does not have any impact on the relative runoff coefficient and the optimal locations for RWH. After calculating the annual average rainfall for each station, it was interpolated using the kriging method in ArcGIS, which was later averaged over each sub-basin. With the mean runoff depth obtained from the SWAT model, and the average rainfall over each sub-basin, a mean runoff coefficient was determined for all 43 sub-basins within the study area.

After preparing the selected criteria and assigning normalized scores to each sub-criterion, the final map for each criterion was generated and classified into three classes. For all the criteria, except soil depth, soil type, and land use, the classification method of Jenks natural breaks was applied [57,58]. The Jenks natural breaks method, which is considered a conventional method of classification, has been frequently used by previous studies [30,59,60]. Natural break classes are based on natural groupings inherent in the data. Class breaks are identified such that similar values are grouped and the differences between the mean values of classes maximized [57]. It should be noted that the classification of soil depth, soil type, and land use was done based on the literature and expert judgment.

2.4.2. Comparison of Approaches

To compare the potential RWH maps resulting from the two approaches, the correlation coefficient (CC) between the two maps was calculated using the band collection statistics tool in ArcGIS [61]. The correlation between two grid datasets was calculated as the covariance of two grids divided by-product of standard deviations of two grids, as follows:

$$\mathrm{CC}=\frac{{\mathrm{Cov}}_{\mathrm{ij}}}{{\mathsf{\sigma }}_{\mathrm{i}}{\mathsf{\sigma }}_{\mathrm{j}}}$$

(1)

where ${\mathsf{\sigma }}_{\mathrm{i}}$ and ${\mathsf{\sigma }}_{\mathrm{j}}$ are the standard deviation of grid i and j, which here can refer to the resulting potential RWH maps of approach 1 and 2, respectively;${\mathrm{Cov}}_{\mathrm{ij}}$ is the covariance between all pairs of cells of the two grids, calculated as follows:

$${\mathrm{Cov}}_{\mathrm{ij}}=\frac{{{\displaystyle \sum }}_{\mathrm{k}=1}^{\mathrm{N}}{(\mathrm{V}}_{\mathrm{ik}}-{\mathsf{\mu }}_{\mathrm{i}}{\left)\right(\mathrm{V}}_{\mathrm{jk}}-{\mathsf{\mu }}_{\mathrm{j}})}{\mathrm{N}-1}$$

(2)

where V is the cell value, µ is the mean value of each grid, N is the total number of cells in each grid, and k denotes a particular cell in the grid. The CC ranges between −1 to 1, such that CC equal to 1 and −1 showing a fully positive and negative correlation between two grids, respectively.

2.5. Multi-Criteria Decision Analysis

2.5.1. Analytic Hierarchy Process (AHP)

To produce the potential RWH map, AHP was selected as an MCDA technique [62]. AHP has been widely used due to its simple interpretation and implementation, as well as the consistency of its results [59,63,64,65]. However, the main challenge with the MCDA method is how to select the criteria and relative weights, considering the judgment of experts. Therefore, through the use of expert knowledge and extensive review of RWH articles, effective criteria were selected and proper weights and scores were assigned [16,25,66]. Each of the mentioned criteria had different feature classes; the appropriate score for each feature class of the layers was determined and then normalized. Proper scores were assigned to selected thematic layers on a scale of 1–9, proposed by Saaty [67]. A pair-wise comparison of the assigned weights matrices was made using the Saaty AHP approach to determining the weights of the criteria. Subsequently, these desired priority vectors were computed using the eigenvector technique [62] and finally the assigned weights were tested for consistency by computing the consistency ratio as follows [67]:

Consistency Ratio (CR) = (λ_max − n)/((n − 1) × RI)

(3)

where λ_max is the principal eigenvalue computed by the eigenvector technique, n is the number of criteria, and RI denotes random index. More details about the eigenvector technique are provided by Saaty (1980) [67].

It is worth mentioning that the AHP method requires independence between each criterion; while the parameters used in approach 2 (i.e., runoff coefficient, soil type, and slope) are interrelated. Nevertheless, the aims of selecting and including criteria into AHP are two-fold. The selected area for RWH should have, first, the potential of runoff generation and, second, the potential of runoff harvesting. Two areas having the same runoff coefficient and the same amount of rainfall do not necessarily provide the same chance of RWH; the area with a gentler slope and less pervious soil provides a better opportunity for RWH. The runoff coefficient is a key parameter in RWH in the second approach that shows the percent of water that can be available for RWH and has the highest weight compared to the other factors. Rainfall, slope, and soil type are, accordingly, other parameters that determine the feasibility of harvesting rainwater, having lower effective weights in the second approach, respectively.

2.5.2. Rainwater Harvesting Potential Index (RWHPI)

In this study, for each of the two approaches mentioned in previous sections, different criteria were used, as listed in

Table 1. The criteria used in RWHPI.

Thematic Layer	Abbreviation	1st Approach	2nd Approach
Runoff Coefficient	RC	-	√
Soil Type	ST	√	√
Slope	S	√	√
Rainfall	R	√	√
Land Use	LU	√	-
Soil Depth	SD	√	-

. These criteria were selected based on their proven effectiveness in locating potential RWH zones in ASARs, as documented in a literature review [16,25]. All thematic layers, along with their normalized weights were integrated using ArcGIS software to map the potential RWH zones for both approaches. The total normalized weights of different features were overlaid in the integrated raster layer using the weighted linear combination (WLC) as Equations (4) and (5), which were used for the first and second approaches, respectively.

$${\mathrm{RWHPI}=\left(\mathrm{R}\right)}_{\mathrm{c}}{\left(\mathrm{R}\right)}_{\mathrm{s}}{+\left(\mathrm{S}\right)}_{\mathrm{c}}{\left(\mathrm{S}\right)}_{\mathrm{s}}{+\left(\mathrm{ST}\right)}_{\mathrm{c}}{\left(\mathrm{ST}\right)}_{\mathrm{s}}{+\left(\mathrm{SD}\right)}_{\mathrm{c}}{\left(\mathrm{SD}\right)}_{\mathrm{s}}{+\left(\mathrm{LU}\right)}_{\mathrm{c}}{\left(\mathrm{LU}\right)}_{\mathrm{s}}$$

(4)

$${\mathrm{RWHPI}=\left(\mathrm{R}\right)}_{\mathrm{c}}{\left(\mathrm{R}\right)}_{\mathrm{s}}{+\left(\mathrm{RC}\right)}_{\mathrm{c}}{\left(\mathrm{RC}\right)}_{\mathrm{s}}{+\left(\mathrm{S}\right)}_{\mathrm{c}}{\left(\mathrm{S}\right)}_{\mathrm{s}}{+\left(\mathrm{ST}\right)}_{\mathrm{c}}{\left(\mathrm{ST}\right)}_{\mathrm{s}}$$

(5)

where RWHPI is rainwater harvesting potential index, subscript c shows the normalized weight of each criterion, and subscript s denotes the normalized score of a feature class of each criterion. RWHPI is a dimensionless indicator that is useful for finding the high-potential RWH zones within the study area.

2.6. Sensitivity Analysis

It is very important to check the reliability and robustness of the methodology used in this study and the results should be validated to evaluate the impact and importance of each of the factors affecting the model. Sensitivity analysis was performed to obtain insights into the dependence of model outputs on certain model variables [68]. These analyses, in addition to determining the impact of each variable on the RWH analysis, reduce the subjective aspects of the various criteria. The most sensitive identified criteria are ranked as the first priority for future analyses, to ensure more accurate measurements due to their high weights in the analysis. Single parameter analysis introduced by Napolitano et al. [69] was used to estimate the contamination vulnerability of aquifers in several studies [70,71]. In addition, single parameter sensitivity analysis manages the over-parameterization that commonly occurs in hydrological modeling [48,49,72]. The sensitivity analysis method replaces the initial criteria weight used for the AHP with “effective weights” calculated using Equation (6):

$${\mathrm{W}}_{\mathrm{i}}={\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}\frac{{\mathrm{F}}_{\mathrm{i}}{\times \mathrm{S}}_{\mathrm{ij}}}{{\mathrm{V}}_{\mathrm{j}}} \times 100/\mathrm{n}$$

(6)

where W_i is the effective weight of each factor, ${\mathrm{F}}_{\mathrm{i}}$ indicates the initial weight of criterion, ${\mathrm{S}}_{\mathrm{ij}}$ is the value score of ith criterion in jth pixel of the raster map, and V_j is the value of the applied index (i.e., RWHPI) and n is the number of pixels in the raster map of the case study. The sensitivity analysis was applied to both approaches.

After performing the sensitivity analysis, the effective weights were applied to obtain the modified RWHPI. The modified RWHPI assumes the same class scoring and parameters as the base RWHPI but using effective weight values obtained from Equation (6), (i.e., W). The modified RWHPI is obtained considering the effective weight of each criterion rather than the initial assigned weight, which reduces the initial subjectivity of the assigned weights selected by using expert judgment.

2.7. Identification of Suitable Sites for Different RWH Techniques

After identifying potential RWH areas and also discarding unsuitable areas for RWH, different RWH techniques for the potential areas were evaluated. For this purpose, the most common techniques used in ASARs for agricultural and domestic applications were considered, including percolation tanks, terracing, and pond and pan [25]. The suitability map for each RWH technique was generated based on common local practice criteria, as well as a comprehensive review of previous studies. The most common parameters applied for the development, planning, and implementation of such techniques are listed in

Table 2. The common techniques and criteria used for RWH site selection in ASARs.

RWH Techniques	Rainfall (mm)	Slope%	Soil Type	Land Use/Cover	Reference
Pond and Pan	>200	<5	Silty loam, sandy clay loam	Planting forests, low-density pasture, semi-dense pasture, woodland and scrubland, high-density pasture, irrigated cropland and pasture, dry farming	[24,31,73,74,75,76]
Terracing	200–1000	5–30	Clay loam, sandy clay, sandy loam	Semi-dense forest, low-density pasture, woodland and scrubland, planting forests, dry farming, semi-dense pasture, rock protrusions, high-density pasture	[37,77,78]
Percolation tank	<1000	<10	Clay loam, silt loam	Semi-dense pasture, woodland and scrubland, low-density pasture, rock protrusions, high-density pasture, dry farming	[74,76,79]

[25]. In the following sections, some explanations about several RWH techniques are provided.

3. Results

3.1. Selected Criteria in MCDA

In the present study, two approaches for locating potential RWH zones were implemented and compared to examine the effect of incorporating SWAT as a hydrological model to the frequently used approach of MCDA-GIS [24,25,31,75]. In the calculation of RWHPI and generating RWH potential maps, the highest harvesting possibilities of rainwater, as well as runoff generation, were considered. The RWHPI for the first and second approaches were calculated using Equations (4) and (5), respectively. To guarantee the consistency of the analyses, the procedure of weighting thematic layers in AHP was set such that the consistency ratio for all of the obtained thematic layers would be less than the limit of 0.10 [67]. The selected RWH criteria, as well as the assigned corresponding weights, are shown in

Table 3. Weights of the thematic layers (criteria) and their features (sub-criteria).

Thematic Layer	Feature Class	Assigned Score	Layers Relative Weights
			Approach 1	Approach 2
Runoff Coefficient	0.031–0.041	2	-	0.387
	0.042–0.058	4
	0.059–0.087	6
	0.088–0.12	8
	0.13–0.18	9
Rainfall (mm/year)	161–202	2	0.380	0.275
	203–272	4
	273–331	6
	332–395	8
	396–466	9
Slope (%)	0–5	9	0.234	0.198
	5–10	8
	10–15	7
	15–30	5
	>30	3
Soil Type	Silt loam or loam	5	0.179	0.140
	Sandy clay loam	7
	Clay loam or clay	9
Soil depth	Water bodies	9	0.107	-
	Very shallow	8
	Shallow	7
	Relatively deep	5
	Deep	3
Land use	Semi-dense forest	2	0.100	-
	Woodland and scrubland	3
	Planting forests	4
	Scatter forest	5
	Dry farming, irrigated cropland	6
	High density pasture	7
	Low and semi-dense pasture	8
	Urban and built-up land, water bodies	9

. Their selection was based on expert knowledge and extensive literature reviews over several RWH publications for ASARs.

Most of the selected biophysical criteria in this study have been widely used in previous studies in ASARs. Adham et al. [25] showed that, in 48 studies regarding the potential RWH evaluation through GIS, land cover/use (75% of studies), soil type (75% of studies), rainfall (56% of studies), and slope (83% of studies), were used as influential criteria. The rainfall data are one of the necessities/prerequisites to identify potential RWH zones and technically more rainfall in an area means a higher-potential RWH [80]. In ASARs, in terms of the practical effectiveness of RWH techniques, a minimum available rainfall of 200 mm/year was recommended [81,82]. Daily global rainfall CFSR data for 30 years (1981 to 2010) were used and mean annual rainfall data over the sub-basins were obtained. According to the rainfall map shown in Figure 3b, the annual rainfall ranges from 161 to 466 mm, and over 60% of rainfall occurs during December, February, March, April, and May. The northern sub-basins receive a higher amount of rainfall (about 400 mm/year), while less rainfall (around 161–202 mm/year) is received in the southern sub-basins. The rainfall received in the middle parts of the study area ranges from around 203 to 395 mm/year. In this study, direct consideration of evaporation was ignored due to a lack of data and simplification. However, in the second approach, evaporation was indirectly considered. This is because of the fact that the runoff coefficient map utilized as a biophysical criterion was obtained from a SWAT model, which considers evaporation in the process.

The slope is another criterion for the mapping of potential RWH zones. As the percentage of slope increases, the runoff generation increases, while the harvesting opportunity of rainwater for most RWH techniques decreases—this results in less suitability for RWH [29,77]. Hilly areas respond to rainfall events in the way that a high amount of rainwater rapidly runs off and becomes unavailable, even if the mean amount of rainfall is fairly noticeable. That is the case, particularly in the peak domestic and agricultural demand periods [77]. On the other hand, milder areas retain water for a longer time and facilitate groundwater recharge. Moreover, a large amount of earthwork along with implementing erosion control measures are required for constructing RWH systems on hilly areas, which undermines the economic aspect of building RWH systems in areas with a high degree of slope [31,42]. Typically, areas with higher slopes (>5%) are more prone to erosion due to the uneven distribution of runoff over the area [83]. A slope map is shown in Figure 3c. All around the study area, except for the southeastern border, is surrounded by mountains (high elevation areas), while central portions of the study area are located at lower elevations. The slope ranges between 0 and 1% in almost 37% of the study area and the slope of 9.1% of the study area is in the range of 1 to 3%. Steep areas (>15%) cover 4.1% of the basin. The rest of the desired region is divided into three almost equal portions (around 16% of the study area) that have slope ranges of 3–5%, 5–10%, and 10–15%.

Another important criterion for RWH planning and selecting optimal locations for implementing RWH techniques is the soil type [84,85]. Soil characteristics determine the infiltration capacity rate, water holding capacity, and runoff generation potential [78,85,86,87]. The textural features of soil are determined by the percentage of sand, silt, and clay. Naturally, a higher sand portion in the soil implies a higher infiltration rate, lower water storage capacity, and lower runoff generation potential, while the reverse of these properties is the case when the soil consists of a higher clay percentage [78,86,88]. A soil map is shown in Figure 3d. Sandy clay loam with low infiltration rates forms about 25.4% of the soil in the study area; besides, silt loam, or loam with a moderate infiltration rate forms about 49.3% of the study area. The rest of the basin (almost 25.3%) consists of clay/clay loam with a low infiltration capacity rate and consequently higher runoff potential. There are several factors (e.g., vegetation coverage, soil texture) involved in the amount of runoff that is generated from rainfall events which are reflected in a runoff coefficient [87,89]. Since most of the influential factors in runoff generation are considered in SWAT’s structure, the resultant runoff coefficient map reasonably represents case study characteristics.

The runoff coefficient map (Figure 3a) indicates the average runoff coefficient of each sub-basin over 30 years (1981–2010). The runoff coefficient map values range from 0.03 to 0.18. The northern parts indicate higher values while lower runoff coefficients are shown in the southern sub-basins. Typically, a higher runoff coefficient denotes a higher runoff potential; nevertheless, it should be mentioned that successful RWH is not exclusively the function of high runoff potential, but the potential of harvesting the generated runoff should also be considered, which itself is a function of other criteria to ascertain potential RWH.

Land use characteristics dictate the quantitative and qualitative relations between captured rainfall by the basin, groundwater recharge, evapotranspiration loss, and the remaining water running off a basin [88]. Infiltration is lower in urban and pasture-covered areas, which results in higher overland flow. While, as the vegetation coverage becomes intense, water abstraction, interception, and infiltration are increased and its consequences in less generated runoff [40,48,84]. The land use map (Figure 3f) shows that dominated land covers of low-density pasture, semi-dense pasture, irrigated agriculture, and dry farming constitute 23%, 24%, 19%, and 25% of the study area, respectively. It is worth mentioning that due to the incorporation of the runoff coefficient map, the land use layer was not directly included in the second approach. This is because the runoff coefficient map is highly dependent on land use and SWAT requires a land-use map as an input. Thus, the impact of different land uses on RWH is intrinsically reflected in a runoff coefficient map.

The classification of the soil depth map (Figure 3e) based on FAO suggestions (Food and Agriculture Organization of the United Nations) [90] includes very shallow (<30 cm), shallow (30–50 cm), moderately deep (50–100 cm), and deep (100–150 cm). The soil depth map (Figure 3e) indicates that central portions of the study area are mainly covered by deep soils (20.9%), while outer areas are mostly shallow soils (46%). The rest of the basin (27.6%) has relatively deep soil. Areas with shallow soils are potentially more suitable for RWH than areas with deep soils (with some exceptions, such as terrace (earthwork)) due to lower permeability and higher water generation [91]. In other words, a runoff coefficient is higher in areas with a shallower soil depth. The soil depth layer was not directly included in the second approach due to the high dependency of runoff coefficient on soil depth. A runoff coefficient is interrelated to the soil depth and assumed as its proper representative in the second approach.

3.2. Rainwater Harvesting Potential Maps

RWH Potential maps were produced by integrating the selected layers (Table 3) of each approach using the WLC technique in ArcGIS Spatial Analyst Tools. Afterward, RWHPI was computed using Equations (4) and (5) for the first and second approaches, respectively, to generate the potential RWH maps, shown in Figure 4a,b. The output maps of both approaches were classified into three categories of “low”, “moderate”, and “high” RWH potential. The same classification method of Jenks natural breaks was implemented to classify both potential RWH maps to facilitate the comparison of the two approaches. The RWHPI for both approaches ranged in almost the same intervals for all RWH potential classes, as shown in Figure 4a,b. The low classes ranged in 0.083–0.181 and 0.102–0.193, moderate classes ranged in 0.182–0.253 and 0.194–0.242, and the interval of high classes were 0.254–0.352 and 0.243–0.337 for the first and second approaches, respectively. According to Figure 4a,b, although there are some discrepancies between the two developed maps, the general patterns of both maps are almost identical. The areas with high potential RWH are located in the northern parts of the case study area, while a lesser potential RWH was predicted by both approaches for the southern parts. The northwestern and central parts of the study area had moderate RWH potential.

Apart from the almost identical potential RWH area distributions of both approaches over the case study area, the coverage areas of each class are also analogous (

Table 4. Three classes of potential RWH using two approaches.

RWH Potential Classes	Approach 1	Approach 2	Percentage of Overlap of the Two Approaches (%)
	The Portion of Study Area (%)	The Portion of Study Area (%)
Low	29.9	34	32.9
Moderate	24.3	30.5	19.1
High	45.8	35.5	48.0

). Moreover, the overlapping map of both approaches is shown in Figure 4c, where navy blue color shows a high overlap and the white-colored areas indicate no overlap. The overlap percentage of each RWH potential class is also provided in Table 4.

According to Table 4, the low potential RWH class comprises almost 30% of the study area in both approaches, such that 80.4% of the area is in common between the two developed maps. This fact shows the conformity of the two approaches in discriminating areas with less potential for RWH. Areas with high RWH potential cover 45.8% and 35.5% of the study area in the first and second approaches, respectively. While there is a noticeable difference in the percentages of the coverage area, 76.5% of the high potential RWH areas are overlapped in both maps. These areas that encompass 3394 km² (Figure 4c) are mainly located in the northern part of the study area and are considered the most promising sites for RWH since they are obtained by both approaches. In addition, high potential overlap areas can be listed as top priority zones for RWH interventions by decision-makers. The models’ agreement on their high potential RWH increases their credibility. The moderate class comprises almost the same portions (24.3% and 30.5% for the first and second approach, respectively) of the basin; nevertheless, there is less agreement between the two models in this class and the overlapping coverage is 57.4%. The moderate class in the first approach is located in a more concentrated way, while there are more scatterings in the second approach. To more thoroughly compare implemented approaches, the correlation between the two obtained potential RWH maps was evaluated utilizing the Spearman method. The calculated correlation coefficient is 0.914, which indicates that the two maps are strongly correlated, and both approaches yield almost the same results. Based on the high correlation coefficient, percentages of overlapped area for each of the classes (Table 4), and the general pattern of the two potential RWH maps (Figure 4a,b), the following conclusion is obtained.

3.3. Sensitivity Analysis

Criteria weights and sub-criteria scores have a vital impact on the evaluation of the obtained results. Given that their determination is generally made through expert interpretation, the selection of potential RWH zones can be sensitive to weight changes in criteria-related decision weights [92,93]. For example, considering a feature in which the distribution of feature classes is rather homogeneous, this would have a fairly equal influence on the whole study area. Solving this problem requires sensitivity analysis in such a way that the weight values constantly change to check the amount of change in the final results. At first, as an initial weight, the runoff coefficient was considered as the dominant parameter. The next considered significant parameters included rainfall, slope, soil type, soil depth, and land use, respectively.

The results of the sensitivity analyses (Figure 5) showed that the least sensitive factor in approach 1 was land use (with the mean value of 0.075), which was expected due to the weight factors assigned (Table 3). Additionally, in approach 2, the least sensitive parameter was slope (with the mean value of 0.207). Performed comparison of the RWHI map produced from the sensitivity analyses of approach 1 and approach 2 indicates that these two approaches are well comparable in general, therefore it would be rational to conclude that the assumed weight was fairly suitable and the results are reliable as both weights and final results were rather similar.

The result of sensitivity analysis for both approaches revealed that the effective influence of the assumed weights for rainfall and runoff coefficient layers would be less in reality (

Table 5. Statistics of the effective weights from sensitivity analysis.

Approach 1
Parameter	Initial Weight (W), (AHP)	Min (Wi)	Max (Wi)	Mean (μ)	SD (σ)
Rainfall	0.380	0.105	0.639	0.329	0.126
Slope	0.234	0.030	0.475	0.243	0.126
Soil Type	0.179	0.134	0.565	0.265	0.083
Soil Depth	0.107	0.033	0.261	0.088	0.035
Land Use	0.100	0.014	0.197	0.075	0.023
Approach 2
Parameter	Initial Weight (W), (AHP)	Min (Wi)	Max (Wi)	Mean (μ)	SD (σ)
Runoff Coefficient	0.387	0.133	0.580	0.358	0.100
Rainfall	0.275	0.071	0.441	0.226	0.073
Slope	0.198	0.023	0.473	0.207	0.119
Soil Type	0.140	0.099	0.555	0.209	0.087

). This can be attributed to the rather homogeneous distribution of rainfall and runoff coefficient layers and the fact that they were averaged for each sub-basin (see Figure 3a,b). Revising the initial weight of criteria with the effective weight obtained from sensitivity analysis, Equation (6), slightly increased the overall overlap of the two approaches with a previous study. The overlap of approaches 1 and 2 after sensitivity analysis were, thus, 47.9 and 37.7%, respectively. To be more specific, using sensitivity analysis decreases the sensitivity of features by minimizing the difference between the estimated weights and their real influence. Examining the results of the sensitivity analysis of approach 1, it is observed that the parameters of rainfall and slope, both had the greatest impact on the model, in other words, have the highest sensitivity. In approach 2, the runoff coefficient parameter had the highest sensitivity. According to the RWH overlap map of the two approaches, about 50% of the study area has a high potential RWH, about 28% moderate, and about 22% low potential RWH (

Table 6. Three classes of potential RWH in the two approaches after sensitivity analysis.

RWH Potential Classes	Approach 1	Approach 2	Percentage of Overlap of the Two Approaches (%)
	The Portion of Study Area (%)	The Portion of Study Area (%)
Low	21.9	30.3	22.1
Moderate	28.4	28.2	27.8
High	39.7	41.4	50.1

). Generally, the northern half of MPB has a high potential RWH; however, the southern half of the study area has a moderate to low potential RWH.

In conclusion, using CFSR data does not result in an identical output as using observed data due to the high sensitivity of locating potential RWH areas according to rainfall data; however, the global availability of CFSR data enables RWH analysis all over the world. Thus, more attention should be paid to data quality and also to the fact that the high uncertainty associated with them should not be neglected. In the next step/section, after identifying the high-potential RWH areas, the feasibility of different RWH structures in the study area is assessed.

3.4. Potential Zones for RWH Techniques

After identifying the most suitable areas for RWH, various RWH techniques were determined. As mentioned earlier, the most common RWH techniques are ponds and pans, terracing, and percolation tanks. Figure 6 shows the suitable locations for various RWH techniques from the perspective of both approaches 1 and 2. The whole study area is divided into low, medium, high, and unsuitable classes in terms of suitability for different techniques. The suitability map of the study area from the perspective of different techniques was prepared using criteria related to each technique for both approaches 1 and 2 in the ArcGIS environment, the results of which are according to

Table 7. The total potential area of the various RWH techniques in the study area as percentages.

RWH Techniques	Class	Approach 1		Approach 2
		Without Sensitivity Analysis	With Sensitivity Analysis	Without Sensitivity Analysis	With Sensitivity Analysis
Terracing	Poor	5.7	2.2	6.8	6.1
	Moderate	3.3	5.6	1.9	2.2
	Very Good	11.1	12.3	11.3	11.8
	Total	20.0	20.0	20.0	20.0
Percolation Tank	Poor	8.9	7.4	9.8	9.5
	Moderate	11.8	12.8	7.1	10.4
	Very Good	9.2	9.7	13.1	10.0
	Total	29.9	29.9	29.9	29.9
Pond and Pan	Poor	10.3	4.0	10.7	10.9
	Moderate	15.3	19.3	9.1	11.4
	Very Good	18.2	20.5	24.0	21.4
	Total	43.7	43.7	43.7	43.7

. As presented in Table 7, the most suitable areas among the various techniques studied are related to ponds and pans, in which the total area with low, moderate, and high potential for this type of technique is about 44% (or about 4295 km²). After ponds and pans, the most suitable areas, respectively, are tanks with about 30% (about 2928 km²) and terracing with about 20% (about 1952 km²).

3.5. Potential of RWH in MBP Water Management

MPB has been struggling with water scarcity and shortage problems for decades. Currently, surface water supplies about 22% of the water demand in Khorasan Razavi Province, encompassing the MPB, and the rest is provided through groundwater resources. Of the total capacity of the large dams in Khorasan Razavi Province (1549 million cubic meters, MCM), most dams operate at only about 34% of their total capacity. On average, about 6404 MCM are extracted annually from groundwater sources, about 85% of which is exploited through wells, about 9% through qanats, and 6% through springs. However, groundwater recharge potential is only 5300 MCM, which indicates an average annual groundwater shortage of 1100 MCM. Groundwater over-extraction and the low operating capacity of existing dams have eventuated the province’s water crisis. Mashhad, as the center of the province, with an average annual water consumption (industrial, agricultural and domestic water consumption) of 26.46 MCM, has the highest over-exploitation of renewable water in the province (Khorasan Razavi Regional Water Authority [55]).

Nevertheless, based on the demonstrated results (Figure 5), the study area has a high RWH potential. An annual rainfall of 2946 MCM falls on the whole study area of which 1505 and 1569 MCM are on the high potential RWH class, and 952 and 780 MCM is received by the moderate class in the first and second approaches, respectively. Moreover, the rainfall received by the overlap of high potential RWH areas in both approaches is 1436 MCM. The resulting surface runoff within a catchment is one of the decent water resources, once it is managed efficiently, and can be utilized to supply a high portion of water (domestic and agricultural) demands. The results indicate the high potential of the case study area in harvesting rainwater and alleviating the existing water crisis. Finding and implementing efficient water management strategies at a basin/sub-basin scale is an urgent requirement of the MBP. Water shortages will continue to be a major problem in the study area unless efforts are made to achieve sustainable and efficient use of potential RWHs.

4. Conclusions

RWHs are recognized as an applicable and favorable method of water supply, as it conserves existing water resources, while contributes to water scarcity alleviation, particularly in ASARs. Nevertheless, identification of potential RWH areas is challenging due to the inherent uncertainty of input data and the subjectivity involved in selection of influential parameters. In the present study, two approaches for locating potential RWH areas were implemented. In the first approach, a frequently used MCDA-GIS method was utilized, while in the second approach, the SWAT model was included in the analysis in order to examine the effect of incorporating SWAT as a hydrological model into an MCDA-GIS method for RWH assessment. In this study, the SWAT model, as a more comprehensive hydrological model compared to the SCS curve number method was employed for generating the runoff coefficient map.

The resultant potential RWH maps of the two approaches indicated a similar pattern for potential RWH areas. In both potential RWH maps, the northern parts of the study area were categorized as the high potential areas, whereas less harvesting potential of rainwater was demonstrated for the southern areas of the case study. In addition, the overlap of the resultant maps of the two approaches indicated that 3394 km² of the study area is considered a high-potential RWH area. These areas, mainly located in the northern parts of the study area receiving an average annual rainfall of 1366 MCM. Meanwhile, the total capacity of the large dams of Khorasan Razavi Province, where the study area is located, is 1549 MCM. Accordingly, there is high RWH potential in the study area that, if be managed and utilized efficiently, could contribute to water supply and ensure long-term water security.

The high value of the calculated correlation coefficient between the two resultant potential RWH maps (i.e., CC = 0.914), along with provided results, demonstrated that both approaches yield almost identical results. Therefore, it can be stated that the inclusion of SWAT as a hydrological model under the described methodology does not necessarily result in different outputs from conventional MCDA-GIS for RWH evaluation purposes, whereas it demands a higher degree of effort to run the hydrological models. On the other hand, the inclusion of hydrological models can be considered as an efficient strategy to reduce uncertainties embedded in RWH assessment using MCDA. Additionally, there would be less uncertainty associated with the overlapping high-potential RWH areas, since these areas are suggested by a combination of the two approaches, rather than a single conventional MCDA-GIS method.

Moreover, the results of the conducted sensitivity analyses indicate that the rainfall and slope criteria (with weights of 0.329 and 0.243, respectively) have the highest impact on the model in approach 1 whereas the criterion of runoff coefficient (with weights of 0.358) has the highest sensitivity in approach 2. The most-sensitive identified criteria could be ranked as the first priority for future data augmentation to ensure more accurate measurements due to their high weight in the analysis.

Based on the discussed results, a suitability map of the study area from the perspective of three different techniques was developed using various criteria related to each technique for both approaches 1 and 2, and for both models, with and without sensitivity analyses. Results indicate that the ponds and pans technique is the most suitable option for MPB, covering more areas with a high potential of RWH than the two other techniques for both approaches (21% of the study area on average). After that, terracing and percolation tank cover 11.6% and 10.5% of the study area with high potential for RWH on average, respectively.

It should be noted that the modeling cost, data availability and accuracy, and efficiency of these methods can vary significantly from basin to basin. In more developed countries, data can be obtained with minimal effort and have a high accuracy. On the other hand, the cost of developing a hydrologic model may defeat the purpose of integrating the hydrologic model with the MCDA-GIS approach. The abovementioned cost-benefit analysis is clearly beyond the scope of this project. It is recommended that the above approach be considered in other case studies where the flow of data is readily available and the SWAT model parameters can be calibrated. Therefore, the findings of this research can be compared against those of other areas and the reliability of these results can be further investigated.

Moreover, applying hydrologic models as comprehensive as SWAT are usually demanding in terms of data, time, and effort and their advantages need to justify their application. Meanwhile, since in this study generated maps from two approaches that were highly correlated which may not justify its application. This may be because of the fact that we intended to evaluate the efficiency of the model for areas suffering from data scarcity and, consequently, we did not calibrate the SWAT model. Considering this, we recommend applying this approach for areas with sufficient in situ data in order to calibrate the SWAT model and to assess the results from two approaches to check the justification for using a hydrologic model. It is possible that other hydrologic models result in a different conclusion than those drawn in this research and justify the use of hydrologic models for RWH. It is recommended that a group of different hydrologic models be considered to compare two approaches and further evaluate the use of hydrologic modeling in RWH. Moreover, the comparison between the two approaches might be impacted by the utilized MCDA method. It is recommended that a similar analysis be considered using other MCDA methods, such as ANP, and assess possible impacts. In the end, the identified high-potential RWH areas are based on the initial assessment of this study area. More specific and comprehensive studies are required to consider impacts of other hydrological processes, including evapotranspiration, to estimate the actual harvestable rainwater in these areas and fine-tune the water resources management of MPB.

Results of this research can potentially be helpful for further studies for populated areas in a developing country, such as Iran. While there are many technical issues for implementing RWH techniques in a populated area, such as MPB, it is an inevitable need to help satisfy the ever-increasing demand for domestic water. Consequently, it is well within expectations that this approach, or any other comparable one, be considered for application in the near future.