Simulation of flood routing using the Muskingum method

Abstract

This paper gives an overview of flood routing and its simulation using the Muskingum method. According to (K Subramaniam, 2008) computation of peak discharge and duration of the high-water level is an important task in flood forecasting and for this paper, I have used the basic Muskingum method to compute outflow hydrograph. The geospatial hydrological analysis module HEC-GeoHMS and HEC-HMS were used for watershed delineation and for the sub-watershed division.  [2]

Keywords: HEC-GeoHMS, HEC-HMS, Muskingum method, flood routing

1. Introduction

The following sections describe the background, purpose of this report:

1.1 Background

To determine the temporal and spatial variations of a flood wave in a channel flow hydraulic and hydrologic flood routing methods are used by engineers [1]. The hydrologic method, on one hand, depends on the solution of the conservation of mass equation and relationship between the discharge and storage in a catchment or reservoir. On the other hand, hydraulic methods generally based on Saint-Venant equations consisting of the continuity and momentum equations. Before the computer, solving these equations was complex but computers have saved time and nowadays all those complex equations are solved using computer software’s which produces more vigorous solutions to the problems of flood routing. [2]

Since the 20th century, lots of simplified hydrological routing methods have been developed which was applied successfully in channel flow and reservoir. Presently, the Muskingum method (McCarthy 1938) is widely used due to its enough level of accuracy and their good relationship between their channel properties and parameters [3]. Here I have focused on the Muskingum method of flood routing. One of the advantages of the Muskingum method is that it can be solved using a minimum number of data required at the same time each flood routing method’s data requirement is different. The selection of method depends on the flow characteristics, data availability such as the topography of floodplain and stream, the roughness of channel as well as hydrological analysis purpose, etc. [2].

1.2 Purpose

The main objective of this paper is to give an overview of the hydrological flood routing method focusing more on the Muskingum method of routing. Next, this paper computes the flow hydrograph of the downstream station in Saint- Lawrence river (Sainte-Anne-De-Bellevue) using ArcGIS Hydro & HEC-GeoHMS tool and using the generated output from ArcGIS to HEC-HMS. The output hydrograph is computed using the Muskingum method in HEC-HMS [3].

2. Flood Routing

Routing is a technique that is used to calculate the changing shape of

a hydrograph as water moves through a river channel or a reservoir. It is used to lowland flooding potential through predicting the hydrograph shape in different sub-catchments of the watershed. Factors such as watershed overall shape, sub-catchment area shapes, the topography of the area, geology or hydrogeology of the area where forests cover and aquifer act like a sponge that absorb the rainfall water and release it later weeks or months, duration of rainfall events, length of the stream reach all pay an important role in changing the shape of outflow hydrograph.

The out-flow hydrograph results sometimes show additive effects for example: when a large flood peak of all sub-catchments arrives at the outlet of the watershed mouth at the same time, as a result of that causing a stacking of the hydrograph peak. [4] [5] [10] The upstream section of a river is denoted by the inflow hydrograph and the outflow hydrograph is determined from the water flow at the downstream end. In the hydrograph, two changes of the flood wave are described i.e., ‘attenuation’ which describes the decrease in the discharge and ‘translation’ which refers to the peak discharge delay moving downstream. [6] Figure-1: Flood routing hydrograph for natural channel (Safa Elbashir, 2011) [6]

Other uses of flood routing include floodplain studies, reservoir, and spillways design, etc. Routing can be categorized into reservoir routing and channel routing. In reservoir routing generally include the flood wave studies entering a reservoir. The effect of flood wave entering the reservoir is studied using volume-elevation characteristics of the reservoir and outflow-elevation characteristics for the outlet of any structures in the reservoir which is used to predict the reservoir elevation and outflow discharge variations with time. This type of reservoir routing is essential in designing the capacity of spillways outlet or any type of reservoir outlet, also used to find the location as well as calculating the size of the reservoir according to the requirement.

Channel routing is generally used to predict flood hydrograph at the various sections of the reach using input hydrograph and channel reach at the upstream end. Duration of high-water level and flood-peak attenuation information is obtained using channel routing is essential in flood-forecasting operations as well as flood protection works [6].

