Integrating WaTEM/SEDEM model and GIS-based FAHP Method for Identifying Ecological Rainwater Harvesting Sites in Ziz upper watershed, SE Morocco

Moroccan southeast areas have limited water ressources, vulnerable to climate change and characterized by a significant spatio-temporal variability. In response, to ensure the availability of water for local comunity, it is advised to develop some alternatives that improve the local water resources management throughout these areas. Rainwater harvesting (RWH) has been used widely as an alternative technique towards water scaricity. However, taking into account socio-economic constraints for identifying ecological sites for RWH remains a complex task for water managers. The present study was conducted using WaTEM/SEDEM model, GIS techniques and Fuzzy Analytical Hierarchy Process (FAHP) method to identify important ecological RWH sites. For this purpose, several data sources were employed to generate needful thematic layers. The soil conservation service curve number (SCS-CN) method was utilized for preparing a yearly runoff potential map. Then, the thematic layers were weighted for generating RWH suitability map. The results show that yearly surface runoff ranges from 136 to 500 mm. Moreover, the spatial distribution map of soil erosion of WaTEM/SEDEM shows that Ziz upper watershed can be classed into four classes: (i) slight (very suitable), (ii) moderate (suitable), (iii) high (less suitable) and (iv) severe (not suitable). Aproximatly 76.1% of the study area falls within slight soil erosion class. The GIS tools were used for generating the important ecological RWH suitability map. The GIS’s intersect tool was used to eliminate insuitable soil erosion classes and the buffer tool was used for integrating the socio-economic factors including the distance from residential areas and agriculltural fields . Prediction accuracy of the resultant map of RWH suitability showed the value of the area under the curve (AUC) equal to 59.6% for FAHP method in this study. Ecological RWH map, mainly, introduced western areas and some central parts of Ziz upper watershed as suitable RWH areas. The present study demonstrated that coupling WaTEM/SEDEM model with FAHP method and GIS tools provide a valuable approach for identifying the ecological RWH sites in large semi-arid areas.

In the last decades, the Moroccan government has adopted the policy of building dams. The number of dams is the most important in North African countries. However, water scarcity has become increasingly alarming as dams are unable, largely due to loss of rainwater through evapotranspiration. During dry periods, water availability is limited due to low water storage capacities, even in areas with high rainfall and low evapotranspiration.
Moroccan south-eastern areas are constantly affected by water scarcity problems. The population especially the farmers are faced with major challenges due to the succession of dry and rainy periods leading to situations of climatic uncertainty and aridity. Thus, water stress in these areas can be alleviated by collecting runoff during rainy periods, and this harvested water can be worn for domestic uses and other activities during drought periods (Rejani et al. 2017).
According to Mekdaschi and Liniger 2013 RWH consists on collecting and managing rainwater runoff to enlarge water availability for domestic and agricultural use as well as ecosystem sustenance and other researchers reported that RWH is really precious, especially in dry areas, to address water shortage, mitigate groundwater extraction, expand crop yields, increase growth and development of pastures, fight erosion, and maintain water resources (Kaposztasova et al. 2014).
Nevertheless, many studies showed that identifying potential RWH sites is a significant stage towards water security and land yield in semi-arid areas (Isioye et al. 2012). According to previous studies, especially in semi-arid areas the most influencing factors in RWH site selection are hydrologic, climatic, agronomic, topographic, soil, and socioeconomic conditions (FAO 2003).
Consequently, several criteria were employed using different techniques and/or methodologies for siting suitable RWH sites are available (Mahmoud and Alazba 2014;Ahmad 2013;Al-Adamat 2008;de Winnaar et al. 2007). Also, in other studies the optimal RWH sites were identified by taking into account many factors including socio-economical and physical aspects of a studied area (Al-Adamat et al. 2010).
Besides the key parameters for siting potential RWH sites, runoff is one the most significant (Rejani et al. 2017;Adham et al. 2016). For simulating rainfall and runoff in the present study, The Soil Conservation Service curve number (SCS-CN) method was used. The SCS-CN method is one of the most used due to its simplicity and little data requirement (de Winnaar et al. 2007). In many studies, the spatial estimates of runoff from the SCS-CN model directly have been adopted for implementing the potential sites for RWH precisely (Rejani et al. 2017;Ajaykumar et al. 2012;Pandey and Sahu 2002;Napoli et al. 2014).
Nevertheless, in areas with high surface runoff, the soil productivity can be decreased by the soil loss and its nutrients, and can then conduct to the decline of water quality by deposing soil particles in watercourses and reservoirs. Consequently, in the riskiest areas of soil erosion, runoff's reduction can increase the ecosystem's capacity to store rainwater instead of releasing surface runoff.
Nowadays, Geographic Information System (GIS) plays an important role for identifying suitable RWH sites, through which several geospatial data can be displayed, processed and analyzed. Consequently, GISbased Multi-Criteria Decision Analysis (MCDA) techniques are very important and valuable approaches to manage enormous spatial and non-spatial data into intended knowledge that along with subjective judgments of decision makers in some critical decisions (Chen et al. 2010).
We therefore assessed the spatial distribution of soil erosion by the WaTEM / SEDEM model and we used geographic information system software (ArcGIS 10.5) to analyze, process and overlay the spatial data and their attributes for slope, LULC, soil type and the soil erosion map to identify optimal sites for harvesting water.
Ziz upper watershed was selected as a typical semi-arid region in southeastern Morocco, to check out potential RWH sites. This area has had periodical water shortages due to dried seasons, particularly with the increase in population and the increasing demand for water from farmers. Then, water neediness has faltered as a major factor restricting local sustainable development, and big amounts of rainwater drain straight without being used or stored.

