Runoff Harvesting Site Suitability Analysis for Wildlife in Sub-Desert Regions

Abstract

Site selection for runoff harvesting at large scales is a very complex task. It requires inclusion and spatial analysis of a multitude of accurately measured parameters in a time-efficient manner. Compared with direct measurements of runoff, which is time consuming and costly, a combination of a Geographic Information System (GIS) and multi-criteria techniques have proven feasible to address this challenge. Although the accuracy of this new approach is lower than the direct method, conducting in-situ measurements over large scales is not feasible due to its financial issues, a lack of sufficient human resources, and time limitations. To achieve this purpose, climatic, topographic, and soil parameters were used to estimate a runoff coefficient and volume for a single event with the 33%-exceedance probability of maximum daily rainfall in the Kavir National Park of Iran. The main challenges ahead of this research have been a) the large area of the park and the inability to directly evaluate site suitability for runoff harvesting, b) the need for a quick and reliable site evaluation to implement water harvesting measures to address water scarcity, and c) the lack of discharge volume data from water streams (as there are no permanent water streams in the site) and the necessity of reliably estimating runoff in different parts of the park to design water harvesting structures which have been addressed by using GIS and a rainfall-runoff model (Soil Conservation Service Curve Number (SCS-CN)). Site suitability was evaluated for the natural territory of two important wildlife species of the park, namely Gazella dorcas and Ovis orientalis, as the main important food sources of an endangered species named Acinonyx jubatus, commonly known as Persian Cheetah. Saving Persian Cheetah from extinction is currently the top priority for the park managers, which is the main factor behind the species chosen for this research. The Analytic Hierarchical Process (AHP) and fuzzy membership functions were employed to assign weights and standardized thematic layers, respectively. The layers were then integrated using the weighted linear combination method (WLC) to obtain the final suitability map. Accordingly, 38% of the area (846 km²) is suitable or highly suitable for runoff harvesting, while 62% (2623 km²) has a very low potential for this purpose. Afterward, 11 suitable locations were identified to collect runoff. The results indicated that suitable catchments are mainly located on the southern slopes of the Mount Siahkouh as the only major elevation in the area. The storage capacity of the earth embankment in each catchment was estimated based on the upstream area of the catchment and runoff volume. Based on the population of the intended wildlife species and their average water requirement, there is a need for 6500 m³ of drinking water annually. In the best-case scenario and under the circumstance of receiving five rainstorm events a year, only 257 m³ is collectible from all runoff harvesting structures, which is only 4% of the total water demand.

1. Introduction

Most of northern Africa and the Middle East is covered by arid lands. The Middle East is believed to be the most water-scarce area on earth [1]. The climate of this region is characterized by low precipitation and high evaporation, with the annual rainfall rarely exceeding 250 mm. Precipitation mostly occurs in the form of intermittent intense rain showers during the wet season of the year, which in most cases does not exceed three to six months. Because of the lack of enough precipitation and high rate of evaporation, a major deficit in surface water resources occurs that adversely affects the flora and fauna of the region [2]. The problem of water scarcity is compounded by the introduction of climate change into the equation [3]. Decreasing precipitation, increasing temperature, increasing intensity of rain showers and frequency of flash floods are the most important signs of climate change [1]. There is growing competition over water resources in the Middle East as the result of increasing water demand [4], which can exacerbate water scarcity for wildlife [5]. Iran as a special case has experienced major population growth in a very short period (over about 40 years), which has put significant stress on its water resources. This increased pressure on water and land resources could ultimately lead to a major threat to wildlife ecosystems and endangered species endemic to this place [6].

Precipitation in arid and semi-arid regions is scarce, unreliable, and unevenly distributed throughout the year. A large proportion of the precipitation is converted into runoff due to various factors such as a lack of vegetation and high intensity [7]. On the other hand, an extremely high evaporation rate in these areas results in the loss of rainfall [8]. Collecting runoff can be very useful in areas without access to other water sources or where digging wells is not feasible for their high costs and/or poor groundwater quality [9,10]. The stored water can be used to supply villages, schools, houses, small gardens, livestock, and wildlife. Therefore water harvesting has been historically practiced in arid and semi-arid regions to collect runoff for later use [11]. In recent years, Iran has experienced the occurrence of severe droughts which have put unprecedented pressure on water resources. Drought and water scarcity are two important factors threatening the survival of wildlife, especially in arid areas. During a drought, supplying water to population centers is prioritized, putting the life of the wildlife at risk. Therefore, augmenting the amount of water resources in natural arid regions can prevent possible losses of animals in dry years [12].

National parks, wildlife sanctuaries, protected areas, etc. are all intended to protect the survival of the wildlife species whose natural ecosystems have been severely destructed by human activities (urbanization, agriculture, animal husbandry etc.). Many of these sanctuaries have very unfavorable climatic conditions, which has naturally prevented the settlement of human populations. The intensification of drought periods, along with the effects of human activities on land and water resources (shrinkage of surface and groundwater resources) have added to the severity of protected areas degradation. For instance, Aryal et al. [13] reported the interactions of wildlife species and the human population as the result of climate change and grazing pressure in the Trans-Himalaya region of Nepal. Supplementing drinking water seems to be indispensable for dealing with water scarcity and drought. With this in mind, employing water collection methods to prevent wildlife losses along with improving the distribution of wildlife species in the region could benefit their reproduction and ensure balanced use of the resources of wildlife sanctuaries. Ogutu et al. [14] also found that the distribution of watering points could favor ungulate distributions in an East African savanna.