2.1 Runoff Generation and Stream flow

Streams are the paths through which almost all the runoff after separated from evaporated and transpiration portion of rainfall is returned to the ocean and completes the global hydrologic cycle. There are three important areas that make the study of this process important. Firstly, the proper management of water resources constitutes by the streams to use it in a sustainable way for human use as well as for stream organisms. Secondly, for flood forecasting which is done through estimating the probability of flood occurrence of various magnitudes. This prediction helps in the design of bridges, levees, and dams as well as in preparing floodplain land-use plans and regulations. Lastly, its study is necessary to regulate the water quality which may affect the suitability of water for use by various natural organisms and humans. [7]

2.1.1 The Watershed and Stream Network

A watershed is defined as an area that topographically appears to contribute all the water that passes through a specified cross-section of the stream outlet (S Lawrence Dingman, 2015) For watershed delineation determining its boundary also known as divide or a ridgeline is an important task. It begins with the selection of watershed outlet which is usually a stream-gauging station because its easy to monitor the streamflow for quantitative studies of water budgets or stream response. On the other hand, for geomorphological analyses of stream networks and landscapes, stream junction or the point where streams enter an ocean, or a lake is usually selected as an outlet. [7] [10]

Drainage of excess precipitation over a watershed which is usually larger than a minimum area occurs via a hierarchical stream network. Strahler (1952) has given the most common approach to designate stream orders. Referring to the below figure, according to this approach (Strahler, 1952) streams with no tributaries are designated as first-order streams, second-order streams are usually said when two first-order streams meet and join, the confluence of two second-order streams adds to form a third-order stream. [7] [11] Figure-2: Stream order designation (S Lawrence Dingman, 2015) [7]

2.2 The Hydrograph

Hydrograph refers to the graph which shows the change of discharge with respect to time of a river. Consider a concentrated rainfall of duration D over a catchment. The excess rainfall reaches the stream after the losses like an infiltration, evapotranspiration, transpiration, etc. Over time a certain amount of storage is build up over the land and channel flow phases.  This storage gradually reduces after rainfall stops. Thus, between the occurrence of the rainfall in the catchment and the time when water passes the outlet of the gauging station there exists a time lag. The runoff measured at the gauging station end is usually given a standard hydrograph shown below. This single-peaked skew distribution hydrograph is also known as flood hydrograph or storm hydrograph (K Subramanya, 2008) [2] [12] Figure-3: Standard Flood Hydrograph basic elements (K Subramanya, 2008) [2]

The typical hydrograph usually has three main characteristics i.e., the rising limb (AB) where A is the starting point and B is the point of inflection, the crest segment (BC) having peak point (P) which indicates the peak discharge at a particular point of time (t) and the depletion curve or falling limb (CD). This limb starts from point C which is the endpoint of the crest segment limb to the start of groundwater flow point D. The start point C indicates the maximum storage condition of the basin. After the cessation of rainfall, depletion of storage takes place and from this point, the shape of the hydrograph depends entirely on the characteristics of basin and independent of storm characteristics. Water is stored in the catchment in various forms like surface storage, interflow storage, groundwater storage also known as base flow storage. (K Subramanya, 2008) [2] [13]

2.2.1 Base Flow Separation

In hydrograph, water that is supplied from groundwater sources is said as base flow and it doesn’t respond quickly to the rainfall. After separating the base flow from the total storm hydrograph, we can determine the surface-flow hydrograph. There are two basic methods to separate base flow i.e., graphical separation method and filtering techniques. The graphical method is based on determining the points of intersection between the base flow and the rising and falling limbs of the quick flow response (Safa Elbashir, 2011). On the other hand, the entire hydrograph is used in determining base flow separation in case of filtering techniques. [2] [6]

2.2.1.1 Graphical Separation Method

This method includes determining the point of intersection of the base flow and the falling limb as shown in the figure below. The base flow until the start of the new hydrograph for the next rainfall event. Four methods are used in this case: Figure-4: Graphical base flow separation method where 1a represent constant discharge method, 1b represent constant slope method and 1c represent concave method (Safa Elbashir, 2011) [6]

