Application of the HEC-HMS hydrological model in the Beht watershed (Morocco)

This work focused on the collection and preparation of the data required for the hydrological modelling of the Beht catchment area, which covers an area of 4560 km with a perimeter of 414 km, by combining the various spatial technologies, in particular geographical information systems (GIS), remote sensing, and digital terrain models (DTM), with hydrological models in order to prepare for spatial hydrological modelling used for flood forecasting. The methodology consists, at first, in the automatic extraction of the sub-basins and the drainage network. Then, edit these data using the HEC-GEO-HMS extension, and the preparation of the land use and land cover data for the elaboration of a Curve Number (CN) map of Beht watershed, then the import of the basin model into the Hydrologic Modeling System (HEC-HMS) to simulate the surface runoff using six extreme daily time series events.


Introduction
Over the last years, Morocco has experienced a number of tragic flood events that have generated flooding in several regions of the country due, on the one hand, to population growth and urban, agricultural, industrial and tourism development which lead to an increasing occupation of vulnerable areas and, on the other hand, to the aggravation of extreme conditions (drought and floods) as a result of climate change [1]. To deal with this flood risk, a set of tools has been developed to understand the hydrological functioning of basins. In this context, hydrological modeling is the most adequate tool to understand the water cycle on small and large scales.
The Soil Conservation Service (SCS) curve number (CN) method is one of the most popular methods for computing the runoff volume from a rainstorm [2]. It is popular because it is simple, easy to understand and apply, and stable, and accounts for most of the runoff producing watershed characteristics, such as soil type, land use, hydrologic condition, and antecedent moisture condition [3,4]. The SCS-CN method was originally developed for its use on small agricultural watersheds and has since been extended and applied to rural, forest and urban watersheds [5]. Due to its low input data requirements and its implementation within GIS, it has been incorporated in many widely used hydrological models. In recent years, the method has received much attention in the hydrologic literature. The SCS-CN method was first published in 1956 in Section-4 of the National Engineering Handbook of Soil Conservation Service (now called the Natural Resources Conservation Service), U. S. Department of Agriculture. The publication has since been revised several times. Despite several limitations of the method and even questionable credibility at times, it has been in continuous use for the simple reason that it works fairly well at the field level. [6,7,8,9,10,11] The Hydrologic Modeling System HEC-HMS, which is a hydrologic modeling software developed by the US Army Corps of Engineers Hydrologic Engineering Center (HEC) is an integrated modeling tool for all hydrologic processes of dendritic watershed systems. It consists of different component processes for rainfall loss, direct runoff, and routing. HEC-HMS has become very popular and been adopted in many hydrological studies because of its ability in the simulation of runoff both in short and longtime events, its simplicity to operate, and use of common methods [12]. Hydrographs developed by HEC-HMS either directly or in conjunction with other software's are used for studies of urban drainage, water availability, future urbanization impact, flow forecasting, flood damage reduction, floodplain regulation, and systems operation [13].
The objective of this study is to apply GIS software and remote sensing to determine Curve Number for Beht watershed to study a rainfall-runoff model based on the HEC-HMS, to calculate runoff volume and peak discharge.

Study area
The Beht catchment area is located in northwestern Morocco and covers an area of approximately 4560 km2 in the southwestern part of the Sebou basin. It is bounded to the north by the Gharb plains and the Meknes shelf, to the south by the Oum-Erbia basin, to the west by the Bouregregreg basin and to the east by the Middle Atlas. Its boundaries are located between the meridians 5° and 6° West and the parallels 33° and 34° North. This basin is located between the Lambert coordinates (X1 = 430347.24; Y1 = 281864.43) and (X2 = 529704.23; Y2 = 386110.82). The Beht watershed has an elongated shape following a SW-NE direction. The Gravelius index of compactness, calculated for this basin, is about 1.86. It is therefore eight times longer than it is large, which allows a rapid collection of water towards the outlet. It can also be assimilated to a rectangle of the same surface area, which is 202.5 km long and 22.5 km in width. The Beht catchment area has a Mediterranean climate (semi-arid to humid). It presents a double gradient of decreasing intensity from South to North and from East to West. This climate is marked by frequent summer droughts and violent stormy rainfall. Rainfall is marked by annual fluctuations. They vary from 550 mm in the North-West of the basin to about 900 mm in the South-East. Temperatures show a clear variation in space and time. High altitudes are characterized by low temperatures, ranging from -0.9°C in winter to 25°C in summer. While low altitude regions record temperatures of around 15°C in the winter and 34°C in the summer.
The hydrological regime is characterized on the one hand by floods recorded mainly during wet periods and which ensures 80% of the annual liquid flow. On the other hand, low water levels during the dry season, when liquid flows are very weak at around 20%.