Study area
Ziz upper watershed is one of the Moroccan southeast regions which is bordered by upper Moulouya on the north, Ghris watershed and Tafilalet plain on the south, Central High Atlas on the west and Guir basin on the east. Administratively, this area consists of one urban commune and some rural communes. The administrative center of this area is Rich center located between two cities Errachidia and Midelt. Ziz upper watershed is bounded by latitude 32 ° 05' 48 '' -32 ° 64 '19'' N and longitude 04 ° 11' 72 ''-05 ° 46'20'' W, covering a total area of 4435 km2 (Fig. 1). The study area's altitudes range from 1023 to 3386 m asl. and it experiences a semi-arid climate, characterized by a harsh winter and a moderate summer with hot temperatures varying depending on the altitude. The average minimum and maximum temperature of Ziz upper watershed are 10.2 °C and 19.2 °C, respectively. The rainfall is characterized by the existence of two rainy seasons, autumn and spring, winter is less rainy and summer is dry. Average annual rainfall ranges from 119 to 377 mm yr -1 (Manaouch et al. 2020). Land use/land cover (LULCs) types of the study area are classified into five major classes such as agricultural fields, very degraded forest, rangelands, settlement and water bodies with rangelands and degraded forest types have higher proportions.

Methodology, data processing and thematic layers
For conducting the present study, the FAHP method and the spatial distribution map of soil erosion risk were integrated in GIS environment. The first stage was the preparation of thematic layers including five criteria such as runoff, slope, hydrologic soil group, drainage density and LULC. For estimating the potential runoff map, SCS-CN model was used. For weighting and integrating all the thematic layers, the AHP method and GIS tools were employed by comparing each factor to each other's as a couple according to the expert's judgment. Then, all these layers were transformed into raster maps using fuzzy set theory. For generating the ecological RWH suitability map, the potential RWH map was overlaid with the spatial distribution of soil erosion produced by WaTEM/SEDEM model (Fig. 2). The last stage was the integration of socio-economic factors. Then, the resulted map was validated using a group of current RWH strictures.