• Empirical Method: This method contains a numerical equation that is used to determine the point on the hydrograph along the falling limb (see figure 4). After finding this point a straight line is drawn which separates the quick flow from the base flow. [6] (i) Where D is the number of days between the crest limb of the storm and the end of quick flow, A is the catchment area (Sq Kms). Here constant value 0.2 depends on the geology and characteristics of the catchment. [6]

• The constant discharge method: This method assumes a constant base flow throughout the storm hydrograph. [6]

• The constant slope method: This method assumes a constant slope between the start of the rising limb and the inflection point of the falling limb.[6]

• The concave method: represents the initial decline in base flow during the rising limb by projecting the decrease in a hydrograph trend prior to the rainfall event to directly under the crest of the flood hydrograph (Safa Elbashir, 2011). From the above four techniques, constant discharge and the constant slope method are most commonly used to determine base flow. [6]

2.2.1.2 Filtering Separation Methods

In this method, an automated index is used related to the response of the base flow of a catchment. The base flow index (BFI) or reliability index is defined as the ratio between base flow and total flow. This ratio is determined from a hydrograph smoothing and separation procedure using daily discharge (Safa Elbashir, 2011). [6] [13]

3. Hydrological Method

3.1 Basic Equation

The flow of flood hydrograph either through a channel or a reservoir is an unsteady flow event. In terms of open channel, it is usually said as gradually varied unsteady flow (K Subramaniam). This paper focuses on the hydrological method of channel routing which uses the principle of continuity equation to solve the mass balance of inflow, outflow and the volume of storage. Methods included in the hydrological approach require a storage-stage-discharge-relation to determine the outflow for each time step. It uses numerical techniques that introduce translation or attenuation to the hydrograph of inflow discharge (Safa Elbashir, 2011) [6].

According to the equation of continuity, the rate of change of storage is equal to the difference between inflow and outflow rates i.e., [2] [6] [13] (ii) Where I(t) is the inflow rate at the upstream section and O(t) is the outflow rate at the downstream section in (m3/s) and S is the storage in m3. Above equation can be written infinite differential form over the finite interval of time t and t+Δt as: (Safa Elbashir, 2011) [2] [6] (iii) The hydrological method consists of two forms one is the Muskingum method which depends on the solution of the continuity equation and the other is the Muskingum-Cunge method which depends on the continuity equation and includes the diffusion form of the momentum equation. (Safa Elbashir, 2011) [2] [6]

3.2 Muskingum Method

3.2.1 The Basic Muskingum Method

This method of routing assumes the storage in a channel reach is a combination of prism and wedge storage (see figure below). The total storage, S in a channel reach can be expressed as: (Safa Elbashir, 2011) [2] [6] (iv) Where K and O are the travel time through the reach and O is the flow through the prism respectively, X refers to the weighted factor which lies in the range of  0 ≤ X ≤ 0.5 and when X is 0, it indicates that inflow has little or no effect on the storage. In natural streams value of X lies between 0 and 0.3 (Safa Elbashir, 2011). Writing above equation at time increments of Δt, the storage S at times jΔt and (j+1) Δt can be written as: [2] [6]           (v)                                (vi)

Subtracting above equations we get the change in storage over the time interval Δt i.e., [2] [6] (Vii) Equation 3 can be written in its discretized form of the continuity equations: [2] [6] (viii) Using Equation (vii) and (viii) and combining we get the routing expression: [2] [6] Where C0, C1, C2 are given by: [2] [6] Where m = 2K(1-X) + Δt C0 + C1 + C2 = 1

3.2.2 Muskingum Parameters

The parameters K and X can be determined using available inflow and outflow data and graphically using the plot of S and XI+ (1-X) O. This will produce a straight line with a slope of K. Trial and error is done to compute the value of X and the value that gives the narrowest loop as shown in the figure below is taken as correct X value and value of K is taken from its slope (Safa Elbashir, 2011) [2] [6] [13] Figure-5: Storage loop (Safa Elbashir, 2011) [6]

4. Study Area