Data processing
The work methodology focuses on the preparation of the data necessary for the spatial hydrological modeling of the basin, from Arc Hydro and HEC-GeoHMS; extensions of a geographic information system (GIS), as well as the elaboration of land use and soil maps and the calculation of the Curve Number grid, then the import of the basin model into HEC-HMS.

Delimitation of the watershed
The traditional method used to delineate a watershed area from the topographic map takes time and is inaccurate and it has been replaced by the automatic extraction from a digital Terrain Model [14].
So, the first step consists in the automatic delimitation of the Beht watershed based on the digital terrain model, derived from the ASTER sensor, which is characterized by its 30 m spatial resolution. Then, the delimitation of sub-basins and the extraction of the drainage network from the basin DTM [15]. Nine operations are carried out to get the schematization of the basin model [16].

Land use map
It is determined through a supervised classification on satellite images "ASTER" using an image processing tool (ENVI: Environment for Visualizing Images). There are six main types of land use in the Beht watershed: Pastures covering almost 1/3 of the area, which is equivalent to an area of 1471 km2. They are geographically dispersed throughout the basin. This natural vegetation develops according to the type of soil conditions and climate. Followed by bare land, which representing 23% of the total surface area, i.e. 1052 km2. They are mainly located upstream. Forests represent 20.6% of the surface area. They are grouped in two lots located respectively on the middle and the southern extremity of the basin. The agricultural lands represent 18% of the land; they are mainly located upstream of the basin. Matorrals appear in the extreme northwest of the study area representing 5.6% of the surface area. (Fig.2)

Fig. 2. Land use map
Because of the specific requirements of the chosen modular combination, specifically the NRCS CN (Natural Resources Conservation Service Curve Number) method as a production function, the elaboration of a land use map over the entire study area was a necessary step.
The information contained in this map should be authentic to the classification of the NRCS (Natural Resources Conservation Service) [17], so we had to make connections between the NRCS classes and the map prepared by the satellite image classification method. Then, we reclassified the thematic classes defined above as shown in the table below:

Soil map
The nature of the soil affects the rate of flood rise and volume, as well as, the infiltration rate, moisture content, storage capacity, initial losses, runoff coefficient (Cr) are all related to the soil type.  The soil cover for the entire watershed shows a significant dominance of poorly developed soils, which can form associations with crude mineral and calcimagnesic soils (33.1%). The poorly developed soils and the association "poorly developed soils and raw mineral soils" are located in the upstream part of the watershed; they are mainly associated with alluvial deposits. The association "poorly developed soils and calcimagnesic soils" is present mainly on the right bank downstream of the basin.
Approximately 30% of the watershed is covered by brown soils and associations of brown soils and hydromorphic or poorly developed soils. These soils are mainly present in the middle of the watershed, although they can also be found further south in the upstream part of the watershed.
In the northeast and southwestern extremities, vertisols and assimilated soils are present, with a percentage around 11%. Isohumic soils and the soil association isohumic and calcimagnesic soils downstream of the watershed, have a proportion of 11.2%.
The basin has other types of soil and soil association, such as: hydromorphic soils, fersiallitic soils, the association of raw mineral soils, poorly evolved soils and hydromorphic soils, the association of fersiallitic and poorly evolved soils, the association of brown soils and raw mineral soils etc..., but in small or even very small proportion.
The soil classification used by the Soil Conservation Service method is the hydrological classification. It is a classification that consists of grouping soils into four hydrological groups (A, B, C, D), [18], based on their estimated infiltration potential. As a result, soils are assigned to one of four groups based on their infiltration rate; A; soil having high infiltration rates, B; soils having moderate infiltration rates, C; soils having slow infiltration rates, and D; soils having very slow infiltration rates [17].
The transition from soil classification to hydrological classification is made by providing information on soil texture according to the composition of sand(S), silt (St), clay(C) and organic matter (O), because Soil texture information is essential to determine the runoff coefficient [19]. The values of these components are given in the following table: From this map it can be concluded that the most dominant class is class C, which shows that the soils have slow infiltration rates, therefore a relatively high runoff.

Calculation of the CN grid
The SCS has developed a soil characterization system based on the hydrology and land use group called the Curve Number (CN). Values range from 0 to 100; A CN value of 0 indicates no runoff potential, while a value of 100 indicates that all precipitation runs off [18]. In other words, a value of 100 is assigned directly to the water surface and 0 for highly permeable soils with a high infiltration potential.