In addition to determining water supply methods, their location and prioritization should also be considered. In this regard, the physical conditions of the area, the type, and distribution of wildlife, and requirements for each method can be very highly involved in determining suitable locations. The technological application of Geographic Information System (GIS) in the literature is gaining momentum as it has simplified the calculation of complicated formulas and the illustration of the results in a user-friendly way [15,16]. More recently, attention has been channeled into the application of multi-criteria decision-making methods along with GIS for site selection. Multi-criteria decision making (MCDM) is used to decide among a set of alternatives based on conflicting criteria [17]. A large and growing body of literature has applied MCDM, Remote Sensing (RS) and GIS technologies to investigate land suitability for water harvesting. In a study which set out to determine rainwater harvesting potential and identify suitable sites for rainwater harvesting, Singh et al. [18] used a GIS-based multi-criteria decision analysis method. They argued that the application of GIS and multi-criteria methods are not only time saving and cost-effective, but they are very efficient for water resources planning and management on larger scales. In another study conducted by Wu et al. [19], GIS was successfully utilized for site selection for rainwater harvesting in Northeastern Guatemala. In the case study provided by Bhagwat et al. [20], GIS and RS were employed to identify suitable sites for runoff harvesting and artificial recharge in India. Besides GIS and MCDM techniques, Analytic Hierarchy Process (AHP) and Soil Conservation Service Curve Number (SCS-CN) have been frequently cited for weighting thematic data layers and mapping runoff potential in site suitability analysis [21,22,23]. Although extensive research has been carried out on-site selection for agricultural and drinking water purposes, up to now, far too little attention has been paid to site selection for water collection for wildlife. Among the few cases available, Jafari Shalamzari et al. [12,24] employed GIS and MCDM for site selection for the installation of small-scale wildlife water guzzlers and small-scale solar ponds in the Kavir National Park of Iran. An overview of the application of RS and GIS for wildlife habitat selection has been provided by Acharyya et al. [25]. Therefore, the main objective of this study is to evaluate the suitability of a semi-desert area in Iran for runoff harvesting for wildlife. The Kavir National Park, in central Iran, was evaluated for its runoff harvesting potential. Lack of weather and hydrometric stations, lack of permanent streams, and difficulty of direct access to the site (covering an area larger than 7200 Km²) necessitate using other methods to evaluate site suitability for runoff harvesting. Application of GIS and hydrological models gives us the possibility to reliably evaluate the area for water harvesting. The results of this study could be very helpful for managers of the park to weigh up their options for water collection and decide on the best option. There is no information available on whether there has been any attempt to harvest water on the site, but as it appears that most of the water deficit is satisfied by carrying water in tankers to watering points in the park. Therefore, this work is unique in evaluating this site’s suitability for water harvesting for the first time. This research is in fact an extension of previous attempts to evaluate site suitability for water harvesting in Kavir National Park [12,24] in which the site suitability was evaluated for small scale rain collectors and solar ponds for water desalination. In this paper, however, we have attempted to evaluate the park’s suitability for large-scale runoff harvesting, and to locate the best candidate catchments for constructing earth embankments.

2. Materials and Methods

2.1. Study Area Delineation

The Central Kavir Area of Iran consists of the Kavir National Park and a Protected Area enclosing 680,000 hectares (Figure 1). The Kavir National park as one of the main remaining territories of Persian Cheetah (which also includes parts of the Tehran, Semnan, Kerman, Golestan, Yazd, and Khorasan Provinces of Iran), was selected as a case study in this work, although the methodology adopted will be applicable to other similar areas in Central Iran. Kavir National Park is bounded from the north by the Varamin, Iwanaki and Garmsar plains, from the east by the Kavir Playa, from the south by the Kavir Salt Pans and flood plains, and from the west by the Salt Lake and the flood plain of Qom. This area was approved as a national park in 1976, as the only park in Iran where no construction of mines, issuance of grazing licenses, and residence of the population is allowed. The study area’s elevation ranges from 684 m to 1922 m above sea level. Much of the Park could be classified as a plain area bounded by the northern and southern heights. Considering the above characteristics, the study area can be divided into three zones consisting of the northern heights, the middle plain to southern heights, and the eastern hills and sand dunes. The southern heights and the middle plain are volcanic with igneous formations. The soil formed on these geological formations is not restricted in terms of salinity and alkalinity, and is seen as a series of interconnected alluvial fans. The northern section is part of a larger playa which contains mainly saline alluvial quaternary formations. The hilly area to the east is mainly composed of moraine formations, with an eroded surface forming a badland. The hills are interspersed with small flat plains. Most of the area, especially in the northern section, is covered by heavy soil texture (loamy to clay) with a low infiltration rate which provides suitable conditions for runoff generation (see Figure 8). The infiltration rate is comparatively higher (measured using the double-ring method) in the eastern section of the area, since it is covered with lighter soil textures (loamy to loamy sand) and sand dunes [25]. The area has cold winters and hot summers. The temperature fluctuates between -5 to 50 degrees Celsius, with an average annual temperature of 19 °C in flat areas declining to 13 °C at higher elevations. Average annual precipitation is 240 mm (50–300 mm), with the rainy season lasting from November to May, with occasional snowfall in January and February [25]. Figure 2 illustrates the Ombrothermic Diagram of Garmsar as the closest weather station to the Kavir National Park, where precipitation is plotted against temperature, indicating that there is a dry season when the temperature line runs above the precipitation line. Precipitation varies considerably between years, with some moderately wet years in between. Therefore, and as the result of the lack of reliable surface water sources, most of the required water in the area comes from the few drilled wells. Groundwater, which reaches the surface in mountainous areas as springs, is mostly saline and cannot be used by wildlife. In flat areas, groundwater is mostly found at very shallow depths with saline to brackish water due to the rapid evaporation and the presence of saline geological formations. Out of seven wells in the area, four are in current use, one is dried and two others are of poor quality. The high temperature of the area causes strong evaporation of water resources, leading to the development of salt marshes around the lake’s margins. Mean annual evaporation based on the evaporation pan measurements is 1541 mm [25], which is exacerbated by hot winds during summer. The presence of saline geological formations has led to the deterioration of groundwater and surface water quality, and acts as a limiting factor for water resources.

