HYDROLOGIC ANALYSIS OF A SEMI-ARID WATERSHED USING KINEMATIC WAVE AND SCS FLOW MODELS By Atiq Ur-Rehman Syed A DISSERTATION Submitted to Michigan State University in Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY Biosystems Engineering 2012 ABSTRACT HYDROLOGIC ANALYSIS OF A SEMI-ARID WATERSHED USING KINEMATIC WAVE AND SCS FLOW MODELS By Atiq Ur-Rehman Syed This research explores the efficiency of kinematic wave and Soil Conservation Service (SCS) flow models at a watershed scale in a semi-arid environment. The scope of this research is based on the hypothesis that flow models based on the simplest approximation of the full dynamic equations (kinematic wave, hydraulic) produce output variables that are representative of the watershed system compared to flow models that rely only on the continuity equation (SCS, hydrologic). The overall objective of this research study is to provide an improved understanding of kinematic wave and SCS flow models and compare their efficiencies to the observe flow data. Physical data such as precipitation, runoff, soils, and topography was derived from the Walnut Gulch Experimental Watershed (WGEW) in the southwest United States. Several important conclusions have emerged from this study that can prove useful to a practicing engineer/hydrologist. First, the kinematic-wave model proved to be a satisfactory tool to predict surface runoff in semi-arid watersheds, where transmission losses are a significant factor besides initial abstraction in the overall water budget computations. Analysis of the “Peak-Weighted Root Mean Squared Error” (PRMSE) values between the computed models (kinematic wave and SCS flow) and observed flow data for the three study watersheds show that the kinematic wave flow model has lower values of objective function compared to SCS flow model. Since PRMSE function is an implicit measure of comparison of the magnitudes of the peaks, volumes, and times of peak of the two hydrographs, it means that the kinematic wave flow model is more accurate than the SCS flow model. Second, the percent difference in peak flows between the observed data and computed flow results indicates that the kinematic wave model is no more likely to over-predict than to under-predict. On the other hand, the majority of the percent difference in peak flows between the observed and the SCS flow model indicates that the SCS model has a strong tendency to under predicted peak flows. Finally, the kinematic wave accuracy is demonstrated with data encompassing a relatively wide range of field conditions, where the kinematic wave flow model proved advantageous in that it can process spatial and/or temporal rainfall and overland and channel roughness variations, which the SCS model, by virtue of it being a lumped model, cannot. DEDICATION I wish to dedicate this work to my wife Rubina and my children Hamza and Zumruck who all helped in their own way, and to my mother who inspired my love of higher education at an early age. iv ACKNOWLEDGEMENTS To my advisor, Dr. Steve Safferman I owe a debt of gratitude for his patience and advise in completion of this work. To the members of my guidance committee, I must thank them as a team that as a whole provided more than guidance. I individually wish to thank my co. advisor, Dr. Pouyan Nejadhashemi for his unrelenting insight of hydrologic analyses and methodology and taking this research work to a new level; Dr. David Lusch for his efforts to provide direction to and synthesis of this work; Dr. John Bartholic for his oversight and support and Dr. Larry Segerlind for his encouragement. I also wish to thank the USDA / Forest Service and the Southwest Watershed Research Center (SCRC), Tucson, AZ, for providing me the resources for completion of this research. v TABLE OF CONTENTS LIST OF TABLES ....................................................................................................................... viii  LIST OF FIGURES ....................................................................................................................... ix  LIST OF SYMBOLS ...................................................................................................................... x  Chapter 1: INTRODUCTION......................................................................................................... 1  1.1.  Overview ..........................................................................................................................1  1.2.  Scope and Objectives .......................................................................................................4  1.3.  Physical Settings ..............................................................................................................7  Chapter 2: REVIEW OF THEORY AND LITERATURE............................................................. 8  2.1.  Hydrologic Modeling .......................................................................................................8  2.1.1.  Hydrologic Models with Distributed Parameters ..................................................12  2.1.2.  Hydrologic Models with Lumped Parameters .......................................................13  2.1.3.  Use of GIS in Hydrologic Models .........................................................................14  2.2.  Flow Routing Models ....................................................................................................16  2.2.1.  Hydraulic Flow Routing ........................................................................................17  2.2.2.  Hydrologic Flow Routing ......................................................................................27  Chapter 3: METHODOLOGY...................................................................................................... 35  3.1.  Synopsis .........................................................................................................................35  3.2.  Flow Model Development .............................................................................................38  3.2.1.  Development of Kinematic Wave Flow Model .....................................................38  3.2.2.  Development of SCS Flow Model .........................................................................43  3.3.  Watershed Characteristics and Development of Model Input Data...............................51  3.3.1.  Drainage Area Calculations ...................................................................................51  3.3.2.  Watershed Lengths.................................................................................................52  3.3.3.  Watershed Slopes ...................................................................................................53  3.3.4.  Soils Survey and Topographic Data (Digital Elevations) ......................................54  3.3.5.  Precipitation Data...................................................................................................54  3.4.  Development of a Grid Based Curve Number Estimator for Runoff Volume Computations.................................................................................................................57  Chapter 4: RESULTS, ANALYSES, AND DISCUSSIONS OF MODEL RESULTS ............... 61  4.1.  Results and Analysis of Precipitation and Runoff Data.................................................61  4.2.  Results and Analysis of Model Output Data..................................................................62  4.2.1.  Analyses of Observed and Computed “Peak-flows” .............................................62  4.2.2.  Analyses of Observed and Computed “time to peak” values ................................63  4.2.3.  Peak Weighted Root Mean Squared Error .............................................................64  4.2.4.  Analyses of Observed and Computed “Center of Mass” values............................66  vi Graphic/visual Comparison of Observed and Computed Hydrographs ................67  4.2.5.  4.3.  Discussions of Model Results ........................................................................................68  Chapter 5: CONCLUSIONS and RECOMMENDATIONS ........................................................ 71  5.1.  Conclusions ....................................................................................................................71  5.2.  Recommendations ..........................................................................................................72  5.3.  Future Applications/Benefits .........................................................................................74  APPENDICES .............................................................................................................................. 76  APPENDIX A: GIS MODEL INPUT DATA .......................................................................... 77  APPENDIX B: RAINFALL AND RUNOFF DATA .............................................................. 78  REFERENCES ............................................................................................................................. 79  vii LIST OF TABLES Table 1: Study Drainage areas. ..................................................................................................... 52  Table 2: Study watersheds lengths................................................................................................ 53  Table 3: Study watersheds slopes. ................................................................................................ 53  Table 4: Definition of hydrologic soils group (adapted from USACE-HEC-HMS, Manual, 2000). ................................................................................................................59  Table 5: Curve number estimates for condition I, II and III (adapted from USACE-HECHMS, Manual, 2000). .....................................................................................................60  viii LIST OF FIGURES Figure 1: Simplified watershed with kinematic wave representation. ...........................................39  Figure 2: Kinematic wave parameters for various channel shapes. (USACE, 1998). ...................43  Figure 3: Figure 3: Conceptual watershed illustrating travel time from the centroid (grey dot) of each band of area to the watershed outlet, (adapted from NRCS, NEH Manual, 2004). ..............................................................................................................50  ix LIST OF SYMBOLS Symbol Definition A cross-sectional area Ac cross-sectional flow area B channel surface water surface width B y t rate of rise CN runoff curve number F actual retention g acceleration due to gravity I inflow rate to the reach Ia initial abstraction i-f rate of excess rainfall mo value of 5/3 from the Manning’s equation N a resistance factor that depends on the cover of the planes O outflow rate from the reach P total precipitation Q discharge rate qo overland flow per unit length (q p) is the peak discharge R hydraulic radius S potential maximum retention t time t change in time interval (t p ) the time to peak (tr ) recession time x Symbol Definition V velocity v average flow velocity of the reach y hydraulic depth yo mean depth of overland flow x distance along the flow path  and m y = x V V = g x 1 V = g t V = x y VB = x A parameters related to flow geometry and roughness pressure gradient convective acceleration local acceleration prism storage wedge storage xi Chapter 1: INTRODUCTION 1.1. Overview Arid and semiarid regions constitute over a third of the world’s landmass yet are under increasing population pressure. In the semiarid southwestern U.S. population is projected to increase over 50% by 2030 in comparison to 5-15% in other U.S. regions (SWRC and WGEW, 2007). This will dramatically increase society’s need to manage its water, soil, and nutrient resources to support people, agriculture and the environment. To effectively manage watershed resources under the stress of increasing population and climatic variations, hydrologic research has been the subject of increasing attention in recent years. The term “watershed hydrology” is defined as the branch of hydrology that integrates hydrologic processes at the watershed scale to determine the watershed response (Singh, 2002). The concept of watershed is basic to all hydrologic designs. Large watersheds are made up of many smaller basins but it is necessary to define the watershed in terms of a point that is usually located at the watershed “outlet”. With respect to the outlet, the watershed consists of all land area that “sheds” water to the outlet during a precipitation event. Using this concept a watershed is defined by all points enclosed within an area from which precipitation falling at these points will contribute water to the outlet. The conceptual representation of hydrologic processes in a watershed can be defined in a systematic manner with the help of a mathematical or physical model. Modeling can be mathematical, if described by mathematical equations or physical, if a scale model is built to represent dimensional similitude to the actual watershed. In either case the model is a conceptualization of the actual watershed. Traditionally, hydrologic analyses have been 1 performed either through stochastic or deterministic modeling techniques. Stochastic models use statistical concepts that links a specific input like rainfall to the model output such as runoff using regression and neural networks (Vieux, 1988). Such a model is often referred to as a “black box” because the fundamental relationships causing the effect are not considered. Deterministic models use mathematical representations of the underlying regularities that are produced by the entities being modeled and generate theoretically perfect data (Chow et.al, 1988). Deterministic hydrologic models are further classified into two categories, lumped and distributed parameters. Lumped parameter models utilize the average of a set of independent variables that represent a sub-basin or an entire watershed, whereas distributed parameter models utilize the spatial location of the independent variables and compute the dependent variable directly at the spatial location of each independent variable (Hann et.al, 1994). The origins of mathematical hydrologic models dates back to the rational method developed by Mulvany (1850) and an “event” model by Imbeau (1892) for relating storm runoff peak to rainfall intensity. In 1932, Sherman introduced the concept of “Unit Hydrograph”, which relates direct runoff response to the rainfall excess and provides a means for manipulating a known volume of runoff from a basin to represent the timing when the volume arrives at the basin outlet. In this method, the ordinates of all graphs of runoff for unit time are directly proportional to the net (effective) depths of rainfall in that same time. About the same time, Horton developed the theory of infiltration and improved hydrograph separation techniques. Horton suggested that infiltration capacity rapidly declines during the early part of a storm and then tends towards an approximately constant value after a couple of hours for the remainder of the event. This results because previously infiltrated water fills the available storage spaces and reduces the capillary forces drawing water into the pores. Clay particles in the soil may swell as 2 they become wet and thereby reduce the size of the pores. In areas where the ground is not protected by a layer of forest litter, raindrops can detach soil particles from the surface and wash fine particles into surface pores where they can impede the infiltration process. In 1945, a month before his death, Horton published his work in The Bulletin of the Geological Society of America, known as Horton’s Laws, which constitutes the foundation of quantitative geomorphology. This work is considered the founding of modern stream chemistry modeling, since it was the first comprehensive set of mathematical models to link basin hydrology with a water pollutant, namely sediment. During the 1950s, hydrology was approached with a more theoretical basis than in the past, facilitated by advances in the physical understanding of hydrological processes and by the advent of computers. Later in 1956, the Soil Conservation Service (SCS), now called the Natural Resources Conservation Service (NRCS) of the U.S. Department of Agriculture developed a hydrologic model known as the SCS-curve number method. Rallison (1980) gives a detailed synopsis of the development of this procedure from its inception to its final form and application to ungauged watersheds. Work on this procedure began in the 1950s in response to the passage of the “Watershed Protection and Flood Prevention Act (P.L. 83-566). Due to work authorized by this act, SCS anticipated the need for a simplified method of hydrologic computation. Based on extensive analyses of gauged, experimental watersheds and infiltrometer studies, a relation between rainfall and runoff was developed (Richard H. McCuen, 1982). A major effort employed the theory of linear systems, which led to the theory of the instantaneous unit hydrograph by Nash in 1957, and then the generalized theory by Dooge in 1959 (Singh, 2000). The decade of the 1960s witnessed the digital revolution that made possible the integration of models of different components of the hydrologic cycle and simulation of virtually the entire 3 watershed (Singh, 2000). Examples of models that became popular are the watershed models of Dawdy and O’Donnell (1965) and HEC-1 (Hydrologic Engineering Center 1968). HEC-1 was developed by the US Army Corps of Engineers to estimate river flows as a result of rainfall. It was written in the FORTRAN language and, until 1984, could only be run on a mainframe computer. The HEC-1 was replaced with the HEC-HMS (Hydrologic Engineering CenterHydrologic Modeling System) beginning in 1992, which has long been considered a standard for hydrologic simulation. The new HEC-HMS provides almost all of the same simulation capabilities, but uses advances in numerical analysis that take advantage of the significantly faster computers available today. Apart from the computational capabilities, hydrologic models require physiographic information such as location of drainage divides, channel lengths, slopes, and sub-basin geometric properties. Historically, most of these parameters were obtained from maps or field surveys but during the last two decades this information has been increasingly derived from Geographic Information System (GIS) data sources (Jenson and Domingue, 1988). 1.2. Scope and Objectives The history of this research dates back to the Red Cedar River project conducted during 2002 on the Michigan State University (MSU) campus. In the Red Cedar River project, an approach was developed to understand how specific land uses affect water quality using one of the basins as a pilot study area. The ultimate goal was to model the transport of chemical and biological pollutants introduced to surface water as a result of different land uses in the study watershed. For the hydrologic model component, the kinematic wave flow routing theory was applied to simulate surface runoff in a storm sewer network (McElmurry, S.P., Aslam, I., Syed, A.U., Syed 4 et.al, 2003). In order to achieve this, it was first necessary to accurately describe the hydrologic characteristics of the watershed system. This was achieved by delineating sub-basins using a micro-watershed approach based on inlet locations of the storm sewer networks, which ultimately discharged into the Red Cedar River. The micro-watershed approach, developed through this research, provided an accurate description of the hydrologic processes in a watershed. In order to model the complex nature of flows generated from surface runoff into the storm sewer network, the application of kinematic wave theory proved useful. The application of kinematic wave theory in the Red Cedar River project achieved a 17percent difference in flow volume (cubic-meters) and a 12-percent difference in peak flow (cubic meters per sec) between the computed and observed flow hydrographs. Although this approach achieved reasonable results compared to the actual flow data, its application was limited to an urban setting dominated by impervious surfaces and flow pathways through storm sewer pipes of known geometric shapes. Consequently, its range of applicability could not be fully explored in this research project. Even though the kinematic wave method has gained wide acceptance both in the US and abroad in solving a variety of hydrologic engineering problems there is a continuing controversy regarding its accuracy and applicability in natural settings (i.e. mixed land use watersheds). Researchers and practitioners alike have reported successes and failures of the model, with papers continuing to appear in the literature describing what it can and cannot do (Hromadka and DeVries, 1988). A current area of discussions focuses on whether the kinematic wave approach can eventually replace other well established methods of surface flow generation such as the SCS unit hydrograph method. In order to answer some of these research questions this research tests the applicability of kinematic wave theory in a mixed watershed setting and examines its 5 efficiency in a semi-arid environment compared to the well-established SCS unit hydrograph method. The scope of this research is based on the hypothesis that flow models based on the simplest approximation of the full dynamic equations (kinematic wave, hydraulic) produces more reliable results compared to the SCS method in mixed land use watersheds. Validation of this hypothesis will be accomplished if the flow model based on the simplified form of the full dynamic equations (kinematic wave) accurately predicts the observed flow hydrograph, compared to the flow model based on the continuity equations (SCS). The specific objectives of the research are:  Objective 1: Develop a deterministic watershed model using spatially variable data.  Objective 2: Develop surface runoff flow models using kinematic wave and SCS methods.  Objective 3: Analyze the accuracy of both flow models by comparing the computed flow hydrographs from the kinematic wave and SCS models to the actual flow hydrographs (observed data) and validate the research hypothesis for event based storms using “goodnessof-fit” criteria between the computed and actual flow hydrographs. The conclusions from this study will lead to a better understanding of both the kinematic wave and SCS flow models and identify their limitations and strengths for effective management of water resources and the environment in semi-arid watersheds.   6 1.3. Physical Settings In this research physical data such as precipitation, runoff, soils, and topography are derived from the Walnut Gulch Experimental Watershed (WGEW), which is shown in Figure 1. The Southwest Watershed Research Center (SWRC) operates the Walnut Gulch Experimental Watershed in southeastern Arizona as an outdoor laboratory for studying semi-arid rangeland hydrologic, ecosystem, climate, and erosion processes. The WGEW encompasses the 150 square kilometers in southeastern Arizona, that surrounds the historical western town of Tombstone is contained within the upper San Pedro River Basin that encompasses 7600 square kilometers in Sonora, Mexico and Arizona. The watershed is representative of approximately 60 million hectares of brush and grass covered rangeland found throughout the semi-arid southwest and is a transition zone between the Chihuahuan and Sonoran Deserts. Elevation in the watershed range from 1250 m to 1585 m MSL. Cattle grazing is the primary land use with mining, limited urbanization, and recreation making up the remaining uses. Hydro-meteorologic and soil erosion/sedimentation data are collected from 125 instrumented locations on WGEW. Precipitation is measured with a network of 88 weighingtype recording rain gauges arranged in a grid throughout the watershed (Figure 2). Various runoff measuring structures are used to monitor small watersheds (< 40 ha) runoff. These structures include a broad-crested V-notch weir, H-flumes, and Santa Rita supercritical flow flumes. The largest flume, at the outlet of the WGEW has a flow capacity of 650-cubic meters/sec. 7 Chapter 2: REVIEW OF THEORY AND LITERATURE The focus of this research is to investigate the accuracy of kinematic wave and SCS flow models. The hypothesis that flow models based on the simplest approximation of the full dynamic equations produces output variables that are representative of the natural system compared to flow models that rely only on the continuity equations. This review expounds the theory in literature on deterministic models of overland flow and hydraulic and hydrologic flow routing methods. 2.1. Hydrologic Modeling The physically correct representation of the surface runoff and flow routing processes in a watershed depends on many factors. To categorize these factors, several distinctions should be made, in general, as to the modeling process that seeks to represent the physical processes. A conceptual model is the first step towards the unknown. The mathematical model then describes those essential processes contained in the conceptual model. We will examine first some mathematical models that seek to model deterministically the physical process of the rainfall event. This process may include overland flow of the rainfall excess and channel routing of the lateral inflow from the overland flow portion of the watershed. The overland flow/surface runoff in deterministic models is generated when the rainfall intensity exceeds the infiltration rate. The overland flow travels over the ground surface to the main channel, from where it moves to the watershed outlet as channel/stream flow. The combined overland/surface runoff and channel flow is the total watershed outflow. The different flow pathways originating as precipitation excess and resulting in watershed outlet flow are typically examined using flow routing models/techniques. In technical terms, flow routing refers 8 to the tracking in time and space of a wave characteristic such as a peak discharge or stage as it moves along the flow path but superimposed on the physical flow itself (Sturm, 2001). In general, routing techniques may be classified into two categories: hydraulic routing and hydrologic routing. Hydraulic routing is a distributed system method, which determines the flow as a function of both space and time (Chaudhry, 1993), and is based on the solution of the partial differential equations of unsteady open channel flow. These equations are often referred to as the St. Venant equations or the dynamic wave equations. Hydrologic routing employs the continuity equation and either an analytical or an empirical relationship between storage within the reach and discharge at the outlet. In hydrologic routing the momentum equation is integrated spatially in the flow direction so it becomes a lumped system spatially, with no variation of parameters within the resulting control volume (Chow, Maidment, and Mays 1988). Smith and Woolhiser (1971) developed a mathematical model that simulated a coupled system of two complex, natural processes on an elemental watershed. The conceptual model included only infiltration and overland flow. Channel routing was not included because the conceptual model was limited to the upland portion of the watershed where channel flow is not well established. The infiltration model provides insight into the process by which rainfall becomes either runoff or subsurface water. The kinematic equations provided insight into the depth and velocity of the runoff as it is accelerated down the watershed. The infiltration and kinematic equations were coupled mathematically such that the rainfall excess, as defined by the infiltration model, was the boundary value for the solution of the kinematic equation. A distributed, deterministic system results when the inherent spatial nature of the processes are preserved in the solution method. A model such as the Smith and Woolhiser’s (1971) provides 9 the opportunity to model the outflow of the watershed, and more importantly the spatial and temporal distribution of the runoff-infiltration processes within that watershed. During the period 1970-1995, several state of the art papers dealing with watershed modeling appeared. Kibler and Woolhiser (1970) made an analysis of the kinematic cascade as a distributed parameter mathematical watershed model. They were able to determine the numerical phenomenon of the kinematic shock. According to the simple wave theory, when there is a change in slope between the planes in the cascade, a shock or wave front is propagated within the system. The shock represents the numeric difficulty in the computation of the hydrograph. C. L. Chen and Ven Te Chow (1971) proposed the hydrodynamic approach and considered watershed hydrology as a distributed continuum, where the hydrodynamic principles of fluid flow apply. The hydrodynamic equations have been derived and solved by various methods. There are two distinct categories of flow in a watershed, 1) overland flow, and 2) channel flow characterized by well-drained channel geometry. The boundary between these two flow changes with time and distance and therefore is difficult to model. Chen and Chow (1971) formulated a comprehensive watershed flow model. They classified watershed hydrology by a molecular, microscopic hydrodynamic and macroscopic hydrodynamic approach. Both of these approaches derive the Navier-Stokes equation of motion for fluid flow with suitable boundary equation. Clark (1973) discussed important issues regarding model identification and diagnosis and parameter estimation and showed that interdependence between model parameters required extensive exploration of error objective function, particularly when the model is used to determine the likely effects of land-use change. 10 Huggins and Barney (1982) observed that hydrologic modeling is most differentiated by the manner in which parameter or input values are handled. They identified that distributed parameter models treat the individual input parameters directly without lumping. Such models avoid the errors caused by averaging of nonlinear variables or threshold values (Barney and Higgins, 1982). El-kady (1989) reviewed numerous watershed models and concluded that the surface water-groundwater linkage needed improvement, while ensuring an integrated treatment of complexity and scale of individual component processes. Goodrich and Woolhiser (1991) reviewed progress in catchment hydrology in the United States and emphasized that a detailed process based understanding of hydrologic response over a range of catchment scales still eluded the hydrologic community. Hornberger and Boyer (1995) identified the need and importance of spatial variability and scaling and the linkages among hydrology, geochemistry, environmental biology, meteorology, and climatology. They discussed the use of digital elevation models (DEMs), and raised the question of subgrid variability and the effects of pixel size on model calibration. Jayawardena et. al. (2006) did a comparative analysis of data-driven and GIS-based conceptual rainfall-runoff model. They investigated the suitability of a conceptual technique along with a data-driven technique, to model the rainfall-runoff process. The conceptual technique used is based on the Xinanjiang model coupled with GIS for runoff routing and the data-driven model is based on genetic programming (GP), which was used for rainfall-runoff modeling in the recent past. They verified that for a small, steep-sloped catchment the conceptual model outperformed the data-driven model and provided a better representation of the rainfallrunoff process in general, and better prediction of peak discharge, in particular. 11 2.1.1. Hydrologic Models with Distributed Parameters The principal advantage of the distributed models is that the geographical variation of data within the watershed is preserved. However, these models are complex, require more computing time, and increase input data. However, today’s computing technologies, allow their use. Huggins and Barney (1982) observed that the hydrologic modeling is most differentiated by the manner in which parameters or input values are handled. Lumping or averaging certain parameters yields a lumped parameter model. Distributed parameter watershed models treat the individual input parameters directly without lumping. Such models avoid the errors caused by averaging of non-linear variables or threshold values (Huggins, L.F. and J.R. Burney, 1982). A.S. Donigian, Jr., B.R. Bicknell and J.C. Imhoff (2005) describes the HSPF (Hydrologic Simulations FORTRAN-Programming) model as a distributed model that can simulate the continuous, dynamic, or steady state behavior of both hydrologic/hydraulic and water quality processes in a watershed, with an integrated linkage of surface, soil, and stream processes. HSPF is commonly recognized as the most complete and defensible process-based watershed model for quantifying runoff and addressing water quality impairments associated with combined point and nonpoint sources (Bicknell et al., 2005). HSPF contains hundreds of process algorithms developed from theory, laboratory experiments, and empirical relations from instrumented watersheds. The model consists of a set of modules arranged in an organized structure, which permit the continuous simulation of a comprehensive range of hydrologic and water quality processes. HSPF's design incorporates a hierarchy of program subroutines, each of which performs a major task during the program's execution. The subroutines are grouped into different levels of operations in a hierarchical structure. The importance of this program structure lies in its modular design. This allows for the addition and/or replacement of individual 12 modules and allows HSPF to be easily adapted for special applications designed by the user. HSPF has been applied to watersheds ranging in size from the Chesapeake Bay, with roughly 99,780 square kilometers of tributary area, down to a few square-meters. HSPF can simulate any period from a few minutes to hundreds of years on an hourly time step. It has been applied to such diverse climatic regimes as the tropical rain forests of the Caribbean, arid conditions of Saudi Arabia and the Southwestern U. S., the humid Eastern U.S. and Europe, and the snow covered regions of the Eastern Canada. 2.1.2. Hydrologic Models with Lumped Parameters The lumped parameter model simplifies the description of the behavior of spatially distributed physical systems in a topology consisting of discrete entities that approximate the behavior of the distributed system under certain assumptions. Mathematically speaking, the simplification reduces the state space of the system to a finite number and partial differential equations of the continuous time and space model of the physical system into ordinary differential equations with finite number of parameters. NRCS developed a physically based lumped parameter model to estimate direct runoff from ungaged watersheds. Rallison (1980) gives a detailed synopsis of the development of this procedure from its inception to its final form and application to ungaged watersheds. Work on this procedure began in the 1950s in response to the passage of the “Watershed Protection and Flood Prevention Act (P.L. 83-566)”. Due to work authorized by this act, SCS anticipated the need for a simplified method of hydrologic computation. Based on extensive analyses of gauged, experimental watersheds and infiltrometer studies, a relation between rainfall and runoff was developed (Richard H. McCuen, 1982). The basic relation was derived by plotting the accumulated natural runoff versus the accumulated rainfall. It was observed that the relation is 13 asymptotic to a line at a 45-degree slope. This shows that the runoff rate approaches the rainfall rate as the accumulation of both continues. In addition, maximum retention, the difference between rainfall and runoff, approaches a constant value. The SCS model estimates precipitation excess as a function of cumulative precipitation, soil cover, land use, and antecedent moisture. Another important element in the SCS model is time forecasting, the time required for the water to flow from the most remote (in time of flow) point of the area to the outlet once the soil has become saturated and minor depressions filled. This is evident from the fact that most hydrologic methods include a time variable as input. It is the time required for the water to flow from the most remote (in time of flow) point of the area to the outlet once the soil has become saturated and minor depressions filled. It is assumed that when the duration of the storm equals the time of concentration, all parts of the watershed are contributing simultaneously to the discharge at the outlet. With respect to the unit hydrographs: the time of concentration is the time from the end of rainfall excess to the point of inflection on the recession (Richard H. McCuen, 1982). 2.1.3. Use of GIS in Hydrologic Models The most recent advances in watershed modeling were achieved from the use of GIS, remotely sensed data, and environmental tracers (Singh, 2000). Apart from the computational capabilities, hydrologic, water-quality, and climatic models require physiographic information such as location of drainage divides, channel lengths, slopes, and sub-basin geometric properties. During the late 90s there was a noticeable shift in hydrologic modeling research from theoretical/conceptual improvements to data driven capabilities using GIS techniques. Historically, most of these parameters were obtained from maps or field surveys, but during the last two decades this information has been increasingly derived from GIS data sources (Jenson 14 and Domingue, 1988). Similarly, the use of Global Positioning System (GPS) provides many new and affordable options for the collection of large number of elevation data sets (Wilson 1999). Digital elevation data sets are usually organized into one of two main data structures, 1) regular grids, and 2) triangulated irregular networks (TINs). Square grid digital elevation models (DEMs) have emerged as the most widely used data structure during the past decade because of its simplicity and ease of computer implementation (Wise 1998). However, an inherent problem with square grid DEM data is the production of nonphysical depressions due to noise in the elevation data that affects interpolating schemes used to describe variation in elevation between raster points (DeVantier et.al., 1993). The results are an unwanted termination of drainage paths in pits during hydrologic analyses. These types of problems can be seen particularly in flat areas. Similarly, DEMs are dependent on the grid size for certain computed topographic parameters and are unable to adjust the grid size to the dimensions of topographic land surface features (Fairfield and Leymarie, 1991). To overcome some of the disadvantages of the grid DEM’s data structures, TIN, and contour-based structures are the preferred choice. A TIN is a vector based representation of the physical land surface, made up of irregularly distributed nodes and lines with three dimensional coordinates (x, y, and z) that are arranged in a network of nonoverlapping triangles (Vivoni et.al., 2004). TIN’s are often derived from the elevation data of a rasterized digital elevation model. Various factors motivate the use of TINs compared to DEMs to represent the watershed topography. Traditionally, terrain data in hydrologic models has been represented in two ways: (1) aggregating or resampling grid-based DEM’s to coarser resolutions, or (2) topographic distribution function that classifies catchment locations according to an elevation index (Vivoni 15 et.al, 2004). Both methods attempt to account for the spatial variability in topography without adding computational burden to hydrologic models that operate over large domains. Neither approach, however, can incorporate all the information on high-resolution topographic data currently available from land surveying, aerial photogrammetry (Gesch et al., 2002), or light detecting and ranging (LIDAR) (Ritchie, 1996). A GIS based hydrologic model requires accurate depiction of terrain features since the surface elevation properties (slope, curvature, aspect) determine the hydrologic response to a meteorologic forcing function. If the model domain increases in size, the resolution and accuracy retained in the terrain representation must decreases to allow efficient model simulation (Vivoni et.al, 2004). As a result, poorly resolved hydrologic models typically have terrain inaccuracies that propagate directly to model predictions of stream flow and soil moisture (Vieux 1993; Kuo et al., 1999). For climate, hydrology, and weather models operating at large spatial scale, inaccurate depiction of the topography and its spatial variability is recognized as an important source of model error (Koster et al., 2000). 2.2. Flow Routing Models Flow routing refers to the tracking in time and space of a wave characteristic such as a peak discharge or stage as it moves along the flow path but superimposed on the physical flow itself (Sturm, 2001). Flow routing procedures range in complexity from simple storage routing methods to relatively complex procedures based on simultaneous solutions to the hydrodynamic equations representing the conservation of momentum and mass. Flow routing methods are classified into a number of ways; one of the most important distinctions is between hydrologic routing and hydraulic routing. In hydrologic routing, the momentum equation is integrated 16 spatially in the flow direction so that it becomes a lumped system spatially, with no variation of parameters within the resulting control volume (Chow et al., 1988). Conversely, hydraulic routing is a distributed system method, which determines the flow as a function of both space and time (Chaudhry, 1993). 2.2.1. Hydraulic Flow Routing Chow (1959), Henderson (1966), Chow et. al (1988), and others discussed hydraulic flow routing based on equation of spatially varied (non-uniform), unsteady flow (Saint Venant equations). These continuity (1) and momentum (2) equation as given below; Q A  0 x t (1) Where: Q = flow rate (cubic-meters per sec.) t = time in secs x = distance measured in downstream flow direction A = cross-sectional area (meters squared) u u h  u  g  g ( So  S f )  0 t x x   --------------- ---------------------------------------------------------------  Kinematic term  Diffusion term  Dynamic term Where: u = x- component of mean velocity 17 (2) t = time in secs x = distance measured in downstream flow direction g = acceleration due to gravity h = mean depth S o = average bottom depth Sf = friction slope defined by manning’s n Dynamic routing includes all the terms in the momentum equation, while diffusion routing neglects the inertia terms (local and convective acceleration), and kinematic routing includes only the gravity and flow resistance terms. With regard to spatial variation, all hydraulic routing methods can be considered distributed models and the errors generated as a result of neglected terms from equation (2) are small relative to a full dynamic equation (Sturm, 2001). Hydraulic routing is usually accomplished by a numerical solution of the governing equations or by the method of characteristics. The relative importance of the various terms in the above equations determines largely the degree of simplification warranted. Hann et al. (1994) concludes that that the use of full momentum equation for routing is limited to situations involving backwater effects, tidal flows, surges, and flow junctions where large tributaries enter the main channel. Fread (1985) put together a large generalized one dimensional routing program known as FLDWAV that use the full momentum equations. Models like FLDWAV can only be used for complex flow routing situations like dam breach analyses and comes at the expense of computer size, speed, and input data requirements. Bradley et. al (1996) performed floodplain mapping using continuous hydrologic and hydraulic simulations. They estimated floodplain limits in their research. A hydraulic model with 18 flood-routing capabilities was coupled with a continuous-simulation hydrologic model to compute river stages. Simulations of historical precipitation data were used to produce long, continuous flow records. Simulations were also made using precipitation data for extreme storms to gain additional information on overbank flooding. They did a statistical analysis of the simulated peak stages to estimate peak stage exceedence probabilities. To illustrate the benefits of this approach, floodplain limits for an Illinois stream were estimated with the new method and compared to estimates based on a design storm approach and conventional frequency analysis was applied to the peak stages. Perumal and Raju (1998) proposed an approach for developing a simplified variableparameter stage-hydrograph routing method from the Saint Venant equations for routing floods in any shape of prismatic channel and flow following a generalized friction law. In this approach, they demonstrated that the parameters of the routing equation could be related to the channel and flow characteristics, which also enables the development of a theoretically based procedure for varying these parameters at every routing time level. They demonstrated simultaneous computation of the discharge hydrographs corresponding to a given input-stage and routed-stage hydrographs. Haktanir and Ozmen (1997) did a comparison of hydrologic and hydraulic routing methods for three long reservoirs in Turkey. They computed outflow hydrographs for three dams using both hydrologic routing (level-pool routing) and hydraulic routing methods. The results were compared with three inflow hydrographs of different peaks with three unregulated ogee spillways of different capacities. In all these cases, the difference between outflow hydrographs was greatest at the peak value, growing larger as the spillway capacity became smaller, relative to 19 the magnitude of the inflow hydrograph. They concluded that the peak outflow by hydraulic routing was smaller than that by hydrologic routing for all the routing combinations. Tseng et al. (2001) did analyses of channel routing with surges. They presented two high-resolution, shock-capturing schemes for the simulation of one dimensional, rapidly varied open-channel flows. The proposed algorithms were assessed using several steady and unsteady problems to verify its accuracy and robustness in capturing strong shocks in open-channel flows. They compared the results of dynamic flood routing and steady routing to demonstrate the risk of using steady routing for flood mitigation. Guo (2004) presented a hydrologic based approach to assess storm water detention basins using new routing schemes in which the storage-outflow curve was approximated from the inflow hydrograph to the basin and the maximum allowable release from the basin. This hydrologic procedure significantly simplified the storm water detention modeling technique. In addition, he rearranged the continuity principle to derive two new reservoir routing functions. Both functions provided a direct solution using all variables at the same time step without iteration. Xiong and Melching (2005) did a comparison of the kinematic wave and non-linear reservoir routing in an urban watershed. They used previously unpublished experimental data from a watershed experimentation station at the University of Illinois, Urbana Champaign, to test the accuracy of the two methods. For the non-linear reservoir method, the U.S. Environmental Protection Agency Storm Water Management Model (SWMM) was used. For the kinematic wave method, the Dynamic Watershed Simulation Model (DWSM) was used. The DWSM (Borah et al. 2002) is a computer program that uses kinematic wave routing theories to simulate water and sediment discharges from storms. They concluded that nonlinear reservoir routing 20 method as applied in SWMM may provide acceptable results for storms with durations longer than the watershed time of concentration (42 experiments) (average model fit efficiency E = 0.88), but for storms with durations less than or equal to the time of concentration (26 experiments), poor results were obtained (average E = 0.07). More accurate results generally were obtained using kinematic-wave routing (average E = 0.928 and 0.807 for storm durations exceeding and not exceeding the watershed time of concentration, respectively). These results fit with the theoretical basis of the kinematic-wave theory that considers the actual physical processes in surface flow generation, while the nonlinear reservoir method does not consider the impact from the time lag needed for the flow depth to grow so that runoff can commence. Borah et al. (2009) did a review of the mathematical-based watershed models to provide insightful facts for the end users, so that an informed decision can be made when selecting a model. They reviewed fourteen watershed models and formulated six summary tables that compile, rank, and compare these models. The dominant procedures used in the models were compared, including rainfall excess or infiltration, overland runoff routing, channel flow routing, and subsurface flow simulation. The rankings were based on the formulations, their relative complexities, and accuracies for routing overland runoff and channel flows. These compilations and rankings are very helpful to managers and modelers in understanding and comparing the models, selecting the most suitable model for a project or an application, and using it to its full potential. 2.2.1.1. Kinematic Wave Theory The kinematic wave theory was developed by Lighthill and Whitham (1955) describes flood movement in long rivers. This method has gained wide acceptance worldwide and is used in solving a variety of applied hydrologic engineering problems. Kinematic waves are often 21 classified as uniform, unsteady flows. It represents the changes in discharge, velocity, and water surface elevation with time at any one location along a stream channel or overland flow (Bedient and Huber, 1992). The Kinematic wave is based on the assumption that inertial and pressure does not affect the flow routing. Instead it is dependent on the weight or gravity force of fluid which is approximately balanced by the resistive forces of the bed friction. Kinematic wave flows in the downstream direction without crest subsidence and without much acceleration. The simplification of the full dynamic equation is achieved by combining the continuity and momentum equations with inertia and pressure terms dropped, from equation (2). u u h u g 0   t x x (3) And: So  S f  u2 (4) 2 c h Where: So = average bottom depth Sf = friction slope defined by manning’s n c = wave celerity The derivative in the above equation can be eliminated by differentiating it w.r.t. x and t separately, 22 2u u u 2 h 0  2 x c h c 2 h 2 x (5) 2u h u 2 h 0  2 t c h c 2 h 2 x (6) This yield, h 2h u  x u x (7) h 2h u  t u t (8) Substituting the above equation into the continuity equation yields, 2h u 2h u u 0 u h u t u x x (9) Or u 3 u  u 0 t 2 x (10) Or 23 u u c 0 t x (11) Where: c 3 u 2 (12) Equation (12) is the kinematic wave equation. For a channel of arbitrary shape, the above equation becomes a a 0  C1 x t (13) Where: a = cross sectional area C1 = roughness factor And C1 is given as: C1  S f a (14) 1  S f q Where: Sf = friction slope defined by manning’s n 24 q = lateral inflow per unit length The kinematic wave theory in the HEC-HMS (Hydrologic Engineering CenterHydrologic Modeling System), model describes the flow as a function of depth only, for all x and t, given that there are no appreciable backwater effects. Q  * ym (15) Where: Q = discharge in cfs  , m  Kinematic wave routing parameters The Kinematic wave equations for overland flow, on a wide plane with shallow flows can be derived from equation (8) and the Manning equation for overland flow is shown below 5 1.49 q S o Yo 3 N (16) An important consideration in using the above equation is to use large Manning’s roughness coefficient values compared to the one used in channels (Bedient and Huber, 1992). Rewriting equations in terms of flow per unit width for overland flow, we have q o   o y o mo (17) Where: 25 o= 1.49 So N = Conveyance factor, mo = 5/3 from the Manning’s equation, So = average overland slope, y o = mean depth of overland flow. The continuity equation is, y o q o  i f x t (18) Where: i-f = rate of excess rainfall (ft/sec), q o = flow rate per unit width (cfs/ft), y o = mean depth of overland flow (ft). 26 By substitution of equation (17) into (18) we have; y o mo 1 y o   o mo y o i f t x (19) Equations (19) form the complete kinematic wave equations for overland flow. The above equation can be solved numerically for y o  f  x, t , i  f  . The solution for y o is then substituted into equation (17) to find a value for q o . Kinematic wave overland flow equations can be used for channel flow. Simple crosssectional shapes, such as circles, trapezoids, and rectangles are used as representative collectors or stream channels. Input data requirements include slope, length, cross-sectional dimension, shape, and Manning’s n value. The basic form of the equation is similar to equation (10) and (12). Ac Qc   qo t x (20) Qc   c Ac mc (21) 2.2.2. Hydrologic Flow Routing Sturm (2001) reviewed storage routing and considered it to be the simplest form of flow routing method. In storage outflow models, the storage in the watershed is approximated by a reservoir 27 whose storage is considered to be a function of the inflow and outflow. The continuity equation in terms of hydrologic routing can be written as: dS =I–O dt (22) In the above equation S = storage in the reach (control volume); I = inflow rate to the reach; and O = outflow rate from the reach. An additional equation is required to solve for the outflow in equation (22) and it is provided by a known functional relationship between storage and the inflow; S = f (I, O). If equation (22) is written over a finite time interval of t , it becomes; I1  I 2 O  O2 t  1 t  S 2  S1 2 2 (23) In equation (23), the subscripts 1 and 2 refer to conditions at the beginning and end of a time interval, respectively. I represent the inflow, O the outflow and S the storage in a channel reach. The time interval is represented by t . The storage in a channel reach depends on the channel geometry and depth of flow. Normally channel routing requires that the channel be divided into several reaches. The outflow from one reach becomes the inflow to the next reach downstream. The time interval t should not exceed one-fifth to one-third of the time to peak of the hydrograph being routed. The routing interval should not exceed the travel time through the reach otherwise it could lead to significant errors in the output hydrograph. 28 Hann et.al, (1994) provide a synopsis of the Muskingum flow routing method, which is based on the modification of the storage routing and considers a linear change in depth along the reach. This routing technique is based on the idea that the flow depth is not constant along the reach because the inflow rate would not be the same as the outflow rate. If the flow were in the rising stage, the inflow would exceed the outflow. Thus, the depth of flow at the upstream end of the reach would exceed the depth at the downstream end. To overcome this non-uniformity, the Muskingum method makes the storage in the reach a linear function of both the inflow and the outflow rate as; S  k xI  1  x O The coefficient k (24) is known as the storage coefficient and is approximately equal to the travel time through the reach. A value of zero for the coefficient x corresponds to reservoir storage routing; and a value of x = 0.5 makes the storage a function of the average flow rate in the reach based on the inflow and outflow. Through manipulation of equations (23) and (24) a linear expression for outflow can be obtained as: O2  C o I 2  C1 I1  C 2 O1 (25) Where: Co   kx  0.5t k (1  x)  0.5t (26) 29 C1  kx  0.5t k (1  x)  0.5t C2  k (1  x)  0.5t k (1  x)  0.5t Curves relating C o , C1 and (27) (28) C 2 to the outflow can be constructed. Singh and Scarlatos (1988) developed a mathematical model for border irrigation using a spatially lumped continuity equation and the Muskingum type storage relation. The storage relation was determined from physical border characteristics. They compared the Muskingum model with the models of Sherman and Singh (1978), and Fork and Bishop (1965) for the “advance phase” and the models developed by Wu (1972) and Sherman and Singh (1982) for the “recession phase”. In the advance phase, the momentum equation is replaced with the assumption of constant surface water depth and an integral equation is derived for advance time in terms of inflow per unit length, mean surface water depth, and cumulative infiltration. Wu (1972) assumed the recession outflow to be linearly related to surface water storage. He derived recession time in terms of border length, slope, roughness, and infiltration constants. The advance phase predicted by the Sherman and Singh (1978) was in agreement with the proposed model. Similarly, the Wu (1972) model and the proposed model predicted the recession phase accurately. They concluded that the Muskingum model performed satisfactorily and was comparable to the other models. 30 Perumal et al. (2001) applied the Muskingum method variable parameter approach to a flood routing scenario of rivers in Australia and the United Kingdom. This study illustrates how to estimate the routing parameters at every routing time interval using limited channel cross section data and the wave speed-discharge relationship developed for the routing reach; which was derived from past observed flood hydrographs or the rating curves available at the inlet and outlet of the study reach. They illustrated that through this approach no information on channel roughness and no calibration were required to estimate the parameters. The ability of the method to reproduce the observed flood hydrographs was evaluated using the Nash-Sutcliffe criterion. Das (2004) did analyses on parameter estimation for the Muskingum model. He was able to minimize the outflow prediction errors subject to the satisfaction of the stream flow-routing equations for all time stages in the routing process. The routing equations incorporated the Muskingum channel storage models. In this research, he developed an algorithm for parameter estimation that iteratively solved the governing equations to identify the Muskingum model parameters. McCuen (1982) presented a channel routing procedure similar to the Muskingum method known as the “Convex routing”. This method only involves inflow-outflow hydrograph. In the Convex method the routing equation is; O2  (1  C )O1  CI1 (29) The parameter C can be estimated using the following equation, 31 C v (1.7  v) Where, (30) v is the average flow velocity of the channel reach. The proper routing interval to use with the convex method is; t  CK Where K (31) is similar to the parameter travel time through the reach. The method as C k of the Muskingum method and may be estimated as the value may also be approximated for the x in Muskingum C  2 x if approximate x is available. In 1983, the SCS replaced the Convex method in their practice standards with the Att-Kin method (Hann et.al, 1994). The Att-kin method combines elements from the Kinematic method with an attenuation procedure based on storage routing. It is solved using the equation below; 2t    2 t  O2   I 1  1  O1   2 K  t   2 K  t  The value of K K (32) is computed by the equation below; L mV (33) 32 In the above equation, V  And L is the reach length and V is the velocity, and V is solved as; q A (34) A is related to q by the rating curve equation; q  xA m (35) Where: A = cross-sectional flow area (sq. meters). q = discharge per unit length (cubic meters per second/meter) If the discharge is derive using Manning’s equation, then 1.49 2 / 3 1 / 2 1.49 1 / 2  A  q R AS  S A  n n P The parameter m in Equation 35 is equal to 5/3, and 33 2/3  1.49 S 1 / 2 nP 2 / 3 A 5 / 3 (36) x 1.49 S 1 / 2 nP (37) 2/3 Which means that m is a function of the velocity-versus-area relationship, and x is a function of the cross-section geometry. 34 Chapter 3: METHODOLOGY This chapter describes 1) research synopsis, 2) development of flow models, 3) watershed characteristics data and development of model input parameters, and 4) development of a grid based curve number estimator for runoff volume computations using ArcMap. 3.1. Synopsis Even though the basic differential equations capable of describing one-dimensional unsteady flows were originally developed over a century ago, it has only been applied during the last thirty years to solve hydrologic engineering problems. This is because of advancements in computer and GIS technologies, as previously it was not possible to efficiently solve these equations without a high-speed computing power. This development has led to the integration of hydrologic cycle components into simulation of comprehensive watershed models. During the last decade, hydrologic research related to flow routing was primarily focused on improving existing methods/models. Very few studies have provided in-depth investigations regarding the choice of flow routing formulations or its relative impacts on watershed model results. For example, researchers like Tseng et al. (2001), Guo (2004), Kim et al. (2009), WenCheng et al. (2009), and Haltas et al. (2009) has mainly contributed to the development and improvement of existing flow routing methods in order to attain computational efficiency. One of the few recent studies that compared the efficiency of different flow routing methods (hydraulic and hydrologic) was completed by Xiong and Melching (2005). In their research, they used previously unpublished experimental data from a watershed experimentation station at the University of Illinois, Urbana Champaign, to test the accuracy of the kinematic wave (simplified full dynamic equation) and non-linear reservoir routing (hydrologic) in an urban 35 watershed. Even though their results fit with the theoretical basis of the kinematic wave theory (hydraulic flow routing), however is not comprehensive. The data used in Xiong and Melching (2005) research was derived from WES, a system “designed to study only the surface runoff with an impervious surface” (Chowand Yen, 1974). Some of the limitations/research gaps in their work follow.  Duration time of the rainfall is short compared to rainstorms in nature (60 sec and 120 sec).  Experimental watershed area is no more than 12.2 meters by 12.2 meters, which is merely as big as a rooftop or driveway.  Roughness coefficient is much lower than that of natural channels or natural/constructed overland surfaces since the surface of the basin used in the research was the smooth side of aluminum plates. This experimental apparatus can only simulate an urban area that has 100% imperviousness, and without any initial abstraction.  Due to the simplicity of the watershed basin, it can only simulate the watershed with a single subarea and a single collecting channel, which is rare in nature.  Because of the nature of the experimental equipment and experimentation, in general, the exact rainfall intensity was unknown. Intensity was computed as peak discharge divided by area. The estimated volume typically was greater than the measured volume, but the full hydrograph was never measured and, thus, the measured volume is less than the true experimental volume. Since the actual rainfall intensities were not known, the average values of the nominal intensities were utilized, which resulted in disagreement as large as 8– 12% between the intensity for the observed and simulated hydrographs (determined on the basis of the variation in the intensities of experiments with plateau discharges). The average value was either higher or lower than the actual intensity of a single event, hence, the peak 36 discharge and/or volume of the simulated hydrograph could be either higher or lower than that of the observed hydrograph. This research proposes to address these gaps/limitations by performing a comparative analysis of flow routing methods in a natural watershed setting. In this research, the physical data such as precipitation, runoff, soils, topographic, and meteorological data is derived from an actual natural system in the WGEW. Below is a summary of the physical data and watershed characteristics of WGEW, which addresses the limitations identified in the research work of Xiong and Melching (2005).  Duration time of the rainfall is based on actual rainstorms in nature.  Experimental watershed area is based on actual watershed size in a natural setting with a 2 2 size range of 395 m to 7.2 km .  Roughness coefficients are based on published data of natural land-cover and use.  Five GIS data layers provide elevations, geology, geomorphology, soils, and land-use type at Walnut Gulch. In my research, the entire watershed is simulated based on the natural conditions on the ground and the initial abstraction is computed using the SCS method.  Watershed is divided in several basins, sub-basins, and streams of varying geometric shapes. The total surface runoff at the watershed outlet is composed of an overland flow of the rainfall excess, and channel flow/stream runoff, which is produced as a result of the overland flow.  Precipitation record is observed via digital gauges and consists of rainfall depths at 1- min intervals during periods of rainfall. Each data logger clock time is checked daily via telemetry and periodically reset to National Institute of Standards and Technology (NIST) 37 standard time. The rainfall intensities in my research are based on actual site-specific data derived from precipitation depths and durations. 3.2. Flow Model Development Modeling can be mathematical, if described by mathematical equations or physical, if a scale model is built to represent dimensional similitude to the actual watershed. In either case the model is a conceptualization of the actual watershed. Mathematical modeling generally seeks to define the mathematical relationship between a set of independent variables and a response or dependent variables. In this research study, we are interested in developing a mathematical model for the actual watershed using the kinematic wave theory and the SCS method. Below is a summary of kinematic wave and SCS flow model development and its application to the study watersheds at the Walnut Gulch Experimental station. 3.2.1. Development of Kinematic Wave Flow Model The kinematic wave theory was originally developed by Lighthill and Whitham (1955) and they have given a full account of the theory for describing flood movement in long rivers. This theory is based on the simplification of the full dynamic equations, which is achieved by combining continuity and momentum equations with inertia and pressure terms dropped. In this research, the kinematic wave flow routing method was conceptualized as shown in Figure (3). This represents the watershed as two plane surfaces over which water runs until it reaches the channel. The water then flows down the channel to the outlet. At a cross section, the system would resemble an open book, with the water running parallel to the text on the page (down the shaded planes) and then into the channel that follows the book’s center binding. The 38 kinematic wave overland flow model represents behavior of overland flow on the plane surfaces. The model was also be used to simulate behavior of flow in the watershed channels. Figure 1: Simplified watershed with kinematic wave representation. The kinematic wave overland flow model is based on the fundamental equations of open channel flow: the momentum equation and the continuity equation. Flow over the plane surfaces is primarily one-dimensional flow. In one dimension, the momentum equation is: S f  So  y V V 1 V   x g x g t (38) 39 The above equation and terms are described in detail in Chow (1959), Chaudhry (1993), and many other texts. The energy gradient can be estimated with Manning's equation given below, which can be written as: Q 2 1 CR 3 Sf 2 N A (39) Where: Q = flow R = hydraulic radius A = cross-sectional area, and N = a resistance factor that depends on the cover of the planes (note that this is not Manning’s n). An important consideration in using the above equation is to use large Manning’s roughness coefficient values as compared to the one used in channels (Bedient and Huber, 1992). For shallow flow, bottom slope and the energy gradient are approximately equal and acceleration effects are negligible, so the momentum equation simplifies to: S f  So (40) The above equation can be simplified to: 40 Q  A m (41) Where  and m are parameters related to flow geometry and surface roughness. The second critical equation, the one-dimensional representation of the continuity equation, is: qA y y V  VB  B x x t (42) Where: q = lateral inflow per unit length of channel B = water surface width A V = prism storage x VB B y = wedge storage x y = rate of rise t The lateral inflow represents the precipitation excess, computed as the difference in precipitation losses. With simplification appropriate for shallow flow over a plane, the continuity equation reduces to equation (43): 41 A t  Q q x (43) Combining the above two equations yields: q A A  a  mA m 1 t x (44) The above equation is a kinematic wave approximation of the equation of motion. HECHMS represents the overland flow element as a wide rectangular channel of unit width with;  =1.486S1/2/N and m=5/3. N is an overland flow roughness factor (Table 1). Kinematic wave overland flow equations can be used for channel flow. Simple crosssectional shapes, such as circles, trapezoids, and rectangles are used as representative collectors or stream channels. Input data requirements include slope, length, cross-sectional dimension, shape, and Manning’s n value, which were derived from GIS data of the study watersheds. Figure 4, shows values of alpha and m for various channel shapes used in HEC-HMS. The kinematic-wave approximation was solved in the same manner for both the overland and channel flow, as described below.  Partial differential equations were approximated with a finite-difference scheme.  Initial and boundary conditions were assigned.  Resulting algebraic equations were solved to find unknown hydrograph ordinates. In this research kinematic wave equations were solved using a finite difference scheme in HEC-HMS. A finite difference method presents “a point wise” approximation to the governing 42 partial differential equations for an array of stationary grid points located in space and time at which the discharge and water surface elevations are computed. Computations advance along the downstream direction for each time step until all the flows and stages are calculated along the entire distance or the reach. Then the computation is advanced ahead in time by one delta and the computation for discharge and water surface elevations are performed once again. Figure 2: Kinematic wave parameters for various channel shapes. (USACE, 1998). 3.2.2. Development of SCS Flow Model The “SCS runoff curve number” method represents the combined hydrologic effect of soil, land use, agricultural land management practices, hydrologic and the antecedent soil moisture 43 conditions (McCuen, 1982). In the SCS method the volume of runoff (Q) depends on the volume of the precipitation (P) and the volume of storage that is available for retention (McCuen, 1982). The actual retention (F) is the difference between the volumes of precipitation and runoff. Furthermore, a certain volume of the precipitation at the beginning of the storm, which is called the initial abstraction, will not appear. The SCS assumed the following rainfallrunoff relation. Q F  S P  Ia (45) S = the potential maximum retention. The actual retention, when the initial abstraction is considered, is: F  P  I a   Q (46) Substituting equation (45) into (46) yields the following: P  I a   Q S  Q P  Ia (47) Solving equation (47) for Q yields: P  I a 2 Q P  I a   S (48) 44 The initial abstraction consists mainly of interception, infiltration during early parts of the storm, and surface depression storage. It can be determined from observed rainfall-runoff events for small watersheds, where lag is minimal, as the rainfall that occurs before runoff begins. Interception and surface depression storage may be estimated from cover and surface conditions, but infiltration during the early part of the storm is highly variable and dependent on such factors as rainfall intensity, soil crusting, and soil moisture. Establishing a relationship for estimating Ia is not easy. Thus, Ia was assumed to be a function of the maximum potential retention, S. An empirical relationship between Ia and S was expressed as: I a  0 .2 S (49) Factors affecting Ia would also affect S, substituting equation (49) in (48): Q P  0.2S 2 (50) P  0.8S S can be estimated as: S 1000  10 CN (51) 45 in which CN = runoff curve number, and it is a function of land use, and antecedent soil moisture and other factors that affect runoff and retention. The curve numbers were estimated at 10-meter resolution from the SSURGO soils and land use datasets using an automated procedure in ArcMap (see next section for details). The initial abstraction (volume of the precipitation at the beginning of the storm that does not appear) was based on curve numbers as; S 1000  10 ; in which CN = runoff curve number, and CN it is a function of soil type, land use data, and antecedent soil moisture condition that affect runoff and retention. The estimated curve numbers and initial abstraction values were then used as an input into the HEC-HMS model to compute the rainfall excess or volume of the inflow hydrographs The overland flow for the hydrologic routing component was determined using the SCS Unit Hydrograph method. In the SCS Unit Hydrograph model the basin outflow results from one unit of direct runoff generated uniformly over the drainage area at a uniform rainfall rate during a specified period of rainfall duration. The underlying concept of the unit hydrograph is that the runoff process is linear, so the runoff from greater or less than one unit is simply a multiple of the unit runoff hydrograph. There are five important concepts in this definition. First, the runoff occurs from precipitation excess, which can be defined as the difference between precipitation and losses, which includes interception, depression storage, and infiltrated water that does not appear as direct runoff. Second, the volume of runoff is one cm, which is the same as the volume of precipitation excess. Third, the excess is applied at a constant rate (uniform rate). 46 Fourth, the excess is applied with a uniform spatial distribution. Fifth, the intensity of the rainfall excess is constant over a specified period of time, which is called duration. The unit time or unit hydrograph duration is the duration for occurrence of precipitation excess. The optimum unit time is less than 20 percent of the time interval between the beginning of runoff from a short duration, high-intensity storm and the peak discharge of the corresponding runoff. The storm duration is the actual duration of the precipitation excess which duration varies with actual storms. There are several types of unit hydrographs. In this research, the SCS dimensionless hydrograph was developed for the study watersheds. The SCS method use dimensional unit hydrographs that are based on an extensive analysis of measured data. Unit hydrographs were evaluated for a large number of actual watersheds and then made dimensionless. An average of these dimensionless unit hydrographs (UH) was developed. The time base of the dimensionless UH was approximately 5 times the time-to-peak and approximately 3/8 of the total volume occurred before the time to peak; the inflection point on the recession limb occurs at approximately 1.7 times the time to peak, and UH had a curvilinear shape. The area under the unit hydrograph equals the volume of the direct runoff Q that was estimated in equation (48):  1 Q  q p t p  tr 2  (52) Solving equation (52) for the peak discharge yields: 47 Q qp  tp   2   1 tr t p     (53) Letting K replace the contents within the brackets yields: qp  KQ tp (54) In order to have the peak discharge in cubic feet per seconds, and the recession time in hours, and Q in inches, it is necessary to divide by the area in square miles and to multiply by the constant 645.3; the recession time is 1.67 times the time-to-peak, equation (52) becomes; qp  484 AQ tp (55) The time of peak (also known as the time rise) is related to the duration of unit excess precipitation as: tp  t  t lag 2 (56) Where: t = is the excess precipitation duration (which is the computational interval in HEC-HMS) t lag = time difference between the center of mass of rainfall excess and the peak of the UH. 48 With qp and tp known UH can be found from the dimensionless form which is included in the HEC-HMS, by multiplication (HEC-HMS uses a computational interval ( t ), less than 29% of t lag (USACE, 1998). The final step in the SCS overland flow model is the transformation of the inflow hydrograph into channel/stream flow using the SCS lag method. The lag method relates the time lag (L), which is defined as the time in hours from the center of mass of rainfall excess to the peak discharge, to the slope (Y) in percent, to the hydraulic length (HL) in feet and the maximum retention (S) (Hydrologic Modeling System, Technical Reference Manual,2000). The data for the SCS lag equation was derived from the 10-m DEMs and the computed curve numbers using ArcMap. The computed lag values were then used to route the overland flow as channel/stream flow. L H l 0.8 S  10.7 1900Y (57) 0.5 In equation (57), S was estimated as; S 1000  10 and CN was estimated as CN described in section (C) of this chapter. The relationship between the time of concentration and lag time was determined using the below equation, which has been derived empirically in the SCS method. tc  5 Lag 3 (58) 49 Figure 3: Conceptual watershed illustrating travel time from the centroid (grey dot) of each band of area to the watershed outlet, (adapted from NRCS, NEH Manual, 2004). For interpretation of references to color in this and all other figures, the reader is referred to the electronic version of this dissertation. 50 3.3. Watershed Characteristics and Development of Model Input Data A total of three separate basins from within the Walnut Gulch watershed were selected and delineated using 10-meters DEM to test the proposed research hypothesis (Figure 2). The map is projected in UTM (Universal Traverse Mercator) zone 11. Figure 9 shows watershed No. 11, which is one of the study watersheds that consist of five sub-basins, while the other two, Watershed No. 121 (Figure 10) and 125 (Figure 11) are composed of a single basin because of total smaller areas. The coordinates generated through the delineation process in ArcMap were used as an input into HEC-HMS to create a background map for the study watersheds. Watershed boundaries and stream networks were both defined in an ASC II file format. Each section begins with a keyword “MapGeo” followed by a colon and either “Boundary Map” or “River Map”. All the coordinates data used for generating the boundary maps and stream network maps for the study watersheds have been provided in the Appendices section of this report. 3.3.1. Drainage Area Calculations The drainage area is an important characteristic for hydrologic analysis. It reflects the volume of water generated from rainfall and is used to indicate the potential for rainfall to provide a volume of water. It is common in hydrologic design to assume a constant depth of rainfall occurring uniformly over the watershed. Under this assumption, the volume of water available for runoff would be the product of rainfall depth, and drainage area. Thus, the drainage area is required as an input to the models ranging from simple linear prediction equations to complex computer models. Table 2 lists the drainage areas for all the study watersheds; these areas were calculated using ArcMap from the delineated watershed boundary maps. Each study watershed was 51 projected in the form of a polygon. Through the use of “list” command, areas in square meters for individual sub-basins were identified, and these values obtained were used as input into HEC-HMS software. Table 1: Study Drainage areas. STUDY WATERSHEDS DRAINAGE AREA (square-meters) Watershed No. 11 8,235,353 Watershed No. 121 54, 228 Watershed No. 125 59,084 3.3.2. Watershed Lengths The length of a watershed is important in hydrologic computation. Watershed length is usually defined as the distance measured along the main channel from the watershed outlet to the basin divide. Thus the length is measured along the principal flow path. Since this length is used for hydrologic calculations, it is commonly known as the hydrologic length. The drainage area and length are both measures of watershed size; they reflect different aspects of size. The length is used in computing a time parameter, which is a measure of the travel time of water through a watershed. Hydrologic lengths for each study watershed were computed using ArcMap (Table 3). These measured lengths were used in the computation of time of concentrations and watershed lag times. 52 Table 2: Study watersheds lengths. STUDY WATERSHEDS WATERSHED LENGTHS (meters) Watershed No. 11 5,383 Watershed No. 121 395 Watershed No. 125 358 3.3.3. Watershed Slopes Watershed slope reflects the rate of change of elevation with respect to distance along the principal flow path. Typically, the principal flow path is delineated, and the watershed slope is calculated as the difference in elevation between the endpoints of the principal flow path divided by the hydrologic length of the flow path. The elevation difference may not necessarily be the maximum elevation difference within the watershed since the point of highest elevation may occur along a side boundary of the watershed rather than at the end of the principal flow path. The slopes of the study watersheds were computed from the 10-meters DEM data using an automated procedure in ArcMap (Table 4). Table 3: Study watersheds slopes. STUDY WATERSHEDS WATERSHED SLOPES (percent) Watershed No. 11 2 Watershed No. 121 4 Watershed No. 125 3 53 3.3.4. Soils Survey and Topographic Data (Digital Elevations) The first soil survey of the WGEW was conducted by the NRCS in the late 1960s (Gelderman, 1970) and contained pedon descriptions and locations of 21 soil map units. Physical and chemical properties of the soil series of the map units became available in 1974 (USDA Soil Conservation Service, 1974). Currently three GIS soil surveys are available for the WGEW, 1) STATSGO, consisting of three soils map units, 2) SSURGO, consisting of 18 soil map units, and 3) a more detailed survey (Breckenfeld, 1994) that is based on the SSURGO data and consists of 25 soil map units. In this research the watershed model was developed based on SSURGO data consisting of 25-soils units in a shape file format (Figure 12). Five GIS data layers provide the geology, geomorphology, soils, potential, and actual vegetation on Walnut Gulch. Standard USGS 10 and 30 m digital elevation model (DEM) data sets cover the study area. In addition, a special mapping effort was undertaken with aerial photography (1:12,000 average photo scales) and corresponding ground control surveys in 1988. This effort resulted in orthorectified 1:5000 map sheets with 5 m contour intervals that, in conjunction with a high-resolution stream map, formed the basis for the creation of a 10 m DEM and as the base maps for subsequent GIS data layer development. These maps meet or exceed national map accuracy standards. In this research topographic data was derived from the standard USGS 10-m digital elevation model (Figure 2). The soils survey and 10-m DEM data along with the land use data was used to compute the SCS curve numbers using an automated procedure in ArcMap (see Appendix A for details). 3.3.5. Precipitation Data The precipitation data was derived from real-time digital precipitation gauges at the study watersheds (Figures, 9, 10, and 11). The precipitation record observed via the digital gauges 54 consists of rainfall depths at 1- min intervals during periods of rainfall. Each data logger clock time is checked daily via telemetry and periodically reset to National Institute of Standards and Technology (NIST) standard time. The time kept by the base station computer is manually set to NIST standard time. Though the base station computer may deviate from NIST time by about ±2 min per month, the network of 88 data logger clocks stay within less than one minute of each other and standard time (Keefer et al., 2008). Inaccuracy in spatially distributed precipitation fields can contribute significantly to the uncertainty of hydrological states and fluxes estimated from land surface models. Garcia, et. al., examined the results of selected interpolation methods for both convective and mixed/stratiform events that occurred during the North American monsoon season over a dense gauge network at the Walnut Gulch Experimental Watershed. Spatial interpolation was performed using both inverse-distance-weighted-squared (IDW) and multiquadric-biharmonic (MQB) methods. Their conclusion was that the order of IDW method is important to the results and under some conditions be just as accurate as the MQB method. In this research the IDW method using HEC-HMS (Hydrologic Engineering CenterHydrologic Modeling System) was used to define the spatial and temporal extents and distribution of the precipitation both for kinematic wave and SCS watershed models. In HECHMS the IDW method relies on the notion of “nodes” that are positioned within a watershed such that they provide adequate spatial resolution of precipitation in the watershed (Figure 18). Watershed No. 11 which because of large surface area (8.2-km2) compared to the other two study watersheds (Watershed No. 121 and 125), was divided into five sub-basins each with three nodes distributed along the centroidal flow path. The precipitation hyetograph was computed for each node using gages near that node. To select these gauges, HEC-HMS constructs 55 hypothetical north-south and east-west axes through each node and finds the nearest gage in each quadrant defined by the axes. Weights were then computed and assigned to these gauges in inverse proportion to the square of the distance from the node to the gage. 3.3.6. Stream Gauge/Observed Flow Hydrographs The WGEW runoff database has the longest period of record of runoff in the world for a semiarid location. The runoff data have been the basis for semi-arid region flood frequency analysis and, in conjunction with rainfall data from the intensive recording rain gauge network, are the basis for understanding rainfall-runoff processes at a range of scales and watershed-scale model development, testing, and validation. The observed flow data was derived from Flume No. 11, 121, and 125 are located at the outlet of each respective watershed (Figures 9, 10, 11). Since 1999, all of these flumes have been operating as digital recorders consisting of potentiometers attached to the stilling well gear mechanism and a Campbell Scientific CR-10 data logger. Prior to 1999, data for these flumes were collected using analog data recorders. For this research, eleven major runoff events were selected for Watershed No. 11 and three events each for Watershed No. 121 and Watershed No. 125 to test the efficiency and performance of kinematic wave and SCS flow routing methods. These runoff events occurred as a result of the precipitation events that happened between 1999 and 2009, which have been discussed in the previous section. These observed storm events produced peak flows between 1and-50 m3/sec, with flow volume between 4000- and 90,046-m3 (Figures 19-21). The runoff data was used to validate the hypothesis in this research. 56 3.4. Development of a Grid Based Curve Number Estimator for Runoff Volume Computations The curve number is a hydrologic parameter that is used to determine the amount of precipitation excess that results from a rainfall event over the basin. It is a function of land use, soil type, and soil moisture (Table 4). Therefore estimation of a curve number requires mapping of the soil and land use within the drainage basin boundaries, and specification of unique soil types and land use category. The manual calculation of curve numbers for large areas or many drainage basins can be cumbersome and time-consuming; therefore an automated procedure using ArcMap GIS was used to develop curve number estimates for the study watersheds. Curve number generation in ArcMap requires three shapefiles: (1) the watershed or drainage basin boundaries for which curve number(s) were calculated, (2) the land use map, and (3) the soil type map. The curve number generator also requires two user-defined look-up tables: (1) the soil group table that provides the conversion from soil types to hydrologic soil groups, and (2) the Curve Number table that defines the land use-soil group categories and curve numbers, similar to Table 4 shown previously. Once the appropriate shapefiles, look-up tables and fields were specified, an error check is performed. Various fields were verified, including the list of drainage basin names in the drainage basin shapefile. Next, the land use and soil type records were compared with the records listed in the look-up tables. Shapefile records that could not be found in the look-up tables were flagged and corrections were made accordingly. After all data specifications and checks were performed, curve number calculations were initiated. These curve numbers were generated based on drainage basin boundaries. For watershed with multiple sub-basins, such as Watershed No. 11, curve numbers were calculated 57 for each individual sub-basin. The program proceeds by clipping the soils and land use shapefiles with the drainage basin boundaries. Soil types were converted to hydrologic soil groups by joining the soil group look-up table to the clipped soil shapefile. The curve number method classifies soils into four hydrologic soil groups (A, B, C, and D), which indicates the amount of infiltration the soil will allow (Table 5). Since, the study watersheds are mainly composed of sandy and gravelly soils, hydrologic soil group (A) dominated. The shapefile were then joined twice, first to the drainage basin shapefile, second to the land use shapefile. This creates a number of smaller polygons inside the drainage boundaries. The curve number look-up table was joined to this compiled shapefile, and a curve number was assigned to each polygon based on the combination of its soil group and land use records. In this way, all the data necessary to determine an area-weighted curve number were formatted into one shapefile, with each polygon having a record for the drainage basin name, soil group, land use, and curve number. An important characteristic of the soil is the “antecedent soil moisture condition”. The curve number method classifies the soil into three antecedent soil moisture conditions.  Condition I: soils are dry but not to wilting point; satisfactory cultivation has taken place.  Condition II: average conditions.  Condition III: heavy rainfall, or light rainfall and low temperatures have occurred within the last five days; saturated soil. For the study watersheds soil moisture condition II was initially used to estimate curve numbers, later these curve numbers were adjusted, depending on soil condition III, or I, Table 6 provides the adjusted values of curves number based on antecedent soil moisture conditions. 58 Table 4: Definition of hydrologic soils group (adapted from USACE-HEC-HMS, Manual, 2000). Hydrologic Soil Group Soil Group Characteristics A Soils having high infiltration rates, even when thoroughly wetted and consisting chiefly of deep, well to excessively-drained sands or gravels. These soils have a high rate of water transmission. B Soils having moderate infiltration rates when thoroughly wetted and consisting chiefly of moderately deep to deep, moderately fine to moderately coarse textures. These soils have a moderate rate of water transmission. C Soils having slow infiltration rates when thoroughly wetted and consisting chiefly of soils with a layer that impedes downward movement of water, or soils with moderately fine to fine texture. These soils have a slow rate of water transmission. D Soils having very slow infiltration rates when thoroughly wetted and consisting chiefly of clay soils with a high swelling potential, soils with a permanent high water table, soils with a claypan or clay layer at or near the surface, and shallow soils over nearly impervious material. These soils have a very slow rate of water transmission. 59 Table 5: Curve number estimates for condition I, II and III (adapted from USACE-HEC-HMS, Manual, 2000). CN for Condition II Corresponding CN for Condition I II 100 95 90 85 80 75 70 65 100 87 78 70 63 57 51 45 100 99 98 97 94 91 87 83 60 55 50 45 40 35 30 25 20 15 10 5 0 40 35 31 27 23 19 15 12 9 7 4 2 0 79 75 70 65 60 55 50 45 39 33 26 17 0 60 Chapter 4: RESULTS, ANALYSES, AND DISCUSSIONS OF MODEL RESULTS Results and analyses are presented in the following order, 1) Results and analysis of precipitation and runoff data, 2) Results and analyses of model output data, and 3) Discussions of model results. 4.1. Results and Analysis of Precipitation and Runoff Data Spreadsheet analyses of the precipitation record show the occurrence of several storm events between 1999 and 2009. In spite of the large number of storm events that occurred during this time frame, only selected storms resulted in runoff events with significant peaks and flow volumes that could be used to test the research hypothesis. These storm events occurred almost exclusively from convective storms during the summer season (Figures 13, 14, and 15). Similarly, data analyses of rainfall intensity and volume indicates that intensity is a dominant factor in the generation of runoff excess compared to the total volume of rainfall. In other words, given the same total rainfall volume, a high intensity event has a higher probability of producing runoff compared to the low intensity event. This is typical of semi-arid regions and many researchers like Dubief, 1953 note that summer rainfalls as low as 0.6 centimeters will yield runoff if the intensity approaches 2.54 centimeters/hr., whereas no flow may result from larger amounts of less intensive winter rain. All of the study watersheds are typical of many semi-arid regions in that the channels are dry for most of the year. Typically, runoff occurs as a result of thunderstorm rainfall, with the flood peak arriving very quickly after the start of runoff, and the duration of runoff being brief. Analysis of runoff data shows that almost all of the annual runoff and all of the largest events occur between July and September due to high-intensity, short-duration thunderstorms of limited 61 aerial extent (Figures 19, 20, and 21). Runoff occurs infrequently in the early fall as a result of tropical cyclones and in the winter as a result of slow moving frontal systems. Both cover large areas and have rainfall of low intensities and long durations, which rarely produces well defined hydrographs and therefore could not be used to test the proposed research hypothesis. On average, there are approximately nine runoff events per year independent of drainage area. The impacts of infiltration of the flood wave into the dry channel bed (transmission losses) and the location of the rainfall producing runoff on runoff peak and volume are discussed in detail by (Renard et al. 2008). 4.2. Results and Analysis of Model Output Data Several different analyses were conducted to explore the accuracy of both the kinematic wave and SCS flow routing methods compared to the observed flow data. 4.2.1. Analyses of Observed and Computed “Peak-flows” First, the peak flows generated by the kinematic wave and SCS models were compared to the observed peak flow hydrographs for selected runoff events. Hydrograph peak flows play a key role in the design, analyses, and performance of hydraulic structures. For Watershed No. 11, the kinematic wave flow model peak flow results agreed within a ± 5 range with the observed hydrograph for ten runoff events (Figure 22). Only a single runoff event shows a negative 25 percent difference to the observed flow hydrograph. In contrast, percent difference in peak flow results of nine runoff events from the SCS flow model are in a range of between 0 and negative 18. Only two runoff events produced a +5 percent difference in peak flows to the observed flow hydrographs. Similarly, for Watershed No. 121 and 125, kinematic wave flow model results agreed within a ± 1 percent with the observed peak flows for all six runoff events (Figure 23). 62 Results from the SCS flow model for these watersheds showed a negative 1 percent difference in peak flows between the observed and computed hydrographs for five runoff events. Analysis of Figures 22 and 23 indicates model bias for both the kinematic wave and SCS flow models as a consequence of the selected methodology. The straight line on the plot represents equality of calculated and observed peak flows. The percent difference in peak flows between the observed data and the kinematic wave flow model results fall near and almost in equal numbers above and below the line. This indicates that the model is no more likely to overpredict than to under-predict. On the other hand, the majority of the percent difference in peak flows between the observed and the SCS flow model results falls below the equality line, which indicates that the model consistently under predicted peak flows. There are smaller differences in observed and computed peak values between Watersheds No. 121 and 125 compared to Watershed No. 11 because of the significant differences in watershed size, topography, land cover, and rainfall distribution patterns (Figures 9, 10, and 11). Less variability in the physical characteristic of a watershed means smaller errors in input data, which ultimately leads to less significant errors in model results. 4.2.2. Analyses of Observed and Computed “time to peak” values A similar type of analysis was performed to explore the differences in “time to peak” values (minutes) between the computed (kinematic wave and SCS flow models) and observed flow hydrographs. “Time to peak” is the time from the rising limb of the hydrograph to the peak flow and is considered an important variable in water quality analysis of natural streams. It can impact the concentration of water quality samples (suspended sediment concentration/rates) if 63 sampling time is determined based on the “time to peak” of the predicted flow hydrograph (Syed, 2005). Figures 24 and 25 shows differences in “time to peak” values in minutes between kinematic wave and SCS flow model results compared to the observed flow hydrographs. Both of the plots indicates that the “time to peak” differences between the kinematic wave flow model results and the observed flow hydrograph fall closely and almost in equal numbers above and below the line meaning that the model is no more likely to over-predict than to under-predict. In contrast, the majority of the difference in “time to peak” values between the SCS flow model results and observed flow hydrograph falls above the equality line, which indicates that the model has a strong tendency to over predict the time to peak flows. There are smaller differences in observed and computed “time to peak” values between Watersheds No. 121 and 125 compared to Watershed No. 11 because of the significant difference in watershed size (Figure 25). Watershed No. 11 is 5,383 m2, while Watersheds No. 121 and125 are 395 and 358 m2, respectively. Smaller basin size means smaller runoff magnitudes and less variability in physical data such as soils, topography, land cover, and rainfall, which ultimately leads to less significant errors in model results. Model results from large basins have more defined errors partially because of lumping/averaging of the physical data and the larger magnitude of runoff events. 4.2.3. Peak Weighted Root Mean Squared Error Another analysis of the “goodness-of-fit” between the computed flow hydrographs (kinematic wave and SCS flow models) to the observed flow hydrographs was performed using the PeakWeighted Root Mean Square Error (PRMSE) as an objective function (USACE, 1998). An 64 objective function measures the degree of variation between computed and observed data. It is equal to zero if the hydrographs are exactly identical. Although several other methods such as “the sum of absolute errors (Stephenson, 1979)” and “sum of squared residuals (Diskin and Simon, 1977)” could be used to compute the “goodness of fit” indices, PRMSE was selected because it provides an implicit measure of comparison of the magnitudes of the peaks, volumes, and times of peak of the two hydrographs. The peak-weighted root mean square (PRMSE) objective function is a modification of the standard root mean square error (RMSE). It compares all ordinates, squaring differences, and it weights the squared differences. The weight assigned to each ordinate is proportional to the magnitude of the ordinate. Ordinates greater than the mean of the observed hydrograph are assigned a weight greater than 1.0, and those smaller, a weight less than 1.0. The peak observed ordinate is assigned the maximum weight. The sum of the weighted squared differences is divided by the number of computed hydrograph ordinates; thus, yielding the mean squared error. Taking the square root yields the root mean squared error. The function is defined as follows: n 2  Qo t   Q s t  Z  t 1 Qo t   Q A 2Q A n 1 n Q A   Qo n t 1 65 (59) (60) Where: Z = objective function Qo t  = observed flow at time t Qs t  = computed flow at time t Q A = average observed flow. The objective function is evaluated for all times t in the objective function window. The PRMSE results for all the three study watersheds show that the kinematic wave outflow hydrographs have lower objective function values compared to SCS outflow hydrographs (Figures 26 and 27). Since this function is an implicit measure of comparison of the magnitudes of the peaks, volumes, and times of peak of the two hydrographs, it means that the kinematic wave flow routing model is more accurate than the SCS flow routing model because of smaller objective function values. 4.2.4. Analyses of Observed and Computed “Center of Mass” values An analysis of the center of mass of the computed hydrographs (kinematic wave and SCS model) to the observed flow hydrographs was performed. Figures 28 and 29 shows the differences in “center of mass” values between the computed and observed flow hydrographs. For Watershed No. 11 ten out of eleven hydrographs computed with the kinematic wave model produced an error margin between zero and +10, meaning that the model is biased and is likely to overpredict the center of mass of computed hydrographs. In contrast, the center of mass of the SCS hydrographs fall in equal numbers above and below the line but has a bigger spread (margin of error, ±15), which means the SCS model does not match the observed data very well. In other 66 words, random errors in the prediction model are large relative to the magnitude of the observed flows. The center of mass of hydrographs for Watersheds No. 121 and 125 produced less significant errors compared to the observed data (Figure 29). Again, this can be attributed mainly due to the relatively smaller basin sizes of Watershed No. 121 and 125 compared to Watershed No. 11. Large basin size can produce significant variability in physical data due to lumping/averaging and can lead to errors in model results. Similarly, because of the differences in basin sizes, there is a significant difference in peak and volume of flows between Watersheds No. 121and125 compared to Watershed No. 11. 4.2.5. Graphic/visual Comparison of Observed and Computed Hydrographs A graphical/visual comparison of the observed and computed flow hydrographs for all three study watersheds and associated runoff events are provided in Figures 30, 31 and 33. These figures indicates that outflow hydrographs computed with the kinematic wave model show a better match with the observed hydrographs compared to the outflow hydrographs computed with the SCS model. It is evident from the observed and computed flows that the peak flow of the SCS model during most events is lower than the observed peak flow. Similarly, there is a forward shift in the “time to peak” values in the SCS hydrographs compared to the observed “time to peak” hydrographs. 67 4.3. Discussions of Model Results Analyses of the SCS model results shows that it produces less accurate results if the aerial distribution of precipitation is irregular and non-uniform over the watershed or if two storm events occur consecutively within a short duration of each other. In nature, these types of precipitation events/conditions typically lead to the formation of double-peaks or irregular shape hydrographs. This is because the SCS model (unit hydrograph) assumes uniform rainfall distribution and intensity over the catchment area during the duration of rainfall excess. In other words, the SCS unit hydrograph theory assumes that watersheds behave as linear systems and that the duration of direct runoff is always the same for uniform-intensity storms of the same duration, regardless of the intensity. In practice, these types of conditions are rarely satisfied. Hydrologic systems are usually nonlinear due to factors such as storm origin and patterns as well as stream channel hydraulic properties. Therefore, if the peak flow produced by a storm of certain intensity is known, the peak corresponding to another storm (of the same duration) with twice the intensity is not necessarily equal to twice the original peak. If all available information indicates that the aerial distribution is inconsistent between different storms or the shape of the watershed and configuration of the drainage network cause multiple peaks for even simple storms, then the SCS model (Unit Hydrograph) should not be used. The alternative to Unit Hydrograph theory is kinematic wave theory and other distributed hydrologic models. The kinematic wave method can describe spatial and/or temporal rainfall and roughness variations, which the SCS method, by virtue of it being lumped, cannot. Likewise, the kinematic wave model is not universally applicable: Ponce (1991) for example, argues that because of numerical properties of the solution algorithms, the method “…is intended primarily for small watersheds [those less than 1 sq mi (2.5 km2)], particularly in the cases in 68 which it is possible to resolve the physical detail without compromising the deterministic nature of the model.” However, this was not the case observed during this research study. For example, analyses of results for Watershed No. 11, which is approximately 8.2-km2, indicates that outflow hydrographs computed with the kinematic wave model show a better match with the observed flow hydrograph compared to the outflow hydrographs computed with the SCS- Unit Hydrograph model. This is primarily because of the high level of discretization of watershed characteristics (physical data) achieved for the study area with the use of GIS data; prior to 1991, use and availability of GIS technology and high-resolution data was limited. In spite of its versatility, certain considerations should be exercised when using the kinematic wave solution to watershed flow problems because it omits or simplifies some terms in the equations to arrive at a solution. These include (but are not limited to) the following:  Backwater effects.  Occurrence of subcritical and supercritical flow.  Channel slope and hydrograph characteristics. Errors introduced due to averaging of physical data in time of concentration (lag) computations can cause a shift in “time to peak” of runoff hydrographs in the SCS unit hydrograph model. Reasons for the variation in lag time may include amount, duration and intensity of rainfall, vegetative growth stage, and available temporary storage. On the contrary, the kinematic wave solution is a distributed parameter and hydraulic data-intensive method (requiring geometric and frictional parameters), in which the flow velocity becomes a function of channel geometry and Manning’s “n” (roughness coefficient). Therefore, the percent difference in peak flows between the observed data and the computed flow results indicates that the 69 kinematic wave model is no more likely to over-predict than to under-predict. On the other hand, the majority of the differences in peak flows between the observed and the SCS flow model results indicate that this model has a strong tendency to under predict peak flows. Other considerations that needs discussion is the effect of curve numbers on watershed model results. Although much care was exercised in the estimation of curve number by developing an automated procedure in ArcMap (see methodology chapter) to eliminate computation errors, it is impossible to avoid the inherent limitations and assumptions associated with the curve number estimation method. It is important to know that the curve number equation does not contain an expression for time and, therefore, does not account for rainfall duration or intensity. Likewise, infiltration rate will approach zero during a storm of long duration, rather than the assumed constant rate. This could be one of the reasons for underpredicted peak flows in the SCS model results. The estimated curve numbers are part of both the SCS unit hydrograph and lag equations (see methodology chapter, section “C” for details of unit hydrograph and lag equations). 70 Chapter 5: CONCLUSIONS and RECOMMENDATIONS The overall goal of this research was to conduct a comparative analysis of the kinematic wave and SCS flow models at a watershed scale in semi-arid environment. Conclusions derived from this research together with recommendations and future applications of the methods developed follow. 5.1. Conclusions From this research the following conclusions were reached.  PRMSE results for all the three study watersheds show that the kinematic wave outflow hydrographs have lower objective function values compared to SCS outflow hydrographs meaning that the kinematic wave flow model is more accurate than the SCS flow model.  Percent difference in peak flows between the observed data and computed flow results indicates that the kinematic wave model is no more likely to over-predict than to underpredict. On the other hand, the majority of the percent difference in peak flows between the observed and the SCS flow model results indicates that the model has a strong tendency to under predict peak flows. In order for the kinematic wave solution to be useful, the discretization must reflect what is actually occurring in the field. If sufficient field data are not available there is a risk that the amount of lumping introduced may interfere with the deterministic character of the method and its ability to simulate overland flows in a distributed context. On the other hand, the SCS method is a spatially lumped conceptual model of runoff generation that is based exclusively on hydrologic data (i.e. streamflow measurements). 71  Since, the overland flow kinematic wave solution is primarily applicable to small catchments, and the SCS method is applicable to mid-size catchments, it seems that there should be little overlap between the two methods. However, the kinematic wave solution has a significant advantage in that it can describe spatial and/or temporal rainfall and roughness variations, which the SCS method, by virtue of it being lumped, cannot do.  Despite the fact that previous researches such as Ponce (1991) only recommended the kinematic wave model for small size watersheds (<2.5 km2), the results of this study showed that the kinematic wave model outperformed the SCS model for watershed greater than 2.5 km2. For Watershed No. 11, which has an area of 8.2-km2, the kinematic wave model produced more accurate results compared to the SCS model. This is primarily because of the high level of discretization of watershed characteristics (physical data) achieved for the study area with the use of GIS; prior to 1991, use and availability of GIS technology and high-resolution data was limited. 5.2. Recommendations Based on the modeling results obtained and taking into account the capability of both the kinematic wave and SCS methods, the following recommendations are made.  Further research should be conducted to better describe the phenomenon of kinematic shock. In the past numerous studies were conducted to analyze the cause and effects of kinematic shock but the subject continue to mystify researchers and practitioners alike. The shock arises due to the nonlinear feature of kinematic wave, which under the right set of circumstances can result in the kinematic wave steepening to the point where it 72 becomes for all practical purposes a wall of water. In the overland flow situation the wall of water would be a small discontinuity in the water surface profile. The shock is a direct consequence of the nonlinear steepening tendency, which is abetted when the 1) wave is kinematic as opposed to diffusion (or dynamic), 2) there is a low base-topeak flow ratio, 3) there is a hydraulically wide and sufficiently long channel, and 4) there is a high Froude number flow.  The initial abstraction method needs further analysis and research. This method is a conceptual model of hydrologic abstraction of storm rainfall. Its objective is to estimate direct runoff depth from storm rainfall depth, based on a parameter referred to as the "curve number." The method does not take into account the spatial and temporal variability of infiltration and other abstractive losses; rather, it aggregates them into a calculation of the total depth loss for a given storm event and drainage area. The method describes average trends, which are based on soil type, land use/treatment, and surface conditions embodied in the concept of antecedent condition. Its main disadvantages are the absence of clear guidance on how to vary antecedent condition and the fixing of the initial abstraction ratio at 0.2, preempting a regionalization based on geologic and climatic setting.    Both the kinematic wave and SCS models should be tested with a digital precipitation model such as NEXRAD (Next-Generation Radar) dataset, to determine if there is an improvement in model results. NEXRAD is a network of 159 high resolution s-band Doppler weather radars operated by the National Weather Service (NWS). It detects 73 precipitation and atmospheric movement or wind and returns data which when processed can be displayed in a mosaic map which shows the pattern of precipitation and its movement. An accurate spatial and temporal extent of a precipitation model is likely to improve watershed model results.    In this research the comparative analyses of kinematic wave and SCS flow models was limited to a semi-arid environment. It is recommended that both of these models are tested in different climatic and geographic regions to explore their capabilities and limitations under a variety of physical and environmental scenarios. Both of these models could produce more accurate results in an environment with moderate weather and uniform land use / cover patterns. 5.3. Future Applications/Benefits The overall objective of this research study was to improve our understanding of mathematical models (kinematic wave and SCS), which are used as tool to quantify surface runoff and water budgets under the stress of increasing population and climatic variation. Several important conclusions have emerged from this study that can prove useful to a practicing engineer/hydrologist. First, the kinematic-wave analysis should be a satisfactory tool to predict surface runoff in semi-arid watersheds, where transmission losses are significant factor besides initial abstraction in the overall water budget. Second, its accuracy is proven and demonstrated with data encompassing a relatively wide range of field conditions in semi-arid environment. And finally, it is proven that the kinematic wave methodology is simple to program and execute. 74 This should enhance our ability to manage watersheds for reliable water supply, water quality, and ecosystem health by improving our ability to quantify semi-arid water budget components, and developing new model components and decision support systems. 75 APPENDICES 76 APPENDIX A: GIS MODEL INPUT DATA In this research physical data such as land-use/land cover, soils, topography, and geomorphological data were derived from the Walnut Gulch Experimental Watershed (WGEW). The Southwest Watershed Research Center (SWRC) operates the Walnut Gulch Experimental Watershed in southeastern Arizona as an outdoor laboratory for studying semi-arid rangeland hydrologic, ecosystem, climate, and erosion processes. Five GIS data layers provide the geology, geomorphology, soils, potential, and actual vegetation on Walnut Gulch. Standard USGS 10 m and 30 m digital elevation model (DEM) data sets cover the study area. In addition, a special mapping effort was undertaken with aerial photography (1:12,000 average photo scale) and corresponding ground control surveys in 1988. This effort resulted in orthorectified 1:5000 map sheets with 5 m contour intervals that, in conjunction with a high-resolution stream map, formed the basis for the creation of a 10 m DEM and as the base maps for subsequent GIS data layer development. These maps meet or exceed national map accuracy standards. 77 APPENDIX B: RAINFALL AND RUNOFF DATA The precipitation data were derived from real-time digital precipitation gauges at the study watersheds. The precipitation record observed via the digital gauges consists of rainfall depths at 1- min intervals during periods of rainfall. Each data logger clock time is checked daily via telemetry and periodically reset to National Institute of Standards and Technology (NIST) standard time. The observed flow data was derived from Flume No. 11, 121, and 125 are located at the outlet of each respective watershed. Since 1999, all of these flumes have been operating as digital recorders consisting of potentiometers attached to the stilling well gear mechanism and a Campbell Scientific CR-10 data logger. Prior to 1999, data for these flumes were collected using analog data recorders. The WGEW runoff database has the longest period of record of runoff in the world for a semi-arid location. The runoff data have been the basis for semi-arid region flood frequency analysis and, in conjunction with rainfall data from the intensive recording rain gauge network, are the basis for understanding rainfall-runoff processes at a range of scales and watershed-scale model development, testing, and validation. 78 REFERENCES 79 REFERENCES Bedient, P.B., and Huber, W.C., 1992, Hydrology and Floodplain Analysis, Addison-Wesley Publishing Company, USA Bennett J.P., 1995, Algorithm for resistance to flow and transport in sand-bed-channels: Journal of Hydraulic Engineering, v. 121, no. 8, p. 578-590. Bennett J.P., 2001, User’s guide for mixed-sized sediment transport model for networks of onedimensional open channels: U.S. Geological Survey Water-Resources Investigations Report 014054, 33 p. Bicknell, B.R., J.C. Imhoff, J.L. Kittle Jr., T.H. Jobes, and A.S. Donigian, Jr. 2005. Hydrological Simulation Program - FORTRAN (HSPF). User's Manual for Release 12.2. U.S. EPA National Exposure Research Laboratory, Athens, GA, in cooperation with U.S. Geological Survey, Water Resources Division, Reston, VA. Boyer, E. W., Hornberger, G. M., Bencala, K. E. and McKnight, D. M. (1997), Response characteristics of DOC flushing in an alpine catchment. Hydrological Processes, 11: 1635–1647 Carpenter, T. M., and Georgakakos, K. P., (2006), “Intercomparison of lumped versus distributed simulations on operational forecasts scales”, “Journal of Hydrology”, volume 329, Issues 1-2, pages 174-185., available at URL: Chapra, S. C., 1997, Surface Water-Quality Modeling, McGraw-Hill Series in Water Resources and Environmental Engineering, McGraw-Hill, New York, 844 p. Chaudhry, M. H., (1993), “Open channel flow”, Englewood Cliffs, N.J., Prentice Hall, 483 p. Chow, V. T., D. R. Maidment, and L. W. Mays, (1988), “Applied Hydrology”. McGraw-Hill, New York. Cosgrove, W. J., and Rijsberman, F. R., “World Water Vision: Making Water Everybody's Business”, 2000, Published with the World Water Commission. 80 C. T. Haan, B. J. Barfield, and J. C. Hayes., “Design Hydrology and Sedimentology”., Academic Press, A Division of Harcourt Brace and Company, 1994. DeVantier, B. A., Feldman, A. D., 1993, “Review of GIS Applications in Hydrologic Modeling.” Journal of Water Resources and Planning. 119(2): 246-261. Fairfield, J., and P. Leymarie (1991), Drainage networks from grid digital elevation models, Water Resour. Res., 27(5), 709–717. Francisco Olivera., “Extracting Hydrologic Information from Spatial Data for HMS Modeling”. Journal of Hydrologic Engineering, 2001, volume 6, pp.524-530 Garcia, M., Peters-Lidard, C., Goodrich, D., 2008., “Spatial interpolation of precipitation in a dense network for monsoon storm events in the southwestern United States.”, Water Resources Research, Vol. 44, W05S13, doi:10.1029/2006WR005788. Gesch, D., and R. Wilson, 2002. Development of a seamless multisource topographic/ bathymetric elevation model of Tampa Bay, Marine Technology Society Journal, 35(4): 58–64. Gleick, P., (2002-2003), “The World’s Water; The biennial report on freshwater resources”. Gleick, P., (1989), “Climate Change, Hydrology, and Water Resources”, Reviews of Geophysics RVGPB4, Vol. 27, No. 3, p 329-344, August 1989. 5 fig, 3 tab, 94 ref. Goodrich, D. C., Woolhiser, D. A., and Keefer, T. O., 1991., “Kinematic routing using finite elements on a triangular irregular network.”, Water Resour. Res., 27(6), 995–1003. Guo-an, Tang, Hui Yang-he, Josef Strobl, and Liu Wang-qing (2001), The impact of resolution on the accuracy of hydrologic data derived from DEMs, J GEOGR SCI, 11(4), 393 Haltas, I., and Kawas, M. L., 2009, “Modeling the Kinematic Wave Parameters with Regression Methods”., J. Hydrol. Eng. 14, 1049; Haktanir, T., and Ozmen, H. 1997, “Comparison of hydraulic and hydrologic routing on three long reservoirs., J. Hydraul. Eng. 123, 153 81 Hann, C. T., Barfield, B. J., and Hayes, J. C., (1994), “Design Hydrology and Sedimentology for Small Catchments”, Academic Press, A Division of Harcourt Brace and Company, 588 p. Hromadaka, T. V., and DVries, J. J. (1988) “Kinematic wave and computational error.” J Hydr. Engrg., ASCE, 114(2), 207-217. Huggins, L.F. and J.R. Burney. “Surface Runoff, Storage and Routing, Chapter 5.” In Hydrologic Modeling of Small Watersheds”, edited by C.T Hann, ASAE Monograph No. 5 (1982) Hubert Chanson., “The Hydraulics of Open Channel Flow”., Arnold Publishing Company, 1999 edition. “Hydrologic Modeling of Small Watershed Hydrology”., American Society of Agricultural Engineering, An ASCE Monograph Number 5 in a series published by the American Society of Agricultural Engineers, 1982 edition. Jayawardena, A. W., Muttil, N., and Lee, J. H. W., 2006, Comparative analyses of data driven and GIS based conceptual rainfall runoff model., J. Hydrol. Eng. 11, 1061; Jenson, S. K., and J. O. Domingue., 1988. “Extracting Topographic Structure from Digital Elevation Data for Geographical Information System Analysis.”, Photogrammetric Engineering and Remote Sensing, 54(11):1593-1600. Kenward, T., Lettenmaier, D. P., Wood, E. F., and Fielding, E., 2000., “Effect of digital elevation model accuracy on hydrologic predictions.”, Remote Sens. Environ., 74, 432–444. Kibler, D. F., and Woolhiser, D. A. 1970, The kinematic cascade as a hydrologic model, Colorado State University Hydrology No. 39., 27pp, 21fig, 4tab, 18ref, 3append., OWRR ProjectB-030-Colo(2). Kuo, W.-L. et al., 1999., “Effect of grid size on runoff and soil moisture for a variable-source area hydrology model.” Water Resour. Res., 35(11), 3419–3428. 82 Lighthill, M. J.; Whitham, G. B. (1955). "On Kinematic Waves. I. Flood Movement in Long Rivers". Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 229 (1178): 281. Leavesley, G. H., (2005), “Modeling the effects of climate change on water resources – a review”, “Earth and environmental science”, Vol. 28, numbers 1-2, Springer Netherlands. Mathsoft Engineering and Education, Inc., 2001, User’s guide Mathcad 2001i: Cambridge, Mass., Mathsoft Engineering and Education, Inc., 505 p. McCuen, R.H., (1982), “A guide to hydrologic analysis using the SCS Methods:” Englewood Cliffs, N.J., Prentice-Hall, 145 p.Moore, I. D., R. B. Grayson, and A. R. Ladson. 1991., “Digital Terrain Modeling: A Review of Hydrological, Geomorphological and Biological Applications.” Hydrological Processes, 5(1): 3-30. McElmurry, S.P., Aslam, I., Syed, A.U., Voice, T.C., Long, D.T., Phanikumar, M.S., Lusch, D.P., Northcott, W.J., Determining the effects of various land uses within Michigan State University's watershed on water quality in the Red Cedar River, 8th International Conference on Environmental Science and Technology, Lemnos Island, Greece, Proceedings CEST2003: pp. 558-565 (2003) Natural Resources Conservation Service, 2004, Estimation of direct runoff from storm rainfall: Part 630 Hydrology National Engineering Handbook, Chap. 10, 51 p. Neitsch, S.L., Arnold, J.G., Kiniry, J.R., and William, J.R., 2001. Soil and Water Assessment Tool User's Manual, Version 2000. U.S. Department of Agriculture, Agricultural Research Service, Temple, TX. Perumal, M., and Raju, K. G. R., 1998, Variable-Parameter Stage-Hydrograph Routing Method. I: Theory J. Hydrol. Eng. 3, 184, 2(109). Philip B. Bedient., Wayne C. Huber., “Hydrology and Flood Plain Analysis”, Addison-Wesley Publishing Company, second edition. Ray K. Linsley, Jr., Max L. Kohler., Joseph L. H. Paulhus., “Hydrology for Engineers”., McGraw-ill Series in Water Resources and Environmental Engineering. 83 Ponce, V.M., and Yevjevich, V. (1978). “Muskingum-Cunge method with variable parameters.” Journal of Hydraulics Division, ASCE, 104(HY12), 1663-1667. Ponce, V.M., (1991). “The Kinematic Wave Controversy.” Journal of Hydraulic Engineering, 1991, Vol. 117, No. 4, pp 511-525. Rallison, R.E. “Origion and Evolution of the SCS Runoff Equations”, Proceedings of the Symposium on Watershed Management, ASCE 1980. R. Edward Beighley., Glenn E. Moglen., “Trend Assessment in Rainfall-Runoff Behavior in Urbanizing Watersheds”., Journal of Hydrologic Engineering, 2002, volume 7, issue 1, pp 27-34 Smith, R.E. and D.A Woolhiser “Mathematical Simulations of Infilterating Watersheds” Hydrology Papers, N0. 47, Colorado State University, Fort Collins, January 1971. Singh, V. P., Scarlatos, P. D., Collins, J. G., Jourdan, M. R., 1988, Breach erosion of earthfill dams (BEED) model., Journal of natural hazards, Volume 1, Number 2, 161-180 Singh, V. P. (2002). “Mathematical models of watershed hydrology” Journal of Hydrologic Engineering, 270-292, Vol. 7, No. 4. Singh, V. P., and Frevert, D. K., (2002), “Mathematical models of large watershed hydrology”, Water Resources Publications LLC (WRP). Southwest Watershed Research Center and Walnut Gulch Experimental Watershed, (2007), “Brochure Southwest Watershed Research Center, Tucson, Arizona”. Sturm, T.W., (2001), “Open channel hydraulics” New York, McGraw-Hill Series in Water Resources and Environmental Engineering, 493 p. Syed, A.U., Bennett, J.P., and Rachol, C.M., 2005, “A pre-dam-removal assessment of sediment transport for four dams on the Kalamazoo River between Plainwell and Allegan, Michigan”: U.S. Geological Survey Scientific Investigations Report 2004-5178, 41 p. at URL: http://pubs.usgs.gov/sir/2004/5178/pdf/sir2004_5178.pdf 84 Syed, A.U., Fogarty, L. R., 2005, “Trends in Surface-Water Quality at Selected National Stream Quality Accounting Network (NASQAN) Stations, in Michigan”: U.S. Geological Survey Scientific Investigations Report 2005-5158, 38 p. at URL: U.S. Army Corps of Engineers, Hydrologic Engineering Center, (2000), “Hydrologic Modeling System, Technical Reference Manual, HEC-HMS”. U.S. Army Corps of Engineers, Hydrologic Engineering Center, (2002), “River Analyses System, User Reference Manual, HEC-RAS”. U.S. Soil Conservation Service, 1985, Hydrology, sect. 4, Soil conservation national engineering handbook, Washington, D.C., U.S. Department of Agriculture. U.S. Geological Survey (USGS). (1998)., “Standards for digital elevation models. National mapping technical instructions.” Dept. of Interior, USGS, Reston, VA Vieux, E. B., (1988), “Finite Element Analyses of Hydrologic Response Areas using Geographic Information Systems”, A Ph. D., dissertation at Michigan State University. Vivoni, E. R., Ivanov, V. Y., Bras, R. L., and Entekhabi, D., 2004., “Generation of Triangulated Irregular Networks Based on Hydrological Similarity.” Journal of Hydrologic Engineering. 9(4): 288-302. Vieux, E. B., (1988), “Finite Element Analyses of Hydrologic Response Areas using Geographic Information Systems”, A Ph. D., dissertation at Michigan State University. Wurbs, R. A. (1998). “Dissemination of generalized water resources models in the United States.” Water Int., 23, 190-198. Xiong, Y., Melching, C. S., (2005), “Comparison of kinematic-wave and non-linear reservoir routing of urban watershed runoff”, “Journal of hydrologic engineering”, Vol -10, No-1, 39-49 p. 85