4. METHODOLOGY

Toggle Video

This chapter describes the methodologies used for the selection and analysis of streamflow data, as well as the flood routing methods applied to selected events. The Muskingum K and X parameters were calibrated using observed inflow and outflow hydrographs employing the Muskingum (M-Cal) and Three-parameter Muskingum matrix (M-Ma) methods. At ungauged sites, the Muskingum-Cunge equation was utilized, with flow variables estimated using both an empirical (MC-E) approach and an assumed cross-section (MC-X) approach. The computed outflow hydrographs were then statistically and graphically compared with the observed hydrographs. The methodology is detailed as follows:

1. For selected sub-catchments, observed events were chosen, and the outflow hydrographs were computed using calibrated Muskingum K and X parameters (M-Cal) as well as the Three-parameter Muskingum method (M-Ma). These computed hydrographs were then compared to the observed hydrographs to evaluate the performance of the Muskingum method with calibrated parameters,

2. The Muskingum K and X parameters were estimated for each reach based on empirically determined flow characteristics, including flow top width (W), wetted perimeter (P), flow depth (y), hydraulic radius (R), and average velocity (V_av) (Section 4.2.1). This method is referred to as MC-E,

3. The Muskingum K and X parameters were estimated for each reach based on assumed channel cross-sectional dimensions observed in the field, including maximum flow depth (y) and maximum flow top width (W) (Section 4.2.2). This method is referred to as MC-X,

4. The computed hydrographs were compared to the observed hydrographs through statistical and visual analyses, considering flood volume, peak flow magnitude and timing, as well as hydrograph shape (Section 4.3), and

5. Sensitivity analyses were conducted on various catchment variables, including routing reach length (ΔL), routing time step (Δt), Manning roughness coefficient (n), channel geometry, river slope (S), and Muskingum K and X parameters (Section 4.3).

The details of the calculation steps are elaborated upon in the following sections.

4.1 Flood Routing Using Observed Inflow and Outflow

The Muskingum K and X parameters were calibrated using observed inflow and outflow hydrographs, employing the Muskingum-Cunge calibrated (M-Cal) and the Three-parameter Muskingum (M-Ma) methods. The procedure details are discussed in Sections 4.1.1 and 4.1.2.

4.1.1 The M-Cal method

In the Muskingum method of flood routing, the river reach may be subdivided into sub-reaches (Section 2.2). According to the US Army Corps of Engineers (1994a), a general rule of thumb suggests that the computation interval should be less than 1/5 of the time of rise of the inflow hydrograph, which was initially used as an estimate for the routing time interval. Viessman et al. (1989) and Fread (1993) noted that the routing time interval Δt is often assigned any convenient value between the limits of (K/3) <= Δt <= K.

The Muskingum parameter (K) in Equation 2.9, equivalent to the wave travel time in the reach, may be estimated by the lag between the peaks of inflow and outflow hydrographs, as illustrated in Figure 4.1.

The number of sub-reaches was determined by dividing the estimated total travel time (K) by the routing time interval ( Δt).

Hence, the celerity may be estimated by substituting the value of K in Equation 2.9 as:

where

VW	=	celerity [m/s],
K	=	wave travel time [s], and
ΔL	=	reach length [m].

Alternatively, celerity may be estimated from the average velocity, using Equation 4.2 and Table 4.1, assuming a parabolic cross-section.

Rearranging Equation 4.2, the average flow velocity for a parabolic cross-section can be calculated as:

where

V av	=	average velocity [m/s].

Table 4.1 Estimation of celerity for various channel shapes (Viessman et al., 1989).

Channel shape	Manning equations	Chezy equation
Wide rectangular	5/3 Vav	3/2 Vav
Triangular	4/3 Vav	5/4 Vav
Parabolic	11/9 Vav	7/6 Vav

The average velocity calculated from Equation 4.3 was used to calculate cross-sectional flow area for a reference flow, computed as shown in Equation 4.4.

From Equation 2.11 (Section 2.2) the reference flow is estimated as:

where

Q0	=	reference flow [m3/s],
Qb	=	minimum discharge [m3/s], and
Qp	=	peak discharge [m3/s].

Using the continuity equation (Chow, 1959; Linsley et al., 1988; KenBohuslay, 2004 ), the cross-sectional area of flow obtained from Equation 4.4 is as follows:

where

A	=	cross-sectional area [m2].