2.2. Criteria Selection for Suitability Assessment

Criteria for site suitability analysis in this paper were selected based on the literature review and with respect to the experience of the experts working in the area. For this purpose, we included those studies using a multi-criteria analysis approach for runoff suitability analysis.

Table 1. Criteria used in the literature for runoff harvesting site selection.

Reference	Criteria
Rishab Mahajan and Khitoliya [26]	Rainfall duration and intensity, vegetation cover, land-use, slope gradient, evapotranspiration, ecological consideration
Mbilinyi et al. [27]	Rainfall, soil depth and texture, land cover, vegetation cover, drainage network
Munyao [28]	Slope, land cover, runoff coefficient, soil texture, distance to road and buildings
Weerasinghe et al. [29]	Soil depth, soil texture, land cover, monthly precipitation
Kadam et al. [30]	Land cover, slope, soil texture, rainfall, drainage density, residential areas
Jha et al. [31]	Land cover, soil texture, daily rainfall data, slope, drainage density, runoff coefficient
Krois and Schulte [32]	Climate, hydrology, topography, agronomy, soils, and socioeconomic aspects
Ochir et al. [33]	Land cover, rainfall, slope, soil type, distance to roads
Tiwari et al. [34]	Drainage network, rainfall data, soil texture, land cover
Wu [19]	Slope, land cover, distance from field and road, runoff coefficient, soil texture

provides a summary of the literature reviewed along with the applied criteria. Accordingly, rainfall, soil texture, and infiltration, slope and drainage network were selected as the basic thematic layers in this research. The choice of layers was basically aimed at satisfying the requirements of the rainfall-runoff model. However, other factors such as forage abundance and quality, location and distance to guarding posts, distance to roads, location of current watering points among others were left out of the analysis for the sake and simplicity of analysis, easier interpretation of the results, and the lack of data. Compared with earlier studies on runoff harvesting site suitability analysis, we also included soil salinity as an important measure controlling collected runoff quality. Moreover, the distribution of wildlife species significantly controls the final suitability of the area and hence was selected as a factor. However, we did not include the role of vegetation cover (or land cover in general) in water harvesting site suitability analysis since the area has a very sparse vegetation which cannot significantly affect runoff generation and water loss and the fact that there are no other limitations in terms of land cover for runoff harvesting such as residential areas, mining activities, water bodies, etc. We also did not include evapotranspiration as it affects the storage of runoff and not its collection.

2.3. Data Collection

We collected data on precipitation, soil (texture, infiltration, salinity), topography (Digital Elevation Model (DEM)) and wildlife distribution to conduct this research. As there is no climatological station inside the area, we obtained the data from 12 climate stations distributed around the Park in central Iran. The slope gradient and drainage network layers were generated by using the ASTER Global Digital Elevation Map at the spatial resolution of 30 m in ArcGIS 10.3 (developed by Environmental Systems Research Institute (ESRI) located in Redlands, California, US) and the ArcHydro Tool, respectively. In order to prepare soil texture, infiltration and salinity maps, paper map sheets of the area at 1: 250,000 scale were obtained from the Semnan Province Environmental Protection Head Office and digitized in ArcGIS environment. Data on the distribution of the two species of wildlife (Gazella dorcas and Ovis orientalis) was also provided by Semnan Province Environmental Protection Head Office. Daily water requirement for wildlife species of Ovis orientalis and Gazella Dorcas was determined from the literature and used to estimate water demand based on the number of animals (data obtained from Semnan Province Environmental Protection Head Office).

2.4. Soil Conservation Service Curve Number (SCS-CN) Method

