Modeling the Hydrological Impacts of Vegetation Cover Changes in the Upper Oum Er-Rbia Watershed (Morocco)

In Morocco, the mountainous areas are often exposed to bulky and vicious flows of water and sediment. This process is exacerbated by the decrease in vegetation cover and the disruption in rainfall-runoff conditions that frequently cause significant flooding. By exploring the main hydrologic elements of these processes, it is possible to understand the behavior and hydrological response of watersheds and thus plan accordingly. In this study, the authors focused on determining the morphometric characteristics of the upper Oum Er-Rbia River basin (UOERRB by assessing/ evaluating the land use and land cover changes for a period of 32 years (1984–2016). Hydrologic Engineering Center’s Hydrologic Modeling System (HEC-HMS) was applied to simulate four daily hydrological events. The concentration time was 7.7 hours. The four storm events examined to calibrate and validate the simulated outflow at the outlet indicated a good agreement between the hydrographs of the measured and simulated flows, with an average Nash-Sutcliffe efficiency (NSE) value ranging from 0.63 to 0.76. Between 2002 and 2016, an average 6.21 percent increase in vegetation cover of with annual rainfall increasing from 690 to 714.1 mm/year was observed. These results can contribute to a better understanding of the hydrologic processes and better estimation of the return flows and thus guiding management decisions and developments in the UOERRB.


INTRODUCTION
For arid and semi-arid countries, there is high anthropogenic pressure on watersheds due to land mis-use, over-exploitation of forests and unsustainable agricultural practices among other causes. This often results in rapid runoff, reducing infiltration and sometimes increasing soil erosion (Dregne, 2002;Ezeaku and Alaci, 2008;Pachauri et al., 2014;Yjjou et al., 2014). Loumagne et al. (1991) noted that land use is one of the key parameters influencing soil erosion in steeply sloping landscapes and often disrupts the rainfall-flow relationships. The soil characteristics and the distribution of surface runoff due to diverse management practices are key determinants of the significance soil loss in a particular watershed (Kosmas et al., 1997;Marchandise, 2007). The relationship between the vegetation cover and hydrological processes is an important factor in maintaining soil against erosion and mitigating the flood risk (Rey et al., 2004).
The Upper Oum Er-Rbia River Basin (UOER-RB), experienced significant changes in land-use and land cover since 1984 to 2016. A significant change in UOERRB climate is causing a disruption of rainfall intensities and amounts (Driouech, 2010;El Ghachi, 2018). The cause-effect relationships between soil types, land cover, land use patterns and hydrological behavior has been studied by several researchers for instance (Chaponnière, 2005; Carlson et Arthur, 2000; Laborde, 2000;Kingumbi, 2006;Rey et al., 2004).
Remote sensing and geographic information systems are among the recent modern tools used for studying complex environmental phenomena at different spatio-temporal scales. Landsat satellites provide the data that is required in mapping and analyzing land cover (Kumar et al., 2020;Obodai et al., 2019). On a global scale, high spatial resolution sensors such as TM, ETM+ and Oli from Landsat satellites have been widely used for land cover mapping of large watersheds with limited accuracy (Hansen and Loeland, 2012). For this study, HEC-HMS (Hydrologic Engineering Center-Hydrologic Modeling System) model was used to simulate the water runoff process in the upper watershed of Oum Er Rbia. In order to calculate the runoff volume, peak flow, unit hydrograph, SCS curves and routing methods were used. The runoff simulation is based on four (4) randomly chosen rainfall events from thirty years of data . Out of the four events, three were selected for model calibration and one for validation. Several criteria were used to judge the effectiveness of the simulated models, including: the coefficient of determination R 2 , ratio of the root mean square error to the standard deviation of measured data (RSR) and the Nash-Sutcliffe Efficiency (NSE). The main aim of this study was to explore the spatio-temporal land use and land cover (LULC) changes in UOERRB and its effect on hydrology for 1984 and 2015 based on remote sensing data. Using the HEC-HMS model to simulate event scenarios and predict basin response allows planning for the events that can result in destructive floods and damages at the watershed level.

Study area
The study area is part of the large watershed of the Oued Oum Er-Rbia which has its source at an altitude of 1800 m at 47 km from the city of Khénifra ( Figure 1). It is located between latitude 33° and 33°5' North and longitude 5°01'and 5°08' West. UOERRB covers an area of 1049 Km², with the outlet located at the Tahrat meteorological station. at the level of the upstream part on the left bank, UOERRB is hydrologically fed by several small tributaries born from the springs which gush out at the foot of the Atlas and some wadis the basins of which develop on small surfaces.
The upstream part of UOERRB is mountainous with altitudes ranging between 864 m and 2400 m, characterized with a diversity of reliefs and geological formations. The main substrates encountered in the basin are Paleozoic schists and quartzites, Triassic basalts and red clays, Lias, limestones and dolomites as well as Quaternary alluvium and colluvium (Bouabdelli and Piqué, 1996). The average annual runoff of UOERRB at upstream and downstream stations is 0.9 m 3 /s (Tamchachat) and 21 m 3 /s (Tahrat), with a moderated weak correlation between the two stations.

Methodology
The morphometric characteristics of the basin were calculated based on 30 m resolution ASTER-DEM (Advanced Space borne Thermal Emission and Reflection Radiometer-Digital Elevation Model) data, and topographic maps of Khenifra, Mrirt, Krrouchen, Timahdite, El Hammam and Itzer at scale 1/50000. These datasets were used to delineate the watershed and to calculate the map of slopes, altitude ranges, KC index, time of concentration, etc.
The nature of the soil influences the rise of floods, the infiltration rate, the retention capacity and the initial losses. The UOERRB soils were classified based on the SCS classification method that allocates soils into four hydrological groups (A, B, C, D) based on the estimates of their infiltration potential (Arlen D, Feldman). The existing soil map was completed in the eastern part, and subsequently each unit was assigned to soil type classes.
The hydrometric data were analyzed using the R software. The categories of flows exceeding 25 m 3 /s and their corresponding rainfall were identified and correlations between these different events were sought.
The rainfall-runoff relationship was simulated using the HEC-HMS model taking into account the SCS (Soil Conservation Service) production function, Curve Number Loss and the SCS (Unitary Hydrograph) transfer function. This method was chosen due to its robustness in reproducing the rainfall-runoff relationship in several contexts, especially in semi-arid areas, besides limited rainfall stations in the UOERRB.
Hydrological and morphometric parameters, such as surface, slope, sub-watershed boundary and hydrographic network, were calculated using ArcGIS 10.3 software and its extensions Arc Hydro and HEC-GeoHMS. Remote sensing image processing, such as geometric and atmospheric corrections, NDVI calculation, supervised image classification was performed by means of ENVI 5.1 (Environment for Visualizing Images) software. The general methodology adopted for this study is shown in (Figure 2).
The method used for atmospheric correction is Dark Object Subtraction (DOS) a simple empirical atmospheric correction method which  assumes that the reflectance of dark objects includes a significant component of atmospheric scattering (Gilmore et al., 2015). Support Vector Machine (SVM) classification was used to distinguish the five land cover classes (bare soil, water surface, dense forest, less dense forest, and sparse forest) of the study area.

Physical characteristics of the basin
The component features of UOERRB were calculated through morphometric analysis (Table 1). The digital terrain model (DTM) shows that the altitudes of the upper Oum Er Rbia basin vary between 800 m and 2400 m in the South-Western and North-Eastern part of the basin respectively. The UOERRB is characterized by a predominance of low and medium slopes, with slope classes below 0.2 covering approximately 82% of the total watershed area. The compactness index (KG) is 2.2, reflecting an elongated basin shape; therefore, the average water concentration time is 7.71 hours. The difference in elevation is 1540 m, which facilitates the circulation of water within the basin.

Analysis of the rainfall-flow relationship
It is necessary to have reliable and precise information on extreme rainfall events in terms of flood risk in the watershed for better management of water resources. The analysis of the rainfall and hydrometric series was carried out with the aim of qualifying the homogeneity of rainfall and flow on a spatio-temporal scale. The two stations used are Tamchachat (upstream station) and Tahrat (downstream station). The hydro-meteorological data are provided by the Agency of the Hydraulic Basin of Oum Er Rbia (ABHOR).
The wettest months at both stations are between November and February. The driest months are June, July and August. Comparing the precipitation in these two locations, it was found that both stations have the same precipitation pattern and a significant difference in flow rates between upstream and downstream; this is explained by the fact that the Tamchachat station is located on a stream upstream with a low flow. The average interannual flow at the Tahrat station is 21m 3 /s, while that of the Tamchachat station does not exceed 0.9 m 3 /s ( Figure 3).
In order to assess the quality of precipitation in the UOERRB on a spatio-temporal scale between the downstream (Tamchachat) and upstream (Tahrat) stations, daily rainfall was analyzed. The correlation between daily rainfall at Tahrat and Tamchachat is moderate ( Figure 4) with a correlation coefficient of 0.67. In addition, Tamchachat recorded more rainfall than Tahrat due to the differences in altitude.

Type of soil
In order to determine the amount of surface runoff in a watershed using the SCS-CN method,  The soil types were classified into hydrological soil groups (A, B, C and D) indicating the infiltration rate for bare soil after prolonged wetting ( Figure 5). Soil groups A and D cover the largest surface area of 476.93 km² and 241.84 km², respectively. Soil groups B and C cover 95.7 Km² and 234.77 Km², respectively. The corresponding runoff potential is low, moderate, medium and high for soil groups A, B, C, and D, respectively.

Land use mapping
The LULC in UOERRB was carried out by classifying the mapping of Landsat images for the selected dates. Figures 6 and 7    The generated LULC maps represent five classes, namely bare soil, water surfaces, dense, less dense, and sparse forest. Figure 8 represents an example of land cover mapping from the July 23, 2015 Landsat8 image. The dominant LULC class category is bare soil, occupying 40% of the basin area with very little urban fabric in the form of dispersed villages. The different classifications were evaluated by the confusion matrix and through the calculation of the Kappa coefficient as well as the overall classification accuracy which is between 70% and 80%. Table 2 represents the confusion matrix calculated for the classification of the Landsat7 ETM image of 23/07/2015. The Kappa coefficient is 0.75, and thus the classification is acceptable.
According to this matrix, there is some confusion, for example, between less dense forest and dense forest as well as between less dense forest and sparse forest.   There is also the following empirical relation Ia = 0.2 S This gives us: Now the maximum retention potential S and the characteristics of the basin are related by the number of curve CN by: The Curve Number (CN) is calculated on the basis of the hydrological soil group and land use. Curve Number values range from 30 to 100. The higher the CN value, the greater the surface runoff.
The weighted CN value over the entire basin was calculated using the following equation: where: CN w : Weighted Curve Number Ai: surface of the sub-basin Aj: total surface of the basin The SCS-CN method takes into account the initial soil moisture conditions that influence the soil capacity for runoff (SCS, 1972). Curve Number values were assigned to each soil unit corresponding to the intersection of land use and soil hydrological classes using the standard look-up table NRCS TRR55.
A great variation of CN values can be noticed in UOERRB (Figure 9), with low values (30-62) for permeable soils where dense forest is frequently found, and high values (64 and 90) for less impermeable soils; on the other hand, for water surfaces the maximum value is 100.

Modeling by HEC-HMS
After the calculation of all the parameters necessary for the modeling (slopes, CN, concentration time, Lag time, etc.), the schema created in HEC-GeoHMS that divided the basin into four sub-basins according to their degree of heterogeneity were exported ( Figure 10).
In order to determine the values of the parameters that give a good simulation performance of the observed flood hydrographs, parameterization of the modules via the HEC-HMS software was performed by calculating the input parameters of each sub-basin (Table 3).

Model calibration
The model calibration was conducted using the objective function on the peak flow. In order to determine the sensitivity of the parameters, CN, the Lag time, the impervious surface and the initial abstraction (Ia) were optimized by manually estimating the initial value of each parameter and change until obtaining a simulated flow close to the observed flow.
The four flood events were selected basing on the availability of datasets as presented in table 4 (three events for calibration and one for validation).

Event 1
For December 2000, the simulated flows correspond to the observed model, presenting well the general appearance of the observed flood. Indeed, the simulated peak flow is very close to the observed one in terms of peak and slightly shifted in terms of time ( Figure 11).

Event 2
For December 01, 2001 to December 30, 2001, it can be noticed that the general flood pattern is respected; the peak flow is well simulated in terms of time with an overestimation of the flow for the second peak ( Figure 12).

Event 3
For the third event, it can be noted that the simulated flow rate corresponds with the observed flow rate, the peak flow rate is well simulated in terms of time and its value is very close to the observed value with minimal spikes, with an over estimation of the flow rate for the second peak, and an underestimation for the first and the third peak ( Figure 13).

Event 4
The validation phase of the model aimed at ensuring that the flood hydrographs are reasonable and close enough to reality. It is a question of finding a peak flow in an independent way from   the optimized calibration parameters. During validation in the 4 th event (Figure 14), the simulated flows match with the observed flows except for some fluctuation at the beginning with an underestimation of the peak flood.
It is generally noticed that the validation parameters tend towards the input parameters of the initial model, except for the Lagtime and the initial losses, there is a significant increase at the level of all the sub-basins, while for the CN, the values are slightly modified and tightened between 73 and 80. For imperviousness, there is a decrease of 22% at BV1, while there is a slight increase for the other sub-basins ( Table 5).

Evaluation of the quality of the model
Once the model is established, performance assessment is essential to determine its relevance and reliability. The evaluation is done based on standardized performance indicators/ parameters (i.e., calculate and judge against a reference value) ( Table 6).
• The Nash-Sutcliffe criterion is a performance indicator, with values ranging between -∞ and 1]. Nash-Sutcliffe and RSR provide the information on the total modeling error. From Table 7, it can be noticed that the model results are good for events 2 and 4, and satisfactory for events 1 and 3. The stronger the correlation, the better the simulation and this is explained by the correlation coefficient R². In general, this parameter is satisfactory for events 1 and 4 with values of 0.711 and 0.752, respectively, and fairly good for events 2 and 3 with values of 0.806 and 0.773, respectively ( Figure 15). Figure 16 illustrates the relationship between peak flood flows and the amount of vegetation cover determined from Landsat images. It was    22.29 noted that the influence of the vegetation variation on the peak flow is clear in most cases. For the years following 2002, the peak flows of the floods are low with a more developed vegetation cover. On the other hand, the floods of the years before 2002 are characterized by higher peak flows and a weak vegetation cover. It can be argued that when the latter is dense, part of the precipitated water is absorbed and therefore the quantity of water that will contribute to direct runoff will be reduced. However, vegetation plays an important role in attenuating flooding. When the vegetation is developed the runoff is delayed and the flood peak is attenuated. Moreover, the flow being longer, the share of infiltration water increases and the volume of flood decreases.

CONCLUSIONS
The use of modeling as a tool to find the link between land use and hydrological process is a better method recognized as both simple and productive. From a digital terrain model, Landsat  The present study focused on the use of a rainfall-runoff model to show the effect of the change of land-use on the hydrological functioning of the Upper Oum Er Rbia basin and to predict the peak flows in case of flood, especially upstream of the city Khénifra (Tahrat Station).
The modeling step was preceded by the elaboration of land cover and soil type maps according to NRCS rules. The classification of Landsat satellite images allowed the identification of five land cover classes (bare soil, water surface, dense forest, less dense and sparse). The soil data for the study area was transformed from a texture-based terminology to a hydrological HGS terminology, which was then used to create the Curve Number (CN) map using the HEC-GeoHMS model. The final Curve Number map shows an average CN of 69; this means that the basin has moderately high runoff. The set of optimized parameters to validate the estimated peak flows at the outlet of the upper Oum Er Rbia basin was deduced based on 4 events.
The results of the land cover mapping show a clear increase in vegetation cover of 6.21%. The application of the rainfall-runoff model in the basin shows a positive correlation between observed and simulated flows. The model is satisfactorily calibrated with a NSE value of 0.79.