The St Lawrence River is an important river that servers over 30 million Americans and almost 15 million Canadians. It consists of America’s third-largest hydrological system after Mississippi and McKenzie rivers. It consists of lake Michigan, Lake Superior, Lake Erie, Lake Ontario, Lake Huron, the St Lawrence River. It covers the surface area of almost 1.6 million km2 with 55% forest area, 22% urban, 20 % of cropland and other land cover comprises of 3%. [8] This study mainly focuses on the watershed covering half Montreal on the river St Lawrence and Lake Saint-Louis. Montreal is located at the confluence of the Ottawa River with St Lawrence.

In my case study, I have used four gauging stations on the St Lawrence River. Beginning from upstream to the downstream these gauging stations are (2) Saint Laurent a Lasalle (02OA016), (4) L’acadie (02OJ026), (3) Chateauguay (02OA054), (1) Outaouais (02OA033). The below figure shows the catchment area along with four gauging stations. [8] (See Appendix-A for input flow graph) Figure-6: Study Area- the river St Lawrence and Lake Saint-Louis along with four stations. (Source: Environment Canada) All these river systems flowing from Canada and enters the United States play a crucial role in boosting the economies of both countries. This study aims to study the flow hydrograph of the year 2017 most importantly month of May 2017 when the flood-hit Montreal. Seeing the importance of the region it is necessary to understand the hydrology of the river nearby to minimize the effect of extreme events. [8]

4.1 Methodology

4.1.1 Preparing Base Map Using ArcGIS Hydro & HEC-GeoHMS Tool

4.1.1.1 Terrain Preprocessing

It is usually done to identify the surface drainage pattern using DEM (Digital elevation model). It is an important task that needs to be done before delineating a watershed and generating stream network. (Source: ESRI) [9] DEM Manipulation: It contains all those useful functions using which one can edit to the original DEM. They are Level DEM, DEM reconditioning, assigning stream slope, burn stream slope, build walls & fill sinks. (Source: ESRI) [9]

• Level DEM: It is used to modify a DEM by setting the cells of the selected polygon features to the associated ‘fillElev’ value. It uses two inputs i.e., raw input and a polygon feature class, for example, a lake. [9]