Figure 2:
Flow chart showing the adopted methodology for ecological rainwater harvesting sites.

Climatic data.
In the present study, Climatic data were derived from the hydrologic watershed agency of Ziz (ABH-GZR 2020). Rainfall data from four meteorological stations in or near the study area were collected and transformed from monthly averages to yearly averages, and data for surface runoff were organized in the same way. Geographic coordinates for the meteorological stations are displayed in Table 1. Accordingly, The rainfall map was developed by using the interpolation method IDW (Inverse Distance Weighted), which allow to estimate the rainfall value of any point. Precipitation is between 128 and 443 mm. year-1. The southern and eastern parts are the least watered. However, the western parts are experiencing high values. The wettest areas are in the western parts, representing the mountains upstream borders of the studied area (Fig. 3A).

Slope
Slope is one main factor that influence the harvesting water. For this study, the slope data were extracted from the SRTM-DEM with a resolution of 30 m. A slope map was produced using the slope tool of ArcGIS 10.5. Gentle slope's areas (0-5°) in the Ziz upper watershed represent 26.8% of the total study area, areas with slopes between 5-15 ° represent 36.4% and areas with between 15 and 25 ° then greater than 25 ° respectively represent 20.3 and 16.8% (Fig. 3B). About 62% of the study area had slopes <15 °.

Land use/ land cover (LULC)
The LULC of a watershed has a major effect on runoff and is also an important criterion for identifying suitable sites for rainwater harvesting. LULC data were extracted by processing a satellite imagery of Landsat 8 OLI acquired in February 2021 (Mohamed et al. 2020). LULC map was classed into 4 categories based on the supervised classification in ArcGIS 10.5 (Fig. 3F). The main part of the area included rangeland followed by degraded forest (Table 3). These areas were inspected as suitable sites for potential RWH.

Soil texture and Hydrological Soil Groups (HSGs)
Soil properties have an important impact on surface runoff by influencing infiltration rate. Data for soil texture were acquired from the prior studies, which is compiled by Manaouch et al. 2020. Soil group was reclassified from the United States Department of Agriculture -Hydrologic Soil Group (USDA 2009) (USDA-HSG) (Fig. 3C), represented by categories A to D, for treatment in the SCS-CN model. HSG-A is for high permeability with low runoff potential, accounts 0% in the Ziz upper watershed, HSG-D represents soil with high runoff potential, accounts 8.5% of the study area and HSG-B and -C are middle classes, account 73 and 18.5%, respectively.

Drainage density
The drainage density parameter is computed using the total stream length per unit catchment area (Horton 1945). It is calculated as: Where Dd is the drainage density, n is the number of streams, l is the stream length (km), and A is the drainage area (km 2 ).
The drainage density was generated from DEM using the ArcGIS 10.5 software. The drainage density map was done based on the length of streams presents in each sub-watershed. An area with high Dd was considered more suitable for RWH structures compared to other areas with low Dd (Fig. 3E). (Prinz et al. 1998)

Runoff map by SCS-CN model
The SCS-CN model, developed by the US Soil Conservation Service in 1972, has been extensively used for assessing runoff and managing water resources due to its simplicity (SCS 1972;He 2003). In this method, some properties such as soil characteristics, LULC, topography and antecedent soil moisture conditions, are taken into consideration (Luxon and Pius 2013). Recently, the SCS-CN method has been developed for application to larger areas by weighting CNs based on the LULC of the studied area (Ramakrishnan et al. 2008). The SCS-CN model was designed for any form of LULC and soil shape, and the use of, and the use of GIS tools have facilitated the runoff prediction more accurately. The runoff estimation using ArcGIS tools is described in details on the methodology flowchart above. The SCS-CN model calculates: Where Q is the runoff depth (mm), P is the total rainfall (mm), S is the potential maximum retention after the beginning of runoff (mm), λ is surface runoff abstraction (dimensionless) and CN is the curve number that ranges from 0 to 100, reflecting the surface runoff.
A standard condition of the SCS value considers the parameter λ to be equal to 0.2 (SCS 1985) and for CN, the greater values means that a lot of the rainfall is transformed into surface runoff and vice versa.