In a channel where top flow width (W) exceeds mean flow depth by a factor of 20, the mean flow depth (d) can approximate the hydraulic radius (R) (Barfield et al., 1981 cited in SAICEHS, 2001). The hydraulic radius approximated from hydraulic mean depth (d) relationships for various cross-sectional shapes is shown in Table 4.2.

Table 4.2 Hydraulic mean depth (after Chow, 1959 and Koegelenberg et al., 1997)

Cross-section	Hydraulic mean depth (d)
Parabolic section	(2/3) y
Rectangular	y
Triangular	0.5y

In this study, based on field observations, a parabolic cross section was assumed for all reaches.

The mean flow depth (d) for a parabolic cross-section is given by the Euqation 4.5a:

where

d	=	mean flow depth [m], and
y	=	flow depth [m].

The hydraulic radius can thus be estimated from Equation 4.5b:

where

R	=	hydraulic radius for parabolic section [m].

The depth of flow was extracted from rating curves shown in Figures 3.15 to 3.17.

The slope of the river reach was estimated from the Digital Elevation Model (DEM) used by Agro hydrology and Climatology Atlas of South Africa (Schulze et al., 1997). The base roughness value (n_b) was estimated from values contained in Table 2.2 (Section 2.9.3) and based on field observations. corrections for the base roughness value for a particular site were taken from Table 2.3 (Section 2.9.3). Assuming a constant river cross-section along the river reach, the area (A) and the hydraulic radius (R) can be used to calculate the wetted perimeter (P) as shown in Equation 4.6.

where

P	=	wetted perimeter [m].

The formulas for geometrical cross-sectional area of various channels are illustrated in Table A.1 in Appendix A. The top flow width (W) for a parabolic river cross-sectional area may be estimated as shown in Equation 4.7 (Chow, 1959; Koegelenberg et al., 1997):

where

A	=	cross-sectional area [m2], and
y	=	depth of flow [m].

The Muskingum X parameter was calibrated by minimizing the error between the peak discharge of the observed and computed hydrographs.

After the Muskingum K (Figure 4.1) and X parameters were calibrated, the routing coefficients were estimated using Equations 2.6a, 2.6b, and 2.6c (Section 2.1).

As suggested by Viessman et al. (1989), negative values of C₁ must be avoided. Negative values of C₂ do not affect the routed hydrographs. The negative values of C₁ can be avoided by satisfying Equation 2.7 (Section 2.1), which is repeated below:

After the X parameter was determined, the routing time interval was adjusted using relationship given in Figure 2.3 (Section 2.1) (Cunge, 1969; cited by NERC, 1975).

Since Q_t, It, and I_t+1 are known for a given time increment, Q_t+1 is computed using Equation 2.13 (Section 2.2) and repeated for successive time increments to estimate the outflow hydrograph. The lateral inflow per unit length (q) for longer river reaches was estimated using Equation 2.15d (Section 2.2) and added to the computed outflow. It was assumed that the inflow was estimated from a single upstream catchment.

From the continuity equation,

where

Q	=	discharge [m3/s],
A	=	area [m2], and
Vw	=	celerity [m/s].

The area term in Equation 2.15d can be substituted by discharge, as shown in the following Equation 4.8 (Fread, 1998).

From the continuity equation,

The value of β is between 0.5-1 (Fread, 1998). Based on the results obtained, a value of 0.7 was used in this study for Reach-II and Reach-III.

Flow characteristics such as the magnitude and timing of peak flow, hydrograph shapes and flow volume, were statistically compared against the observed events. The K and X parameters were calibrated for each event analysed Using this M-Cal method.

The Three-parameter Muskingum estimation method (M-Ma), outlined in Section 2.3, was also applied.

4.1.2 The M-Ma method

The K, X, and α parameters were estimated using the M-Ma method by treating the entire reach as a single reach. The three coefficients (d₁, d₂, and d₃) in Equation 2.20 were estimated directly using the matrix method. After the primary flow data were reformatted into fixed time steps, matrix inversion was performed for each selected event using the SPSS (Version 11.0.1) statistical software package (SPSS11, 2003).

Using the routing coefficients (d₁, d₂, and d₃) obtained directly from the matrix computation, and since I_t and I_t+1 are known for every time increment, routing is accomplished by solving Equation 2.5 (Section 2.1) for successive time increments.

The three C_i coefficients in Equation 2.5 (Section 2.1) can be derived from the d_i coefficients of Equations 2.20a, 2.20b, and 2.20c, as indicated in Section 2.3. Computations were used to calculate the K, X, and α parameters numerically from Equations 2.21a, 2.21b, and 2.21c (Section 2.3).

