Flood routing in ungauged catchments using Muskingum methods

River stage or flow rates are required for the design and evaluation of hydraulic structures. Most river reaches are ungauged and a methodology is needed to estimate the rates of flow, at specific locations in streams where no measurements are available. Flood-routing techniques are utilised to estimate the stages, or rates of flow, in order to predict flood wave propagation along river reaches. Models can be developed for gauged catchments and their parameters related to physical characteristics such as slope, reach width, reach length so that the approach can be applied to ungauged catchments within the same region. The objective of this study is to assess the Muskingum-Cunge method for flow routing in ungauged river reaches, both with and without lateral inflows. The Muskingum-Cunge method was assessed using catchment-derived parameters for use in ungauged river reaches. Three sub-catchments in the Thukela catchment in KwaZulu-Natal, South Africa were selected for analyses, with river lengths of 4, 21 and 54 km. The slopes of the river reaches and reach lengths were derived from a digital elevation model. Manning’s roughness coefficients were estimated from field observations. Flow variables such as velocity, hydraulic radius, wetted perimeters and flow depth were determined both from empirical equations and assumed cross-sections of the reaches. Lateral inflows to long river reaches were estimated from the Saint-Venant equation. The performance of the methods was evaluated by comparing both graphically and statistically the simulated and observed hydrographs. The results obtained show that the computed outflow hydrographs generated using the Muskingum-Cunge method, with variables estimated using both the empirical relationships or assumed cross-sectional shapes, resulted in reasonably accurate computed outflow hydrographs with respect to volume, peak discharge, timing of peak flow and shape of the hydrograph. From this study, it is concluded that the Muskingum-Cunge method can be applied to route floods in ungauged catchments using derived variables in the Thukela catchment and it is postulated that the method can be used to route floods in other ungauged rivers in South Africa.


Introduction
As defined by Fread (1981) and Linsley et al. (1982), flood routing is a mathematical method for predicting the changing magnitude and celerity of a flood wave as it propagates down rivers or through reservoirs.Numerous flood-routing techniques, such as the Muskingum flood-routing methods, have been developed and successfully applied to a wide range of rivers and reservoirs (France, 1985).Generally, flood-routing methods are categorised into two broad, but somewhat related applications, namely reservoir routing and open channel routing (Lawler, 1964).These methods are frequently used to estimate inflow or outflow hydrographs and peak flow rates in reservoirs, river reaches, farm ponds, tanks, swamps and lakes (NRCS, 1972;Viessman et al., 1989;Smithers and Caldecott, 1995).
Flood routing is important in the design of flood protection measures in order to estimate how the proposed measures will affect the behaviour of flood waves in rivers so that adequate protection and economic solutions can be found (Wilson, 1990).In practical applications, two steps are involved in the prediction and assessment of flood level inundation.A flood-routing model is used to estimate the outflow hydrograph by routing a flood event from an upstream flow-gauging station to a downstream location.Then the flood hydrograph is input to a hydraulic model in order to estimate the flood levels at the downstream site (Blackburn and Hicks, 2001).
Flood-routing procedures may be classified as either hydrological or hydraulic (Choudhury et al., 2002).Hydrological methods use the principle of continuity and a relationship between discharge and the temporary storage of excess volumes of water during the flood period (Shaw, 1994).Hydraulic methods of routing involve the numerical solutions of either the convective diffusion equations or the one-dimensional Saint-Venant equations of gradually varied unsteady flow in open channels (France, 1985).
As noted by the US Army Corps of Engineers (1994), several factors should be considered when evaluating the most appropriate routing method for a given situation.The factors that should be considered in the selection process include, inter alia, backwater effects, floodplains, channel slope, hydrograph characteristics, flow network, subcritical and supercritical flow (US Army Corps of Engineers, 1994).The selection of a routing model is also influenced by other factors such as the required accuracy, the type and availability of data, the available computational facilities, the computational costs, the extent of flood wave information desired, and the familiarity of the user with a given model (NERC, 1975;Fread, 1981).
The hydraulic methods generally describe the flood-wave profile adequately when compared to hydrological techniques, but practical application of hydraulic methods is restricted because of their high demand on computing technology, as well as on quantity and quality of input data (Singh, 1988).Even when simplifying assumptions and approximations are introduced, the hydraulic techniques are complex and often difficult to implement (France, 1985).Studies have shown that the simulated outflow hydrographs from the hydrological routing methods always have peak discharges higher than those of the hydraulic routing methods (Haktanir and Ozmen, 1997).However, in practical applications, the hydrological routing methods are relatively simple to implement and reasonably accurate (Haktanir and Ozmen, 1997).An example of a simple hydrological flood-routing technique used in natural channels is the Muskingum flood-routing method (Shaw, 1994).
Among the many models used for flood routing in rivers, the Muskingum model has been one of the most frequently used tools, because of its simplicity (Tung, 1985).As noted by Kundzewics and Strupczewski (1982), the Muskingum method of flood routing has been extensively applied in river engineering practices since its introduction in the 1930s.The modification and the interpretation of the Muskingum model parameters in terms of the physical characteristics, extends the applicability of the method to ungauged rivers (Kundzewicz and Strupczewski, 1982).Most catchments are ungauged and thus a methodology to compute the flood-wave propagation down a river reach or through a reservoir is required.One option is to calibrate floodrouting models on gauged catchments and relate their parameters to physical characteristics (Kundzewicz, 2002).The floodrouting models with derived parameters can then be applied to ungauged catchments in the region (Kundzewicz, 2002).
In this study, the Muskingum-Cunge method is adopted to estimate the model parameters because of its simplicity as well as its ability to perform flood routing in ungauged catchments by estimating the model parameters from flow and channel characteristics.The Muskingum-Cunge parameter estimation method utilises catchment variables such as flow top width (W), slope (S), average velocity (V av ), reference discharge (Q 0 ), celerity (V w ), and reach length (L) to estimate the parameters of the Muskingum method.
When performing flood routing in ungauged catchments, the model parameters have to be estimated without observed hydrographs.The inflow hydrograph could be generated using a hydrological model such as the ACRU Model (Schulze, 1995).For this study, observed inflow hydrographs are used as input to simulate outflow hydrographs.
The objective of this study was to assess the performance of the Muskingum-Cunge method in ungauged catchments using parameters derived using both variables estimated from empirical equations developed for different river reaches and using variables estimated from selected (assumed) cross-sections within the river reach.The study was performed using data from the Thukela catchment in South Africa.