Fuzzy set theory and fuzzy maps
In the present study, the well-known Zadeh's fuzzy logic theory has been adopted for giving to each parameter its value. Then, each selected parameter has been converted into a fuzzy raster map by using the linear membership function (LMF) tool of ArcGIS 10.5. In the fuzzy set theory, the member's value 0 means that this member is not a party of a set, whereas the memder's value 1 means that this element is a full member of a set (Zadeh 1965). An illustration of fuzzy set is given in the following equation (Mcbratney and Odeh 1997).

= {x, μA(x) for each x}
Where, μA is the MF (membership of x in fuzzy set A) so that: If x is not a member of A then μA = 0. If x is a full member of A then μA = 1. If x belongs in a certain level to A then 0 < ( ) < 1 As introduced in the equation above, LMF was used for runoff, hydrologic soil group, drainage density and LULC (Feizizadeh and Blaschke 2013).
Where x is the considering parameter and a, b are the minimum and maximum values, respectively.
For assessing the slope's effect on RWH suitability, the following LMF was used (Feizizadeh and Blaschke 2013) Where x is the slope's value and a, b are the maximum and minimum values, respectively.
In this work, five criteria for RWH suitability site were used: runoff, slope, drainage density, hydrologic soil group and LULC. Each pixel in Fuzzy map of a factor has its value, and it ranges between the minimum 0 and the maximum 1.

Analytic Hierarchy Process (AHP) and criteria's weights
AHP is a well-known method in multi-criteria decision analysis. In this work, AHP and GIS tools were employed for weighting and integrating the thematic layers. The five factors including runoff, slope, HSG, drainage density and LULC were used in this study. After preparing these layers in the GIS environment, each of them is prioritized to one another layer using expert's judgment. Then, to determine the relative value of the factors, each factor is compared to each other's as a couple and computing the pairwise comparison consistency index until obtaining acceptable results. The following equation shows the consistency ratio computation in AHP method (Saaty 1990).

= − − 1
Where CR is consistency ratio, is the average value of consistency vector, and n is the number of criteria or factors. Based on this relationship, if CR is less than 0.1, acceptable level of conflict has been paired comparisons, and otherwise, the rates that reflect the judgment are inconsistent (Saaty 1990). The CR indicates that the cumulative judgments derived from the pairwise matrix are satisfactory and consistent.

WaTEM/SEDEM is a spatially distributed soil erosion and sediment delivery model that combines the Water and Tillage Erosion Model (WaTEM) and Sediment Delivery Model (SEDEM). It is based on Revised Universal Soil Loss Equation (RUSLE) and a Transport Capacity equation (TC) to predict sediment delivery in drainage network
Error! Reference source not found. Verstraeten et al. 2002). WaTEM/SEDEM comprises three modules, the first calculates annual soil loss using the RUSLE model (Renard et al. 1991).

= 2
Where A is the mean annual soil loss (kg m −2 yr−1 ), R is a rainfall erosivity factor (MJ mm m-2 h-1year-1), K is a soil erodibility factor (kg h MJ-1 mm-1), C is a dimensionless cover and management factor, P is a dimensionless erosion control practice factor, LS2D is a slope-length factor (Desmet and Govers 1996) The second estimates the sediment flux from slopes to the stream network (Van Rompaey et al 2001).

