Extreme Rainfall-Runoff Events Modeling by HEC-HMS Model for Koudiet Rosfa Watershed, Algeria

GeoScience Engineering Vol. 67 (2021), No. 4 http://gse.vsb.cz pp. 144–155, ISSN 1802-5420 DOI 10.35180/gse-2021-0060 EXTREME RAINFALL-RUNOFF EVENTS MODELING BY HEC-HMS MODEL FOR KOUDIET ROSFA WATERSHED, ALGERIA Ali HADDAD , Boualem REMINI Blida 1 University, Faculty of Technology, Department of Waters Science and Environment, 9000 Blida Algeria E-mail: [email protected], [email protected] ABSTRACT Changing of precipitation regime and intensification of extreme storms in semi-arid regions because of climate change requires the use of numerical models to forecast the outlet hydrographs. In this paper, HEC-HMS software was applied using a loss method the curve number CN to estimate the precipitation excess and a parametric unit hydrograph model to compute the transformation of precipitation excess into direct runoff over the watershed. The Muskingum-Cunge routing model was used to simulate the propagation of direct runoff through the main streams of Koudiet Rosfa watershed. The curve number CN and lag time parameters were used to calibrate the model towards several storms. The Nash–Sutcliffe efficiency coefficient (NSE) was adopted to assess the performance of the model to reproduce the observed hydrographs. Volume of the storms, peak discharges, times of peak and times of center of mass between the simulated and observed discharges were used to validate the model. The simulated discharges reproduce the hydrological events. The calibrated model was used to simulate the different hypothetical storms that could occur in the future in order to ensure the safety of Koudiet Rosfa dam towards extreme rainfall-runoff events. Keywords: HEC-HMS model; Hydrologic model; Rainfall-runoff; Semi-arid region. 1 INTRODUCTION Populations of semi-arid regions are facing serious challenges in the form of climate events such as perturbation of rainfall regime and intensification of extreme storms that generate spectacular disasters. Protection against flooding has taken the attention of researchers who develop flood forecasting models. Development of flood model prediction varies with application conditions and watershed size. It aims to compute the discharge hydrograph of the watershed and according to flow values, these may lead to either taking emergency measures or preventive operations on hydraulic structures like dams, or even to the evacuation of population during crisis situations. Flow modeling could be lumped or distributed. The first consists of computing the hydrograph at the downstream of the watershed. The second, also known as hydraulic routing, model consist of the different flow characteristics that are computed for every cross section along the channel [1]. Hydraulic models often use numerical methods to resolve the Saint-Venant equations [2]. Rainfall-runoff events modeling has become absolutely necessary for the ungagged watershed to fill the gaps of hydrological data. The most popular models used to simulate extreme floods are models that compute runoff volume, direct runoff or channel flow [3]. The runoff volume model assesses the volume of precipitations that fall over the watershed, the water infiltration rate, the runoff volume resulting from impervious surfaces and time of runoff beginning. The direct runoff models contain empirical models such as the unit hydrograph (Clark’s UH, Snyder’s UH, CSC UH, ModClark) and Kinematic wave which is a conceptual model. It describes the becoming of water that has not infiltrated or stored in the watershed. The models of channel flow (so-called hydraulic models), based on the resolution of the fundamental equations of open channel flow and given an upstream hydrograph as a boundary condition, simulate the flow along the channel by computing the downstream hydrograph [4–6]. The HEC-HMS is Hydrologic Modeling Software developed by the US Army Corps of Engineers of the Hydrologic Engineering 145 GeoScience Engineering Vol. 67 (2021), No. 4 http://gse.vsb.cz pp. 144–155, ISSN 1802-5420 DOI 10.35180/gse-2021-0060 Center (HEC) [6], it contains integrated tools of lumped and distributed models for modeling hydrologic processes. It consists of several components for calculating rainfall loss, direct runoff and routing [7]. The HEC-HMS model has been widely used to simulate the rainfall-runoff process over large watersheds because of its simplicity and capability to be used in common methods [7]. Tassew et al. (2019) [8] used HEC-HMS to perform a rainfall-runoff simulation of the Lake Tana Basin of the Gilgel Abay watershed in the upper Blue Nile basin in Ethiopia where the coefficient of determination and the Nash Sutcliffe Efficiency (NSE) were used to evaluate the model performance between simulated and observed hydrographs. Zhuohang et al. (2019) [9] have studied the applicability of HEC-HMS for Flash Flood Simulation in fourteen typical small catchments in hilly areas across China. The results show that the HEC-HMS distributed hydrological model is suitable to simulate the flash floods caused by intense rainfall. The main objective of this paper is to apply the HEC-HMS using its integrated components of coupled lumped-distributed models to compute the inflow hydrographs at Koudiet Rosfa dam by simulating the rainfall-runoff processes over the watershed. Results of simulation could be used to improve the safety of Koudiet Rosfa dam towards extreme storms. 2 METHODOLOGY 2.1 Area study Koudiet Rosfa watershed is located in semi-arid region of Algeria. It is a part of the hydrologic watershed of Cheliff (Fig.1). It is also considered as a subbasin of Oued El Fodda watershed. Its surface extends over 440 km2, with 88 km of perimeter, and it is discharged by a main stream of 32 km length until the dam that has taken the name of the watershed. Watershed elevations range from 600 m to 1786 m. The mean elevation is 904 m and the mean slope is 2.5 %. It is limited by UTM geographical coordinates (X1: 383055.47 m; Y1: 3943781.61 m; and X2: 406838.79 m; Y2: 3976609.04 m). Figure 1. Koudiet Rosfa watershed 2.2 Land use of Koudiet Rosfa watershed Only the downstream part of Koudiet Rosfa watershed that represents 15.6 % of the global surface is covered by clear forest. The rest of area is used for cereal agriculture and pasture. Fig. 2 represents the map of the land use. 146 GeoScience Engineering Vol. 67 (2021), No. 4 http://gse.vsb.cz pp. 144–155, ISSN 1802-5420 DOI 10.35180/gse-2021-0060 Figure 2. Land use of Koudiet Rosfa watershed 2.3 Characteristics of Koudiet Rosfa Dam Koudiet Rosfa dam is an embankment dam. It was built during the period of 1998 to 2004. With 57 m of height, the dam has as initial storage capacity 73 million m3. , Figure 3. Photo of Koudiet Rosfa Dam 2.4 Watershed modeling by HEC-HMS The area study was divided into 11 subbasins in order to conceive the watershed model by HEC-HMS as shown in Fig. 4. 147 GeoScience Engineering Vol. 67 (2021), No. 4 http://gse.vsb.cz pp. 144–155, ISSN 1802-5420 DOI 10.35180/gse-2021-0060 Figure 4. HEC-HMS model for Koudiet Rosfa watershed 2.5 Direct runoff model equations 2.5.1 Loss method The excess of precipitation was estimated by using the curve number CN method, which is a function of the ground using, the cover, the cumulative precipitation and past moisture, by the equation (1) [6]: 𝑃𝑒=(𝑃−𝐼𝑎)2 (𝑃−𝐼𝑎+𝑆) (1) where Pe = excess of accumulated precipitation at the time t; P = depth of accumulated rainfall at the time t; Ie = initial loss; S = potential maximum retention. If the accumulated rainfall does not exceed the initial loss, the runoff will remain zero. By many analyses for small experimental watersheds, a relationship between Ia and S was developed by the SCS [6] as follows: 𝐼𝑎=0.2 𝑆 (2) The accumulated precipitation excess will be: 𝑃𝑒=(𝑃−0.2𝑆)2 (𝑃+0.8𝑆) (3) The relationship between the maximum retention S and the curve number CN is given as following [6]: 𝑆=25400−254𝐶𝑁 𝐶𝑁 (4) 148 GeoScience Engineering Vol. 67 (2021), No. 4 http://gse.vsb.cz pp. 144–155, ISSN 1802-5420 DOI 10.35180/gse-2021-0060 Mccuen (1982) [10] discusses the use of the SCS runoff model in detail. The SCS published the values of CN that consider the soil type and land use for arid and semi-arid rangelands [5]. To consider the different ground types and land use of the watershed, a composite CN will be computed as follows [6]: 𝐶𝑁𝑐𝑜𝑚𝑝𝑜𝑠𝑖𝑡𝑒 =∑𝐶𝑁𝑖𝐴𝑖 ∑𝐴𝑖 (5) with CNcomposite = the composite number for runoff computations; i = the index of watershed subdivision for uniform land use and ground type; CNi = CN of subdivision i; Ai = area of subdivision i. 2.5.2 Transform method The soil conservation service (SCS) proposed a parametric unit hydrograph model to compute the transformation of excess precipitation into direct runoff [6]. The SCS unit hydrograph requires only a time to reach the peak which is related to the duration of the unit of excess precipitation as [6]: 𝑇𝑝=∆𝑡 2+𝑡𝑙𝑎𝑔 (6) For ungagged watersheds, the SCS relates the lag time to time of concentration as [6]: 𝑡𝑙𝑎𝑔 =0.6𝑡𝑐 (7) Time of concentration is defined as the time needed for water to flow from the most remote point in a watershed to the outlet watershed, it is calculated by the GIANDOTTI formula which is expressed as [11]: 𝑡𝑐=4√𝐴+1.5𝐿 0.8√𝐻𝑚𝑒𝑑−√𝐻𝑚𝑖𝑛 (8) where ∆𝑡 is the excess precipitation duration and tlag the difference of time between the center of mass of rainfall excess and the peak of the UH. tc is the concentration time, A the sub watershed surface, L the main stream length, Hmed the medium elevation and, Hmin the minimum elevation. Tab.1 shows the concentration and lag times for every subbasin. Table 1. Times of concentration and lag times N° Sub Basin Surface (km2) Time of concentration (min) Lag time (min) 1 12.603 162 97 2 5.783 205 123 3 26.445 260 156 4 64.608 243 146 5 8.518 362 217 6 113.220 315 189 7 13.950 382 229 8 51.573 537 322 9 22.300 237 142 10 22.470 265 159 11 81.140 462 277 149 GeoScience Engineering Vol. 67 (2021), No. 4 http://gse.vsb.cz pp. 144–155, ISSN 1802-5420 DOI 10.35180/gse-2021-0060 2.6 Channel flow modeling The Saint Venant equations composed by the continuity and momentum equations represent the fundamental equations for open channel flow modeling [6]. The first represents the volume of water in the reach including the outflow of the reach and that stored. The second represents the forces that act on the body of water for open channel [5]. In one dimension, these equations are written as [6]: 𝜕𝑄 𝜕𝑥+𝜕𝐴 𝜕𝑡 =𝑞𝐿 (9) 𝑆𝑓=𝑆0−𝜕𝑦 𝜕𝑥−𝑉 𝑔𝜕𝑉 𝜕𝑥 −1 𝑔𝜕𝑉 𝜕𝑡 (10) where A = wetted surface; qL = lateral inflow per unit length of channel; Sf = friction slope; S0 = bottom slope; y = hydraulic depth; x = distance along the flow path; V = velocity; t = time; g = acceleration due to gravity; 𝜕𝑦/𝜕𝑥 = pressure gradient; 𝑉 𝑔𝜕𝑉 𝜕𝑥 = convective acceleration and 1 𝑔𝜕𝑉 𝜕𝑡 = local acceleration [6]. The continuity and momentum equations are written by taking accounts of these assumptions: - horizontal water surface and constant Velocity for any channel section, - gradually and varied flow regime, - neglected vertical acceleration, - a trapezoidal channel was assumed, - Strickler roughness coefficient was assumed constant throughout the main stream. The Muskingum-Cunge model is based upon the diffusion form of the momentum equation [12]: 𝑆𝑓=𝑆0−𝜕𝑦 𝜕𝑥 (11) Combining this equation with continuity equation, we get the convective equation [12]: 𝜕𝑄 𝜕𝑡 +𝑐𝜕𝑄 𝜕𝑥 =𝜇𝜕2𝑄 𝜕𝑥2+𝑐𝑞𝐿 (12) where c = celerity wave; µ = hydraulic diffusivity. 𝑐=𝑑𝑄 𝑑𝐴 (13) 𝜇= 𝑄 2𝐵𝑆0 (14) where c = celerity wave; µ = hydraulic diffusivity; B = max width of water surface. By taking account of these assumptions, the Muskingum-Cunge hydraulic model is based upon the approximation of the continuity equation by using a simple finite difference [6]: 𝑄𝑡=𝐶1𝐼𝑡−1 +𝐶2𝐼𝑡+𝐶3𝑂𝑡−1 +𝐶4(𝑞𝐿∆𝑥) (15) The coefficients: 𝐶1=𝛥𝑡 𝐾+2𝑋 𝛥𝑡 𝐾+2(1−𝑋) (16) 𝐶2=𝛥𝑡 𝐾−2𝑋 𝛥𝑡 𝐾+2(1−𝑋) (17) 𝐶3=2(1−𝑋)−𝛥𝑡 𝐾 𝛥𝑡 𝐾+2(1−𝑋) (18) 150 GeoScience Engineering Vol. 67 (2021), No. 4 http://gse.vsb.cz pp. 144–155, ISSN 1802-5420 DOI 10.35180/gse-2021-0060 𝐶4=2(𝛥𝑡 𝐾) 𝛥𝑡 𝐾+2(1−𝑋) (19) K and X = coefficients expressed in terms of flow, channel and finite difference cell parameters as [13–16]: 𝐾=𝛥𝑥 𝑐 (20) 𝑋=1 2(1− 𝑄 𝐵𝑆0𝑐𝛥𝑥) (21) where Δx = space increment; c = wave celerity; Q = flow and S0 = channel bed slope. 2.7 Boundary conditions The necessary boundary conditions for the direct runoff models are the precipitations, which will produce runoff over the watershed. For the routing model, the boundary condition is the calculated hydrograph at the upstream of the reach. At the downstream of the reach, it is mostly recommended to use stage series or rating curve as downstream boundary condition to calibrate the model. In our area study, the main stream of Koudiet Rosfa Watershed does not contain control station with available flow data. Therefore, to perform the routing calculations, we assume in this paper that flow is normal for any channel section and, consequentially, the downstream boundary condition is the normal depth. 2.8 Initial conditions Initial conditions represent the state of the system before beginning the computation. There is necessary to calculate the discharge of every sub-basin that contributes to the channel flow. Water elevation in each cross section of each channel must be known. 2.9 Meteorological model In order to apply the HEC-HMS software on Koudiet Rosfa watershed, we have used the data of precipitations and dam exploitation during the period of 2004 to 2020 (NADT,2021) [16]. We have taken account of all storms that cumulated precipitation exceeding 40 mm. Using the criteria of matching between hyetographs and their corresponding observed hydrographs has enabled us to take account of only 4 events for the meteorological model. In fact, the storms of February 24/2014, February 20/2015, March 16/2018 and April 24/2018 (Fig. 5) were selected to calibrate the model where the cumulated precipitation ranging from 40 to 80 mm and the peak inflow from 12 to 100 m3/s. 151 GeoScience Engineering Vol. 67 (2021), No. 4 http://gse.vsb.cz pp. 144–155, ISSN 1802-5420 DOI 10.35180/gse-2021-0060 Figure 5. Hydrological events for model calibration 3 RESULTS AND DISCUSSION 3.1 Calibration and validation To calibrate the model, the storms of February 24/2014, February 20/2015, March 16/2018 and April 24/2018 were used as input data of precipitation. By assumption, the rainfall is considered uniformly distributed over the all sab-basins of the watershed. The computed hydrographs will be compared to the observed hydrographs. The curve number CN for loss method and the lag time for transform method were used to calibrate the model towards all events. Tab. 2 shows the parameters of calibration for each sub-basin. Table 2. Parameters of calibration Parameters N° of Sub-basin Loss method Transform method Curve Number (-) Lag time (min) 1 82 97 2 82 123 3 82 156 4 79 146 5 82 217 6 82 189 7 79 229 8 79 322 9 79 142 10 79 159 11 79 277 152 GeoScience Engineering Vol. 67 (2021), No. 4 http://gse.vsb.cz pp. 144–155, ISSN 1802-5420 DOI 10.35180/gse-2021-0060 The performance of the model to reproduce the observed hydrographs was evaluated by using the Nash–Sutcliffe efficiency coefficient (NSE) expressed by the Eq. (22) [17]. This parameter varies from negative values to 1. Several attempts were made in order to decrease the difference between simulated and observed values of both volume and peak flow for all considered events. After calibrating the model and performing the simulation, we have noted that the model reproduces the observed hydrographs with a NSE ranging from 0.60 to 0.87 (Tab.3). According to Moriasi et al. (2007) [18], an NSE near of 1 means a good accordance between observed and computed hydrographs. 𝑁𝑆𝐸=1−∑(𝑄𝑖,𝑜𝑏𝑠−𝑄𝑖,𝑠𝑖𝑚)2 𝑁 𝑖=1 ∑(𝑄𝑖,𝑜𝑏𝑠−𝑄𝑜𝑏𝑠        )2 𝑁 𝑖=1 (22) where Qi,obs = observed discharge; N = the number of observations; Qi,sim = simulated discharge; Qobs = the average observed discharge. Figure 6. Resulting hydrographs Table 3 gives the values of the objective functions for optimization.