Methodology
This section describes the flood-routing methods applied to selected events.At ungauged sites the Muskingum-Cunge equation was used, with flow variables estimated using both an empirical approach (MC-E) and using variables computed from a user-defined river reach cross-section (MC-X).The computed outflow hydrographs were computed and compared statistically and graphically with the observed hydrographs.The methodology is detailed as follows: • The Muskingum K and X parameters were estimated for each reach on the basis of empirically determined flow characteristics such as flow top width (W), wetted perimeter (P), flow depth (y), hydraulic radius (R) and average velocity (V av ).The method is referred to as MC-E.• The Muskingum K and X parameters were estimated for each reach on the basis of assumed (selected) channel crosssectional shape and dimensions, as observed in the field, such as maximum flow depth (y) and maximum flow top width (W).The method is referred to as MC-X.• Outflow hydrographs computed using the estimated parameters were compared to the observed hydrographs both statistically as well as visually with regard to flood volume, the magnitude and timing of peak discharge as well as the shape of the hydrograph.
The details of the calculation steps are explained in the following sections.

The MC-E approach
The depth of flow and discharge in a reach can be derived from empirical relationships recommended by the US Reclamation Service (Chow, 1959).Hence, the Muskingum K and X parameters can be estimated in ungauged river reaches using inflow hydrographs and channel dimensions estimated from empirical relationships.In this study, the observed inflow hydrographs obtained from the Department of Water Affairs and Forestry ( 2003) and empirically estimated channel variables are used to estimate the Muskingum K and X parameters.

381
where: Wilson and Ruffini (1988) defined a reference flow as: (4) where: For a given hydrograph the reference flow (Q 0 ) can be estimated using Eq. ( 4).The roughness coefficient (n) and slope (S) were estimated from field visit and topographic maps, respectively.Thus, the section factor in Eq. (3b) has to be estimated using either empirical relationships which relate discharge, depth and slope or by using charts that relate section factor to the depth of flow (Chow, 1959).Equation ( 5) relates the wetted perimeter with discharge for natural rivers as follows Lacey (1930;1947, cited by Punmia andPande, 1981): for stable river channels (5) where: The hydraulic radius (R) and hydraulic mean depth (d) can be assumed to be equal when the top flow width exceeds the mean flow depth by a factor of 20 m (Chow, 1959;Barfield et al., 1981, cited in SAICEHS, 2001).Since the flow depth in most river reaches is small relative to the top width of a reach (W), the wetted perimeter (P) may be assumed to be equal to W. The area of parabolic section is computed as (Chow, 1959;Koegelenberg et al., 1997): (6a) where: For wide parabolic channels, the hydraulic radius may be estimated from Table 2 as: (6b) where: Equations ( 6a) and (6b) can be substituted in Eq. (3b) and for a parabolic section results in: (6c) But P is assumed to be approximately equal to W and hence P can be substituted in Eq. (6c) to solve for y as shown in Eq. ( 7): (7) For the MC-E approach, Eq. ( 7) is used to estimate the depth of flow in ungauged catchments in order to estimate the Muskingum K and X parameters, as shown below.