Preparation of the CN-Lookup table
The look-up   The final map of curve number of Beht watershed (Fig.5) shows an average CN of 78 which means that the basin have a moderate high runoff, and this is due to the clay type soil dominated by poorly developed soils, brown soils and vertisols assimilated soils and also, the vegetation cover that is marked by significant presence of pasture lands. These results are nearly similar to the study realized by [20] who found that the average curve number in the Sebou basin is 82.
The choice of the Curve Number depends, in addition to the soil type and land use, on the antecedent soil moisture conditions (AMC). These can be dry (I), moderate (II) or wet (III) [21]. The values provided in the Attribute Table (tab.3) are representative of average initial moisture conditions (CNII) (Fig.5) and the Curve Number CNI and CNIII are calculated directly using the [22] equations below: Soil moisture status is determined based on precipitation in the watershed during the last five days before the event in question, and by season (low and rainy seasons).
The curve number in the conditions I and III are 51 and 89 respectively.

HMS model
The HEC-HMS model is physically based and conceptually semi-distributed model designed to simulate rainfall-runoff processes in a wide range of geographic areas, from large river basin to small urban and natural watershed runoffs. It makes it easy to perform huge tasks related to hydrological studies, including losses, runoff transform, open channel routing, weather data analysis, rainfall-runoff simulation and parameter estimation [23,24]. Moreover, the HEC-HMS uses separate models that compute runoff volume, models of direct runoff, and models of baseflow. It has nine different loss methods; some of it is designed for event simulations, whereas others are for continuous simulation. It also has seven different transformation methods, Clark Unit Hydrograph Banitt [25] methods that have been applied successfully to simulate long-term stream flows.
In this study, the Soil Conservation Service Curve Number loss method was selected to estimate direct runoff from a specific or design rainfall [26]. This method normally calculate the runoff volume by computing the volume of water that is intercepted infiltrated, stored, evaporated, or transpired and subtracting it from the precipitation. It was chosen because it is a simple method for the estimation of the direct runoff from a storm rainfall event, and it relies only on the curve number, which is a function of the soil type and land use/cover that are the major runoffproducing watershed characteristics and also it is commonly used in different environments and provides better results compared to initial and constant loss rate method [27]. The SCS-CN method is based on the fact that the accumulated rainfall-excess depends on the cumulative precipitation, soil type, land use and the previous moisture conditions as estimated in the following relationship [21]: = (P−Ia) 2 P−Ia+S (3) Where Pe is the accumulated precipitation excess at time t (mm); P is the accumulated rainfall depth at time t (mm); Ia is the initial abstraction (mm) = 0.2S; and S is the potential maximum retention (mm). The maximum retention, S, and watershed characteristics are related through an intermediate dimensionless parameter, the curve number (CN) as: Where CN is the SCS curve number used to represent the combined effects of the primary characteristics of the catchment area.
Regarding the transform method, the Soil Conservation Service Unit Hydrograph model was chosen to transform excess precipitation into runoff. The lag time (Tlag) is the only input for this method. It is the time from the center of mass of excess rainfall to the hydrograph peak and is calculated based on the time of concentration Tc, as: = 0.6 (5) Where Tlag and Tc are in minute.
The routing method selected is the Muskingum method, which was developed by McCarthy [28]. It allows calculating the outflow hydrograph at the downstream end of the channel reach from the inflow hydrograph at the upstream end. Two parameters are needed; travel time (K) of the flood wave through routing reach; and dimensionless weight (X) which corresponds to the attenuation of the flood wave as it moves through the reach. The routing parameters in the models are usually derived through calibration using measured discharge hydrographs [29].

Model Calibration and Validation
In order to identify the key parameters and parameter precision required for calibration, a sensitivity analysis is usually used in most modelling studies [30,31]. The sensitivity parameters were selected based on their effect on peak discharge and total volume. The model calibration was done with the Univariate Gradient optimization package and Peak-Weighted Root Mean Square Error (PWRMS) objective function because of their simplicity and performance [26]. In this study, six flood events (four for calibration and two for validation) in the period from 1996-2014 were selected for calibration and validation. Watershed parameters such as curve number, initial abstraction, Tlag and baseflow possibly will need modification to produce the best fit between simulations and observations. For validation, the simulated data must be computed and compared with the observed data. Statistical assessment criteria as relative bias error functions proposed by Najim [32], Nash-Sutcliffe Efficiency (NSE) by Nash and Sutcliffe [33], Ratio of standard deviation of observations to root mean square error (RSR) by Moriasi et al. [34] and coefficient of determination (R 2 ) as described in Neter et al. [35] were applied to evaluate the performance of the model and the selected loss and transform methods.