A broad overview of rainfall-runoff models and their inter-comparison is provided in reference [35]. As the author argues (also see reference [36]), in data-scarce regions, using fully distributed hydrological models (such as SWAT (Soil and Water Assessment Tool) [37,38], SHE (Système Hydrologique Européen (SHE) [38] or tRIBS (TIN-based Real-time Integrated Basin Simulator) [39]) is not feasible mainly because of the lack of data and human resources familiar with complex modeling techniques. On the other hand, for ungaged data-scarce catchments and for the time the discharge values need to be estimated for long term regulatory purposes, lumped and semi-distributed models (such as IHACRES (identification of unit hydrographs and component flows from rainfall, evapotranspiration and streamflow) [40] or Soil Conservation Service Curve Number (SCS-CN) [22] among others) have proven to be reliable tools for runoff estimation, given their relatively simple logic and computation process [41]. The SCS-CN model has a long history of application in arid and semi-arid parts of the Middle East and Iran [42,43,44] in ungaged watersheds and therefore we also adopted this model for the purpose of rainfall-runoff simulation in this research. The SCS-CN method is a simple and commonly used technique for the calculation of direct runoff for a given duration of rainfall events [45] using commonly available data inputs such as topography, land cover, climate data, and soil properties. The surface runoff is estimated as:

$$\mathrm{Q}=\frac{{\left(\mathrm{P}-\lambda \mathrm{S}\right)}^2}{\mathrm{P}+\left(1-\lambda \right)\mathrm{S}}$$

(1)

where Q is the direct surface runoff depth from an area (mm); P is the total precipitation (mm); S is the potential maximum retention after runoff initiation (mm); and, $\lambda$ is surface runoff abstraction (dimensionless).

The parameter $\lambda$ is assumed to be equal to 0.2 for practical applications [22]. The maximum retention parameter S can be calculated from the curve number as:

$$\mathrm{S}=\frac{25400}{\mathrm{CN}}-254$$

(2)

where CN is the curve number ranging from 0 to100, indicating the runoff response to a given rainfall.

CN values in this study were obtained using raster calculator in ArcMap 10.3 software according to soil hydrological groups, land cover, and antecedent soil moisture condition [46] based on the standard tables provided by the US Soil Conservation Services (USCS) [47]. Soils were assigned to different hydrological groups based on infiltration capacity and soil texture [47]. Land cover of the whole park was considered as subtropical steppe vegetation. Based on soil hydrological groups and land cover, suitable CN values were identified. As soil moisture content due to antecedent precipitation events could affect the amount of runoff, the curve numbers were then adjusted according to Antecedent Moisture Condition I (AMC I) [48] for dry conditions. Higher CN values indicate that a greater proportion of the rainfall is converted into surface runoff [47]. The runoff coefficient was then calculated as the ratio between total runoff depth and the total rainfall depth.

2.5. Fuzzification of Criteria Maps

In order to combine the thematic layers, they must be associated with a common scale of measurement and hence normalized. There are a number of common methods for data normalization such as Min–Max, Z-Score, Median Normalization, and Fuzzy Transform (see reference [49]). However, Fuzzy Transform is less prone to outliers and give the user more flexibility as to how to reclassify values in the range interested using fuzzy membership functions [50]. In this study, we used the fuzzy transform method for converting data layers into a common 0–1 scale. The fuzzy set theory was first introduced by Zadeh [51] to deal with the conversion of words into mathematical terms. Accordingly, the fuzzy set X has its fuzzy subset A defined by a membership function $f_A\left(X\right)$ which maps each element of x onto a real number in the range of zero and one [52]. A larger $f_A\left(x\right)$ will result in stronger membership for $x$. Among all fuzzy membership functions of Gaussian (assigning membership value of one to the midpoint decreasing to zero at both tails), Near (giving membership value of one to the range defined by the user), Linear (linearly transforming data to membership values at defined ranges), Small (where small values of the data have high membership in the set), Large (where large values of the data have high membership in the set), MSLarge (where M and S stand for mean and standard deviation multiplication respectively), and MSSmall (giving higher importance to lower and higher ranges of the data compared with small and large functions), the small and large functions were chosen and applied in ArcGIS. The Near and Gaussian methods were not selected as we were not interested in the mid-range of the values and the MSSmall and MSLarge functions were not applied as they exaggerated the role of low and high values in the analysis. Large fuzzy membership function is defined as below:

$$f\left(x\right)=\frac{1}{1+{\left(\frac{x}{f^2}\right)}^{-f^1 }}$$

(3)

Similarly, using the small fuzzy membership function, smaller input values have a greater possibility of becoming a member of the fuzzy set. The small fuzzy membership function is hence defined as:

$$f\left(x\right)=\frac{1}{1+{\left(\frac{x}{f^2}\right)}^{f^1 }}$$

(4)

where $f^1$ is the spread parameter defining the shape and character of the transition zone and $f^2$ is the midpoint, after which numbers have a higher possibility of becoming a member of the set.

Accordingly, slope, runoff potential, and drainage density layers were standardized using small, large, and small functions, respectively.

2.6. Runoff Harvesting Potential Assessment

For identifying areas suitable for runoff harvesting, three major factor layers of slope, runoff coefficient, and drainage density maps were combined. Slope map was prepared from the digital elevation map (DEM) and drainage lines were extracted from DEM using ArcHyrdo extension in ArcGIS. The drainage line map was then converted into a drainage density map using the line density function in ArcGIS. The Analytic Hierarchical Process (AHP) and the Analytic Network Process (ANP) developed by Saaty [53] are the most commonly used methods in MCDM for identifying the importance of each data set to a final decision [54]. The ANP compared with AHP is more complicated for implementation as a standard tool for practical decision making or for describing the problem by means of networks rather than hierarchies as in AHP [55]. Therefore, for introducing the importance of each layer into the final suitability analysis, we used the Analytic Hierarchy Process. The layers along with their classes were assigned a score between 1 and 9 (

Table 2. The scoring scheme used for generating pairwise comparison matrices in the analytic hierarchy process (AHP) [53].

Scores	Importance
1	Equally important
3	Moderately more important
5	Strongly more important
7	Very strongly more important
9	Extremely more important
2,4,6,8	Intermediate values between two levels of importance

) by five experts familiar with the topic and the Park from wildlife ecology, hydrology, GIS, RS, and climatology backgrounds. There is no hard rule on the sample size needed for implementing AHP, but as long as the consistency ratio (CR) is below 10%, the results are regarded as unbiased and acceptable [56]. The weights were calculated in Expert Choice Software 11.0 (developed by ExpertChoice located in Arlington, Texas, US).

Over the last few decades, a number of multi-criteria evaluation methods have been introduced into GIS among which one could mention Weighted Linear Combination (WLC), and Boolean Overlay (by combining layers based on intersection (AND) or union (OR)) as the most straightforward methods. The Boolean method, however, does not allow trade-offs between input variables, which is not the case in nature. However, WLC allows for a full trade-off between factors so that a low score on one factor could be compensated for by the high score on the other factor, which is a more flexible method for suitability analysis. Therefore, the WLC method was selected for conducting the analysis in this paper. Eventually, the three layers were combined by applying their corresponding weights to obtain the final suitability map for runoff harvesting in ArcGIS. The three thematic layers were integrated using the Weighted Linear Combination (WLC) in ArcGIS Software as follows [31]:

$$RHSI=R{C}_{wt}+S{L}_{wt}+D{D}_{wt}$$

(5)

where RHSI is the runoff harvesting suitability index; $R{C}_{wt}$ is the normalized weight of the runoff coefficient (by multiplying the thematic raster layer with its corresponding weight in Raster Calculator); $S{L}_{wt}$ is the normalized weight of the slope map; $D{D}_{wt}$ is the normalized weight of the drainage density map.

Two constraining factors were also included in our analysis. Soil salinity was used to exclude those areas that could affect runoff quality. On the other hand, areas outside the wildlife species’ territories were also masked out of the final suitability map. In order to apply the constraining factors, their corresponding layers were converted into binary formats of zero and one where the former indicates areas to be masked out of the final map and the latter indicates areas suitable for water harvesting. The final suitability map was again combined with the constraining factor layers as follows:

$$SI=RHSI\times EC\times WD$$

(6)

where SI is the final suitability index, RHSI is runoff harvesting suitability index, $EC$ is the soil salinity constraint map, and WD is the wildlife distribution map.

Multiplication of the two constraint layers (in binary format) with the RHSI layer (in raster format with values ranging between 0 and 255) using Raster Calculator of ArcMap software results in the identification of those areas that are suitable for runoff harvesting and removal of those areas limited in terms of the constraints introduced. The process of demarcating suitable areas for runoff harvesting is illustrated in Figure 3.

3. Results

3.1. Preparation of Precipitation Map

Table 3. Summary of the climatic parameters measured at the selected climatic stations surrounding the Kavir National Park.

Rain Gauges	Elevation	Distance (km)	$\mathbf{Mean}\mathbf{Annual}\mathbf{Temp}.\left(℃\right)$	Annual Prec. (mm)	Aridity Index	Mean Seasonal Prec. (mm)	$\mathbf{Mean}\mathbf{Seasonal}\mathbf{Temp}.\left(℃\right)$	Maximum Daily Rainfall (mm)	EPMDR *	Data Record (years)
Abardej	900	73	19.6	124.3	2.4	118.6	13.8	36	26.5	13
Aminabad	1000	93	16.2	194.8	4.0	175.2	11.4	50	26	44
Badroud	1056	143	18	110.2	2.2	102.9	12.1	70	25	34
Hamand	1800	81	16	337.4	6.9	287.6	5.7	55	35.5	43
Mamazan	1021	73	18.3	147.2	2.9	135.3	12.7	38	23.5	26
Abali	2465.2	96	10.3	537.9	12.4	473.8	3.5	91	46	27
Firouzkouh Pol	2985.7	95	10.4	399.7	9.2	315.4	0	38	33	16
Firouzkouh	1975.6	124	10.2	289.3	6.7	225.6	5.1	40	33	16
Garmsar	825.2	34	19.7	121.7	2.3	112.4	14.3	38	20	21
Kashan	982.3	67	19.7	128.4	2.4	118.0	13.9	58	24	48
Natanz	1684.9	88	16.4	198.1	4.0	189.1	10.9	50	29	14
Semnan	1684.9	93	18.7	143.4	2.8	127.5	12.4	41	22	45

* 33% Exceedance Probability of Maximum Daily Rainfall (EPMDR).

contains a summary of the climatic variables from the selected weather stations used in this study. Average annual precipitation ranges between 110 and 537 mm, corresponding to the elevations of 1056 m and 2465 m, respectively. Approximately 93% of precipitation falls between the months of October and May. Maximum daily rainfall ranges from 33 mm at Abardej station to 91 mm at Abali Station. In this table, a simple aridity index as mean annual precipitation to mean annual temperature known as the Köpen index (see reference [57]) is also provided for better comparison of the stations as:

$$AI=Mean Annual\frac{Mean Annual Precipitation \left(mm\right)}{Mean annual Temperature \left(℃\right)+33}$$

(7)

where AI is the aridity index (dimensionless), which is higher in wet stations and lower for dry stations.

The estimate of rainfall depth or intensity in terms of a probability for an especial duration (such as 24 hrs. in this case) is required for designing drainage projects or runoff control structures. The probability of exceedance refers to the probability that the actual rainfall during the period of interest exceeds an estimated depth [58]. As there is no data on the maximum daily rainfall for the study area, similar to the methodology of Krois and Schulte [32], the correlation between altitude and rainfall for the selected stations was investigated. Based on our knowledge of the site, we assumed that using the relationship between the climate stations recorded precipitation values and altitude would result in a more accurate picture of the spatial variation of rainfall in the park than using Remote Sensing Precipitation products such as TRMM (Tropical Rainfall Measuring Mission), PERSIANN (Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks), GSMaP (Global Satellite Mapping of Precipitation), CMORPH (Climate Prediction Center’s morphing technique) (a more in-depth comparison of the products is provided in the work of Zeng et al. [59]) or using a single value of the nearest weather station for the whole area. The coefficients of determination for the regression between altitude and 33% exceedance probability of maximum daily rainfall, average annual precipitation and mean seasonal precipitation were 0.76, 0.72 and 0.56, respectively. As for the 33% Exceedance Probability of Maximum Daily Rainfall (EPMDR) [60], Semnan and Abali had the highest values at 22 and 46 mm, respectively. The exponential relationship between the 33% exceedance probability of maximum daily rainfall and altitude was used to produce the distributed map of precipitation for the study area (i.e. precipitation (mm) = 18.847e^{0.0003×Altitude}) (Figure 4). This map was used as the design storm in the model, which is believed to be the most probable scenario for designing runoff harvesting structures [32].

3.2. Preparation of Slope Map

The slope map of the area was prepared using the ASTER DEM 30 m, after filling the sinks in the original DEM image. As also indicated in Figure 5, much of the area falls within the flat to moderately gentle classes (81%). Nearly flat or gently sloping areas (3%–5%) are suitable for the construction of small earthen embankments or even small ponds [18].

3.3. Preparation of drainage network map

Drainage network of the area was prepared using DEM in ArcHydro Extension of ArcGIS. As illustrated in Figure 6, the area has a third-order drainage network better formed in the northern section of the park. According to the Strahler’s stream order classification, the first, the second and the third-order streams had a length of 929 km (65.1% of the total length), 374 km (26.2%) and 125 km (8.8%), respectively. Drainage density map of the study area is illustrated in Figure 7. Accordingly, the whole area was divided into four classes of very poor drainage density (0–0.3 km.km^-2), poor drainage density (0.6–0.9 km.km^-2), moderate drainage density (0.9–1.2 km.km^-2) and good drainage density (0.9–1.2 km.km^-2). The classification of surface area of different categories showed that more than 85% (3394 km²) of the total area falls within the very poor to poor drainage categories, which implies the suitability of the majority of the Kavir National Park for runoff harvesting. The good drainage category, which has simultaneously the lowest potential for water harvesting, encompasses only 0.8% (30 km²) of the total area.

3.4. Soil Map and Soil Hydrological Groups

Soil map of the study area was classified based on soil texture and infiltration into four different hydrological categories of Group A, Group B, Group C, and Group D [46]. Soils falling in the categories A and B have high infiltration rate and low runoff potential while groups C and D have low infiltration rate and high runoff potential. As illustrated in Figure 8, six major soil texture classes exist in the area: Loamy Sand, Loam, Silty Loam, Silty Clay Loam, Silty Clay, and Clay. It is apparent from Figure 8 that the Silty Loam class is the dominant soil texture in the area, covering approximately 1184 km² (26.6%). The next four classes of soil texture in order of coverage were Loam at 928.8 km² (22%), silty clay loam at 856.4 km² (19.2%), silty clay at 673.2 km² (15.1%), and loamy sand at 514.99 km² (11.6%). The very heavy soil texture (clay soils) is found as a strip extending in the northern section of the park by encompassing an area of smaller than 245.6 km² (5.5%). Figure 9a and Figure 9b provide overviews of the distribution of different soil hydrological groups and curve numbers over the study area. As shown, C and D groups each cover approximately 796.9 km² and 1886.1 km², contributing to more than 67.4% of the total area. Hydrological group A covers 1318.2 km² (33.1%) and group B has the smallest coverage with 231.9 km² (5.8%).

3.5. Constraining Factors

Since not all of Kavir National Park is intended for water harvesting, two major limiting factors were introduced into the framework. Soil salinity could affect the quality of the collected runoff and its suitability for drinking. Although different species have different tolerance for water salinity, saline water could severely affect animals’ health and reproduction [60]. Data on soil salinity was obtained from the soil map at the scale of 1:25000 and areas with the salinity of higher than 8 ds.m⁻¹ were identified to be discarded from the final suitability map. Figure 10 provides the results obtained for soil salinity classification. As illustrated, the highest level of salinity exists in the northern and eastern portions of the area spreading over 700 km² (15%). Less soil salinity occurs in the southern part of the area covering 172 km² (4%). Soil salinity map was converted into a binary layer with zero and one indicating unsuitable and suitable areas, respectively.

Another major factor for the assessment of land suitability for runoff harvesting in the Kavir National Park is the distribution of wildlife species. According to the direct interview with the experts of the Semnan Province Environmental Protection Head Office, conservation of two species, namely Ovis orientalis and Gazella Dorcas, is of high importance. According to the latest statistics of the Semnan Department of Environment (direct interviews), Kavir National Park is home to 13 cheetahs, 900 Gazella Dorcas, and 2100 Ovis orientalis. Acinonyx jubatus, commonly known as Persian Cheetah, mainly preys upon these two species which make them important targets for water harvesting. Figure 11 illustrates the distribution of Ovis orientalis and Gazella Dorcas in the Kavir National Park. The wildlife distribution map was also converted into a binary layer where zero indicates the areas outside the territory of the intended wildlife species (and hence unsuitable) and one indicates areas suitable for runoff collection. The total area inside the territory of these two species amounts to 2980 km² (67%).

3.6. Runoff Coefficient and Runoff Potential Map

Runoff potential of the study area for a single event was simulated using the 33%-exceedance probability of maximum daily rainfall and the result is shown in Figure 12a. The obtained map was classified into five categories of equal intervals. Accordingly, most of the central and northern parts of the study area fell under the very poor runoff generation potential class, with an area of 1790 km² (41%). Only 10% of the area (450 km²) fell under very good runoff generation potential class, which indicates a very limited opportunity for runoff harvesting. In terms of runoff volume, the ‘very poor’ class produces 30 m³.ha^-1 runoff on each event reaching up to a maximum of 135 m³.ha^-1 under a ‘very good’ condition. A runoff coefficient map was obtained by dividing the total runoff generation potential by the design rainfall (Figure 12b). The final map was classified into five categories representing very low to very high runoff coefficient values. Like the runoff potential map, most of the northern parts of the area fell under the very poor runoff coefficient class, covering 1450 km² (33%).

3.7. Land Suitability for Runoff Harvesting

Figure 13 illustrates the suitability of Kavir National Park for runoff harvesting. The final map was classified into five categories of land suitability. Accordingly, 38% of the area (846 km²) was suitable to highly suitable for runoff harvesting. The lower three classes from not suitable to moderately suitable enclosed 62% of the total area (2623 km²). Therefore, a large proportion of the park is not suitable for runoff harvesting. As also indicated in Figure 13, only a narrow strip extending in the southern part of the park fell under the optimum suitability class. Most of the central and northern portions of the park had a less favorable condition.

3.8. Runoff Harvesting

Small embankments and check dams are the best ways to collect surface runoff. For the construction of these structures, local natural materials can be used and there is no need for heavy types of machinery. By integrating five thematic maps, suitable locations for runoff harvesting were identified and the results are provided in Figure 14. Locations of these small structures were identified according to final suitability map, stream order (second and third orders are most suitable), and the location of small tributaries. By considering the arrangement of tributaries, those originating from surfaces covered with saline soils were discarded for the sake of the final quality of collected water. In total, 11 suitable locations were identified by visual interpretation of the suitability map. The upstream catchments were separated and the total amount of runoff collectible from the surface was estimated for a single event. Based on the total amount of collectible runoff, the identified options were prioritized. According to the results provided in

Table 4. Estimated amount of collectible water at each runoff harvesting structure based on the quantity of generated runoff from upstream catchments.

Check Dam Code	Collectible Runoff (m3)	Upstream Area (km2)	Ranking
A	0.26	72.8	11
B	1.53	32.7	9
C	1.59	37.9	8
D	18.34	225.1	1
E	5.27	118.6	2
F	7.92	96.9	3
G	3.29	45.7	6
H	4.22	41.8	5
I	1.09	45.5	10
J	5.09	77.1	4
K	2.98	38.4	7

, catchments D, E, and F obtained the highest priority for runoff harvesting. These catchments are located on the southern slope of Mount Siahkouh (Figure 1). The catchments selected in the north and north-western part of the park had the lowest priority for runoff collection.

4. Discussion

In this study, we combined multi-criteria analysis and GIS to identify suitable sites for runoff harvesting in Kavir National Park of Iran. Soft computing techniques have shown to hold a promising capability of simulating runoff [61,62]. The first step in combining GIS and MCDM is the definition of a set of decision factors. These factors were selected based on the literature review, our knowledge of the field as well as the experts familiar with water harvesting and wildlife management. The factors included in site suitability analysis were slope, rainfall, soil texture and infiltration, and drainage network. The topography is an important requisite for site suitability for runoff harvesting. As the slope increases, the suitability of the site for runoff collection and storage decreases. Soil texture and infiltration are critical for runoff generation and storage. Light soil textures generate less runoff and are also less suitable for runoff collection. Therefore, gently sloping areas with heavy soil textures (to prevent water infiltration) are potential locations for runoff harvesting [23,27,63]. However, except for two minor elevations to the east and west of the area, the remaining part of the park is flat and does not impose any limitation on runoff harvesting. As for soil texture, heavy and very heavy soil textures, which are more likely to produce runoff, are distributed as small packets nearly evenly across the whole park with a greater density in the northern and southern sections. Drainage density mainly governs time of concentration and high drainage density reduces the chance for runoff harvesting. As also argued by Singh et al. [64] and Jha et al. [31], areas with lower drainage network density and shorter time of concentration are more preferable for runoff harvesting. The park has a third-order drainage network with its majority being first and second-order streamlines, indicating park’s suitable condition in this regard.

Evaporation plays a major role in water resources management in arid and semi-arid areas. Evaporation is rarely measured in developing countries including Iran and those stations measuring this parameter are not well distributed in the country [65]. There have been efforts though to use soft computing and machine learning techniques for forecasting evaporation for data-scarce regions with acceptable accuracy such such as Ghorbani et al. [66] who used artificial neural network for forecasting pan evaporation in northern Iran, however, data scarcity is still one of the main challenges in the way of using evaporation as a site suitability factor for runoff harvesting. Evaporation causes water depletion and deterioration of water quality by increasing the concentration of soluble solids in water. Low relative moisture and strong winds in the study increase the intensity of evaporation causing a potential evapotranspiraion of higher than 500 mm per month in some parts [67]. The evaporation rate is high during mid-spring to early autumn period and low during the rainy season in late autumn until early spring [67]. Runoff harvesting is mainly possible during the rainy season when evaporation rate is at its lowest, however, later on the collected water could be lost through evaporation and therefore methods to store or conserve water (such as storing water in tankers, or installing shades) should be pursued. As data on evaporation is not available for the site, and as it can only affect the storage of water and given that it’s at its lowest during the water harvesting period, therefore evaporation (and evpotranspiration in general given scarce vegetation cover in desert environments) was not included in our analysis.

Rainfall is the primary factor for runoff generation. However, as there is no weather station inside the study area, and given that precipitation in this part of the country is considered homogeneous and without considerable spatial variations [68], a correlation was established between precipitation from 12 stations around the park and their corresponding altitude. Since water harvesting structures design is dependent on the design rainfall and considering annual rainfall does not provide reliable estimates [32], the correlation between 33%-exceedance probability of maximum daily rainfall and altitude was evaluated and the obtained equation was then used to convert the digital elevation model of the site into a design rainfall layer. In order to measure the total amount of runoff, we used the SCS-CN method, which is a simple and reliable tool for the estimation of runoff especially in small to medium-sized watersheds [21,69]. The SCS-CN method is designed for single storm events and the 33%-exceedance probability of maximum daily rainfall was used to simulate a single event.

Three thematic layers of potential runoff, slope and drainage network density were standardized, weighted and combined using the fuzzy membership function, AHP, and WLC methods. As also noted by several authors [32,70,71,72,73,74,75] the integration of GIS and multicriteria decision making provides an objective and reliable tool for site suitability analysis. The wildlife territory and soil salinity maps were applied as the constraining factors for runoff harvesting. The main objective of this research was to provide wildlife with water to preserve the remaining population of Persian Cheetah from extinction. Therefore, the territory of Ovis orientalis and Gazella Dorcas, as the main preys of cheetah, was identified through field surveys and by updating the maps provided by Semnan Province Environmental Protection Head Office and areas outside their territory were masked out from final suitability map.

Quality of runoff is as important as its volume for water harvesting. Runoff water quality is highly influenced by its origin. Existence of chemical or biological contamination could severely deteriorate the quality of water. Water pollution not only lowers its appeal to the local wildlife but it could also increase fatality in case of high concentration of elements or the presence of pathogens. Rosenstock et al. [76] have carried out important research on the possible negative impacts of water development on wildlife by stating that water quality at these facilities could be deleterious to animal health, although they did not find significant chemical and biological water contamination in their case study in Sonoran Desert of the US. Much of the research regarding water quality in runoff harvesting up to now has been concerned with urban watersheds and green roofs [77,78,79], total collectible runoff and overall suitability of the site for runoff harvesting. Far too little attention has been paid to water quality in natural environments in site suitability analysis. Such expositions are unsatisfactory because of excluding the quality of collectible water mostly because of data limitation. Therefore, considering water quality at the origin could be one of the unique features and strengths of this research. The most important factor affecting runoff quality in Kavir National Park is salinity. As there are multiple patches of evaporative saline geological formations in the area (Figure 10), originated runoff from these patches might not be suitable for collection and hence they were discarded from the map for final suitability assessment.

The results of site suitability analysis finally showed that more than 62% of the total area is not potentially suitable for runoff harvesting, and this technique is only applicable to less than 40% of the park. Eventually, 11 catchments were identified by visual interpretation of the suitability map and total collectible runoff from an event with the 33%-exceedance probability of maximum daily runoff. The catchments were prioritized based on the total amount of collectible runoff [7]. We believe that small embankments and check dams are the best methods to collect surface runoff in our study area as they require simple construction material available at the site, light pieces of machinery and low maintenance. The storage capacity of the earth embankment in each catchment was estimated based on the upstream area of the catchment and runoff volume. Ovis orientalis and Gazella Dorcas, which were the target wildlife species in this research, each requires 1.4–6.9 liters per day [80] and 4.5–9 liters per day [81] of water, respectively. Water requirement of animals depends on numerous factors such as season, sex, age, lactation, having young animals, climate, topography, water security, vegetation composition among other things [82]. Given the total number of these two species in the park and their water requirement (i.e. 4.15 liters per day for Gazella Dorcas and 6.7 liters per day for Gazella dorcas) there is a need for 6500 m³ water each year. In the best-case scenario and under the assumption of receiving five rain storm events a year which is the most logical assumption based on our knowledge of the site, only 257 m³ is collectable from all runoff harvesting structures, which is merely 4% of the total demanded water. The rest of the water requirement should be provided from available wells, most of which do not have suitable water quality. Another option is to mix saline water with the collected runoff. Embry et al. [83] believe that the range between 10 and 4000 ppm is the logical standard for water salinity control. Other methods such as rain guzzlers [12] or establishment of small scale solar ponds [24] could also be piloted in the park. Therefore, we believe that runoff harvesting technique should be pursued in combination with other water harvesting and conservation methods to provide drinking water for wildlife.

5. Conclusions

In this paper, we adopted a simple and yet robust method for identifying areas with adequate runoff volume and suitable condition for water collection using GIS and multi-criteria analysis. This study has shown that 38% of the area (846 km²) is suitable to highly suitable for runoff harvesting. Much of the Kavir National Park equalling to 62% (2623 km²) has a very low potential for runoff collection. Except for a narrow strip, extending from the western side to southern and then eastern side of the park, the rest of the area doesn’t have a suitable condition for runoff collection. Eleven catchments were recognized suitable to collect water, totally amounting to 52 m³ after each rainstorm. The amount of annual collected runoff volume only satisfies 4% of total water demand for Ovis orientalis and Gazella Dorcas. In general, we believe there is only a limited possibility for using runoff in Kavir National Park and other methods of water harvesting such as desalination of brine groundwaters should also be pursued. Studies like ours suffer from data limitation, especially involving climatic variables which require the installation of measuring devices in the park for future modeling and simulations. On the other hand, as wildlife species visit familiar places for watering, the way water development might affect wildlife and how the animals behave towards new water sources should also be studied. Possible sources of pollution should be identified and water quality should also be constantly monitored to prevent the spread of disease among the animals. At the same time that these new watering sources attract animals, they might also be good spots for hunters, which requires increasing patrols and establishing new guarding posts as well.