=
( 2 − 4,1 0,8 ) where TC is the annual transport capacity (kg m-2 year-1) and R, K, and LS factors are the same as in RUSLE's equation, s is the local slope (m m-1) and ktc is a transport capacity parameter (m).
The last module is used for assessing the tillage erosion that represent the net soil flux caused by tillage on a hillslope of infinitesimal length and unit width is proportional to the local slope gradient (Govers et al. 1994): Where Qs, t is the net downslope flux due to tillage translocation (kg m-1 yr-1), ktill is the tillage transport coefficient (kg m-1 yr-1), s is the local slope gradient (m m-1), h is the height at a given point of the hillslope (m) and x is the horizontal distance (m).
A spatial distribution map of soil erosion was prepared based on the previous work (Mohamed et al. 2020). For more details on the WaTEM/SEDEM application and validation, we refer to the work of Manaouch et al. 2020 in Ziz upper watershed.

Integration of socio-economic factors in RWH suitability map
The success of any project is related not only to the technical aspects, but also to how can improve the local community's situation. Many studies illustrated that one of the primary reasons we do not use RWH enough in semi-arid areas is the lack of data on socio-economic contexts. In the present work based on the previous studies and expert judgment's, we are taking into account the distance to the residential farmers and their agricultural fields such as socio-economic criteria affecting the suitability of the installation of RWH sites. Since farmers and rural residents are targeted first, proximity to RWH is an influential parameter in the selection and assessment of suitability of RWH. The distance between residential's home and RWH strictures should be as short as possible. As rural farmers have scattered agricultural fields and each field occupies less than 1 ha and distributed at a distance of approximately 5 to 1000 m from their residence. Therefore, it is very important that the closer the RWH sites are to the fields, the easier and cheaper the RWH control and surveillance operations will be, especially in mountainous areas where transport is not easy. The distance between sites and fields was assessed using Google Earth images and GIS's buffer tool of ArcGIS 10.5 software.

Runoff potential
The Runoff parameter is the most influencing factor for identifying suitable RWH sites. According to the estimated average annual runoffs for different HSG and LULC classes provided in Tables 6 and 7. The yearly runoff average was about 248.5 mm with a standard deviation of 54.7 mm. The maximum runoff and average annual runoff potential in the western parts were 342.30 and 277.65 mm, respectively. The eastern and southern parts of the Ziz upper watershed are on the periphery of the Tafilalet oasis and the Guir's arid watershed, where the potential runoff was the lowest and its average annual was 136 mm. The low runoff yield in these areas was primarily due to lower rainfall. The average annual runoff was also divided into four classes from low to high (Fig. 3D). The very suitable zones of potential runoff accounts for 7.3 % of the study area.

Spatial distribution of soil erosion risk
The risk of soil erosion is narrowly linked to precipitation, soil properties, topography, LULC and human interventions. The transport of soil particles from one place to another can have several impacts, on site such as the decrease of soil productivity. While, the deposition of eroded soil particles can cause imbalances in waterway ecosystems such as deposits in streams and reservoirs which can reduce the quality of their water. The average annual soil loss (A) in the upper Ziz watershed was estimated by exploiting WaTEM/SEDEM model (Fig. 4F). To assess the potential severity of soil loss, the spatial distribution map of soil erosion has been classified into four classes of erosion intensity (Table 10).
According to the classified map produced by WaTEM/SEDEM model, about 76.1% of the area had a slight soil erosion class, only 0.2% had a severe soil erosion class and the rest of the area had a moderate to high risk of erosion. The spatial distribution of soil erosion shows that western and northern mountainous areas in the Ziz upper watershed receive an important rainfall amount and are largely occupied by steep slopes, so they have a quite high soil erosion risk. Whilst southern plains and eastern parts have a rather low rainfall and thus a low soil erosion risk.