Calibration
The results of the hydrological model in this study showed a reasonable fit between the simulated and observed hydrographs after optimization; the shape of the hydrograph and the time of peak had a good match ( fig.6). In the majority of events, the hydrograph shape was accurately reproduced in the model output. However, the volume of runoff was overestimated in events 1, 2 and 3 and underestimated for the event 4. The modelling results of peak discharge, total volume, and their relative errors with respect to the observed data, the Nash Sutcliffe Efficiency, the coefficient of determination and ratio of standard deviation of the observation to the root mean square error (RSR) values during calibration are mentioned in Table below.    The calculated values of the percent error in total volume and peak flow between simulated and observed values in all simulations before optimization was approximately high, with mean values of 6.14% for the peak flow and 14.15% for the total volume. According to this result, a sensitivity analysis was conducted to determine the most sensitive parameter. It was found that the initial abstraction and curve number were more sensitive, travel time K less sensitive and lag time was insensitive.
After optimization, the values in this study were reduced to 0.1 to 6.1% for the peak flow and 0.86 to 17% for the total volume, with mean values of 2.5% and 4.8%, respectively ( Table 4). The result is very good according to Najim et al. [32] and Cheng et al. [36] who recommended that the acceptable ranges of relative percent errors between the observed and simulated values should be below ±20%. The results showed also a relatively close agreement between the observed and simulated peak flow values at the period of calibration (R 2 = 0.842) (Fig.7). Based on the classification mentioned in Zou et al. [37] the mean correlation coefficient obtained in this study can be considered as strong (>0.8). Regarding the Nash-Sutcliffe Efficiency (NSE) criteria, better results were obtained between the simulated and observed values, with a mean NSE value of 81.5% (Table 4), which means that the model performs very well. The model simulation can be judged as satisfactory if Nash-Sutcliffe Efficiency is greater than 50%, good if it is greater than 65%, and very good if it is greater than 75% [34]. The mean value of Ratio of standard deviation of observations to root mean square error (RSR) obtained was 0.1, so according to Moriasi et al., [34], the model can be said as satisfactory if RSR<=0.7. The four statistical evaluation criteria with mean values of REP = 2.5%, REV = 4.8%, NSE = 0.815, R2 =0.842, RSR=0.1 showed good simulation between the estimated and observed values.

Validation
The table below shows the validation of the simulated results ( Table 5) of peak discharge, total volume, and their relative errors, the Nash-Sutcliffe Efficiency (NSE), the coefficient of correlation (R 2 ) and the (RSR) values. Figure 8 and 9 showed the simulated and observed hydrographs and their correlation respectively (Sample validation events are presented).   To increase the performance of the model in simulating the runoff, It is recommended to set up more rain gauge stations in the basin because the use of 3 stations is not sufficient to reduce the effect of spatiotemporal heterogeneity in precipitation, also for the flow data which are not enough to estimate perfectly the flows at the outlet, it is necessary to set up more hydrometric stations in the upper areas of the catchment, and to identify regional valid unit hydrographs and curve number as Zelelew [38] and Hawkins [39] recommend.

Conclusions
Our project is a perspective study which is inserted on the preparation of the data necessary as a first step for the hydrologic simulation using HEC-HMS model. The delimitation of sub-basins, the extraction of the hydrographic network and the development of the soil and land use databases was a very important step in this study. These were used with ArcGIS and HEC-GeoHMS to estimate the curve number in three states (CNI (dry), CNII (medium) and CNIII (wet)), and to determine a curve number map that has been used successfully for estimating surface runoff from the Beht watershed. The results demonstrate that the watershed is characterized by clay soil type dominated by poorly developed soils, brown soils and vertisols assimilated soils and a vegetation cover that is marked by significant presence of pasture lands which occupy more than 32%, followed by forests, agricultural land and bare land. The CN of the Beht watershed is medium to high, with an average value of 78, which means that the basin have a moderate runoff potential, the most runoff-producing areas have a high runoff coefficient. The model calibration was used to optimize the parameters and the most sensitive parameter in the simulation was the initial abstraction followed by the curve number. After optimization the peak flow and total volume of all events are very close to the observations. The performance of the HEC-HMS model was very good in terms of relative error functions, Nash-Sutcliffe Efficiency, coefficient of determination and Ratio of standard deviation of observations to root mean square error RSR based on the selected loss, transform and flow routing methods.
Finally, the results obtained showed that the model can be judged as valid and good in terms of evaluation criteria, but we suggest to set-up other meteorological and hydrometric stations in order to generate more details and to increase the performance of the model in the simulations.