The MC-X approach
For the MC-X approach the depth of flow in ungauged catchment is estimated by developing a rating curve for a selected (assumed) section in the reaches (Smithers and Caldecott, 1995).
For the selected section a rating curve is developed based on the maximum width and depth relationships.Assuming a linear relationship between width and depth of the river section, for each given depth, a corresponding top width can be proportioned from the observed maximum depth and width, as shown in Eq. ( 8): (8) where: The wetted perimeter for each sub-section is calculated from geometrical equations (Chow, 1959;Koegelenberg et al., 1997).For example, in the case of a parabolic section: (9) where: The flow area can be estimated from the continuity equation or geometrical properties for a parabolic shape (Chow, 1959;Koegelenberg et al., 1997). where: The hydraulic radius is computed from the cumulative area and wetted perimeter as shown in Eq. ( 10): (10) Then, the corresponding cumulative discharge is calculated from the Manning equation as shown in Eq. ( 11): (11) where: Since the roughness coefficients (n), hydraulic radius (R), flow area (A) and slope (S) of each reach are known, the discharge can be calculated from Eq. ( 11).Thus, from the derived rating curve, the depth of flow can be estimated for a given discharge and the corresponding width and wave celerity can be computed using Eq. ( 8) and Table 1 respectively in the MC-X approach.Eqs.
(1) and ( 13) are then used to estimate the K and X parameters respectively.
The average velocity calculated using Eq. ( 3a) is also used to calculate flow area of the reference flow from the continuity equation in the MC-E approach (Chow, 1959;Linsley et al., 1988;KenBohuslay, 2004).Alternatively, the area of flow can be estimated using the geometrical parameters as shown in the following equation for a parabolic section: (12) The wetted perimeter (P) and top flow width (W) are assumed to be equal.
The three coefficients computed in Eq. ( 14) are used in Eq. ( 15) to estimate the discharge in a successive time step (NERC, 1975): (15) where: = outflow at time [t] of the j th sub-reach = inflow at time [t] of the j th sub-reach C 3 = lateral inflow term, computed as shown in Eq. ( 17) As suggested by Viessman et al. (1989), negative values of C 1 must be avoided.Negative values of C 2 do not affect the routed hydrographs.The negative values of C 1 can be avoided by satisfying Eq. ( 16): (16) After the X parameter has been determined, the routing time interval can be adjusted using the relationship shown in Fig. 1 (Cunge, 1969;NERC, 1975).
Since Q t , I t and I t+1 are known for a given time increment, Q t+1 is computed using Eq. ( 15) and repeated for successive time increments to estimate the outflow hydrograph.In large catchments, where there is a lateral inflow to the main stream, the volume of the outflow hydrograph may be larger than the inflow volume.Hence, lateral inflow should be considered and added to the main flow as follows (NERC, 1975).
The C 3 coefficient, calculated using Eq. ( 17), is added as a lateral inflow term in Eq. ( 15) (NERC, 1975): (17 where: The lateral inflow per unit length at a specified time (q) is estimated from Saint-Venant equation.The Saint-Venant equation for gradually varying flow in open channels is given as follows (NERC, 1975): (18) Fread (1998) approximated the terms in Eq. ( 18) as: Hence, lateral flow per unit length at a specified time (q) can be calculated from Eq. (19c) and substituted in Eq. ( 17).

(19c)
When the ratio of lateral inflow to the main flow is too large, numerical difficulties in solving Eq. (19b) may arise.Increasing the routing length (∆L j ) of the specified reach may solve the numerical difficulties (Fread, 1993).

Figure 1
The Cunge curve (Cunge, 1969, cited by NERC, 1975)    According to Fread (1998), the value of β is between 0.5 and 1.0.A value of 0.7 was used in this study.
Flow characteristics from the routed hydrograph such as the magnitude and timing of peak flow, hydrograph shapes and flow volume were compared to the characteristics of the observed events.The statistics used to assess the model performance are detailed in the following section.

Model performance
As suggested by Green and Stephenson (1985), in order to compare a model output to the observed data, criteria for making such a comparison must first be identified.Visual comparison by plotting simulated and observed hydrographs provides a valuable means of assessing the accuracy of the model output.However, visual comparisons usually tend to be subjective and need additional statistical analysis.To overcome these difficulties, as well as to highlight certain model peculiarities, statistical goodness-of-fit procedures were employed.
The difference in the observed and computed hydrograph were analysed by means of the root-mean-square error (RMSE) and other goodness-of-fit statistics.A statistical goodness-of-fit procedure implies a procedure employed to measure the deviation of simulated output from the observed input data set (Green and Stephenson, 1985).
Even though numerous goodness-of-fit criteria for assessing the accuracy of simulated output have been proposed, particular aspects may give more weight to certain output interests (Green and Stephenson, 1985).Hence, different goodness-of-fit statistics should be applied to assess different hydrograph components such as flood volume, hydrograph shape, peak flow magnitude and timing.Since an objective of this research is to compute hydrographs in ungauged reaches, the criteria for assessments were selected as described below.
Equation ( 22) was used to estimate the actual errors in the computed hydrographs.RMSE computes the magnitude of error in the computed hydrographs (O'Donnell, 1985;Schulze et al., 1995) As peak outflow is important in a single event model, a comparison of computed and observed peak flow rates, peak timing and volume were computed as shown in Eqs.(23a, b and c) (Green and Stephenson, 1985). where: Even though the RMSE, E peak , E time and E volume statistics may appear to be reasonable, the shapes of the respective hydrographs may be different.Nash andSutclife (1970, cited by Green andStephenson, 1985) proposed a dimensionless coefficient of model efficiency (E).The computed hydrograph is a better fit to the observed hydrograph when the coefficient of model efficiency (E) approaches 1 (Green and Stephenson, 1985).Hence, the hydrograph shape comparison was estimated using Eq. ( 24): ( where: The methodology outlined above was assessed on selected subcatchments in the Thukela catchment.

Study area
Three sub-catchments in the Thukela catchment were selected for analyses, with river lengths of 4, 21 and 54 km.The locations of the gauging weirs are shown in Fig. 2 (on next page).The characteristics of the selected reaches and variables calculated using both methods are contained in Tables 3 and 4 for a selected event in each reach.The computed and observed hydrographs from the applications of the MC-E and MC-X methods for the selected event in the three reaches are shown as Figs. 3, 4 and 5 and the corresponding performance statistics are summarised in Tables 7 and 8.

�L
As shown in Tables 7 and 8, the event in Reach-II has a relatively small RMSE error and volume error with the coefficient of

Results
The results from only one event in each reach with the hydrographs are reported.Additional results reported by Tewolde ( 2005) are indicated for statistical comparisons at the end of this section.The parameters computed using the MC-E and MC-X methodologies are contained in Tables 5 and 6 respectively.efficiency (E) nearly equal to 1.The other events generally have acceptable statistical results.These results indicate a high degree of correlation between the computed and observed hydrographs.
The result obtained from all three events show that the computed hydrographs are similar to the observed hydrographs, with errors of less than 12% for the statistics considered.Similar results were obtained by Tewolde (2005) for additional events from the three reaches.The summary of results of other events computed in the three Thukela catchments by Tewolde (2005) are contained in Tables 9 and 10 for statistical evaluation.
The addition of lateral inflow in the simulated hydrographs was not adequate when compared to the observed outflow.As the length of a reach increases the possibilities of tributary inflows also increases.Hence, in large catchments it is recommended that the tributary flow should be added separately to the outflow hydrograph.

Discussion and conclusions
Hydrological flood-routing techniques are widely accepted and are extensively used in engineering practice.The ability to predict the changing magnitude and celerity of a flood wave as it propagates along rivers or through reservoirs makes flood routing important in designing hydraulic structures and in assessing the adequacy of measures for flood protection.However, in practice, only a limited number of gauging stations are available, and even measured or gauged runoff data are frequently unreliable.To establish gauging stations is an expensive task and ongoing maintenance and service costs are also significant.Other studies such as continuous daily hydrograph simulation using duration curves of a precipitation index (Smakhtin and Masse, 2000) and ACRU model rainfallrunoff generating (Schulze, 1995) may be important when the peak flow and volume of the hydrograph are required for local

387
or destination river site analyses.But in routing a flood the other important additional information obtained includes the travel times, the change in magnitude and celerity of flood wave as it propagates along the river.Understanding the hydraulic behaviour of the incoming flood helps to take protective measures prior to any danger that may be caused by the arrival of excess flood.
From the results obtained, it was generally found that the estimated values of the flow variables in the MC-E and MC-X methods resulted in computed hydrographs with acceptable errors for the peak-flow magnitude, peak timing, volume and small RMSE values.In addition, the coefficients of model efficiency (E) were generally close to 1 in most cases, which indicates that the hydrograph shapes are very similar to the observed hydrographs.Hence, it can be concluded that the MC-E and MC-X approaches can be applied in ungauged reaches in the Thukela catchment where observed data sets are unavailable and it is postulated that the methodology would result in similar performance in other catchments.The MC-E and MC-X method also has application in hydrological modelling where inflow hydrographs are simulated by a model and the hydrographs can be simulated in a reach using either the MC-E or MC-X methodology of parameter estimation. 383 Figure 2Selected gauging stations in the Thukela catchment Figure 4 Observed and computed hydrographs inReach-II