4.2 Flood Routing in Ungauged River Reaches

When routing in ungauged catchments, the parameters of the model have to be estimated without observed hydrographs. Inflow hydrographs to downstream catchments can be simulated using a hydrological model such as the ACRU model (Schulze, 1995). However, for this study, the observed inflow hydrographs to the reaches are used.

4.2.1 The MC-E method

In ungauged catchments, where observed inflow and outflow hydrographs are not available, a methodology to estimate these hydrographs has to be derived. As noted by Angus (1987), runoff from rainfall can be estimated using the ACRU distributed model. Therefore, inflow hydrographs for ungauged catchments may be simulated using a hydrological model. The depth of flow and discharge can be derived from empirical relationships recommended by the US Reclamation Service (Chow, 1959), or other flow regime theories suited to rivers with similar characteristics (Section 2.8). Consequently, the Muskingum K and X parameters can be estimated in ungauged river reaches using inflow hydrographs and empirically estimated channel dimensions. In this method, the observed primary flow data obtained from the Department of Water Affairs and Forestry (DWAF) and empirically estimated channel variables are used for the estimation of the Muskingum K and X parameters.

The Muskingum K parameter is estimated from Equation 2.9 as follows(Section 2.2):

For a parabolic cross-section, the celerity (V_w) may be estimated from Table 4.1 as follows:

The average velocity (V_w) is calculated from the Manning Equation 2.32b as follows (Chow, 1959):

The hydraulic radius (R) is estimated from Table 4.2 once the flow depth (y) has been determined using empirical relationships shown in Equation 2.31b.

From Manning's Equation, a section factor can be calculated as shown in Equation 2.32a (Section 2.8) (Chow, 1959):

From Equation 2.11 (Section 2.2) the reference flow is estimated as:

Since the reference flow (Q₀), roughness coefficient (n), and slope (S) are known for the river reach and events under consideration, the unknown variables in Equation 2.32a (Section 2.8) must be estimated using known variables. This can be achieved by either using empirical rules that relate discharge and depth, slope, or by using charts that relate the section factor to the depth of flow (Chow, 1959). The chart is shown in Figure A.1 in Appendix A.

Equation 2.31b relates the wetted perimeter with discharge for natural rivers as follows (Section 2.8):

where

P	=	wetted perimeter [m], and
Q0	=	reference flow [m3/s].

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 (Section 4.1). According to the flow area equations contained in Table A.1 (Appendix A), for a parabolic section, the wetted perimeter (P) can be assumed to be equal to the top flow width (W). Therefore, the top flow width (W) is substituted by the wetted perimeter (P) as shown in Equation 4.9 for a parabolic channel.

When the hydraulic radius is equated with the mean hydraulic depth, then from Table A.1 (Appendix A), width (W) and perimeter (P) can be verified to be equal for a parabolic equation as follows:

Area of parabolic section is computed as:

and for wide parabolic channels, the hydraulic radius may be estimated as

The wetted perimeter is computed as:

Substituting Equations 4.9a and 4.9b into Equation 4.9c result in:

Equation 4.9d can then be substituted in the following equation as:

Then, solving for y in Equation 4.9e, the flow depth (y) may be estimated as:

For the MC-E method, Equation 4.10 was used to estimate the depth of flow in ungauged catchments in order to estimate the Muskingum K and X parameters.

A representative slope of the river reach may be estimated from the river profile as shown in Section 3.4.

The roughness coefficient (n), slope (S), and hydraulic radius (R) can be substituted in the Manning Equation 2.32b to estimate the average velocity (V_av).

The average velocity (V_av) is substituted in Equation 4.2 to calculate celerity (V_w) and then the celerity term is substituted into Equation 2.9 (Section 2.2) to calculate the Muskingum K parameter.

The average velocity is calculated from Equation 2.32b. Manning's equation is then employed to calculate the celerity and flow area of the reference flow (Chow, 1959; Linsley et al., 1988; KenBohuslay, 2004). Alternatively, the area of flow can be estimated using geometrical parameters, as shown in the following equation for a parabolic section:

where

A	=	cross-sectional area [m2],
P	=	wetted perimeter (From Equation 2.29) [m], and
y	=	depth of flow [m].

From Equation 4.9d, the wetted perimeter (P) and top flow width (W) are assumed to be equal.

