ABSTRACT Title of Thesis: UNCERTAINTY AND SENSITIVITY ANALYSES OF THE MUSKINGUM-CUNGE CHANNEL ROUTING MODEL Stephen Scott Reid, Master of Science, 2008 Thesis directed by: Professor Richard H. McCuen Department of Civil and Environmental Engineering Two widely used channel routing models are the Muskingum and Muskingum- Cunge (M-C) models. However, no significant research has been conducted on how input parameter uncertainty and model sensitivity affect the accuracy of the downstream hydrograph. Therefore, the goal of this research was to perform uncertainty and sensitivity analyses of the Muskingum and M-C models. The results of the sensitivity analyses show that neither model is particularly sensitive to variations in channel routing input parameters, although the Muskingum method is slightly more sensitive. The results of the uncertainty analyses indicate that in all models uncertainty in the channel routing input parameters introduce only a moderate amount of error in the downstream hydrograph; whereas uncertainty in the upstream hydrograph could introduce significant errors. Finally, the M-C method gave negative outflows in the first few time steps, in agreement with other widely reported results in the literature; however, while others have claimed that negative outflows do not have a significant effect on the downstream hydrograph, the results presented herein suggest that this is not always true. UNCERTAINTY AND SENSITIVITY ANALYSES OF THE MUSKINGUM-CUNGE CHANNEL ROUTING MODEL by Stephen Scott Reid Thesis submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial fulfillment of the requirements for the degree of Master of Science 2008 Advisory Committee: Professor Richard H. McCuen, Chair Associate Professor Kaye Brubaker Professor Allen Davis ?Copyright by Stephen Scott Reid 2008 ii Acknowledgments First and foremost, I wanted to thank Professor Richard McCuen for his inspiration in helping me accomplish my educational goals. I am deeply indebted to him for the numerous hours he spent providing close oversight throughout this research effort. His detailed guidance and comments on earlier drafts of this thesis were invaluable and greatly appreciated. I want to thank Professor Kaye Brubaker for the many contributions she made towards my education at the University of Maryland. I would also like to thank Professor Allen Davis for agreeing to be a member of my thesis committee. I would also like to thank Professor Todd Cooke for his continued guidance and support throughout this process. Finally, this accomplishment would not be possible without the continued support of my entire family, but especially my wife and children. I will forever be grateful for their love and support throughout this process. iii TABLE OF CONTENTS ? CHAPTER I INTRODUCTION ................................................................................. 1? 1.1? IMPORTANCE OF ENVIRONMENTAL MODELS AND UNCERTAINTY ........... 1? 1.2? CHANNEL ROUTING MODELS ......................................................................... 2? 1.3? GOAL AND OBJECTIVES .................................................................................. 3? CHAPTER II LITERATURE REVIEW ....................................................................... 5? 2.1? INTRODUCTION ................................................................................................ 5? 2.2? CHANNEL ROUTING MODEL DEVELOPMENT ............................................... 5? 2.3? THE MUSKINGUM METHOD ............................................................................. 8? 2.4? THE MUSKINGUM-CUNGE METHOD ............................................................. 10? 2.5? ESTIMATING THE ROUTING PARAMETERS ................................................ 12? 2.5.1? Muskingum Routing Parameters ............................................................... 12? 2.5.2? Muskingum-Cunge Routing Parameters.................................................... 15? 2.6? VOLUME CONSERVATION, NEGATIVE PARAMETERS, AND THE ?DIP? EFFECT ............................................................................................................ 19? 2.7? SENSITIVITY AND UNCERTAINTY ANALYSIS .............................................. 21? 2.7.1? Sensitivity Analysis .................................................................................... 21? 2.7.2? Uncertainty Analysis .................................................................................. 22? 2.7.3? Monte Carlo Method .................................................................................. 26? CHAPTER III METHODOLOGY .............................................................................. 28? 3.1? INTRODUCTION .............................................................................................. 28? 3.2? DEVELOPMENT OF UPSTREAM HYDROGRAPH ......................................... 28? 3.2.1? Development of a Design Storm ................................................................ 29? 3.2.2? Generation of Unit Hydrograph .................................................................. 32? 3.2.3? Convolution to Compute Direct Runoff Hydrograph .................................. 35? 3.2.4? Final Upstream Hydrograph ...................................................................... 35? 3.3? DEVELOPMENT OF CHANNEL ROUTING MODELS ..................................... 36? 3.3.1? Selection of Channel Characteristics ......................................................... 36? 3.3.2? Development of the Stage-Discharge Relationship ................................... 38? iv 3.3.3? Muskingum Model Routing Coefficients .................................................... 39? 3.3.4? Selection of Time Step .............................................................................. 40? 3.4? SELECTION OF PARAMETERS FOR EVALUATING MODELING RESULTS 41? 3.5? METHOD FOR EVALUATING PARAMETER SENSITIVITY ............................ 41? 3.6? METHOD FOR EVALUATING PARAMETER UNCERTAINTY ........................ 43? CHAPTER IV RESULTS AND DISCUSSION .......................................................... 46? 4.1? INTRODUCTION .............................................................................................. 46? 4.2? RESULTS OF SENSITIVITY ANALYSIS .......................................................... 46? 4.2.1? Muskingum Model ..................................................................................... 47? 4.2.2? 3-Point Variable Parameter Muskingum-Cunge Model ............................. 52? 4.2.3? 4-Point Variable Parameter Muskingum-Cunge Model ............................. 57? 4.3? RESULTS OF UNCERTAINTY ANALYSES ..................................................... 59? 4.3.1? Muskingum Model ..................................................................................... 60? 4.3.2? 3-Point Variable-Parameter Muskingum-Cunge Model ............................. 66? 4.3.3? 4-Point Variable-Parameter Muskingum-Cunge Model ............................. 72? 4.4? DISCUSSION AND OBSERVATIONS .............................................................. 76? 4.4.1? Comparison of Muskingum and VPMC Models ......................................... 77? 4.4.2? VPMC Method and Negative Outflows ...................................................... 79? CHAPTER V CONCLUSIONS AND RECOMMENDATIONS.................................. 82? 5.1? INTRODUCTION .............................................................................................. 82? 5.2? SUMMARY OF RESULTS ................................................................................ 83? 5.3? CONCLUSIONS ............................................................................................... 84? 5.4? RECOMMENDATIONS AND FUTURE RESEARCH ....................................... 85? CHAPTER VI REFERENCES .................................................................................. 89? 1 CHAPTER I INTRODUCTION 1.1 IMPORTANCE OF ENVIRONMENTAL MODELS AND UNCERTAINTY The use of environmental models to represent natural processes has been widely accepted. Such models can be used to support a wide range of activities, from assisting policy-makers in the development of regulations to providing banks and investors information about land transactions and property liabilities. The accuracy of such models is very important in that it ensures that the results are reliable. The accuracy of a model and hence the reliability of the results depend on many factors, including: (1) the input parameters selected, (2) the uncertainty associated with each of these parameters, (3) the model sensitivity to each of these parameters, and (4) the model calibration. Uncertainty and sensitivity analyses play an important role in developing environmental models because they allow environmental professionals, regulatory reviewers, and policy-makers to better understand the results obtained by use of such models and consequently make more effective decisions. However, these analyses are often omitted due to limited resources or time, and the fact that such analyses are typically not explicitly required and/or emphasized by regulatory agencies. Recently, however, increasing attention has been given to uncertainty analyses, but the importance of this concept is still not firmly rooted in practice. For example, the key U.S. document for water and related land resources development (officially known as the Economic and Environmental Principles and Guidelines for Water and Related Land Resources Implementation Studies, U.S. Water Resources Council, 1983) until recently gave little attention to risk and uncertainty analyses. In spite of the fact that a 2007 act of Congress 2 required that the this key document be revised to include special emphasis on risk and uncertainty analyses, the resulting document still fell short of that goal as evident by the comments made by the American Society of Civil Engineers (ASCE, 2008). 1.2 CHANNEL ROUTING MODELS Hydrologic models are an important class of environmental simulation methods. They are valuable tools for engineering design and analysis of stormwater drainage systems, reservoir design, environmental health assessments, etc. For hydrologic design and analyses on small watersheds and low-order streams, in-channel processes such as channel storage are often unimportant with respect to design. However, as the size of a watershed and the length of the stream channel increase, channel processes become more important and must be taken into consideration in order to ensure accurate results. Channel routing is the method used to account for channel processes such as the resistance to flow and temporary storage because channel characteristics such as the roughness, length, and slope, as well as channel geometry, determine the amount and duration of water storage in the reach when routing a storm hydrograph. These channel characteristics are used in channel routing models to calculate routing parameters. Ultimately, the magnitude and timing of a downstream hydrograph depends on the routing parameters in conjunction with an upstream hydrograph. However, the extent to which parameters realistically reflect channel processes is not well understood. Numerous channel routing techniques have been developed and applied to a wide range of rivers and streams. Two commonly used techniques are the Muskingum and the 3 Muskingum-Cunge (M-C) methods. The M-C channel routing method is a modified version of the Muskingum method and is often used as an alternative to the Muskingum method. Several studies have developed methodologies for estimating the routing parameters for characterizing in-channel processes in the Muskingum and M-C models. However, no significant research has been conducted to date to help understand how uncertainty associated with each input parameter affects the magnitude of the routing parameters. Similarly, significant research is also lacking on the sensitivity of the computed downstream hydrograph to the input parameters. 1.3 GOAL AND OBJECTIVES The goal of this research was to perform uncertainty and sensitivity analyses of the Muskingum and M-C models, with the aim of allowing end-users to better understand the effect of errors in the input parameters on the computed downstream hydrograph. In the selection of a channel routing model, engineers could then use the results of this research to compare the advantages and disadvantages of each model and compare how errors in the input parameters affect the downstream hydrographs of each case. For example, an engineer who designs a levee system for a stream reach should know the effect of uncertainty in the design inputs on the selected design height of the levees. The objectives of this research were: ? Develop a representative upstream hydrograph and channel routing model using commonly applied engineering methods, and quantify the uncertainty associated with each of the input parameters used in that development. 4 ? Conduct sensitivity analyses to understand the relative importance of each of the channel routing input parameters in the model. ? Evaluate the effects of input parameter uncertainty on the routing parameters, as well as on the downstream hydrograph. ? Conduct a sensitivity and uncertainty analysis on the Muskingum and M-C models and compare the results to determine which model is more sensitive to input parameter uncertainty. 5 CHAPTER II LITERATURE REVIEW 2.1 INTRODUCTION This chapter discusses the Muskingum and Muskingum-Cunge (M-C) channel routing methods, both of which rely on the channel routing equation to obtain outflow and storage. The potential uncertainties associated with estimating the various hydrologic parameters needed to calculate the downstream hydrograph are also discussed. 2.2 CHANNEL ROUTING MODEL DEVELOPMENT Considering the number of stream-miles in the United States, it would be impractical to collect discharge information for all of the streams. Furthermore, even if such information were collected, it would be of limited use because the location of the gages involved in such an effort may not always be the location of interest for a particular design project. Therefore, hydrologic models have been developed to simulate water movement through streams. Most, if not all, of these models are based on the St. Venant equations for gradually varied, unsteady open channel flow. The St. Venant equations consist of the continuity and momentum equations. The continuity equation is expressed as: nullnull nullnull null nullnull nullnull null0 (2-1) 6 where A is the cross-sectional area of flow (L 2 ), Q is the discharge (L 3 /T), t is the time (T), and x is the distance along the channel (L). The momentum equation is given by: 1 null nullnull nullnull null null null nullnull nullnull null nullnull nullnull nullnullnull null nullnull null nullnull0 (2-2) in which g is acceleration due to gravity (L/T 2 ), V is the velocity (L/T), h is the depth of flow (L), and S f and S o are friction and bed slopes (L/L), respectively. Channel routing models are distinguished from each other based on the number of terms considered in the momentum equation. Dynamic wave models consider all of the terms in the momentum equation. Models that neglect inertial terms are referred to as diffusion wave models, while models that neglect both inertial and pressure forces are called kinematic wave models. The kinematic wave model is expressed as: nullnull nullnull nullnullnull nullnull nullnull (2-3) in which c is the wave celerity (L/T), and the other terms are as previously defined. The continuity equation expresses the idea that, if at any given point in time, the area of flow at a particular point in space is increasing with respect to time, then the discharge passing this point at that time must be decreasing spatially to account for the water that is stored in the increasing cross section (Chapra, 1997). The momentum 7 equation simply states that the accelerating effect of slope is offset by the effect of friction (Chapra, 1997). The conservation of mass equation is the fundamental equation used in channel routing models. It is expressed as: nullnullnullnull nullnull nullnull (2-4) where I is the inflow or upstream hydrograph (L 3 /T), O is the outflow or downstream hydrograph (L 3 /T), S is the storage (L 3 ), and t is the time (T). This equation is derived by integrating the continuity equation (Eq. 2-1) between the two ends of the channel reach. It is noted that using an averaging approach, the conservation of mass equation can be approximated in numerical form as: nullnull null nullnull null null 2 null nullnull null nullnull null null 2 null null null nullnull null ?null (2-5) where I, O, and S are previously defined, the subscripts 1 and 2 represent successive time steps, and ?t represents the difference in time or the time step (also represented as dt). For a given inflow hydrograph, I 1 , I 2 , S 1 , and O 1 are known or assumed; however, because Eq. 2-5 includes two unknowns (O 2 and S 2 ), a second relationship is needed. The Muskingum and Muskingum-Cunge methods represent two approaches to completing Eq. 2-5. 8 2.3 THE MUSKINGUM METHOD The Muskingum method was developed by McCarthy (1938) for channel control studies in the Muskingum River basin in Ohio. The Muskingum channel routing method is based on the conservation of mass equation (Eq. 2-4) and the storage equation, which is defined as: nullnullnullnullnullnull null null1nullnullnullnullnull (2-6) where x and K are the called the ?routing parameters?, and the other variable are as previously defined. More specifically, x (unitless) represents the spatial weighting factor and reflects the importance of storage between two cross sections. The K parameter (T) is referred to as the storage constant and represents a temporal weighting factor. The Muskingum method solves the channel routing equation by substituting Eq. 2-6 into Eq. 2-5, to obtain: nullnull null nullnull null null 2 null nullnull null nullnull null null 2 null nullnullnullnull null null null1nullnullnullnull null null nullnullnullnullnull null null null1nullnullnullnull null null ?null (2-7) This is simplified by solving for O 2 and combining terms to yield the Muskingum routing equation as follows: 9 null null nullnull null null null nullnull null null null nullnull null null null (2-8) where C 0 , C 1 , and C 2 are dimensionless parameters defined in terms of K and x: null null null nullnullnullnull null 0.5?nullnull nullnullnullnullnull0.5?null nullnull 0.5?null nullnullnull null (2-9) null null null nullnull null 0.5?null nullnullnullnullnull0.5?null null nullnull null 0.5?null null (2-10) null null null nullnullnullnullnull0.5?null nullnullnullnullnull0.5?null null nullnullnullnullnull0.5?null null (2-11) Note that the sum of C 0 , C 1 , and C 2 is 1.0 and outflow at time 2 is thus a weighted average of inflows at times 1 and 2, and outflow at time 1. The routing parameters can be obtained by analysis of measured upstream and downstream hydrographs using the conservation of mass equation (Eq. 2-4). McCuen (1998) also discusses a graphical approach to analyzing the upstream and downstream hydrographs to obtain the routing parameters. This is a trial-and-error approach in which storage is plotted against the weighted flow (xI + (1-x)O) for several x values. The resulting graph will include a series of loops, and the optimal x value is the loop that has the least deviation from a straight line. The K parameter is then estimated as the slope of that line. Given the method of model parameter estimation, it is easy to understand why the derived model will be the most accurate for the given reach and flood event used to calibrate the model, and it may lack predictive capabilities for other reaches and events (Barry and Bajracharya, 1995). 10 While the Muskingum method with calibrated values of x and K can produce reliable results for gaged stream reaches with good datasets, most streams are not gaged and other methods must be employed to estimate the routing parameters. One option would be to calibrate the routing model using data from a nearby gaged catchment and apply the fitted model to the ungaged reach of interest (Tewolde and Smithers, 2006). In this case, synthesis of the downstream hydrograph and the accuracy of the routing parameters will depend on the similarity between the ungaged stream and the characteristics of the stream reach(es) used in estimating the routing parameters (McCuen, 1998). 2.4 THE MUSKINGUM-CUNGE METHOD The Muskingum-Cunge (M-C) method is an improvement on the Muskingum model in that the routing parameters are based on measurable hydrologic data (e.g., channel slope, channel width) for the reach of interest instead of being based on historic discharge data (Ponce et al., 1996). The M-C method uses the same routing equations as the Muskingum method (Eq. 2-8 through Eq. 2-11) to determine the downstream hydrograph; however the routing equation is usually expressed in the finite-difference form: null nullnullnull,nullnullnull nullnull null null null,nullnullnull nullnull null null null,null nullnull null null null,nullnullnull (2-12) where i is the station indicator and t is the time indicator. The M-C eliminates the time- consuming and somewhat arbitrary calibration step (as discussed in Section 2.3) and 11 estimates the routing parameters based on the physical characteristics of the reach of interest. The routing parameters for the M-C method can be calculated and held constant throughout the simulated time period, which is referred to as the Constant Parameter M-C method. Alternatively, the routing parameters can be recalculated at each time step with changes in the wave celerity (c) and discharge (q), which is referred to as the Variable Parameter M-C (VPMC) method. If c and q are held constant, then the K and x parameters will be constant and the computed downstream hydrograph will be the same as that using the Muskingum method for equal routing parameters (McCuen, 1998). Both methods are widely used and have been incorporated into common hydrologic models. The popular TR-20 model developed by the Natural Resources Conservation Service (NRCS) incorporates the Constant Parameter M-C method (Merkel, 2002). Conversely, the watershed model HEC-1 offers several channel routing methods, including the VPMC method (HEC, 1998). The Ponce and Yevjevich (1978) procedure for estimating the routing parameters for the VPMC is presented below. This method uses the reach length, slope, kinematic wave celerity, and unit discharge to calculate the routing parameters: nullnull nullnull nullnull? null null (2-13) nullnull0.5null1null null null nullnull null null null (2-14) nullnull null null (2-15) 12 where c is the kinematic wave celerity (L/T), dq/dh is the stage-discharge relationship, W t is the top width of flow (L), q o is the unit-width discharge (L 2 /T), S o is the channel bed slope (L/L), and L is the length of channel reach (L). In this approach, c is a variable parameter that depends on the unit discharge at each time step. The values of x and K are then used with Eqs. 2-9, 2-10, and 2-11 to compute the three routing parameters (C 0 , C 1 , and C 2 ). . 2.5 ESTIMATING THE ROUTING PARAMETERS As discussed above, routing parameter estimation has been the subject of much research since the development of the Muskingum method. The graphical method for estimating the Muskingum routing parameters and the Ponce and Yevjevich (1978) method for estimating the M-C routing parameters were briefly discussed above. The following discussion describes the Ponce and Yevjevich (1978) method in more detail, as well as presents several other approaches for deriving and optimizing the routing parameters for both the Muskingum and M-C methods. 2.5.1 Muskingum Routing Parameters Gill (1978) recognized that the storage equation used in the Muskingum method [S = k(Ix + (1-x)O)] implies that storage versus weighted flow is a linear relationship. However, this is not always the case, and use of the Muskingum equation may lead to considerable error in situations where this relationship is not linear. Gill (1978) presented two channel routing models based on the Muskingum model to address such non-linear scenarios. Gill?s first model is an exponential model, in which x is determined by trial 13 and error, and then the other parameters are subsequently determined by solving a series of simultaneous equations. Gill?s second method, usually referred to as the Segmented- Curve method, essentially divides the data into a number of segments and then uses a straight-line method to solve for K and x for each segment, thus making the second method simpler to use and the results approximate a curve. Neither one of the two methods focuses on optimizing the routing coefficients, and no quantitative error analyses (only graphical analyses) are performed. Wu et al. (1985) presents a simple approach that improves the accuracy of the Muskingum routing coefficients. Traditionally, the Muskingum routing coefficients were identified graphically by plotting known upstream and downstream hydrographs through the Muskingum equation for several x values. The Muskingum approach created looped plots, and the loop that best approximates a straight line was selected as the x parameter. The K parameter could then be estimated as the slope of the line. Wu et al. (1985) presents an approach in which least square regression is used on the forward and backward paths of the looped plots to determine the slope of each path. A t-test is then used to determine if the respective slopes are statistically different. This procedure is repeated for several x values. The x values are plotted against the t-statistic, and the optimal value of x is selected at the minimum t-value. Wu et al. applied this approach to three previously published examples in which the traditional graphical approach had been used to identify the routing coefficients. The approached presented by Wu et al. (1985) resulted in a 5 to 6-fold reduction in the squared errors compared to the traditional method. 14 Gelegenis and Serrano (2000) present a fully implicit and semi-implicit finite- difference scheme for evaluating the Muskingum routing parameters. The semi-implicit scheme uses linear regression. The storage constant (K) is the same for both methods. The weighting factor for the fully implicit (x f ) method is derived from: null null nullnull null null ?null nullnull 2null (2-16) where x s is the semi-implicit weighting factor and ?t in is the time step used to derive routing parameters. If the time step used in the direct problem (?t di ) is the same as ?t in , the semi-implicit and fully implicit scheme yield identical results. If the two time steps are not equal, the semi-implicit scheme is more accurate and, therefore, preferable. The semi-implicit method also leads to the same results as those obtained from the graphical method. However, a stability analysis shows that the fully implicit scheme is more stable in cases when the time steps for river routing are smaller than the time step used to estimate the routing parameters. It should be noted here that the Gelegenis-Serrano approach is not applicable to the M-C method because it estimates constant routing parameters based on discharge data rather than measurable hydrologic data. Recently, Das (2004) developed a more effective methodology to estimate the Muskingum routing parameters. Traditionally, the Muskingum parameters were estimated by calibrating the model by means of historical inflow-outflow hydrograph data for a known reach; and that calibrated model was then used in the routing equation to predict the outflow at an ungaged reach. Several sources of error were associated with this traditional approach, the most important of which was the fact that the model was 15 calibrated using a different reach than the reach of interest. Das?s (2004) methodology minimizes the outflow prediction errors by using observed flow data, subject to the satisfaction of the routing equations. In his constant parameter-estimation process, he combined the storage models within the routing equation and then minimized the outflow prediction error. To do this, a series of constrained nonlinear optimization models were transformed to unconstrained models using Lagrange multipliers. The governing equations for parameter estimation are then derived using the unconstrained models and the necessary condition equations. Finally, an iterative approach is used to solve the governing equations for parameter estimation. While this methodology results in lower prediction errors, it is not applicable to the M-C method because the M-C method requires information on stream characteristics. 2.5.2 Muskingum-Cunge Routing Parameters Ponce and Yevjevich (1978) present the VPMC method in which the routing parameters K and x are allowed to vary in time and space according to flow variability. Prior to this research, K and x were held constant, which made the results dependent on the value of the reference discharge used to derive the constant parameters. In the VPMC method, the wave celerity (c) and unit-width discharge (q o ) vary with inflow, and they are in turn used to calculate K and x for each time step and cross section. Three methods are employed for calculating c and q o : (1) using a two-point average of the values at grid points (j,n); (2) using a three-point average of the values at grid points (j,n), (j+1,n), and (j,n+1); and, (3) iteration using a four-point average calculation in which the value at (j+1,n+1) are initially obtained by the three-point average method for the first iteration. 16 Results of Ponce and Yevjevich (1978) indicate that the two-point approach did not yield sufficiently accurate results due to a low peak, a slow rate of travel, and significant loss of volume. The three-point average and four-point average approaches both yielded accurate results; however, it is noted that the volume of the outflow hydrograph was not fully conserved in either of these methods compared to the observed inflow. Both of the latter two approaches are widely used and yield similar results. For example, McCuen (1998) conducted an analysis on both methods and indicated that the slight increase in accuracy obtained by using the four-point approach generally did not warrant the additional effort associated with it. Ponce and Changanti (1994) revisited the VPMC method presented by Ponce and Yevjevich (1978). Recognizing that the VPMC method is associated with a small but noticeable loss of volume, Ponce and Changanti attempt to better quantify this volume loss over a wide range of peak inflows using five computational methods: (1) Constant Parameter M-C (CPMC), (2) 3-Pt VPMC (VPMC3), (3) 4-Pt VPMC (VPMC4), (4) a modified 3-Pt VPMC (MVPMC3), and (5) a modified 4-Pt VPMC (MVPMC4). In the MVPMC3 and MVPMC4 methods, the unit-width discharge for each cell is calculated in the same way as in Ponce and Yevjevich (1978), but the average celerity for each cell is calculated by means of the celerity equation and the average unit discharge calculated for that cell. Results of this analysis reveal that the MVPMC3 and MVPMC4 perform slightly better than the VPMC3 and VPMC4 methods, but as much as 3.3 percent volume loss was still observed. Quantitative error analyses were not performed, and optimization of the routing coefficients to minimize volume loss was not discussed. 17 Bravo et al. (1994) used the Schaefer and Stevens technique to determine the M-C routing parameters. The Schaefer and Stevens approach routed a series of inflow hydrographs through a hypothetical channel. The M-C parameters were then estimated using a least squares technique. Physical characteristics involved in this approach include: channel bed slope, channel reach length, channel flow resistance (Manning?s n), baseflow or initial discharge, and inflow peak discharge. The M-C parameters were then used to predict outflow hydrographs, which were compared to the hydrographs obtained for the same stream using the Muskingum Routing option in the HEC-1 program. The results of this research indicated that the Schaefer and Stevens approach for estimating the Muskingum-Cunge parameters yields accurate hydrographs for ungaged streams. It should be noted, however, that although this method is accurate for use on long reach lengths and for large drainage areas, Bravo et al. readily admit that additional research is still needed to test the reliability of this approach for use on short reaches and small drainage areas, which are actually the most common conditions where channel routing methods are applied. It should also be noted that in this approach, K and x are held constant for all timesteps (i.e., this is a Constant Parameter M-C method). Barry and Bajracharya (1995) developed an M-C model that solves the diffusion wave model while eliminating numerical dispersion. The M-C model is generally based on the kinematic wave model, which is the simplest approximation of St. Venant equations, because more complex solutions introduce too much numerical dispersion. The approach developed optimizes the solution based on the Courant number, defined as: 18 null null null null?null ?null (2-17) where c is the celerity, t is the time, and x is the channel length. The Courant number was selected as the criterion for evaluation because the main variables that will affect the accuracy of the results are the spatial and temporal step sizes?the ratio of these parameters is essentially the Courant number. Using Taylor series expansions about the spatial and temporal grid points, Barry and Bajracharya developed a third-order accuracy equation that related the Courant number, the spatial weighting factor (x), and the temporal weighting factor (K), and showing that the optimal Courant number is 1/2, which can then be used to fix the spatial and temporal steps when solving channel routing problems by means of the diffusive wave model. It should be noted that a Courant number of 1 also results in a third-order accurate solution, but under this condition, numerical diffusion is eliminated resulting in a kinematic wave model. Bajracharya and Barry (1997) built on their 1995 research, with the main objective of developing accuracy criteria for the two routing parameters. Using a Taylor Series analysis on the terms in the M-C method, they obtained conditions leading to second-order, third-order, and fourth-order accurate solutions of the diffusion wave equation. The optimal second-order solution occurred when the following two conditions were met: ?null null 2null null (2-18) ?null null null?null null2null ?null (2-19) 19 where D is the diffusion coefficient defined as: nullnull null 2null null null null (2-20) It should be noted that Barry and Bajracharya (1995) and Bajracharya and Barry (1997) criteria are valid only for optimizing routing parameters for the Constant Parameter M-C method. Applying these criteria to the VPMC method is not practical as the celerity changes with each time step resulting in a corresponding different Courant number. 2.6 VOLUME CONSERVATION, NEGATIVE PARAMETERS, AND THE ?DIP? EFFECT As discussed above, it is well documented that the VPMC methods results in a small but perceptible loss of volume. The loss of volume is determined as the difference between the inflow and outflow hydrographs: ?null nullnull nullnullnullnull (2-21) where ?V is the change in volume, and I is the volume of the inflow and O is the volume of the outflow. The inflow and outflow volume is determined by integrating I and O over the hydrograph duration (T) as follows: 20 ?null null null null null null nullnull null null null null null nullnull (2-22) A model is volume conservative when Eq. 2-22 equals zero. The Muskingum model and the Constant Parameter M-C model are volume conservative; however the VPMC is not volume conservative. In addition, and likely related to this finding, is the fact that the VPMC method produces negative spatial weighting parameters (x) and negative outflows during the rising limb of the outflow hydrograph. It is noted that while negative x routing coefficients are irrational in the Muskingum method (McCuen, 1998), negative x values are possible in the Muskingum-Cunge method as x is considered to be a diffusion- matching factor (Ponce, 1994). In the physical sense, a negative x value in the storage equation (Eq. 2-6) would indicate that the outflow contributes more to channel storage than the inflow, which wouldn?t occur unless there is significant backwater?an event not simulated with the Muskingum-Cunge model. Thus, while negative x coefficients are theoretically possible, they are physically not logical. Negative x coefficients result in the occurrence of negative outflows during the first several time steps, often referred to as the ?dip? effect. The dip effect is inherent in the VPMC model, and results in an increased peak discharge in the outflow hydrograph (Tang et al., 1999). Many researchers have investigated the causes of these phenomena and some have developed alternate approaches (e.g., Tang et al., 1999, Todini, 2007); however, wide-spread consensus related to the causes of and solutions to these problems has not been reached. In fact, there is not even agreement related to the significance of the volume losses, negative coefficients, and dip effect?for example, many researchers 21 indicate that the effects are insignificant and do not recommend controlling them (e.g., Perumal and Sahoo, 2008). In addition, the effect of the negative parameters on the accuracy of the computed hydrograph has not been explored. 2.7 SENSITIVITY AND UNCERTAINTY ANALYSIS Environmental models are intended to represent or simulate natural environmental processes, often in mathematical or statistical terms. When developing an environmental model, many questions must be answered including: How will natural processes be simulated (i.e., which model will be used)? How will the input parameters be estimated? And what data will be used? With each step in the model development process, many sources of uncertainty exist, and model sensitivity can affect that accuracy and reliability of the results. Therefore, conducting an uncertainty and sensitivity analyses are an important step in modeling. As discussed above, many studies have focused on optimizing the routing parameters using a variety of complex routines; however, few (if any) studies have focused on the uncertainty associated with each of the elements used to calculate the routing parameters, and how sensitive the results are to changes in each of the input parameters. The following sections provide a background on uncertainty and sensitivity analyses. In addition, variability associated with each of the input parameters used in the channel routing equations are presented. 2.7.1 Sensitivity Analysis Sensitivity may be defined as the degree to which the model outputs are affected by changes in selected input parameters. By investigating the relative sensitivity of the 22 each of the input parameters, a user can become knowledgeable about the relative importance of each of the parameters in the model (USEPA, 2003). The greater the parameter sensitivity, the greater the effect of error in that parameter will have on the computed results (McCuen and Knight, 2006). As such, sensitivity analyses allow users to evaluate the effects of input parameter error on modeling results. After significant literature reviews, it does not appear that any sensitivity analyses have been conducted on the Muskingum and Muskingum-Cunge methods. 2.7.2 Uncertainty Analysis Uncertainty is the term used to describe a lack of knowledge about models, parameters, constants, data, and beliefs (USEPA, 2003). Uncertainties that affect model quality include: (1) uncertainty in the underlying science and algorithms of a model, (2) uncertainty in model parameters and data, and (3) uncertainty associated with the appropriate use and application of a model (USEPA, 2003). Model uncertainty is often referred to as ?reducible uncertainty? because this uncertainty can be reduced with further study and characterization of a given system. With any measured quantity, several sources of random error exist and natural variability ultimately increases model uncertainty. Random errors may be human induced (e.g., measurement error) and/or systematic equipment errors. Natural system variability is also a source of model uncertainty. Variability is inherent in a system and is the term used to describe observed differences attributable to true heterogeneity (USEPA, 2003). Therefore, variability is usually not reducible by further measurement and study, and will always exist in any model. 23 The measured parameters in the M-C method and the VPMC method that are associated with random error and/or uncertainty are: inflow (upstream) discharge (q), length of channel reach (L), top width of flow (W t ), depth of flow (h), and channel slope (S o ). Note the kinematic wave celerity (c) and the unit-width discharge (q o ) are not included in this list because they are calculated using other measured parameters (e.g., discharge, depth, and top width of flow). The statistical measure used to describe the uncertainty of each parameter is the coefficient of variation (C v ): where s is the sample standard deviation and nullnull is the sample mean. The C v is a measure of the degree of variation in a data set. 2.7.2.1 Length of Channel Reach (L) The length of the channel reach is the distance along a stream channel between two points of interest. For channel routing, the two points of interest are the locations of a known upstream hydrograph and the location of the desired downstream hydrograph that will become the model output. Many methods are available for calculating channel length including measurements from USGS topographic maps, computer-aided methods using geographic information systems (GIS) software, and field surveys. The reproducibility of measurement will depend on the method used, as well as user error, and the accuracy of the maps and/or equipment. The U.S. Water Resources Council (USWRC, 1981) attempted to describe the variability in channel length measurements for null null null null nullnull (2-23) 24 use in two methods used to estimate peak flow frequencies in ungaged watersheds. The C v for the mean channel length (averaged over both methods) was about 0.14, with an average maximum C v (averaged over both methods) of about 0.42. Thus channel length has low to medium measurement variability (USWRC, 1981). 2.7.2.2 Channel Slope (So) The channel bed slope is the difference in elevation between the two points of the channel divided by the channel length. It describes the steepness of the reach. The methods for determining channel slope are similar to those used for determining channel length and include using USGS topographic maps, GIS software, and field surveys. The U.S. Water Resources Council (USWRC, 1981) described the variability in channel slope measurements for use in four methods used to estimate peak flow frequencies in ungaged watersheds. The mean C v for channel slope (averaged over the four methods) was about 0.22, with an average maximum C v of about 0.7. Based on these results, measured channel slopes have medium variability (USWRC, 1981). 2.7.2.3 Stream Discharge:Upstream Hydrograph (Q or q) Channel routing procedures require that a discharge hydrograph be obtained for the upstream point of the stream channel (referred to as the ?inflow?). Many methods for obtaining stream discharge hydrographs are available including stream gages, computer simulation modeling, and obtaining data from nearby, similar-order streams. Many sources of uncertainty and variability are present with all of these methods. Several variables determine the shape of a hydrograph including regional geology (e.g., 25 mountains, rural plains) and the land use such as the amount of urban development. The Soil Conservation Service (SCS) approach to developing hydrographs uses the peak rate factor (PRF) to control the volume of water on the rising and receding limbs of the hydrograph. The SCS default PRF is 484 (SCS, 1992), which assumes that about 37.5 percent of the runoff volume is under the rising limb of the hydrograph and 62.5 percent is under the receding limb. It is noted that the default PRF of 484 represents typical or average conditions as determined by the SCS, but mountainous regions can have higher PRFs while flat rural areas can have much lower PRFs (NWS, 2005). Thus uncertainty and variability are associated with the selection of the PRF used to generate the inflow hydrograph. Very little information is available on the uncertainty of the PRF. McCuen (1989) presents a range of PRFs for each county in Maryland. Data for Fredrick County were used in this study, and the standard deviation for the PRF is approximately 86. Thus, with a mean of 449, the C v is 0.19. 2.7.2.4 Manning?s n Manning?s roughness coefficient (n) is very difficult to accurately measure and values are usually obtained from tables. McCuen (1998) presents the range of n values for several types of natural stream channels. Assuming the mid-point of the range approximates the average, the C v varied from about 0.11 to 0.33, with an average of about 0.20. 26 2.7.2.5 Width of Flow (W t ) Similar to depth of flow, the width of flow may be estimated in a variety of ways. It may be physically measured, and a relationship can be developed for width of flow and discharge. Width may also be calculated using the discharge continuity equation (Q = AV) for a given discharge, velocity, and height of flow for an assumed channel profile (e.g., rectangular, trapezoidal). The U.S. Fish & Wildlife Service (USFWS, 2003) estimated a C v of 0.104 for width of flow at USGS gage stations in the Maryland Coastal Plain region. 2.7.3 Monte Carlo Method Monte Carlo analysis uses statistical sampling techniques to obtain a probabilistic approximation to the solution of a mathematical equation or model (USEPA, 1997). A goal of Monte Carlo analysis is to quantitatively characterize the uncertainty and variability in model estimates (Sobol, 1994). A single method for conducting a Monte Carlo analysis is not recommended. Instead the term refers to a widely-used probabilistic approach for conducting uncertainty and sensitivity analyses that follow a particular pattern. To conduct a Monte Carlo analysis, input parameters are assigned a probability density function?or statistical distribution (e.g., normal, log-normal, triangular)?with the distribution and standard deviation of the probability density function based on the uncertainty associated with the parameter. Selecting the distribution is an important part of the analysis as the shape of the distribution can greatly affect the outcome (USEPA, 1997). Once each input parameter has been assigned a probability density function, a computer algorithm is used to repeatedly run the model with randomly selected input 27 values based on the defined probability density functions. The model is generally run hundreds or thousands of times using different combinations of input parameters based on the defined distributions of each input parameter. Each time the model is run, the output value is saved. After all of the computer simulations are finished, the output values are analyzed and descriptive statistics and probability plots can be created to describe the likelihood of a particular outcome occurring. 28 CHAPTER III METHODOLOGY 3.1 INTRODUCTION The following sections present the approach to developing the inflow hydrograph, the selection of the model input parameters, and the scenarios evaluated using the channel routing models. The process for evaluating model sensitivity and uncertainty are also discussed. 3.2 DEVELOPMENT OF UPSTREAM HYDROGRAPH A synthetic hydrograph, which was used as the upstream hydrograph, was generated rather than using real streamflow data so that the hydrograph shape and peak flow could be controlled to evaluate the effects on the channel routing coefficients and characteristics of the downstream hydrograph. Also, real streamflow data are usually not available to the design engineer for a site of interest. Therefore, design engineers most often use a synthetic hydrograph based on watershed and channel properties. As discussed in Chapter 1, one of the objectives of this research was to use commonly applied engineering methods to develop a model and evaluate the uncertainty associated with all of the input parameters. The upstream hydrograph was developed using the U.S. Soil Conservation Service (SCS) unit hydrograph method. Prior to selecting the unit hydrograph approach, other methods for developing a synthetic hydrograph were considered. For example, the use of a gamma distribution hydrograph was considered because the gamma probability distribution function nicely mimics the shape of a real stream hydrograph. In addition, 29 variations of the gamma function hydrograph used in previous channel routing research (e.g., Tang et al., 1999) were also considered. However, the unit hydrograph approach was preferred because it is widely used by hydrologic engineers and is a component of commonly used watershed models (e.g., TR-20, HEC-1). In that context, it meets the first objective of this research. Following the procedure outlined in McCuen (1998), the SCS unit hydrograph was used in a four-step process to develop a direct runoff hydrograph. First, a design storm was computed using the Baltimore, MD, intensity-duration-frequency (IDF) curve. Second, a unit hydrograph was generated using the SCS triangular unit hydrograph method. Third, the precipitation excess hyetograph and the unit hydrograph were transformed into a direct runoff hydrograph using a process called convolution. Finally, baseflow was added into the direct runoff hydrograph to create the final upstream hydrograph which was then used in the research analyses. Each of these steps is discussed in detail below. 3.2.1 Development of a Design Storm To develop a design storm, the duration, volume, and distribution of precipitation for a given storm event are needed. To estimate the volume of precipitation (or depth of rainfall), the intensity of a given precipitation event was estimated using the Baltimore, MD, rainfall IDF curve (McCuen, 1998). For this analysis, the precipitation event was assumed to be 10-year, 6-hour storm, which has an average intensity of 0.6 inches/hour. The total depth of precipitation was then calculated as: 30 nullnullnullnullnull (3-1) where Q is the depth of precipitation (L), i is the precipitation intensity (L/T), and D is the duration of the storm event (T). The distribution of rainfall throughout the event was assumed to fit a triangular distribution with the peak precipitation intensity calculated as: null null null2null null null null2nullnull (3-2) where i p is the peak precipitation intensity (L/T). The precipitation hyetograph was then determined as follows: nullnullnullnull null nullnullnull null null 2 null for 0 ? t ? D/2 (3-3a) and nullnullnullnull null2nullnull null null1null null null null for D/2 < t ? D (3-3b) where t is the time since the onset of the precipitation event (T). The rainfall hyetograph was then separated using a loss function to estimate the amount of precipitation that was retained in the watershed and the amount of precipitation excess that became direct runoff. The constant percentage method was selected to estimate the amount of precipitation that was converted to direct runoff. This method uses a constant rate (or percent) and assumes that losses are proportional to the 31 rainfall intensity (McCuen, 1998). This approach is similar to the Rational method for estimating peak discharge, which assumes the peak discharge rate is equal to a fraction of the rainfall intensity. In the Rational method, the fraction of rainfall intensity converted to runoff is called the runoff coefficient (C). A Rational method C value of 0.4 represents runoff generated from residential lots (McCuen, 1998). Since the approach for calculating runoff using the constant percentage method is similar to the Rational method, the assumption that 40 percent of rainfall was converted to runoff (i.e., 60 percent of the rainfall was lost to watershed storage) is reasonable and represents runoff from typical residential lots. Using the constant percentage method, direct runoff was calculated by multiplying each of the precipitation ordinates by 0.4. The precipitation hyetograph and the precipitation excess function are shown in Figure 3-1. Figure 3-1. Rainfall Hyetograph Used to Generate Runoff Hydrograph. 1.40 1.20 1.00 0.80 0.60 0.40 0.20 0.00 0123456 Precipitation ? Intensity ? (in/hr) Duration?(hrs) Precipitation Precipitation?Excess 32 3.2.2 Generation of Unit Hydrograph The next step was to develop a generic unit hydrograph using the SCS unit hydrograph method. The unit hydrograph method is a practical approach to developing a storm hydrograph for ungaged streams, and it is based on evaluations of a large number of watersheds. This approach develops a generic hydrograph that results from one inch of precipitation excess evenly distributed over the entire watershed. Two parameters must be determined to develop the unit hydrograph: the peak discharge and time to peak. The SCS equation for calculating the peak discharge for the unit hydrograph is: null null null 484nullnull null null (3-4) where q p is the peak discharge rate in cubic feet per second (ft 3 /sec), A is the drainage area in square miles (mi 2 ), Q is the runoff depth (assumed to be one inch for unit hydrograph), and t p is the time to peak in hours (hrs). The unit hydrograph was developed for a drainage area of 30 mi 2 . The unit hydrograph historically has been used on much larger watersheds (e.g., USACE, 1994), and it is well suited for drainage areas less than 100 mi 2 (Dunne and Leopold, 1978). The 30 mi 2 area allows for a slightly larger watershed and the use of a longer channel length to better evaluate the significance of in-channel processes. The time to peak was calculated as: null null null 2 3 null null (3-5) 33 where t c is the time of concentration (T), defined as the time for runoff to travel from the most hydrologically distant point in the watershed to the point of interest (NRCS, 2007a). The t c was calculated as follows. null null null null null (3-6) where null null 209null null.null (3-7) nullnullnullnull null.null (3-8) In Eqs. 3-6 through 3-8, t c is the time of concentration, L is the watershed length (ft), A is the drainage area size (acres), V is the velocity of runoff (ft/sec), k is constant representing land cover (unitless), and S is the slope of the watershed (ft/ft). The equation for relating the length of the watershed to the drainage area was obtained from the Soil Conservation Service (1973). To calculate the velocity, the velocity method was used. A value of 7 was used for k, which represents short grass pasture (McCuen, 1998), and the watershed slope was assumed to be 0.05 ft/ft. Using these values, the time to peak was calculated to be 11 hrs. Figure 3-2 shows the curvilinear unit hydrograph and the associated triangular hydrograph. The triangular unit hydrograph method was used to simplify computations and because the research objectives will not be influenced by the shape of the upstream hydrograph. Additionally, the convolution process will smooth the triangular shape to a more conventionally shaped hydrograph. 34 Figure 3-2. SCS Curvilinear Unit Hydrograph and Equivalent Triangular Unit Hydrograph (NRCS, 2007a). Using the calculated time to peak and peak discharge rate, the triangular unit hydrograph was generated using the following equation: nullnullnullnull null nullnullnull null null null for 0 ? t ? t p (3-9a) and null nullnullnull nullnull null null1 null nullnullnull null null null null for t p < t ? t p +t r (3-9b) where q(t) is the discharge rate (ft 3 /s) at time t (hrs), t r is the time to recede, and q p and t p are as previously defined. 35 3.2.3 Convolution to Compute Direct Runoff Hydrograph The process of converting the unit hydrograph into a direct runoff hydrograph that can be used as the design storm in the numerical experiments is called convolution. It is a process of multiplication, translation with time, and addition that is represented as: nullnullnullnull nullnullnullnullnullnull null null nullnullnullnullnullnullnullnull (3-10) where U(t) is the time distributed unit hydrograph, y(t) is the time distributed direct runoff, x(?) is the time distributed rainfall excess, and ? is the lag time between the beginning of rainfall excess and the unit hydrograph. 3.2.4 Final Upstream Hydrograph The final step to developing the upstream hydrograph was to add baseflow. A constant baseflow of 100 ft 3 /sec was added to each of the hydrograph ordinates. As previously discussed, channel routing is used to evaluate how in-channel processes affect the downstream hydrograph during a storm event. In the experiments conducted in this research, a steady-state baseflow was assumed and thus baseflow does not affect the channel routing model. Therefore, a relatively small baseflow compared to the peak discharge rate was used. The final upstream hydrograph used in the Muskingum and Muskingum-Cunge channel routing experiments is presented in Figure 3-3. 36 Figure 3-3. Final Upstream Hydrograph Used in Channel Routing Experiments. 3.3 DEVELOPMENT OF CHANNEL ROUTING MODELS Three commonly used and well document models were evaluated?the Muskingum model, the 3-Point VPMC model, and the 4-Point VPMC model. As discussed in Chapter 2, the basic input parameters for the three routing methods are the same. The difference between the methods is how the x and K routing coefficients are calculated. The following subsections present the development of the channel routing models. 3.3.1 Selection of Channel Characteristics The basic characteristics used in the routing models are the length of the reach, the channel slope, Manning?s roughness coefficient (n), and the channel geometry. The values selected for each of these input parameters are shown in Table 3-1. It should be noted that some of the values selected for the channel properties were selected using an 0 500 1000 1500 2000 2500 0 5 10 15 20 25 Discharge ? (cfs) Time?(hrs) 37 iterative process. The objective was to select values that did not result in significant volume loss due to negative coefficients (refer to Section 2.6) while still resulting in perceptible channel storage (as indicated by a reduced downstream peak discharge rate). For example, a channel reach that is too long results in a significant loss of volume; however, if the channel reach is too short the difference between the peaks of the upstream and downstream hydrographs will be minimal. A channel length of 9000 ft enabled both of these problems to be avoided. Table 3-1. Channel Properties Used in Channel Routing Models Channel Property Value Selected Length of Reach 9000 ft Channel Bed Slope 0.05% Channel Geometry Rectangular channel?100 ft bottom width; trapezoidal channel?100 ft bottom width and 3:1 (H:V) side slope Manning?s n 0.05 One task of the research was to model conditions similar to past research so that comparisons could be made. The channel bed slope and Manning?s n were selected to be consistent with previous research (e.g.,Tang et al., 1999; Perumal and Sahoo, 2008) and to represent a natural channel with some weeds/brush. These values enable sufficient attenuation to be realized without a significant loss of volume. Two channel geometries were evaluated?a rectangular channel and a trapezoidal channel. These two channel shapes are widely modeled in engineering applications. For large watersheds, a rectangular channel is often assumed as flow on the bank area is minimal except during extreme floods. The channel bottom width was assumed to be 38 100 ft based on the magnitude of the inflow hydrograph and the assumption that the ratio width of flow:depth of flow should be about 20:1 or greater for most time steps of the inflow hydrograph. This ratio is the approximate mid-point of the width:depth range used for classifying streams (NRCS, 2007b). For the trapezoidal channel, a side slope (z) of 3:1 (horizontal:vertical) was used. Using both a rectangular and a trapezoidal channel for the three routing methods resulted in six different scenarios. 3.3.2 Development of the Stage-Discharge Relationship The VPMC method uses the stage-discharge relationship to calculate the celerity (c) at each time step as shown in Eq. 2-13 and re-presented in Eq. 3-11: nullnull nullnull nullnull? null null (3-11) where c is the kinematic wave celerity (L/T), dq/dh is the derivative of the stage- discharge relationship, and W t is the top width of flow (L). To develop the stage- discharge relationship, Manning?s equation was used. Manning?s equation for discharge in a rectangular channel: nullnull 1.49 null nullnullnullnull null null null null null null nullnull null null null null null 1.49 null null null null null null null nullnull null nullnull null null null nullnull null 2nullnull null null null (3-12) where n is Manning?s coefficient (unitless), s is the channel bed slope (L/L), R h is the hydraulic radius (L), w is the channel width (L), and h is the depth of flow (L). 39 The stage-discharge relationship was obtained by taking the derivative of the discharge (q) with respect to depth of flow (h) as shown in Eq. 3-13 and 3-14: nullnull nullnull null 1.49null null null null null null null null ? nullnullnull (3-13) where nullnullnullnull null nullnull null nullnull null null null nullnullnull2nullnull null null null (3-14) Taking the derivative of Eq. 3-13 yields the following stage-discharge relationship: nullnull nullnull null 1.49null null null null null null null 5 3 nullnullnullnull2nullnull null null null null nullnullnullnullnull null null null null 4 3 nullnullnull2nullnull null null null null null nullnullnullnullnull null null null nullnullnull2nullnull null null null (3-15) Similarly, the stage-discharge relationship for a trapezoidal channel is: nullnull nullnull null 1.49 null null null null null null null 5 3 nullnull null 2nullnull1nullnull null null null null null null null null null nullnullnull null nullnull null null null null null nullnullnull2nullnullnull null 4 3 nullnullnull null nullnull null null null null null nullnull null2nullnull1nullnull null null null null null null null null null null null1nullnull null null null null null nullnull null 2nullnull1nullnull null null null null null null null null null (3-16) 3.3.3 Muskingum Model Routing Coefficients The Muskingum model routing coefficients (x and K) are usually obtained by analysis of measured upstream and downstream hydrographs. In the absence of upstream and downstream data to calibrate values, the x coefficient is usually assumed. A value of 0.2 was selected for x, which is based on the midpoint of the range of x values for most 40 streams (Dunne and Leopold, 1978). The K coefficient was estimated as the average travel time through the reach. 3.3.4 Selection of Time Step Selecting the proper time step (dt or ?t) to route the flood wave to a downstream reach is a critical aspect of developing a channel routing model. Using a very small time step does not necessarily result in more accurate results (Bajracharya and Barry, 1997). In fact, a small time step can result in negative coefficients and a loss of volume. Section 2.5 presents previous research that has shown that the Courant number can be used to fix the time step relative to the channel length. However, previous research has focused on optimizing channel routing for the constant parameter M-C method?not the VPMC method as used in this study. It would be difficult to apply these criteria to the VPMC method because as the celerity changes with each time step, the Courant number also changes. Ponce (1994) recommends that at a minimum, dt meet the criteria in Eq. 3-17 to ensure that a sufficient number of data points are represented on the hydrograph: null null nullnull null5 (3-17) The t p of 11 hours was selected for the watershed. Based on the above recommendation, a dt of 1 hr was selected for the channel routing models because this allows for a sufficient number of time steps without causing negative coefficients or a significant loss of volume. In addition, a dt of 1 hr is consistent with the dt used to generate the excess rainfall hyetograph. 41 3.4 SELECTION OF PARAMETERS FOR EVALUATING MODELING RESULTS Prior to conducting the sensitivity and uncertainty analyses, it was necessary to determine which parameters would be used to evaluate the effect of uncertainty in the routing models. Two parameters selected?the outflow peak discharge (Q po ) and the outflow time to peak (t po ). These are the main parameters that define a hydrograph shape, and these two parameters were used in previous research to evaluate channel routing models. 3.5 METHOD FOR EVALUATING PARAMETER SENSITIVITY Model sensitivity was determined by varying individual input parameters and evaluating the effect on the outflow peak discharge and time to peak. The Muskingum and VPMC input parameters were adjusted by ?20 percent as shown in Table 3-2. It should be noted only Muskingum and VPMC direct input parameters were assessed in the sensitivity analysis. Secondary parameters?such as those that were used to develop the inflow hydrograph (e.g., drainage area, watershed slope)?not directly used in the channel routing models were not assessed in the sensitivity analysis; however, they were included in the uncertainty analysis. 42 Table 3-2. Parameters Evaluated in the Sensitivity Analysis Channel Property Input Value Values Used in Sensitivity Analysis Length of Reach 9000 ft 7,200 ft, 10,800 ft Channel Bed Slope 0.0005 ft/ft 0.0004 ft/ft, 0.0006 ft/ft Channel Bottom Width 100 ft 80 ft, 120 ft Manning?s n 0.05 0.04, 0.06 Channel Side Slope (trapezoidal channels only) 3 2.4, 3.6 X routing coefficient (Muskingum method only) 0.2 0.16, 0.24 After each parameter was varied, the outflow peak discharge and time to peak was recorded and the relative sensitivity (R s ) was calculated. Relative sensitivity is a useful metric because it is dimensionless and therefore it can be used to compare parameter sensitivities (McCuen, 2003). Relative sensitivities were calculated as follows: null null null ?null nullnull null nullnull null ?null null null null null (3-18) where ?Q po is the change in the peak outflows, and ?F i is the change in the input parameter being evaluated. The R s is essentially the ratio of change in the output relative to the change in the input parameter. A negative R s indicates an inverse relationship between the input parameter and the Q po ?i.e., as one increases the other decreases. 43 3.6 METHOD FOR EVALUATING PARAMETER UNCERTAINTY Monte Carlo simulation techniques were used to evaluate the effect of uncertainty in input parameters on the downstream hydrograph. As discussed in Section 2.7.3, input parameters are assigned a probability density function based on the uncertainty associated with that parameter. Once each input parameter has been assigned a probability density function, a computer algorithm is used to repeatedly run the model with randomly selected input values based on the defined probability density functions and each time the model is run the output value is saved. Ten thousand Monte Carlo simulations were run for each model scenario evaluated in the study. Ten thousand simulations was selected for the Monte Carlo analyses because the histogram was more defined compared to the resulting histogram using less simulations (e.g., 1000 or 5000); however, there was no significant change in the histogram when using more simulations (e.g., 15,000). Note that in the Monte Carlo simulations, all the input parameters were simultaneously varied during each of the 10,000 simulations for each scenario. Two uncertainty scenarios were evaluated for each model. In the first scenario, only the parameters that were used in the channel routing models were varied to evaluate potential outflow hydrographs due to input parameter uncertainty. This is an extension of the sensitivity analysis. In the second scenario, all of the user-defined parameters, including those used to develop the upstream hydrograph, were varied because the goal of the uncertainty analysis was to evaluate how uncertainty in each of the parameters used to develop the inflow hydrograph and channel routing model potentially magnify the error in the outflow hydrograph?i.e., uncertainty propagation. The distribution 44 functions used in the Monte Carlo simulations for each of user-defined input parameters are shown in Table 3-3. Table 3-3. Input Parameter Distribution Functions Used in Uncertainty Analysis Channel Property Mean Coefficient of Variation Distribution Channel Routing Direct Input Parameters Length of Reach (ft) 9000 0.14 Log-normal Channel Bed Slope (ft/ft) 0.0005 0.22 Log-normal Channel Bottom Width (ft) 100 0.104 Normal Manning?s n 0.05 0.2 Normal Channel Side Slope (trapezoidal channels only) 3 0.5 Normal x routing coefficient (Muskingum method only) 0.2 0.5 Normal Secondary User-Defined Input Parameters Watershed Size (mi 2 ) 30 0.025 Log-normal Watershed Slope (ft/ft) 0.05 0.56 Normal Peak Rate Factor 484 0.19 Normal Rational C 0.4 0.2 Normal Velocity Method K 7 0.2 Normal Rainfall Intensity (in/hr) 0.6 0.2 Normal The derivation of the C v s for channel length, slope, Manning?s n, width of flow (or bottom width), and the PRF were presented in Section 2.7.2. The derivation of the other input parameters are briefly discussed below: ? x routing coefficient: The C v of 0.5 is based on the average value of 0.2 and the realistic range for x of 0.1 to 0.3. This assumes the realistic range is within one standard deviation of the mean. 45 ? Channel side slope: The C v of 0.5 is based on the average value of 2 and assumes that the standard deviation for side slope is 1. ? Watershed area: The C v of 0.025 was obtained from USWRC (1981) which indicates that there is very low uncertainty associated with estimating watershed area. ? Watershed slope: The C v of 0.56 was also obtained from USWRC (1981). These results indicate there is moderate variability associated with estimating watershed slope. ? Rational C: The C v of 0.2 was based on the results from USWRC (1981). Thus Rational C has low to medium estimation variability. ? Rainfall Intensity: The C v of 0.2 was based on the results from USWRC (1981). Rainfall intensity also has a low to medium measurement variability (USWRC, 1981). ? Velocity K: Velocity K is very difficult to accurately measure as it requires accurate measurement of both Manning?s n and the hydraulic radius (R h ). Velocity K is calculated as: nullnull 1.486 null null null null null null (3-19) The K value of 7 used in this research was based on an R h of 0.04 and an n of 0.025 (McCuen, 1998). A mini-sensitivity analysis was conducted by varying R h and n at known intervals and recording the change on K. The standard deviation was then estimated to be approximately 1.5, resulting in a C v of 0.2. 46 CHAPTER IV RESULTS AND DISCUSSION 4.1 INTRODUCTION In this chapter, the sensitivity and uncertainty associated with each input parameter used in the channel routing models were investigated using the methods and procedures presented in Chapter III. The results of the sensitivity and uncertainty analyses are presented for the Muskingum model, the 3-Pt VPMC model, and the 4-Pt VPMC model. Finally, a general discussion is included in which the results of the different models are compared, the negative outflows are further evaluated, and the ability to reduce outflow uncertainty is examined. 4.2 RESULTS OF SENSITIVITY ANALYSIS Table 4-1 includes the relative sensitivities (R s ) for each of the channel routing models evaluated based on changing one input parameter and evaluating the effect on the outflow peak discharge (Q po ). Note that this table does not present the sensitivity analysis results for the outflow time to peak (t po ) as this parameter was found not to be sensitive in any of the models because the changes in t po were small relative to ?t. The sensitivity of t po is discussed in more detail in each of the subsections below. The results of the sensitivity analyses indicate that both the Muskingum and VPMC models are not particularly sensitive to errors in the input parameters. While the Muskingum model is slightly more sensitive than the VPMC, all of the R s values were less than 0.1, meaning that Q po is not very sensitive to changes in any of the input parameters. 47 Table 4-1. Summary of Sensitivity Analyses as Determined by the Relative Sensitivity 1 (P = model parameter; S = bed slope; W = bottom width; L = channel length; z = side slope) P Musk. (Rect.) Musk. (Trap.) 3-Pt VPMC (Rect.) 3-Pt VPMC (Trap.) 4-Pt VPMC (Rect.) 4-Pt VPMC (Trap.) S 0.034 0.033 0.029 0.032 0.029 0.032 n -0.064 -0.065 -0.027 -0.031 -0.027 -0.031 W -0.030 -0.021 0.0055 0.0028 0.0055 0.0028 L -0.088 -0.094 -0.014 -0.019 -0.014 -0.019 z -- -0.013 -- -0.0035 -- -0.0035 x 0.028 0.031 -- -- -- -- 1 The R s presented is the average R s when the given parameter was decreased and increased by 20 percent, respectively. Refer to Tables 4-2 through 4-7 for individual results. 4.2.1 Muskingum Model Tables 4-2 and 4-3 present the sensitivity analyses results for the rectangular and trapezoidal channels, respectively. For Q po , the results indicate that the rankings of the parameter importance for the Muskingum rectangular channel model are (R s ): channel length (-0.088), Manning?s n (-0.064), channel bed slope (0.034), bottom width (-0.030), and x coefficient (0.028). The relative parameter importance for the Muskingum trapezoidal channel model are (R s ): channel length (-0.094), Manning?s n (-0.065), channel bed slope (0.033), x coefficient (0.031), bottom width (-0.021), and side slope (- 0.013). The R s for t po is discussed below. 48 Table 4-2. Sensitivity Analysis Results for Muskingum Method?Rectangular Channel (S = bed slope; W = bottom width; L = channel length) Run S (ft/ft) n W (ft) L (ft) x Q po (cfs) ?Q po R s t po (hrs) ?t po R s 1 0.0005 0.05 100 9000 0.2 2131 NA NA 13 NA NA 2 0.0005 0.05 100 7200 0.2 2169 1.8% -0.088 13 0.0% 0.0 3 0.0005 0.05 100 10800 0.2 2094 -1.8% -0.088 14 7.7% 0.38 4 0.0005 0.05 80 9000 0.2 2144 0.6% -0.030 13 0.0% 0.0 5 0.0005 0.05 120 9000 0.2 2118 -0.6% -0.031 13 0.0% 0.0 6 0.0004 0.05 100 9000 0.2 2114 -0.8% 0.040 13 0.0% 0.0 7 0.0006 0.05 100 9000 0.2 2143 0.6% 0.028 13 0.0% 0.0 8 0.0005 0.04 100 9000 0.2 2157 1.2% -0.062 13 0.0% 0.0 9 0.0005 0.06 100 9000 0.2 2103 -1.3% -0.066 14 7.7% 0.38 10 0.0005 0.05 100 9000 0.16 2120 -0.5% 0.027 13 0.0% 0.0 11 0.0005 0.05 100 9000 0.24 2143 0.6% 0.028 13 0.0% 0.0 Table 4-3. Sensitivity Analysis Results for Muskingum Method?Trapezoidal Channel (S = bed slope; W = bottom width; L = channel length; z = side slope) Run S (ft/ft) n W (ft) L (ft) z x Q po (cfs) ?Q po R s t po (hrs) ?t po R s 1 0.0005 0.05 100 9000 3 0.2 2108 NA NA 13 NA NA 2 0.0005 0.05 100 7200 3 0.2 2155 2.2% -0.11 13 0.0% 0.0 3 0.0005 0.05 100 10800 3 0.2 2075 -1.5% -0.077 14 7.7% 0.38 4 0.0005 0.05 80 9000 3 0.2 2118 0.5% -0.024 13 0.0% 0.0 5 0.0005 0.05 120 9000 3 0.2 2101 -0.3% -0.017 14 7.7% 0.38 6 0.0004 0.05 100 9000 3 0.2 2095 -0.6% 0.030 14 7.7% -0.38 7 0.0006 0.05 100 9000 3 0.2 2123 0.7% 0.037 13 0.0% 0.0 8 0.0005 0.04 100 9000 3 0.2 2142 1.6% -0.082 13 0.0% 0.0 9 0.0005 0.06 100 9000 3 0.2 2087 -1.0% -0.049 14 7.7% 0.38 10 0.0005 0.05 100 9000 2.4 0.2 2114 0.3% -0.014 13 0.0% 0.0 11 0.0005 0.05 100 9000 3.6 0.2 2103 -0.2% -0.011 14 7.7% 0.38 12 0.0005 0.05 100 9000 3 0.16 2095 -0.6% 0.030 13 0.0% 0.0 13 0.0005 0.05 100 9000 3 0.24 2121 0.6% 0.032 13 0.0% 0.0 49 It is interesting that the channel length is the most sensitive parameter as this parameter is the only unbounded user-defined parameter. The other parameters are generally constrained and must be selected within a realistic range to simulate natural systems. For example, Manning?s n for natural streams generally ranges from 0.02 to 0.07, while the channel bottom width must be proportional to the depth of flow and watershed size. The channel length can essentially be as small or as large as the user desires based on the goals of the project. This underscores the importance of selecting the point of interest to compute the outflow hydrograph. If the simulated channel length is too short, there will be little to no difference between the Q po and the Q pin (i.e., little to no storage) and the outflow hydrograph will approximate the inflow hydrograph. However, a simulated channel reach that is too long will result in negative C 0 parameters and negative outflows, which is an undesirable and irrational phenomenon for the Muskingum method. Previous research has found that negative outflows can be avoided in the Muskingum model if x and K are selected such that (McCuen, 1998): nullnull 0.5?null null null1nullnull (4-1) and nullnull0.5 (4-2) 50 While channel geometry did affect the outflow hydrograph, the Muskingum method was not very sensitive to channel geometry parameters (i.e., Manning?s n, bottom width, channel bed slope, and side slopes for trapezoidal channels). These parameters are not directly input into the Muskingum model. Instead, they were used in Manning?s equation to calculate the average velocity through the reach, which in turn was used to calculate the routing coefficient K. Thus, the R s of each of these parameters can be assessed by examining Manning?s equation. For power models, the R s of each variable is equal to the exponent of that variable (McCuen and Knight, 2006). Based on this, it is reasonable that Manning?s n is the most important channel geometry parameter as the exponent for this parameter is -1, and the bed slope was comparatively half as sensitive because its exponent is 0.5. For embedded parameters such as channel width and side slopes, R s is not equal to the exponent of the hydraulic radius (0.67), and they need to be determined numerically (McCuen and Knight, 2006). In the Muskingum rectangular channel model, the R s for the channel width was approximately equal to the R s for the bed slope, while in the trapezoidal model the R s of bottom width plus the side slope were approximately equal to the R s for the bed slope. The sensitivities shown in Tables 4-2 and 4-3 indicate that the Muskingum model is not particularly sensitive to the x coefficient. Similar observations were made by Wang et al. (2006), who used a finite-difference Muskingum method and reported only a 2.8 percent change in the Q po when x was varied between 0 and 0.5. Wu et al. (1985) used looped plots to estimate the routing coefficients. While the overall error of the outflow hydrograph was reduced, the Q po did not significantly change (less than 2 percent 51 on average) as the routing parameters were modified. Therefore, the default x value of 0.2 likely doesn?t introduce significant error to the computed downstream hydrograph. The outflow time to peak had a maximum R s of 0.38, which might suggest that this output parameter is potentially more sensitive than Q po . However, upon further analysis, t po does not appear to be significantly affected by perturbations in the input parameters. The t po increased from 13 hours to 14 hours when the channel length and Manning?s n were increased by 20 percent in the rectangular channel model, and when the channel length, Manning?s n, bottom width, side slope were increased or when the bed slope was decreased in the trapezoidal channel model. However, this was a one- directional change, which means that t po did not equally decrease when the given parameter was adjusted by the same magnitude in the opposite direction. The t po probably did not decrease to 12 hours because of the coarse time step selected for routing. The time to peak is not a continuous variable as its value can only be reported as an integer-multiple of the time step (delta t). Delta t was selected as 1 hour for these simulations, so the minimum that t po could change is 1 hour. The true t po is likely slightly larger than 13 and it was rounded down to the nearest time step (13 hours). When input parameters were changed, t po increased enough to round up to 14 hours. However, when the same parameters were changed in the opposite direction the effect on t po was not significant enough to round down to 12 hours. It should also be noted that a change in 1 hour for the t po results in a 7 percent change from the baseline run and a R s of 0.38. Thus, a R s of 0.38 is the minimum non- zero value for t po , and therefore it is potentially not that different from zero. 52 Finally, as Table 4-1 shows, the sensitivity of Q po for the Muskingum and the VPMC method were generally comparable. As discussed below, the R s for t po was zero for all of the VPMC scenarios, similarly implying that t po is not a very sensitive output parameter in the Muskingum method. 4.2.2 3-Point Variable Parameter Muskingum-Cunge Model The results of the 3-Pt VPMC model sensitivity analyses are presented in Tables 4-4 and 4-5 for the rectangular and trapezoidal channels, respectively. For Q po , the results indicate that relative parameter importance for the 3-Pt VPMC rectangular channel model are (R s ): channel bed slope (0.029), Manning?s n (-0.027), channel length (-0.014), and bottom width (0.0055). The relative parameter importance for the Muskingum trapezoidal channel model are (R s ): channel bed slope (0.032), Manning?s n (-0.031), channel length (-0.019), side slope (-0.0035), and bottom width (0.0028). These results indicate that the 3-Pt VPMC method is not particularly sensitive to any of these input parameters and less sensitive than the Muskingum method. The 3-Pt VPMC method is less sensitive to changes in the input parameters than the Muskingum method for several reasons. First, many of the input parameters are correlated to other input parameters, and a change in one parameter is often offset by a change in a different correlated parameter. Second, the VPMC method uses several intermediate calculations prior to calculating the x and K coefficients. For example, the stage-discharge relationship is calculated using Manning?s equation, which is used to calculate the wave celerity, which is then used to calculate x and K. The effect of the change in a parameter used to calculate the stage-storage relationship is muted through this iterative process. 53 Table 4-4. Sensitivity Analysis Results for 3-Point Variable Parameter Muskingum- Cunge Method?Rectangular Channel (S = bed slope; W = bottom width; L = channel length) Run S (ft/ft) n W (ft) L (ft) Q po (cfs) ?Q po R s t po (hrs) ?t po R s 1 0.0005 0.05 100 9000 2216 NA NA 12 NA NA 2 0.0005 0.05 100 7200 2222 0.2% -0.012 12 0.0% 0.0 3 0.0005 0.05 100 10800 2209 -0.3% -0.016 12 0.0% 0.0 4 0.0005 0.05 80 9000 2213 -0.1% 0.0075 12 0.0% 0.0 5 0.0005 0.05 120 9000 2218 0.1% 0.0035 12 0.0% 0.0 6 0.0004 0.05 100 9000 2200 -0.7% 0.037 12 0.0% 0.0 7 0.0006 0.05 100 9000 2226 0.4% 0.021 12 0.0% 0.0 8 0.0005 0.04 100 9000 2227 0.5% -0.025 12 0.0% 0.0 9 0.0005 0.06 100 9000 2203 -0.6% -0.030 12 0.0% 0.0 Table 4-5. Sensitivity Analysis Results for 3-Point Variable Parameter Muskingum- Cunge Method?Trapezoidal Channel (S = bed slope; W = bottom width; L = channel length; z = side slope) Run S (ft/ft) n W (ft) L (ft) z Q po (cfs) ?Q po R s t po (hrs) ?t po R s 1 0.0005 0.05 100 9000 3 2214 NA NA 12 NA NA 2 0.0005 0.05 100 7200 3 2221 0.3% -0.016 12 0.0% 0.0 3 0.0005 0.05 100 10800 3 2204 -0.4% -0.022 12 0.0% 0.0 4 0.0005 0.05 80 9000 3 2212 -0.1% 0.0034 12 0.0% 0.0 5 0.0005 0.05 120 9000 3 2214 0.0% 0.0021 12 0.0% 0.0 6 0.0004 0.05 100 9000 3 2196 -0.8% 0.040 12 0.0% 0.0 7 0.0006 0.05 100 9000 3 2224 0.5% 0.024 12 0.0% 0.0 8 0.0005 0.04 100 9000 3 2226 0.6% -0.028 12 0.0% 0.0 9 0.0005 0.06 100 9000 3 2198 -0.7% -0.034 12 0.0% 0.0 10 0.0005 0.05 100 9000 2.4 2215 0.1% -0.0036 12 0.0% 0.0 11 0.0005 0.05 100 9000 3.6 2212 -0.1% -0.0035 12 0.0% 0.0 54 By comparison, in the Muskingum method the x coefficient is usually obtained by hydrograph analysis or it is assumed, and the K coefficient is calculated in a one-step process based on average reach characteristics. Therefore, these coefficients are impacted more directly by changes in the input parameters. Examples of the VPMC input parameter intercorrelation are provided below. The bed slope appears to be the most important parameter in the VPMC method, while it was only the third-most important parameter in the Muskingum method. This increased sensitivity may be attributed to the number of variables impacted by the bed slope. In the VPMC method, the bed slope directly and indirectly effects the calculation of the K and x routing coefficients. The bed slope is used directly to calculate x (refer to Eq. 2-14). In addition, the bed slope is also used in Manning?s equation to develop the stage-discharge relationship. This relationship is used to calculate the wave celerity, which is used in both the K and x routing coefficients equations. Comparatively, in the Muskingum method, the bed slope is used only for estimating the average velocity through the reach for calculating K. Manning?s n is the second-most important parameter in both the 3-Pt VPMC and Muskingum method. As discussed in Section 4.2.1 above, Manning?s equation is most sensitive to Manning?s n. Therefore, while Manning?s n isn?t used in Eqs. 2-14 or 2-15 to directly calculate K and x, it does impact the unit discharge and wave celerity that are used to calculate the routing coefficients. While channel length was the most important input parameter in the Muskingum method, it was only the third-most important parameter in the 3-Pt VPMC method. This result is not expected as channel length is directly used to calculate both the K and x 55 routing coefficients and it is the most important parameter in the Muskingum method. This may be attributed to how channel length is used to calculate K and x. Channel length is positively correlated and proportional to both K and x; however, K is negatively correlated to Q po and x is positively correctly to Q po . Therefore, the changes in the channel length are almost offset by the opposite effect that K and x have on Q po . Similar to the sensitivity results for the Muskingum method, the 3-Pt VPMC method is not sensitive to changes in either the side slope or the bottom width. However, it should be mentioned that the 3-Pt VPMC model had a positive correlation between the channel bottom width and Q po , while a negative correlation was observed in the Muskingum model. One would expect the bottom width and Q po to be negatively correlated. As the channel becomes narrower, the depth of flow increases and the wetted perimeter decreases. Thus, turbulence is reduced and the velocity increases, resulting in an increased Q po . The positive correlation between bottom width and Q po in the 3-Pt VPMC model may be attributed to intercorrelations with other variables. The bottom width probably affects more routing variables than any other input parameter, which makes it more noteworthy that the VPMC model is not sensitive to this parameter. The relationships between bottom width and other VPMC parameters are shown in Table 4-6. Table 4-6. Correlation Between Bottom Width and Other VPMC Parameters (c = celerity; q o = unit-width discharge) c q o x K C 0 C 1 C 2 Neg. Neg. Pos. Pos. Neg. Pos. Neg. 56 The negative correlation between bottom width and c and and q o was expected. Based on Eq. 2-15, K decreases as c increases and K is unaffected by q o . Based on Eq. 2- 14, x is affected differently by c and q o ?as c increases, x increases, and as q o increases, x decreases. Thus, Table 4-6 shows the expected relationship between bottom width and K?i.e., as the bottom width increases, K increases. However, when both c and q o decrease, the response to x is difficult to determine because x is proportional to the ratio of q o /c. As the bottom width increases, the ratio of q o /c decreases, meaning that q o is not as large compared to c, and therefore x increases when bottom width increases. (Refer to Section 4.4.1 for more of a discussion on the ratio of q o /c.) There is still one more aspect that needs to be discussed to understand why bottom width has a positive correlation with Q po . Ordinarily, Q po is positively correlated to x and negatively correlated to K. However, it was determined that bottom width is positively correlated to both x and K, thus x appears to be the driving force behind the positive correlation between bottom width and Q po . To understand why Q po increases as bottom width increases, one must look at how the relative change in x and K changes with bottom width. Table 4-7 shows the x and K values for different bottom widths at the peak discharge in the outflow hydrograph. The x coefficient is much more sensitive to bottom width (R s is approximately 4) compared to K (R s is approximately 0.2). Thus, as bottom width changes, the larger change in x compared to K results in the positive correlation with Q po . The t po did not change under any of the scenarios evaluated. This is related to the insensitivity of the VPMC model and the coarse time step as discussed in the Muskingum method sensitivity analysis results. 57 Table 4-7. Routing Coefficients at the Time to Peak for Different Bottom Widths Bottom Width (ft) x K 80 -0.24 2187 100 -0.11 2273 120 -0.033 2371 4.2.3 4-Point Variable Parameter Muskingum-Cunge Model The results of the sensitivity analyses for the 4-Pt VPMC model are presented in Tables 4-8 and 4-9 for the rectangular and trapezoidal channels, respectively. These results show that the 4-Pt VPMC method has the same sensitivity to changes in the input parameters as those of the 3-pt method. This is because the 3-Pt and 4-Pt VPMC had the same Q po and t po for each of the scenarios evaluated. It is noted that the entire outflow hydrographs were also very similar, with only minor differences in outflow discharge at a couple ordinates. The method is very insensitive to the computation of the fourth point. The similarities between the 3-Pt and 4-Pt VPMC models are well documented. Ponce et al. (1978), Ponce and Changanti (1994), McCuen (1998), and Tang et al. (1999) have compared the outflow results of these two models under various conditions. The Q po was approximately equal in all conditions (less than 0.5 percent difference), and the t po was equal for all but one of the scenarios evaluated. In one of the Ponce and Changanti (1994) scenarios, the t po for the 4-Pt VPMC method was less than the t po from 58 Table 4-8. Sensitivity Analysis Results for 4-Point Variable Parameter Muskingum-Cunge Method?Rectangular Channel (S = bed slope; W = bottom width; L = channel length) Run S (ft/ft) n W (ft) L (ft) Q po (cfs) ?Q po R s t po (hrs) ?t po R s 1 0.0005 0.05 100 9000 2216 NA NA 12 NA NA 2 0.0005 0.05 100 7200 2222 0.2% -0.012 12 0.0% 0.0 3 0.0005 0.05 100 10800 2209 -0.3% -0.016 12 0.0% 0.0 4 0.0005 0.05 80 9000 2213 -0.1% 0.0075 12 0.0% 0.0 5 0.0005 0.05 120 9000 2218 0.1% 0.0035 12 0.0% 0.0 6 0.0004 0.05 100 9000 2200 -0.7% 0.037 12 0.0% 0.0 7 0.0006 0.05 100 9000 2226 0.4% 0.021 12 0.0% 0.0 8 0.0005 0.04 100 9000 2227 0.5% -0.024 12 0.0% 0.0 9 0.0005 0.06 100 9000 2203 -0.6% -0.030 12 0.0% 0.0 Table 4-9. Sensitivity Analysis Results for 4-Point Variable Parameter Muskingum-Cunge Method?Trapezoidal Channel (S = bed slope; W = bottom width; L = channel length; z = side slope) Run S (ft/ft) n W (ft) L (ft) z Q po (cfs) ?Q po R s t po (hrs) ?t po R s 1 0.0005 0.05 100 9000 3 2214 NA NA 12 NA NA 2 0.0005 0.05 100 7200 3 2221 0.3% -0.016 12 0.0% 0.0 3 0.0005 0.05 100 10800 3 2204 -0.4% -0.022 12 0.0% 0.0 4 0.0005 0.05 80 9000 3 2212 -0.1% 0.0034 12 0.0% 0.0 5 0.0005 0.05 120 9000 3 2214 0.0% 0.0022 12 0.0% 0.0 6 0.0004 0.05 100 9000 3 2196 -0.8% 0.040 12 0.0% 0.0 7 0.0006 0.05 100 9000 3 2224 0.5% 0.024 12 0.0% 0.0 8 0.0005 0.04 100 9000 3 2226 0.6% -0.028 12 0.0% 0.0 9 0.0005 0.06 100 9000 3 2198 -0.7% -0.034 12 0.0% 0.0 10 0.0005 0.05 100 9000 2.4 2215 0.1% -0.0036 12 0.0% 0.0 11 0.0005 0.05 100 9000 3.6 2212 -0.1% -0.0035 12 0.0% 0.0 59 the 3-Pt VPMC method; however, an explanation was not offered for this anomaly. Many of these authors noted that the 4-Pt VPMC is slighter more volume conservative compared to the 3-Pt VPMC; however, the improvement was usually less than one percent. This slight increase in accuracy generally does not justify the additional effort associated with using the 4-Pt VPMC. 4.3 RESULTS OF UNCERTAINTY ANALYSES As discussed in Section 3.6, Monte Carlo simulations with 10,000 runs per simulations were used to evaluate the uncertainty of the channel routing models. Two Monte Carlo simulations were conducted for each model. In the first scenario, only the input parameters used in the channel routing model were allowed to vary. In the second scenario, all of the user-defined input parameters were allowed to vary. The first scenario emphasizes the uncertainty of the routing parameters. The second scenario was intended to estimate the total variability introduced to the downstream hydrograph from all sources of uncertainty using commonly applied techniques. The results of the Monte Carlo simulation were graphically summarized using histograms. A histogram is a graphical display of tabulated data. Bars are used to group results into ranges or intervals, and the length of the bar represents the frequency of individual results that fell into the specified range or interval. The main advantage of histograms is they allow for the analysis of very large data sets by reducing them into a single graph. 60 4.3.1 Muskingum Model 4.3.1.1 Uncertainty Results for the Muskingum Rectangular Channel The results of the uncertainty analyses for the Muskingum rectangular channel are presented in Figures 4-1a and 4-1b and Table 4-10. Figures 4-1a and 4-1b show the histograms for the Q po when only the Muskingum input parameters were varied in the Monte Carlo simulations (Scenario 1) and when all in the user-defined input parameters were varied (Scenario 2), respectively. Table 4-10 presents the statistical results for Q pin , Q po , t pin , and t po . Figure 4-1a shows that, when only the Muskingum input parameters were varied, the distribution of the Q po in the downstream hydrograph was relatively confined?90 percent of the simulated outflows were between 2048 cfs and 2199 cfs. The mean Q po was 2132 cfs, which is equal to the original simulated outflow (2131 cfs). The C v was 0.022, which indicates that there is very low variability in the downstream hydrograph when considering uncertainty in the channel routing input parameters. The maximum t po was 15 hours; however, the mean t po (13.3) and the low C v (0.034) indicate very little variation in t po . Table 4-10 presents the statistical results for the rectangular channel Q pin , t pin , Q po , and t po . When only the Muskingum model parameters were varied, the Q pin and t pin did not vary as expected as the parameters used to calculate the upstream hydrograph were not changed. When all of the user-defined input parameters were varied (Scenario 2), the mean peak discharge at the downstream section (2045 cfs) decreased compared to the original results (2131 cfs). The lower peak discharge may be explained by looking at the shape of 61 Figure 4-1a. Muskingum Rectangular Channel Model Uncertainty Analysis Results? Muskingum Model Input Parameters Figure 4-1b. Muskingum Rectangular Channel Model Uncertainty Analysis?All Input Parameters 62 Table 4-10. Uncertainty Analysis Results for the Muskingum Rectangular Channel Model Muskingum Rectangular Channel?Model Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 2240 2240 2240 0 0 Time to Peak Inflow (t pin ) 11 11 11.0 0 0 Peak Outflow (Q po ) 1895 2237 2132 46 0.022 Time to Peak Outflow (t po ) 12 15 13.3 0.5 0.034 Muskingum Rectangular Channel?All User-Defined Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 100 7891 2148 978 0.46 Time to Peak Inflow (t pin ) 6 80 13.4 7.0 0.52 Peak Outflow (Q po ) 100 7071 2039 922 0.45 Time to Peak Outflow (t po ) 7 80 15.6 7.2 0.46 the histograms for two different scenarios (Figures 4-1a and 4-1b). Both distributions appear to approximate a log-normal population. However, when only the Muskingum input parameters were varied (Figure 4-1a), the distribution of Q po was negatively skewed because Q po was bound by the peak inflow (2240 cfs). Conversely, when all of the user-defined input parameters were varied (Figure 4- 1b) the distribution of Q po was positively skewed because 0 cfs was the lower bound for Q po . Therefore, there was a higher frequency Q po s less than 2100 cfs compared to when only the model input parameters were varied. The standard deviation for Q po increased from 46 cfs (Scenario 1) to 922 cfs (Scenario 2), resulting in a C v of 0.45. This increase is expected as it is a response to the variability in the upstream hydrograph. When all of the parameters were varied, the 63 standard deviation for Q pin was 978 cfs due to the uncertainty in the parameters used to develop the upstream hydrograph. Thus, variability in the downstream hydrograph increases with increased uncertainty in the upstream hydrograph. The t po ranged from a minimum value of 7 hours to a maximum value of 80 hours. The mean was 15.6 hours, which is about 20 percent greater than the baseline simulation (13 hours). Again, the t po is directly related to the t p in the inflow hydrograph. When all of the parameters were varied in the Monte Carlo simulations, t pin would be expected to have more variability. Parameters that affect t pin ,and therefore t po , are: drainage area, watershed slope, velocity K factor, Manning?s n, channel slope, and channel length. However, the results of the sensitivity analyses indicate that t po is insensitive to Manning?s n, channel slope, and channel length. Further, the C v for watershed size was very low (0.025), indicating that there is little uncertainty associated with this parameter. Thus, the increased variability for t po is a result of the increased variability in t pin due to uncertainty in the watershed slope and/or velocity K. 4.3.1.2 Uncertainty Results for the Muskingum Trapezoidal Channel The results of the uncertainty analyses for the Muskingum trapezoidal channel are presented in Table 4-11 and Figures 4-2a and 4-2b. Figures 4-2a and 4-2b show the histograms for the Q po when only the Muskingum input parameters were varied and when all of the user-defined input parameters were varied, respectively. Table 4-11 presents the statistical results of the Monte Carlo simulations. When only the Muskingum input parameters were varied, there was little variability of Q po and t po based on the C v of 0.025 and 0.039, respectively (Table 4-11). 64 Figure 4-2a. Muskingum Trapezoidal Channel Model Uncertainty Analysis Results? Muskingum Model Input Parameters Figure 4-2b. Muskingum Trapezoidal Channel Model Uncertainty Analysis Results?All Input Parameters 65 Table 4-11. Uncertainty Analysis Results for the Muskingum Trapezoidal Channel Model Muskingum Trapezoidal Channel?Model Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 2240 2240 2240 0 0 Time to Peak Inflow (t pin ) 11 11 11 0 0 Peak Outflow (Q po ) 1854 2240 2114 54 0.025 Time to Peak Outflow (t po ) 12 15 13.5 0.52 0.039 Muskingum Trapezoidal Channel?All User-Defined Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 104 7610 2150 979 0.46 Time to Peak Inflow (t pin ) 6 80 13.4 7.1 0.53 Peak Outflow (Q po ) 103 6869 2023 912 0.45 Time to Peak Outflow (t po ) 7 80 15.7 7.2 0.46 These results are very similar to the uncertainty results for the Muskingum rectangular channel. The sensitivity analyses indicate that the Muskingum method is insensitive to changes in channel geometry; therefore, the variability in the trapezoidal channel model would be expected to be similar to the variability in the rectangular model. When all of the input parameters were varied, the C v for Q po and t po increased to 0.45 and 0.46, respectively. As discussed in Section 4.3.1.1, the increased variability for t po is a result of the increased variability in t pin due to uncertainty in the watershed slope and/or velocity K. It is noted that the maximum Q po (2240 cfs) is identical to the Q pin (2240 cfs). The input values used in the Monte Carlo simulation were verified, and the similar peak discharges are a result of negative outflows during the first few time steps. The 66 occurrence of negative outflows is mainly restricted to the VPMC method and is discussed in Section 2.6. Negative outflows increase the Q po and under extreme conditions it is possible for Q po to be greater than Q pin in the Muskingum model. This is a result of numeric dispersion in the channel routing model and in reality Q po is limited to Q pin if lateral inflows are ignored. In this scenario, the negative outflow resulted mainly from a shortened channel length (7488 ft) and an increased x coefficient of (0.48). While an x value of 0.48 is within the 0 to 0.5 range for the Muskingum method, it is greater than two standard deviations from the simulated mean (mean = 0.2, standard deviation = 0.1), and few if any experts recommend such high values for natural channels. Thus, it is very unlikely that the x coefficient would be 0.48 for any natural channel; however, this result indicates that under extreme conditions the Muskingum method can produce negative coefficients and increased Q po s. If this were to occur, it is recommended that the channel length be increased such that the conditions of Eq. 4-1 are satisfied. 4.3.2 3-Point Variable-Parameter Muskingum-Cunge Model 4.3.2.1 Uncertainty Results for the 3-Pt VPMC Rectangular Channel The results of the uncertainty analyses for the 3-Pt VPMC rectangular channel are presented in Table 4-12 and Figures 4-3a and 4-3b. Figures 4-3a and 4-3b show the histograms for Q po when only the 3-Pt VPMC input parameters were varied and when all of the user-defined input parameters were varied, respectively. The VPMC results are consistent with the Muskingum rectangular channel uncertainty results, i.e., there was significantly more variability in the upstream and downstream results when all of the user-defined parameters were varied in the Monte 67 Figure 4-3a. 3-Pt Variable Parameter Muskingum-Cunge Rectangular Channel Model Uncertainty Analysis Results?3-Pt VPMC Model Input Parameters Figure 4-3b. 3-Pt Variable Parameter Muskingum-Cunge Rectangular Channel Model Uncertainty Analysis Results?All Input Parameters 68 Table 4-12. Uncertainty Analysis Results for the 3-Pt VPMC Rectangular Channel Model 3-Pt VPMC Rectangular Channel?Model Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 2240 2240 2240 0 0 Time to Peak Inflow (t pin ) 11 11 11 0 0 Peak Outflow (Q po ) 2032 2244 2211 22 0.010 Time to Peak Outflow (t po ) 12 13 12.0 0.05 0.004 3-Pt VPMC Rectangular Channel?All User-Defined Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 102 7497 2149 987 0.46 Time to Peak Inflow (t pin ) 6 80 13.6 7.9 0.58 Peak Outflow (Q po ) 102 7254 2111 953 0.45 Time to Peak Outflow (t po ) 6 80 14.4 7.9 0.55 Carlo simulations. However, a few distinctions should be made between the 3-Pt VPMC rectangular channel model the Muskingum model counterpart. The variability of Q po and t po for the 3-Pt VPMC method was actually lower than the variability of Q po and t po in the Muskingum method (refer to Table 4-10). This result is expected based on the sensitivity analysis results that revealed that the VPMC method is not as sensitive to input parameter uncertainty as the Muskingum method. When only the model input parameters were varied, the maximum Q po was approximately equal to the Q pin . The Monte Carlo input values for the VPMC simulation were verified, and the slightly increased outflow was not only a result of negative initial outflows but also the effect of the combination of the input parameters on the x and K routing coefficients. In this case, negative outflows did occur during the first few time 69 steps; however, the negative outflows were not significantly different than the negative outflows observed during most of the VPMC simulations including the baseline run. In addition to negative outflows, the increased Q po was likely the result of using a higher bed slope (0.0007 ft/ft), a lower Manning?s n (0.02), and a longer channel length (13,491 ft). This combination of input parameters resulted in a larger x coefficient and a smaller K coefficient, which produced a slightly steeper rising limb of the hydrograph and a greater Q po . 4.3.2.2 Uncertainty Results for the 3-Pt VPMC Trapezoidal Channel The results of the uncertainty analyses for the 3-Pt VPMC trapezoidal channel are presented in Table 4-13 and Figures 4-4a and 4-4b. Figures 4-4a and 4-4b show the histograms for the Q po when only the 3-Pt VPMC input parameters were varied and when all of the user-defined input parameters were varied, respectively. The variability of Q po and t po for the 3-Pt VPMC trapezoidal model are similar to the uncertainty results for the 3-Pt VPMC rectangular channel. Again, these results are expected based on the sensitivity analysis results, which indicate these methods are insensitive to changes in channel geometry. The 3-Pt VPMC trapezoidal channel model had lower variability for Q po compared to the Muskingum trapezoidal channel model (C v = 0.010 compared to 0.025) 70 Figure 4-4a. 3-Pt Variable Parameter Muskingum-Cunge Trapezoidal Channel Model Uncertainty Analysis Results?3-Pt VPMC Model Input Parameters Figure 4-4b. 3-Pt Variable Parameter Muskingum-Cunge Trapezoidal Channel Model Uncertainty Analysis Results?All Input Parameters 71 Table 4-13. Uncertainty Analysis Results for the 3-Pt VPMC Trapezoidal Channel Model 3-Pt VPMC Trapezoidal Channel?Model Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 2240 2240 2240 0 0 Time to Peak Inflow (t pin ) 11 11 11 0 0 Peak Outflow (Q po ) 2061 2244 2208 23.1 0.010 Time to Peak Outflow (t po ) 12 13 12.0 0.18 0.015 3-Pt VPMC Trapezoidal Channel?All User-Defined Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 105 9159 2140 982 0.46 Time to Peak Inflow (t pin ) 6 80 13.5 7.0 0.52 Peak Outflow (Q po ) 105 8775 2101 950 0.45 Time to Peak Outflow (t po ) 6 80 14.4 7.0 0.49 and t po (C v = 0.015 compared to 0.039). Lower variability was also observed in the 3-Pt VPMC rectangular model compared to the Muskingum rectangular model, which was attributed to the VPMC method being less sensitive to changes in input parameters (refer to Section 4.3.2.1). Similarly, the VPMC trapezoidal model is less sensitive to changes in the input parameters than the Muskingum rectangular model based on the results of the sensitivity analyses. 72 4.3.3 4-Point Variable-Parameter Muskingum-Cunge Model 4.3.3.1 Uncertainty Results for the 4-Pt VPMC Rectangular Channel The results of the uncertainty analyses for the 4-Pt VPMC rectangular channel are presented in Table 4-14 and Figures 4-5a and 4-5b. Figures 4-5a and 4-5b show the histograms for the Q po when only the 4-Pt VPMC input parameters were varied and when all of the user-defined input parameters were varied, respectively. The results of the uncertainty analyses for the 4-Pt VPMC rectangular channel are very similar to the uncertainty results for the 3-Pt VPMC rectangular channel. The C v for the 3-Pt and 4-Pt models were exactly the same when only the VPMC model input parameters were varied. The C v for Q po (0.45) was also the same as the C v for the 3-Pt VPMC rectangular channel model when all of the input parameters were varied. However, the C v for t po (0.48) was slightly lower than the C v for t po for the 3-Pt rectangular model (0.55) when all of the input parameters were varied. This is due to the difference in variability in the t pin between the 3-Pt and 4-Pt models. The C v for t pin in the 4-Pt model was 0.51, while the C v for t pin in the 3-Pt was 0.58. The higher variability in the 3-Pt model is likely the result of the random values selected during the Monte Carlo simulations and the sensitivity of t pin to watershed slope and velocity K (refer to Section 4.3.1.1). 73 Figure 4-5a. 4-Pt Variable Parameter Muskingum-Cunge Rectangular Channel Model Uncertainty Analysis Results?4-Pt VPMC Model Input Parameters Figure 4-5b. 4-Pt Variable Parameter Muskingum-Cunge Rectangular Channel Model Uncertainty Analysis Results?All Input Parameters 74 Table 4-14. Uncertainty Analysis Results for the 4-Pt VPMC Rectangular Channel Model 4-Pt VPMC Rectangular Channel?Model Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 2240 2240 2240 0 0 Time to Peak Inflow (t pin ) 11 11 11 0 0 Peak Outflow (Q po ) 2015 2243 2211 21.4 0.010 Time to Peak Outflow (t po ) 12 13 12.0 0.04 0.004 4-Pt VPMC Rectangular Channel?All User-Defined Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 125 7180 2174 983 0.45 Time to Peak Inflow (t pin ) 6 80 13.3 6.8 0.51 Peak Outflow (Q po ) 124 6951 2136 951 0.45 Time to Peak Outflow (t po ) 6 80 14.1 6.8 0.48 4.3.3.2 Uncertainty Results for the 4-Pt VPMC Trapezoidal Channel The results of the uncertainty analyses for the 4-Pt VPMC trapezoidal channel are presented in Table 4-15 and Figures 4-6a and 4-6b. Figures 4-6a and 4-6b show the histograms for the Q po when only the 4-Pt VPMC input parameters were varied and when all of the user-defined input parameters were varied, respectively. The results of the uncertainty analyses for the 4-Pt VPMC trapezoidal channel are very similar to the uncertainty results for the 3-Pt VPMC trapezoidal channel. These results are expected based on the results of the sensitivity analyses, and the 3-Pt and 4-Pt VPMC rectangular channel uncertainty analyses. This means that uncertainty in the input 75 Figure 4-6a. 4-Pt Variable Parameter Muskingum-Cunge Trapezoidal Channel Model Uncertainty Analysis Results?4-Pt VPMC Model Input Parameters Figure 4-6b. 4-Pt Variable Parameter Muskingum-Cunge Trapezoidal Channel Model Uncertainty Analysis Results?All Input Parameters 76 Table 4-15. Uncertainty Analysis Results for the 4-Pt VPMC Trapezoidal Channel Model 4-Pt VPMC Trapezoidal Channel?Model Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 2240 2240 2240 0 0 Time to Peak Inflow (t pin ) 11 11 11 0 0 Peak Outflow (Q po ) 2047 2244 2207 24 0.011 Time to Peak Outflow (t po ) 12 13 12.0 0.2 0.016 4-Pt VPMC Trapezoidal Channel?All User-Defined Input Parameters Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 105 7245 2134 976 0.46 Time to Peak Inflow (t pin ) 6 80 13.5 7.1 0.53 Peak Outflow (Q po ) 105 7073 2095 945 0.45 Time to Peak Outflow (t po ) 7 80 14.4 7.1 0.50 parameters similarly affect the computed downstream hydrograph using the 3-Pt and 4-Pt methods. Therefore, users should not preferentially select one of these models over the other assuming that it will result in a more accurate downstream hydrograph. 4.4 DISCUSSION AND OBSERVATIONS A general discussion of the uncertainty and sensitivity analyses results is provided in the following sections. The Muskingum and VPMC models are compared, and methods to reduce uncertainty are presented. Finally, the negative outflows observed in the VPMC method are explored and methods to reduce them are proposed. 77 4.4.1 Comparison of Muskingum and VPMC Models The results of the uncertainty and sensitivity analyses indicate that the Muskingum model is slightly more sensitive to model input parameters than the VPMC model. While neither channel routing method was particularly sensitive to any of the input parameters, channel length, bed slope, and Manning?s n were the three most sensitive parameters in both models. Channel length was the most sensitive parameter in the Muskingum method, while bed slope was the most sensitive parameter in the VPMC method. The increased sensitivity to bed slope may be attributed to the number of variables in the VPMC method impacted by the bed slope compared to the Muskingum method. Finally, both methods were insensitive to channel geometry parameters such as bottom width and side slopes. The results of the sensitivity and uncertainty analyses for rectangular and trapezoidal channels were very similar, with both scenarios resulting in similar R s values for the three most sensitive input parameters. The results of the uncertainty analysis showed that the C v and distribution of Q po were similar for the rectangular and trapezoidal channels. These led to the conclusion that both methods are relatively insensitive to channel geometry. However, it is important to ensure that the stage- discharge relationship is accurate for the geometry selected as the amount of storage in the channel will be affected by the channel geometry. The uncertainty in the inflow hydrograph appears to be the major factor contributing to error in the downstream hydrograph. When only the channel routing input parameters were varied in the Monte Carlo simulations, there was only a moderate amount of variation in the downstream hydrograph. However, when all sources of 78 uncertainty are considered, the variation in the downstream hydrograph significantly increased. Thus, reducing the uncertainty in the parameters used to develop the inflow hydrograph could significantly reduce the uncertainty in the computed downstream hydrograph. The inflow hydrograph was most sensitive to the rainfall intensity, watershed slope, velocity K, and PRF. A modest reduction in the uncertainty in these parameters, as well as the three most sensitive channel routing parameters, could lead to significant improvement. The C v s for the parameters listed above were reduced by 50 percent and Monte Carlo simulations were re-run for the Muskingum rectangular channel model. The results are presented in Figure 4-7 and Table 4-16. The standard deviation for Q po and t po were reduced by almost 50 percent (compared to Table 4-10), and the range for 90 percent of the peak discharge rates (1231 cfs to 3070 cfs) was reduced by about 40 percent (compared to Figure 4-1b). While improving the accuracy of these parameters does decrease the variation in the downstream hydrograph, uncertainty will always be present and needs to be considered when evaluating results. Table 4-16. Uncertainty Analysis Results for the Muskingum Rectangular Channel Model, Reduced C v Muskingum Rectangular Channel?All User-Defined Input Parameters, Reduced Cv Parameter Min. Max. Mean Stand. Dev. C v Peak Inflow (Q pin ) 318 5452 2170 604 0.28 Time to Peak Inflow (t pin ) 7 75 12.1 3.1 0.26 Peak Outflow (Q po ) 302 5151 2067 568 0.27 Time to Peak Outflow (t po ) 8 79 14.3 3.3 0.23 79 Figure 4-7. Muskingum Rectangular Channel Model Uncertainty Analysis?All Input Parameters, Reduced C v 4.4.2 VPMC Method and Negative Outflows The occurrence of negative outflows in the VPMC model is a common phenomenon that has been widely documented and discussed (refer to Section 2.6). While it is acknowledged that this is undesirable, many experts dismiss them by stating that they have an insignificant effect on the computed outflow hydrograph. The research presented herein suggests that the negative outflows could have a significant effect on the computed outflow hydrograph, especially when all sources of uncertainty are considered. When the same upstream hydrograph and channel properties were routed using the Muskingum and VPMC methods, one would expect similar peak outflows. However, the peak outflow estimated using the VPMC method was, on average, almost 5 percent 80 higher than the peak outflow from the Muskingum method, and the t po was on average 1 hour sooner. This is likely a consequence of the negative outflows, which result in a steeper rising limb and overestimation of Q po and shorter t po . More specifically, the negative outflows are the result of the routing coefficients violating the constraints, which is common with actual data. The errant coefficients force the initial ordinates downward, and negative, and then force the following coefficients upward because of the other values (C 1 and C 2 ). A more in-depth analysis of the VPMC method reveals that several factors may contribute to the undesired negative outflows. First, negative outflows were always accompanied by negative C 0 parameters. Negative C 0 parameters occurred during the first few time steps on the rising limb of the hydrograph. As previously discussed, the first few time steps on the rising limb of the hydrograph is the region where the negative outflows occur. Modeling results were inspected for the presence of negative C 0 parameters and negative outflows, and it was observed that negative outflows occurred for approximately half the number of time steps that negative C 0 parameters occurred. Also, the more negative the C 0 ordinate, the more negative the outflow. Thus, minimizing the number and magnitude of negative C 0 ordinates will reduce negative outflows. The negative outflows also relate to the variability in the wave celerity (c) and the unit-width discharge (q o ). The unit-width discharge varied by greater than an order of magnitude throughout the duration of the simulation, while the wave celerity generally varied by a factor of 5 or less. As q o and c increase during the rising limb of the hydrograph, x and K decrease and C 0 increases. The rate of decrease in K appears to have 81 a bigger impact on negative C 0 ; however, the variability in x also affects the C 0 . The equation for x (Eq. 2-14) shows that x is related to q o , c, S o , and L. Since bed slope and channel length are fixed for all time steps, the magnitude of x is proportional to the ratio of q o /c. As the upstream discharge increases, the stage increases and c and q o increase. However, q o increases at a faster rate than c, and as a result, the ratio of q o /c increases. During the first few ordinates when negative outflows occur, the ratio of q o /c was less than or approximately equal to one. As this ratio increased and became larger than one, negative outflows were reduced. An attempt was made to control the ratio of q o /c by using a reference q o and/or reference c based on the inflow characteristics. When a reference value was used for either parameter and for both parameters, the negative outflows are eliminated; however, some or all of the wave diffusion is eliminated, defeating the purpose of the VPMC method. These analyses show that while negative outflows can be minimized, they are very difficult to eliminate and still simulate wave diffusion. 82 CHAPTER V CONCLUSIONS AND RECOMMENDATIONS 5.1 INTRODUCTION The goal of this research was to perform uncertainty and sensitivity analyses of the Muskingum and Variable Parameter Muskingum-Cunge channel routing methods. While several studies have developed approaches to estimate the routing parameters to characterize in-channel processes, to date significant research to help understand how uncertainty associated with each input parameter affects the computed downstream hydrograph has not been conducted. Similarly, significant research is also lacking on the sensitivity of the computed downstream hydrograph to the input parameters. To accomplish the research goal, an upstream hydrograph was developed using commonly applied engineering methods, and the hydrograph was routed through the channel routing models to obtain a downstream hydrograph. Two channel cross-section geometries were assumed, rectangular and trapezoidal, to evaluate the effect of channel geometry on model uncertainty. Finally, the channel routing input parameters were selected based on previous research and published guidance. The sensitivity analysis for each channel routing model was conducted by changing one input parameter and observing the change in the downstream peak discharge and time to peak. Monte Carlo analysis was used to evaluate how input parameter uncertainty affects the downstream hydrograph. All parameters were varied together for the uncertainty analyses. 83 5.2 SUMMARY OF RESULTS The results of the sensitivity analysis indicate that the Muskingum method is more sensitive to changes in the input parameters than the 3-Pt and 4-Pt VPMC methods. This is likely due to the large number of computations necessary to calculate the VPMC method x and K coefficients, which reduces the sensitivity of x and K to small changes in one parameter, and intercorrelation between the input parameters. In the Muskingum method, x and K are either assumed or calibrated using existing hydrographs; therefore, when a parameter is changed, it more directly impacts x and K. Table 5-1 summarizes the results of the sensitivity analysis. Table 5-1. Sensitivity Rankings of the Inputs to the Muskingum and VMPC Methods Order of Importance Muskingum Method 3-Pt VPMC Method 4-Pt VPMC Method 1 L S S 2 n n n 3 S L L 4 W z z 5 x W W 6 z -- -- The results of the uncertainty analysis were consistent with the sensitivity results. When only the channel routing input parameters were varied in the Monte Carlo simulations, the Muskingum method had higher variability in the downstream hydrograph peak discharge and time to peak compared to the VPMC method. This is expected as the Muskingum method is more sensitive to changes in input parameters. 84 The computations of the VPMC method tend to reduce the flexibility inherent to the Muskingum method. However, when all of the user-defined input parameters were varied, the variability in the downstream hydrographs were very similar for the Muskingum and VPMC methods. The results indicate that error in the input parameters used to develop the upstream hydrograph is the greatest source of error in the downstream hydrograph. Reducing error in these parameters can reduce variability in the downstream hydrograph. 5.3 CONCLUSIONS Several conclusions can be made related to the sensitivity and uncertainty associated with the Muskingum and VPMC channel routing models. Perhaps the most important conclusion of this research is that all sources of error must be considered when using these models, even parameters not directly used in the models. The Muskingum model generally reacted to changes in input parameters in a predictable and rational way. However, it was almost impossible to predict how the VPMC model would respond to changes in certain input parameters. This is illustrated by the positive correlation observed between channel bottom width and Q po (refer to Section 4.2.2), which is attributed to the bottom width being correlated to many other variables that effect Q po . Bottom width would be expected to have a negative correlation with Q po , because as bottom width increases, the wetted perimeter increases and therefore the velocity decreases and Q po should decrease. This was observed in the Muskingum method. 85 Further, the VPMC and the Muskingum methods are based on the storage equation. While the VPMC method introduces wave diffusion to better simulate natural processes, one would expect similar output from both of these models. However, the Muskingum peak outflow was almost 5 percent lower and the t po was 1 hour shorter than the VPMC method. This may not initially seem to be that significant; however, this really is a big difference in Q po considering that in all of the scenarios evaluated in the sensitivity analysis the maximum change in the peak outflow was only 2.2 percent. It is noted that the same watershed characteristics, inflow hydrograph, and channel dimensions were used for both models. Thus, the difference is attributed to how the different models compute the x and K parameters, and the negative outflows introduced by numeric diffusion in the VPMC method. Negative outflows were observed in almost all of the VPMC scenarios evaluated. This is an undesirable and an irrational phenomenon that appears to be inherent in the VPMC method due to numeric dispersion associated with simulated wave diffusion. Many researchers have suggested that negative outflows have an insignificant effect on the outflow hydrograph and they do not need to be controlled. However, the research presented herein suggests that the negative outflows could have a significant effect on the computed outflow hydrograph, especially when all sources of uncertainty are considered. 5.4 RECOMMENDATIONS AND FUTURE RESEARCH Based on the results of this research, the following recommendations should be considered when using these models: 86 ? This research highlights the importance of conducting sensitivity and uncertainty analyses on all channel routing models. When developing an environmental model to simulate natural processes, users should conduct sensitivity and uncertainty analyses so that the reliability of the computed output is better understood and in order to make more effective decisions. ? These results indicate that variability in the downstream hydrograph is significantly affected by uncertainty in the input parameters used to develop the upstream hydrograph. There was a high degree of uncertainty for several of the input parameters, and the parameters that affected the upstream hydrograph the most were: rainfall intensity, watershed slope, velocity K factor, and the peak rate factor. When using the SCS unit hydrograph approach to generate an upstream hydrograph, a reasonable effort should be made to reduce the user error in these parameters. ? Both of these models should not be treated as a ?black box?. All of the time steps (especially those on the rising limb of the hydrograph) should be checked to ensure that the input parameters are resulting in rational x and K coefficients, and C 0 parameters were not negative. The x coefficient and C 0 parameter should not be negative in the Muskingum model. While negative coefficients are hydraulically rational in the VPMC model, effort should be made to reduce their magnitude and the number of time steps in which they occur in order to minimize the negative outflows. ? When selecting a channel routing model, consideration should also be given to the constant parameter M-C method. While not evaluated in this research, the 87 constant parameter method has the advantage of calculating the routing parameters based on the physical properties of the channel without the disadvantage of the loss of volume. ? Regulatory agencies should provide better guidance on selecting input parameters. A consistent approach for selecting input parameters such as the peak rate factor and the velocity K parameter should be developed in order to reduce the user error associated with selecting the proper value. (Also see recommendations for future research below.) This research provides some insight into how the Muskingum and VPMC models route flood waves to downstream reaches in larger watershed. However, there are still many areas of uncertainty to investigate and potential areas to improve these models. The following list presents some recommendations for future research: ? Previous research conducted on the Muskingum and VPMC models used generic inflow hydrographs to generate the downstream hydrograph. In order to assess the accuracy of the model output, the outflow hydrograph was compared to hydrographs generated using other numeric techniques?for example, Wang et al. (2006) compared results to the results obtained from the Crank-Nicolson equation. Little, if any, research has been done using real upstream and downstream hydrographs to evaluate these models. Real hydrograph data should be used in future research for comparison to simulated results. 88 ? The constant parameter M-C method should be further studied and compared to both the VPMC and Muskingum models. ? More research should be conducted on the parameters used to generate the upstream hydrograph. This research should focus on developing better criteria for selecting input parameters in order to reduce user error (see Recommendations above). For example, the SCS unit hydrograph method uses a default peak rate factor of 484. However, there is a significant amount of variability in this parameter, and blindly using this value could introduce significant error to the inflow hydrograph. Similarly, uncertainty in the velocity K parameter used in the velocity method to calculate the time of concentration should be further investigated. 89 CHAPTER VI REFERENCES American Society of Civil Engineers (ASCE). 2008. Comments of the American Society of Civil Engineers on the Proposed Revisions to the Economic and Environmental Principles and Guidelines for Water and Related Land Resources Implementation Studies. American Society of Civil Engineers, Washington, DC. June 5. Bajracharya, K., and Barry, D.A. 1997. ?Accuracy Criteria for Linearised Diffusion Wave Flood Routing.? J. of Hydrol., 195(1-4): 200-217. Barry, D.A., and Bajracharya, K. 1995. ?On the Muskingum-Cunge Flood Routing Method.? Environment International, 21(5): 485-490. Bravo, R., Dow, D.A., and Rogers, J.R. 1994. ?Parameter Determination for the Muskingum-Cunge Flood Routing Method.? J. of the AWRA, 30(5): 891-899. Chapra, S.C. 1997. Surface Water-Quality Modeling. McGraw-Hill, New York. Das, A. 2004. ?Parameter Estimation for Muskingum Models.? J. or Irr. and Drainage Eng., 130: 140-147. Dunne, T. and Leopold, L.B. 1978. Water in Environmental Planning, W.H. Freeman and Company, New York. Gelegenis, J.J., and Serrano, S.E. 2000. ?Analysis of Muskingum Equation Based Flood Routing Schemes.? J. of Hydrol. Eng., 5(1): 102-105. Gill, M.A. 1978. ?Flood Routing by the Muskingum Method.? J. of Hydrolo., 36: 353- 363. Hydrologic Engineering Center (HEC). 1998. HEC-1 Flood Hydrograph Package-User?s Manual. U.S.Army Corps of Engineers, Davis, California. McCarthy, G.T. 1938. ?The Unit Hydrograph and Flood Routing.? Presented at the Conference of North Atlantic Division, US Army Corps of Engineers. McCuen, R.H. 1989. ?Calibration of Hydrologic Model for Maryland?. FHWA/MD- 91/03. MDOT, SHA, Brooklandville, MD. McCuen, R.H. 1998. Hydrologic Analysis and Design, Second Edition, Prentice Hall, Upper Saddle River, New Jersey. 90 McCuen, R.H. 2003. Modeling Hydrologic Change: Statistical Methods, CRC Press, Boca Raton, Florida. McCuen, R.H. and Z. Knight. 2006. ?Fuzzy Analysis of Slope-Area Discharge Estimates.? J. of Irr and Drainage Eng., 132(1): 64-69. Merkel, W.H. 2002. ?Muskingum-Cunge Flood Routing Procedure in NRCS Hydrologic Models.? Presented at the Second Federal Interagency Hydrologic Modeling Conference. July. National Weather Service (NWS). 2005. ?Unit Hydrograph (UHG) Technical Manual.? National Weather Service?Office of Hydrology, Hydrologic Research Laboratory and the National Operational Hydrologic Remote Sensing Center. Web address: http://www.nohrsc.noaa.gov/technology/gis/uhg_manual.html. Last updated October 12, 2005. Natural Resources Conservation Service (NRCS). 2007a. National Engineering Handbook, Part 630-Hydrology ,Chapter 16-Hydrographs. U.S. Department of Agriculture, Natural Resources Conservation Service. March. Natural Resources Conservation Service (NRCS). 2007b. National Engineering Handbook, Part 654-Stream Restoration Design. U.S. Department of Agriculture, Natural Resources Conservation Service. August. Perumal, M. and Sahoo, B. 2008. ?Volume Conservation Controversy of the Variable Parameter Muskingum-Cunge Method.? J. of Hydraul. Eng., 134(4): 475-485. Ponce, V.M. 1994. Engineering Hydrology, Principles and Practices, Prentice Hall, Upper Saddle River, New Jersey. Ponce, V.M., and Changanti, P.V. 1994. ?Variable-Parameter Muskingum-Cunge Method Revisited.? J. of Hydrol, 162(3-4): 433-439. Ponce, V.M., Lohani, A.K., and Scheyhing, C. 1996. ?Analytical Verification of Muskingum-Cunge Routing.? J. of Hydrol, 174: 235-241. Ponce, V.M., and Yevjevich, M. 1978. ?Muskingum-Cunge Method with Variable Parameters.? J. Hydr. Div., ASCE, 104(HY12): 1663-1667. Sobol, I.M. 1994. A Primer on the Monte Carol Method, CRC Press, Boca Raton, Florida. 91 Soil Conservation Service (SCS). 1973. A Method for Estimating Volume and Rate of Runoff in Small Watersheds. SCS-TP-149. U.S. Department of Agriculture, Soil Conservation Service. April. Soil Conservation Service (SCS). 1992. TR-20 Computer Program for Project Formulation Hydrology. Revised by the Hydrology Unit and the Technology Development Support Staff. February. Tang, X.N., Knight, D.W., and Samuels, P.G. 1999. ?Volume Conservation in Variable Parameter Muskingum-Cunge Method.? J. of Hydraul. Eng., 125(6): 610-620. Tewolde, M.H., and Smithers, J.C. 2006. ?Flood Routing in Ungauged Catchments Using Muskingum-Methods.? Water SA, 32(3): 379-388. Todini, E. 2007. ?A Mass Conservative and Water Storage Consistent Variable Parameter Muskingum- Cunge Approach.? Hydrol. Earth Syst. Sci., 11: 1645-1659. U.S. Army Corps of Engineers (USACE). 1994. Engineering and Design Flood-Runoff Analysis. Engineer Manual 1110-2-1417. Department of Army, U.S. Army Corps of Engineers, Washington, DC. August. U.S. Environmental Protection Agency (USEPA). 1997. ?Guiding Principles for Monte Carlo Analysis.? Prepared by the Risk Assessment Forum. EPA/630/R-97/001. March. U.S. Environmental Protection Agency (USEPA). 2003. ?Draft Guidance on the Development, Evaluation, and Applications of Regulatory Environmental Models.? Prepared by the Council for Regulatory Environmental Modeling. November. U.S. Fish and Wildlife Service (USFWS). 2003. ?Maryland Stream Survey: Bankfull Discharge and Channel Characteristics of Streams in the Coastal Plain Hydrologic Region.? Chesapeake Bay Field Office, Annapolis, MD. CBFO-S03-02. U.S. Water Resources Council (USWRC). 1981. Estimating Peak Flow Frequencies for Natural Ungaged Watersheds - A Proposed Nationwide Test: Hydrology Subcommittee, U.S. Water Resources Council, Washington, DC, 346 p. U.S. Water Resources Council (USWRC). 1983. Economic and Environmental Principles and Guidelines for Water and Related Land Resources Implementation Studies. U.S. Water Resources Council, Washington, DC. Wang, G.T., Yao, C., Okoren, C., and Chen, S. 2006. ?4-Point FDF of Muskingum Method Based on the Complete St Venant Equations.? J. of Hydrol, 324: 339-349. 92 Wu, J.S., King, E.L., and Wang, M. 1985. ?Optimal Identification of Muskingum Routing Coefficients.? J. of the AWRA, 21(3): 417-421.