• DEM Reconditioning: This function is used to modify a DEM. Three input values need to be put in this function that are: Stream buffer (number cells for which smoothing need to be done), smooth drop raise value (It is used to interpolate the DEM and buffered area) and sharp drop/raise value (It is used to put additional fencing on top of the smooth buffer interpolation. (Source: ESRI) [9]

• Assign Stream slope: It is used to assign stream slope to the input stream feature class. In this function we need to assign two input values i.e., maximum start elevation that will be assigned to the ‘from node’ of the most upstream feature and secondly, the elevation drops between two nodes. [9]

• Burn stream slope: This tool ensures the water to flow in the digitized stream direction until and unless it reaches a stream. [9]

• Build walls: It is used to build walls in the input grid. Generally, two types of walls are created i.e., outer walls and inner walls. The first one depends on the input polygon feature class while the latter one is based on the input line or point, polygon feature class. After successfully completing this process, the ‘walledDEM’ layer will be added to the map. [9]

• Fill Sinks: It is used to fill the sinks in a grid. If a cell is surrounded by cells with higher elevation values, then the flowing water is trapped and cannot flow further. Fill sink function fills all those sinks which have depth lower than the specified threshold. [9] Terrain Processing: After completing DEM manipulation and getting the required edited DEM, the next thing which is required to be done is the delineation process using suitable functions like Flow Direction Grid, Catchment, and Adjoint Catchment.

• Flow Direction: This function is used to calculate the direction of flow for a given grid. Each cell of a flow direction grid has a specific value which indicates the direction of the steepest slope from that cell.  This function has ‘Fdr Filled’ input which has 8 distinct values, each of which indicates the water flowing direction. [9]

• Flow Direction with Sinks: This function is used to calculate the direction of flow for a grid with sinks to make sure that all cells within a sink watershed flow into watershed’s sink point. This will create a flow direction grid by creating a sink point for each processed sink to make sure that each cell in a sink flows towards its sink point. Both the inputs i.e., Sink watershed grid and Sink link grid under this function are used to mask these areas (area in a grid that drains into each sink) to make sure that no stream links may generate within the sink watershed or sinks while creating stream links using the function Stream Segmentation. [9]

• Flow Accumulation: This tool helps to calculate the flow in each cell flowing from upstream to downstream. Higher flow accumulation values are assigned to lower elevation areas like a valley. Flow accumulation tool allows locating cells through which stream flows. Each cell indicates the pour point through which stream flows out. [9]

• Stream Definition: Based on the flow accumulation grid and a user-specified threshold this function helps to compute a stream grid. Output generated from flow accumulation function is used as an input to this function where the cell that has a value greater than the threshold are assigned stream grid value of 1 and other cells are assigned no data. [9]

• Stream Segmentation: This function help in creating a stream segments grid that has a unique identification. A segment may be defined as a segment between two segment junctions or ahead segment. Cells of a segment are assigned grid codes like that segment. [9]

• Flow direction with Streams: Using the existing stream layer this function creates the drainage line from the input flow direction grid. This function can be used instead of the Stream definition to match with the input stream. Input flow direction is edited using this function to generate an output stream sloped flow direction grid. Further, this function makes sure that the water remains within the stream without jumping between streams.

• Catchment Grid Delineation: This function assigns a specific value (grid code) to a cell and creates a grid. This helps to categorize each cell depending upon its location in the catchment. [9]

• Catchment Polygon Processing: It is used to convert the catchment grid generated using the previous function into a catchment polygon feature class. [9]

• Drainage Line Processing: This function uses the output of stream segmentation function i.e., the input stream link grid into drainage line feature class. [9]

• Adjoint Catchment Processing: This function speeds up the delineation process. This helps to create a polygon constituting the whole upstream area draining to its inlet point.

• Drainage Point Processing: This function is used to create the drainage points on each sub-catchment. [9] After completed all the processes delineated watershed is formed as shown in figure-7 below: Figure- 7: Study Area- Watershed along the St. Laurent River with stations in green dot, red line indicates the inlet to outlet flow connection & blue line indicates stream flow path (Source: Created in Arc Hydro & ArcGeo-HMS)

4.1.2 Flood Routing in HEC-HMS and Result

After completing all necessary processing in HEC-GeoHMS the generated base map is transported to HEC-HMS for further computation of outflow. Four stations and outlet/reach are set up as shown in the figure-8 below: Figure-8: HEC-HMS model showing four stations (Source: Created using HEC-HMS)

4.1.2.1 Result

Using the value of coefficient K as 48 hrs and X as 0.3 hydrograph of outflow is generated (figure-9). Figure-9: Hydrograph of combined inflow (dash line) and result outflow (Source: HEC-HMS tool)

5. Conclusion

We can conclude from the result that the peak inflow was 15283.8 m3/s on May 07, 2017, and peak discharge was 15043.8 m3/s on May 11, 2017. Since real flooding occurred in the same time period on the selected watershed, therefore, it can be said that simulation results were satisfactory. Also, the use of geospatial analysis technique (HEC-GeoHMS) and hydrological model (HEC-HMS) help in determining basin topography parameters. Further, it can be said that the Muskingum model shows high simulation accuracy for generating outflow hydrograph and hence doing flood forecasting. [4] [6] 6. Recommendations and Future Works:

• Use of computer programming: When it comes to simulating a large number of inputs and stations it became time-consuming to process each tool. So, to reduce the time computer programming like python (Arcpy) can be used for watershed delineation. Its main aim is to perform geographical data analysis, data management and map automation with Python. (Safa Elbashir, 2012) [6]

• Muskingum Parameters: To get accurately simulated hydrograph of the outflow discharge coefficient of K and X is important and should be accurately computed. To check its accuracy regression analysis can be done. (Safa Elbashir, 2012) [6]

• Use of Artificial Intelligence: Artificial neural networks (ANN) can also be used to determine the channel routing sensitivity to past time steps discharge which are effective in modeling accuracy. (Safa Elbashir, 2012) [6]