The top flow width (W), reference flow (Q₀), river slope (S), celerity (V_w) and sub-reach length (ΔL) are substituted in Equations 2.9 and 2.10 (Section 2.2) to estimate the Muskingum K and X parameters as follows:

In this method, Equation 4.10 was used to estimate the channel flow depth.

4.2.2 The MC-X method

Another method of estimating flow depth in ungauged catchments is by developing a rating curve for a selected section in the reaches (Smithers and Caldecott, 1995). Hence, flood routing can be applied using selected channel cross-sections. For a selected channel cross-section observed in the field, a rating curve was developed based on the maximum width and depth relationships.

The assumed cross-section is divided into incremental cross-sectional depths, and the corresponding cumulative area is calculated from geometrical properties of a parabolic channel shape, as shown in Table A.1 (Appendix A). The parabolic shape is used to illustrate the calculation steps.

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 ratio contained in Table 3.2 (Section 3.4) as follows:

where

yi	=	given depth [m],
Wi	=	top width [m],
WMax	=	maximum top width [m] from Table 3.2, and
yMax	=	maximum depth [m] from Table 3.2.

The wetted perimeter for each sub-section is calculated from geometrical equations as contained in Table A.1 (Appendix A). The parabolic section is used as example for showing the calculations:

where

Pi	=	wetted perimeter for given sub-section [m],
Wi	=	top width for sub-sections [m], and
yi	=	given depth for sub-section [m].

The flow area can be estimated from the continuity equation or geometrical properties for a parabolic shape illustrated in Table A.1 (Appendix A).

where

A	=	flow area [m2],
W	=	top flow width [m], and
y	=	flow depth [m].

The hydraulic radius is computed from the cumulative area and wetted perimeter as follows:

Then, the corresponding cumulative discharge is calculated from the Manning equation as follows:

where

A	=	area [m2], and
Q	=	discharge [m3/s].

Since the roughness coefficients (n), hydraulic radius (R), flow area (A) and slope (S) of each reach is known, the discharge can be calculated from Equation 4.15. 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 Equation 4.12 and Table 4.1 respectively. Equations 2.9 and 2.10 are then used to estimate the K and X parameters.

4.3 Model Performance and Sensitivity

It is generally accepted that the output of a hydrological simulation model will not be identical in every aspect to the real system it proposes to represent. However, it is required that the output be sufficiently close to the real system so that the model may be considered to be an acceptable model (Green and Stephenson, 1985).

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 can be employed.

Although the reliability of a hydrological simulation model depends on the quality of input data provided, the accuracy of simulated hydrographs should be assessed by comparing the computed hydrographs against the observed hydrographs (Caldecott, 1989). In addition, the reliability of a runoff estimate made for an ungauged catchment is a function of the reliability of the flood runoff model, the type of the predictive equations and their parameters and coefficients as well as the wisdom and experience of the analyst (US Army Corps of Engineers, 1994a). Hence, the difference in the observed and computed hydrograph are analysed statistically by means of Root-Mean-Square Error (RMSE) and 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 4.16 was used to estimate the actual errors in the computed hydrographs. RMSE computes the magnitude of error in the computed hydrographs (Schulze et al., 1995).

where

RMSE	=	Root-Mean-Square Error [m3/s per event],
Qcomp	=	computed outflow [m3/s], and
Qobs	=	observed outflow [m3/s].

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 Equation 4.17 (Green and Stephenson, 1985):

where

Epeak	=	peak flow error [%],
Qp-comp	=	computed peak flows[m3/s], and
Qp-obs	=	observed peak flows [m3/s].

where

Etime	=	peak time error,
tp-comp	=	time when Qcomp occurs [s], and
tp-obs	=	time when Qobs occurs [s].

where

Evolume	=	peak volume error [%],
Vcomp	=	computed total volume [m3], and
Vobs	=	observed total volume [m3].

Even though the RMSE, Epeak, Etime and Evolume statistics may appear to be reasonable, the shapes of the respective hydrographs may be different. Nash and Sutcliffe (1970, cited by Green and Stephenson, 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 as follows:

where

Qm	=	mean of the observed flows [m3/s], and
n	=	number of data values.

In addition to the limitations of the Muskingum-Cunge method, uncertainties in estimation of catchment parameters such as reach slope, roughness coefficients and geometry of the channel are expected to affect the simulated hydrographs. Therefore, in addition to the results from the application of the methodology discussed in this chapter, sensitivity analysis of catchment parameters such as roughness coefficients, reach slope and reach geometry are contained in Chapter 5.