Potential RWH suitability map
Harvesting rainwater involves collecting, storing and managing rainwater for avoiding water scarcity, and for the maintenance of ecosystems (Studer et al. 2013). In the present study, water deficiency in the study area can be resolute sustainably by improving the collection of rainwater runoff and taking into account soil and water conservation.
According to the two maps of potential runoff and spatial distribution of soil erosion, the study area's ecosystems are largely dominated by very low to feeble runoff and by slight to moderate risk of soil erosion. The potential map of RWH suitability was categorized into five categories i-v: C-i represents high runoff and slight erosion class (extremely suitable), C-ii represents moderate runoff and slight erosion class (very suitable), C-iii represents low runoff and slight erosion class (suitable), C-iv represents very low runoff and slight erosion risk (less suitable) and C-v as unsuitable sites for RWH.
The extremely suitable RWH sites covered an area of 4.3%, while the residual of the study area was classified as suitable and unsuitable for harvesting rainwater. Optimal C-i covered a small area (19072 ha) and sub-optimal Group ii and iii RWH sites covered a larger area, accounting for 63% of the study area.
In the study area, rangeland is the dominant LULC category, so a high portion of the RWH area was rangeland, and a lower part was agricultural field (Table 11). The western part of the Ziz upper watershed had a large RWH area, the southern and eastern parts of the study area belonging to Tafilalet plain and desert, respectively, had a smaller RWH area and the northern area had small and scattered RWH areas (Fig.  6).

Fig. 6
Potential RWH suitability map (left) and classified potential RWH suitability map (right).

Integration of soil erosion and socio-economic factors
The result of integrating the socio-economic criteria layers of distance to settlement and distance to agricultural field shows that after removing these constraints from the total suitability area, the net area of RWH suitability map decrease from 443,532 ha to 147,864 ha. The excluded areas have 0 value (white color) which means not suitable area and the rest of areas are classed for RWH (Fig. 7). After the constraint's removal from the classified RWH suitability map, the net area coverage for each class has decreased.

Validation of the suitability of identified potential RWH sites
In this study, the prediction accuracy's of FAHP method for identifying RWH sites was evaluated using the receiver operating characteristics (ROC) by comparing the final RWH suitability map with the scattered existing RWH sites in the study area. According to the field visits, 30 current sites were found for harvesting rainwater. They are all earthen dams. These 30 RWH sites were then used to judge the predicted results accuracy of the FAHP model. Next, each pixel value in the resultant RWH map was classified into four groups as unsuitable, less suitable, suitable, more suitable and extremely suitable. Then, using the ArcSDTM tool in ArcGIS 10.5 software, the AUC of the ROC curve was calculated. This AUC represents the accuracy of the model to predict suitable RWH sites (Chen et al. 2016). A good model prediction has AUC values between 0.5 and 1, while values less than 0.5 represent a random fit. The AUC in our study is 0.596, which corresponds to a prediction accuracy of 59.6%. This acceptable value of AUC shows that FAHP method was a useful tool for predicting suitable RWH sites in this case study.

Fig. 8
Classified RWH suitability map with current RWH sites.

Fig. 9
Prediction rate curve for the potential RWH suitability map produced by FAHP method

Conclusion
Arid and semi-arid areas around the world are characterized by dry seasons. For ensuring the water availability, rainwater harvesting is an alternative source. In the present study, Ziz upper watershed was selected as a representative semi-arid region in southeastern Morocco. For identifying ecological suitable sites for harvesting rainwater, the methodology adopted integrates three components: a potential runoff map based on the SCS-CN model with a spatial distribution map of soil erosion based on the WaTEM/SEDEM model and a constraints layers of socio-economic factors. The studied area had very low to low runoff potential and moderate to slight risk of erosion.
The results showed that 6.9% of the study area was extremely suitable for potential RWH. Optimal RWH sites were scattered in the western parts of the study area. The predicted results verification shows that the accuracy of the method for predicting suitable sites for harvesting rainwater was 59.6% by comparing the existing RWH sites with the classified map of potential RWH sites. The methodology used in this work demonstrates the utility of integrating a spatially distributed model for soil erosion and GIS based FAHP approach for identifying ecological rainwater harvesting areas in large semi-arid areas. This study has provided a roadmap for siting RWH strictures and alleviating water shortage problems in Ziz upper watershed.