CRWR Online Report 99-6 Flood Forecasting for the Buffalo Bayou Using CRWR-PrePro and HEC-HMS by Seth R. Ahrens, M.S.E. Graduate Research Assistant and David R. Maidment, PhD Principal Investigator December 1999 CENTER FOR RESEARCH IN WATER RESOURCES Bureau of Engineering Research ? The University of Texas at Austin J.J. Pickle Research Campus ? Austin, TX 78712-4497 This document is available online via World Wide Web at http://www.crwr.utexas.edu/online.html ii Flood Forecasting for the Buffalo Bayou Using CRWR-PrePro and HEC-HMS by Seth Robert Ahrens, M.S.E. The University of Texas at Austin, 1999 SUPERVISOR: David R. Maidment The Center for Research in Water Resources (CRWR) at the University of Texas at Austin over the past few years has spent significant effort developing CRWR-PrePro. This package, which runs off of a Geographic Information Systems (GIS) platform, provides a link between GIS data and the Hydrologic Engineering Center?s (HEC) Hydrologic Modeling System (HMS). This project focuses on the application of CRWR-PrePro and HMS to generate a hydrologic model of the Buffalo Bayou watershed, which is west of Houston, Texas along Interstate 10. In October, 1994, a severe storm event occurred over the Houston metropolitan area, dropping rainfall totals ranging from 5 to 40+ inches. Since the National Weather Service had recently installed Next Generation Weather Radar (NEXRAD) equipment in Houston prior to this event, high-resolution (1 km 2 ) rainfall data were available. These rainfall data were processed using a new Avenue-based version of the UNIX code (GridParm) the HEC had used previously to create the necessary rainfall data file for HMS. In addition, 30-meter DEM data were acquired via the internet and analyzed by CRWR-PrePro to complete a sub-watershed delineation of the Buffalo Bayou in the area of interest. Using the results from the new GridParm procedure and the output from CRWR-PrePro, a HMS model of the Buffalo Bayou was created. Results from the model indicate that the use of 30-meter DEM data is inadequate when dealing with a place as flat as Houston. Higher resolution terrain data, such as LIDAR, would presumably work much better. Nonetheless, this project revealed the benefits of using CRWR-PrePro to link GIS with HMS. iii ACKNOWLEDGMENTS This research was supported in part by the Hydrologic Engineering Center of the US Army Corps of Engineers under grant DACW 05-97-P-0973. First, for giving me the background and the support to even make it to the graduate level of education, I would like to thank both of parents, Sally and Chris, and my sister, Meg. My family has lent unmitigated support to me my entire life and for that I am eternally grateful. While attending Rice University as an undergraduate, Phil Bedient helped spark my interests in water resources. I would like to thank him for the excellent opportunities he offered me at Rice to help prepare me for graduate school. During my time here at The University of Texas at Austin, I have been absolutely fortunate in my dealings with the faculty. First and foremost, I owe everything to David Maidment, my advisor and rock steady pillar of support. There were days while at the University when I did not know if I would finish or not, yet Dr. Maidment always exhibited confidence in me to get the job done. His willingness to look to the future, no matter what the past may have held, is unparalleled. Next, I would like to thank Edward Holley. The great pains he went to help me get my thesis finished deserve a round of applause. In addition, many thanks go out to Howard Liljestrand. Always with a smile, he was never deterred when it came to helping me get my act together here at the University. Finally, Ty Lehman deserves recognition iv for always being there to offer assistance when I needed it. To all four of you, I thank so very much. Many students at the University also helped to provide me with a unique experience. Fellow students in the department such as Pablo Valdez, Brad Hudgens, Brian Adams, Jerry Perales, Caroline Gerwe, Patrice Melancon, Matt Russell, and Eric Tate all managed to keep my mind sane. I?d also like to thank Seann Reed for the wealth of knowledge always offered when I needed it. Finally, I?d like to thank the U.S. Army Corps of Engineers Galveston Office and the Hydrologic Engineering Center for their contribution to this project. In particular, I?d like to thank James Doan, Tom Evans, and Arlen Feldman at the HEC for their support and guidance during my time working on this project. Thanks also are in order for Mary Lynn Baeck and Jim Smith at the Princeton Environmental Institute and Fred Liscum at the USGS for providing data for this project free of charge. v TABLE OF CONTENTS Section Page ABSTRACT....................................................................................??????. ii ACKNOWLEDGMENTS....................................................................................... iii TABLE OF CONTENTS ....................................................................................... v LIST OF TABLES ................................................................................................ vi LIST OF FIGURES ................................................................................................ viii 1 INTRODUCTION .......................................................................................... 1 2 LITERATURE REVIEW ................................................................................. 6 3 DATA COLLECTION AND PREPROCESSING .......................................... 10 3.1 Preparation of DEM Data for Use in CRWR-PrePro ........................... 13 3.1.1 Photographs in and around Addicks Reservoir ........................... 20 3.1.2 Creation of a Triangular Irregular Network (TIN) ..................... 23 3.2 Preparation of NEXRAD Data ............................................................ 25 3.3 Preparation of Digital Line Graph Data for the Stream Network ......... 32 3.4 Preparation of USGS Flow Data ......................................................... 37 3.5 Acquisition of Anderson Land Use/Land Cover Data ........................ 43 3.6 Acquisition of the STATSGO Data Set ................................................ 45 3.7 Preparation of the Curve Number Grid ................................................ 46 3.8 Acquisition of Digital Orthophoto Quadrangels for Study Area ......... 47 4 TERRAIN PROCESSING USING CRWR-PREPRO .................................... 52 4.1 Descriptions of Functions within CRWR-PrePro ................................. 52 4.2 Use of CRWR-PrePro to Create the HMS Basin Model ..................... 58 5 NEXRAD PROCESSING .............................................................................. 76 5.1 Preparing the Rainfall Grids for the HMS .......................................... 77 5.2 The ModClark Parameter File with the HMS .................................... 78 5.2.1 The SHG Grid .............................................................................. 78 5.2.2 Use of GridParm to Create the ModClark Parameter File ......... 79 5.2.3 Creation of ModClark File in ArcView Using Avenue ......... 81 6 MODEL CALIBRATION WITHIN HEC-HMS ............................................. 86 6.1 Additional Data for the HMS Model ................................................... 87 6.2 Execution of the HMS Model ............................................................ 92 7 CONCLUSIONS AND FUTURE WORK ...................................................... 101 7.1 Modeling Conclusions ........................................................................ 101 7.2 Recommendations for Future Work ................................................... 105 APPENDICES ...................................................................................................... 107 A COMPUTER CODES OF INTEREST ......................................................... 107 B SUPPLEMENTAL DATA .............................................................................. 134 REFERENCES ...................................................................................................... 137 VITA ..................................................................................................................... 139 vi LIST OF TABLES Table Page 3.1 Names and locations of pertinent USGS flow and storage gauges ........... 37 3.2 Partial variable shift table for Buffalo Bayou near Katy, TX .................... 39 3.3 Legend for land use data ............................................................................. 44 3.4 Percentages of land-use types in modeling area ......................................... 44 3.5 Breakdown of hydrologic groups in study area ......................................... 46 4.1 Area comparisons using natural stream network ...................................... 63 4.2 Area comparisons using a combination of natural streams and man-made channels ..................................................................................................... 65 4.3 Area comparisons using stream network derived from TIN ....................... 68 6.1 Summary of data at flow gauges in the HMS model ................................ 100 7.1 Average stream length per unit area ........................................................ 104 vii LIST OF FIGURES Figure Page 1.1 Average daily flow for Buffalo Bayou at Houston flow gauge ............... 3 3.1 Area covered by the GIS database ......................................................... 11 3.2 DEM names of interest on the TNRIS server .......................................... 14 3.3 Gaps existent in the DEM data ............................................................... 16 3.4 Interpolated DEM .................................................................................... 18 3.5 View of berm from inside Addicks Reservoir ....................................... 21 3.6 Blocked culverts in the upstream part of Addicks Reservoir .................. 22 3.7 Exit point of Addicks Reservoir ............................................................ 22 3.8 DEM (left) and TIN (right) of the same area .......................................... 24 3.9 Area covered by NEXRAD rainfall data ................................................ 26 3.10 Inadequacy of 1000 meter cell size in ASCII files ................................. 28 3.11 Upper-right-corner of rainfall grid using cell size of 997 meters ............ 30 3.12 Demonstration of inadequacies in original DEM .................................... 33 3.13 DLG hydrography features of the study area .......................................... 36 3.14 Graphical relationship among elevation vs. time, flow vs. elevation, and flow vs. time (bottom) ..................................................................... 42 3.15 Digital orthophoto of Addicks Reservoir ................................................ 51 4.1 Linking GIS Data and modeling software ............................................. 52 4.2 CRWR-PrePro menu .............................................................................. 53 4.3 A stream network broken into links ......................................................... 56 4.4 Burning a stream network into a DEM ................................................... 60 4.5 Natural stream network used when burning in steams ........................... 62 4.6 Watershed delineation using natural stream network .............................. 63 4.7 Delineated network using a combination of natural streams and man-made channels ................................................................................. 64 4.8 DLG data draped over TIN to show which parts of the DLG data belong in the network ..................................................................... 66 4.9 Delineated system using stream network based on TIN and DLG data 68 4.10 Basin used by CRWR-PrePro to create the HMS basin model ............... 71 4.11 Correlation to determine average watershed velocity ........................... 73 5.1 Rainfall totals (inches) over the Buffalo Bayou .................................... 76 5.2 Visual representation of ModClark code ................................................ 84 6.1 Basin model as imported into HEC-HMS ............................................. 86 6.2 Land use over the modeled region ......................................................... 89 6.3 T c vs. area for the modeling region ...................................................... 91 6.4 Hydrograph of Buffalo Bayou near Katy, TX gauge using original parameters ........................................................................ 92 6.5 Hydrograph of Buffalo Bayou near Katy, TX using T c and R equations ................................................................................................ 95 viii LIST OF FIGURES (continued) Figure Page 6.6 Hydrograph of Buffalo Bayou near Katy, TX using 3*Bedient R values ...................................................................................................... 96 6.7 Hydrograph of the Bear Creek gauge using 3*Bedient R values ............ 98 6.8 Hydrograph of the Langham Creek gauge using 3*Bedient R values ... 99 7.1 Distribution of rainfall depths in inches ................................................ 102 1 Chapter 1: Introduction The greater Houston metropolitan area receives approximately forty-five inches of rain annually. All too often in Houston, sudden rainstorms cause severe flooding, which, as Houston continues to grow, will only get worse. As the risk of significant flooding increases, concerns about flood damage within the city rise. Considering that the U.S. continental record for rainfall in a twenty-four period is forty-three inches in nearby Alvin, TX, these concerns are warranted. If a rainstorm comparable in magnitude struck Houston, it is vital to understand the rainfall/runoff aspects of the terrain. The ability to predict the velocity and magnitude of the flood peaks as they move through the city is extremely valuable. Many of the channel systems in Houston are called bayous instead of rivers. The Buffalo Bayou lies mostly to the west of downtown Houston and drains approximately 360 mi 2 . Within the bayou are two reservoirs?the Addicks and Barker reservoirs. Though the bayou forms in a rural area near the western edge of Harris County, it eventually flows through the heart of Houston on its way to the Houston Ship Channel. Because of this, the bayou is of particular interest since a severe storm event occurring in the western area of the bayou watershed could cause significant damage as the flood wave passes through Houston. In response to these concerns, the Galveston, TX, District of the United States Army Corps of Engineers (the USACE) and the Hydrologic Engineering Center (HEC) in Davis, CA, have teamed up to sponsor the work done to complete this thesis. The overall goal of this 2 study is to develop and calibrate a surface-water model of the Addicks/Barker reservoir system in the Buffalo Bayou upstream of Houston proper using a geographic information system (GIS) in connection with HEC?s surface-water model, the Hydrologic Modeling System (HMS). The process of creating such a model can be broken down into several major steps. Briefly, these include collecting the necessary GIS data, processing them for calculations within the GIS environment, applying a new CRWR-PrePro for connecting GIS systems with other modeling packages, and finally performing the actual modeling using the HMS. The first piece of the study involved the construction of a geographic information system database of the area in and around Houston, which covers about 3420 mi 2 . A GIS provides an environment for displaying and manipulating spatial data. Typical GIS data sets included in the collection are thirty-meter digital elevation maps (DEM), digital line graphs (DLG), Anderson Land Use/Land Coverage data, and STATSGO soils data. These data were organized within an ArcView environment and represent the terrain-based aspects of the model. In addition to terrain data, rainfall data were required for creation of the model. The rainfall data used for this model were NEXRAD-based rainfall data. NEXRAD, from NEXt generation RADar, uses Doppler radar to record both the intensity of rainfall as well as average raindrop size to deduce rainfall rates. For this model, the rainfall data used have a spatial resolution of one km 2 and a time resolution of approximately six minutes and are from a rainfall event which occurred from 16-18 October 1994. These data were generously supplied by the Princeton 3 Environmental Institute. During these few days, a convergence of several weather systems caused a severe, widespread rainfall event in the Houston metropolitan area. While some areas reported only a few inches of rain, typical values of rainfall within Houston were between ten and twenty inches of rain with some small areas reporting values of over forty inches! This is the type of event which flood forecasters in Houston fear the most. As such, this storm serves as a benchmark event in the study of flooding in Houston. Buffalo Bayou at Houston Average Daily Flow 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10/1/94 10/8/94 10/15/94 10/22/94 10/29/94 11/5/94 11/12/94 Date Flow (c fs ) Figure 1.1: Average daily flow for Buffalo Bayou at Houston flow gauge With data assembly complete, the next step of the modeling involved processing the acquired information in a GIS environment. In recent years, the Center for Research in Water Resources (CRWR) at the University of Texas at Austin (UT-Austin) has devoted substantial resources to the development of CRWR-PrePro, 4 which employs the geospatial data collected in phase one to create an input basin model for HMS. While the actual modeling of Buffalo Bayou occurs in HMS, CRWR-PrePro greatly simplifies the creation of this model. Once CRWR-PrePro had been used to generate an HMS basin model, all remaining work occurred within the HMS environment. Most of the time spent working with HMS focused on model calibration. During this flood event, the United States Geological Survey (USGS) maintained five flow gauges within the Buffalo Bayou. These gauges provided the field results to which the model results can be compared. In addition to the flow data, the USGS also collected water level storage data for both of the reservoirs. These data are important for determining flow values downstream of the reservoir system. The remaining sections of this thesis, in addition to providing a pertinent literature review, address the details in carrying out the process previously outlined. This approach involves first collecting data to create a GIS database of not just the Buffalo Bayou but the greater Houston metropolitan area as well. In Chapter 3, both the acquisition and preprocessing of these data sets are addressed. From here, terrain processing in CRWR-PrePro is discussed. The end result of the CRWR-PrePro analysis is a usable basin model for HEC-HMS. Chapter 4 describes how this file was created using CRWR-PrePro. The next section reviews the techniques used to prepare the NEXRAD rainfall data for use in the HMS environment in conjunction with the basin model from CRWR-PrePro. From there, Chapter 6 relates how the HMS model was used to analyze the October, 1994 storm event in the Buffalo 5 Bayou. This chapter includes the results of the modeling process. Finally, conclusions and recommendations for future work are included in the last chapter of this thesis. Overall, this thesis presents an original contribution to knowledge in that it is one of the first large-scale applications of the CRWR-PrePro software package. The benefits of this package, such as a simplified approach to creating a basin model for HMS and a steady link between GIS data and non-GIS modeling software, are presented. In addition, the research done for this project included updating the process created by the HEC to generate the rainfall grid file HMS needs to link gridded rainfall data with one of its basin files. This process, dubbed ModClark by the HEC, while first written using AML files in a UNIX environment, has been rewritten for direct use with CRWR-PrePro. Finally, an end result of this work is a rainfall/runoff model of the Buffalo Bayou for the major rainfall event of October, 1994. As a group, these subtopics work together to present a different take on modeling sizable areas using GIS and HEC-HMS. 6 Chapter 2: Literature Review This chapter summarizes literature reviewed as part of this study. Previous research efforts of interest include the original development of the CRWR-PrePro code, the feasibility of the use of NEXRAD in rainfall/runoff models, past applications of grid-based rainfall data in modeling studies, and previous studies of the Buffalo Bayou region. The research performed for this thesis extends two primary points presented in previous research. First, NEXRAD data have recently become more available over a larger portion of the country in the past few years. HEC (1996) performed an early study of the application of this type of data on the Salt River Basin in Missouri. The HEC, in its documentation of its research, indicates that a primary purpose of the work was to provide an example application of the technology that other field offices of theirs could follow. This report concludes that NEXRAD rainfall data "has the potential to make major improvements to the modeling of spatially varied rainfall events." Furthermore, Reed (1995) concludes that NEXRAD data provide "high quality rainfall estimates." Individual rain gauges do not always have the ability to capture localized convective events, while NEXRAD radar data excel at quantifying the details of events such as these. In addition, using NEXRAD data allows the modeler to break up the terrain surface into a finer resolution, which leads to better model results. The application of NEXRAD data in this study provides an extension 7 of the research presented by HEC in 1996. Further examples of the value of this type of data set will hopefully lead to greater availability and application. HEC itself clearly understands the value of this type of rainfall data. However, to be used effectively, the data must be managed properly. To incorporate NEXRAD data into its HEC-HMS modeling software, HEC has developed the ModClark algorithm for converting NEXRAD data into a format which HMS can handle. When considering grid-based rainfall, HMS requires the modeler to break the terrain into a grid of a pre-specified resolution (Chapter 5) . ModClark, using this resolution, assigns rainfall quantities to each cell in the grid as well as average distance from each location in this cell to the sub-watershed outlet. In addition to discussing this procedure in the aforementioned Salt River Basin study, HEC contributed an entire study on the ModClark procedure in a 1996 study of the Muskingum River Basin in Ohio. While this study employs rain gauge data over NEXRAD data, the conclusions stress the importance of the ModClark procedure in preparing rainfall data for use in HMS by stating that "the ModClark method had significant potential for improving forecasting capability." Since it appears that NEXRAD data are the superior type of rainfall data and that the ModClark procedure is the best algorithm at preparing rainfall data for use in HMS, this study incorporates both of these issues in its scope. To create a more integrating GIS package, the ModClark procedure was rewritten for this project in the Avenue scripting language. By doing this, a seamless connection between CRWR-PrePro and rainfall data preparation could be achieved. 8 At the heart of this research, however, rests CRWR-PrePro. Originally developed in the Arc/Info Macro Language (AML) by Hellweger and Maidment (1997), this procedure has been adapted in the Avenue programming language and provides a bridge between GIS data and application of these data in HMS. Currently, offices such as the Texas Natural Resource Information System offer a wealth of free GIS data of the state of Texas. In addition, groups such as HEC provide effective rainfall/runoff models. However, in the past, packages linking these two primary concepts have not existed. CRWR-PrePro, however, offers a methodology for bridging this gap. This research project represents one of the first applications of CRWR-PrePro on a large-scale scope. In fact, the primary objective of this research is to demonstrate the effectiveness of CRWR-PrePro in manipulating GIS data for use in modeling software packages such as HMS. A final important element of this research presents an updated study of the Buffalo Bayou region. Two groups have performed hydrologic studies of the Buffalo Bayou. The first study, executed by HEC in 1977, addresses the urbanization in the vicinity of the Addicks and Barker Reservoirs in the upstream portion of the bayou on the effectiveness of these reservoirs at flood control. While this study does not directly discuss this issue, the modeling results should provide a modern view of the additional impact of urbanization in the watershed since 1977. In addition to this HEC study, the Harris County Flood Control District (HCFCD) contracted the consulting firm Bernard Johnson, Inc. to perform a management study (1992) of the Addicks Reservoir. As with the HEC report, this research should also supplement 9 this study. The use of NEXRAD and GIS data in the creation of a rainfall/runoff model of the Buffalo Bayou should provide further insight into the management needs of the both the reservoirs and the bayou. 10 Chapter 3: Data Collection and Preprocessing While the Galveston District of the Corps of Engineers requested that the model created for this study specifically address the Buffalo Bayou region upstream of Houston proper, it also requested that the GIS database assembled cover all of the Houston metropolitan area. Figure 3.1 shows the area for which GIS data were collected, a range from 95? W to 96? W in longitude and from 29? 30? N to 30? 15? N in latitude. A fundamental aspect of working within a GIS environment is the map projection one wishes to use for the analysis. Since the earth is a spheroid, to work with GIS data, one must first choose a map projection such that the three-dimensional data may be represented in two-dimensional space. The map projection for this project was specified by the USACE. HEC recognizes two basic map projections for NEXRAD rainfall data, the Hydrologic Rainfall Analysis Project (HRAP) grid and the Standard Hydrologic Grid (SHG). For this project, the District specified that the SHG be used. The parameters for the SHG map projection are as follows: STANDARD HYDROLOGIC GRID MAP PROJECTION PARAMETERS Projection: Albers Datum: North American Datum of 1983 (NAD83) Ellipsoid: Geodetic Reference System of 1980 (GRS80) Map units: meters Central Meridian: 96? W (-96.00) Reference Latitude: 23? N (23.00) 11 Figure 3.1: Area covered by the GIS database 30? 15? N, 96? W 29? 30? N, 96? W 30? 15? N, 95? W 30? 15? N, 95? W 12 Standard Parallel 1: 29? 30? N (29.50) Standard Parallel 2: 45? 30? N (45.50) False Easting: 0.0 False Northing: 0.0 As part of the GIS data preprocessing, all data were converted into this map projection. An additional basic concept within a GIS is the difference between raster and vector data. A raster data set divides a rectangular domain into a mesh of square cells. In a GIS, this raster set is referred to as a grid. Two examples of a grid are a DEM and NEXRAD rainfall data. A vector data set, on the other hand, represents points, lines, or polygons and is referred to as a coverage. A line is an open sequence of points, and a polygon is a closed sequence of lines. Examples of point, line, and polygon coverages could be gauging stations, a stream network, and a watershed, respectively. Other examples of vector data that were assembled for this project include the Anderson Land Use/Land Cover set and the STATSGO soils data. It is important to differentiate between raster and vector data in a GIS. For instance, one may multiply grids together using a process called map algebra in which mathematical operations are done on a cell-by-cell basis to cells in corresponding locations in the various grids defining the variables required. Each grid cell can have only a single value. Moreover, one must convert a coverage to a grid in order to multiply data in a coverage by data in a grid. Conversion operations along these lines are included in the ArcView software package. 13 The following sections describe the acquisition and preprocessing procedures performed on each data set. These descriptions include information such as data source, original map projection, and algorithms used for preprocessing. A full data dictionary may be found in the Appendix. 3.1 Preparation of DEM Data for Use in CRWR-PrePro The Texas Natural Resources Information System (TNRIS) offers thirty-meter DEMs over the internet free of charge at www.tnris.state.tx.us. As a raster data set, this DEM contains elevation points of the earth?s surface in a grid spaced at thirty- meter intervals. These data are interpolated from the contour elevation data contained in 1:24,000 USGS topographic map sheets. The TNRIS organizes the DEM data by USGS 7.5 minute quad names with elevation points measured in feet. For the data collection area of interest, this includes 48 7.5' quadrangles, as shown in Figure 3.2. Note that the names are properly oriented spatially to each other. As taken from the TNRIS, the DEM data are not immediately usable in a GIS. To save storage space and expedite the download time, the TNRIS stores DEM data in text form on its web site. Fortunately, the USGS has specified a particular method for storing DEM data in text form so conversion into a useful product is simple. Included in the Spatial Analyst Extension of ArcView is the Import Data Source command, which allows a user to load these text-based DEMs into the ArcView environment. One simply indicates that the DEM is in USGS DEM format and then 14 Figure 3.2: DEM names of interest on the TNRIS server 29? 30? N, 96? W 29? 30? N, 95? W 30? 15? N, 96? W 30? 15? N, 95? W 15 specifies the location of the file of interest. Once the import process is complete, the user may view the DEM data in ArcView. However, one should pay attention to whether the TNRIS has labeled the data sets properly since some data sets downloaded for this project were incorrectly named while others were in vertical units of meters instead of feet. Fortunately, it was obvious when errors such as these were present since, for example, mislabeled DEM sets were spatially correct even if named incorrectly. Since the DEM data were available from the TNRIS as forty-eight separate data sets, these sets first had to be merged together into a single DEM before any further analysis was performed. Using the Merge command within the Grid subsection of Arc/Info conveniently and quickly accomplished this task. For instance, to merge two DEMs arbitrarily named DEM1 and DEM2 requires Arc: grid Grid: merged_dem = merge(DEM1,DEM2) The merge command looks for overlap between the two DEM sets and appropriately matches them together. Unfortunately, the merged DEM was marred by a severed problem, as may be seen in Figure 3.3. 16 Figure 3.3: Gaps existent in the DEM data While the black border on the left edge represents the physical edge of the data set (i.e., NODATA cells), the black lines present in Figure 3.3 demonstrate the presence of gaps between some of the forty-eight individual DEMs sets when they were combined into a single DEM. For whatever reason, the gap was always six cell values in height while the width of the gaps varied. This, of course, posed a problem since ArcView considers these areas as locations with no data values. 17 After consideration, the hurdle was overcome as follows. First, using the ArcView Export Data Source command, the entire DEM was converted into an ArcView ASCII-based grid. A file of this type has a header at the top which describes the number of rows, the number of columns, the location of the lower-left hand corner of the grid, the grid?s cell size, and the value a cell would possess for ArcView to consider it a no-data cell. After this header, the file contains the values of each individual cell. After removing the header portion of the DEM ASCII file, the file was loaded into a spreadsheet. Since the area of consideration is flat and since the height of the data gaps was only one hundred eighty meters, which is small considering the size of the data set, a simple interpolation technique was used to generate values to fill in the data gaps. After replacing the header in the file, the text- based file could be re-imported back into ArcView. Figure 3.4 shows the interpolated DEM. 18 Figure 3.4: Interpolated DEM The next aspect of the DEM to consider is its map projection. As its standard map projection, the TNRIS employs the Texas State Mapping System, whose parameters follow. TEXAS STATE MAPPING SYSTEM PARAMETERS Projection: Lambert Conformal Conic Datum: North American Datum of 1983 (NAD83) Ellipsoid: Geodetic Reference System of 1980 (GRS80) Map units: meters Central Meridian: 100?W (-100.0000) Reference Latitude: 31? 10? N (31.166667) Standard Parallel 1: 27? 25? N (27.416667) Standard Parallel 2: 34? 55? N (34.916667) False Easting: 1000000 False Northing: 1000000 19 As this is different from the SHG map projection, a conversion process must be employed to switch between the two map projections. To accomplish this, Arc/Info contains a reprojection command which has the following form: Arc: project grid in_name out_name projection_file.txt The projection file in this command line contains the projection parameters of the original and newly-projected grids. The following projection file was used in this situation: input projection LAMBERT units METERS datum NAD83 spheroid GRS80 parameters 27 25 00 34 55 00 -100 00 00 31 10 00 1000000.0 1000000.0 output projection ALBERS datum NAD83 units METERS spheroid GRS80 Parameters 29 30 00 45 30 00 -96 00 00 23 00 00 0.00000 0.00000 end 20 Note that switching between the Lambert and Albers map projections alters the resolution of the DEM from 30 meters to 30.8 meters. Before moving on to other concepts concerning the area terrain, a final aspect of the DEM preprocessing must be considered. As taken from the TNRIS, DEM elevation points are measured in feet above the employed datum, which in this case is the North American Datum of 1983 (NAD83). The CRWR-PrePro system used to analyze the GIS data in ArcView requires that the DEM elevation points be measured in meters. By dividing the reprojected DEM by 3.281 using floating point division within ArcView, values were converted from feet to meters. 3.1.1 Photographs in and around Addicks Reservoir In order to understand better the topography of the terrain, photographs of the Addicks reservoir area were taken. While previous studies had indicated that slopes in the region averaged around five feet per mile, which is very flat, it was decided that pictures would provide an even better idea of the lack of terrain relief in the area. In addition, the photography trip provided additional insight about the area, such as how the reservoir itself is designed, errors in maps about the reservoir, and drainage problems present. The Addicks reservoir, like the Barker reservoir, possesses berms which were constructed using soil. The picture below, taken from inside the east edge of the Addicks reservoir, indicates the height of these berms. 21 Figure 3.5: View of berm from inside Addicks Reservoir A startling revelation gleaned for acquiring these pictures is the general poor shape in which the drainage system upstream of the Addicks reservoir is maintained. Figure 3.6, which is of a major set of culverts in the soccer fields in the northwest corner of the reservoir, typifies the condition of many of the culverts witnessed while taking pictures of the drainage features of the reservoir.. Finally, the last photo, Figure 3.7, is of the exit point of Addicks Reservoir. It is included to give the reader an idea as to the size of the gates available to discharge water from the reservoir. Notice the man on top of the blue metal structure. 22 Figure 3.6: Blocked culverts in the upstream part of Addicks Reservoir Figure 3.7: Exit point of Addicks Reservoir 23 3.1.2 Creation of a Triangular Irregular Network (TIN) A TIN describes a three-dimensional surface while occupying much less memory than a DEM. Essentially, a TIN is comprised of a set of contiguous, non- overlapping triangles. Each node where triangles meet is assigned an elevation value. Lines connecting the nodes (i.e. edges of triangles) are straight which means that elevations between nodes may be linearly interpolated. The triangles within the TIN are relatively small in areas where the terrain possesses fine details and large in areas which are flat. A TIN is pertinent to this project in that it provides an additional method for viewing the land?s surface throughout the modeling area. Initial inspection of the DEM data set indicated that the slope in the region averaged 0.0010, or slightly over five feet per mile. Because of this low slope, it was a concern that the DEM itself was not going to be able to provide completely accurate information concerning the shape of the stream network. A TIN, however, does provide an additional source of information for deriving the stream network and, in fact, proved invaluable in determining the orientation of the stream network. The ArcView extension 3-D Analyst adds functionality for creating a TIN from a DEM. In this case, with the 3-D Analyst activated, one merely has to "Convert Grid to TIN" under the "Theme" menu list to create a TIN. To show the difference between a TIN and a DEM, the figure below shows a small area within the modeling region using a DEM and using a TIN. Note that around the stream network, the TIN provides a better three-dimensional visualization of the ground 24 surface. This improvement primarily resides in that the finer details of the terrain relief may be gleaned from the TIN since the TIN does a better job of shading small differences in elevation. If a DEM is broken down such that each color represents a ten-meter block of elevation values, all values, say, between ten and nineteen meters would be the same color. In a TIN, however, one could discern slight differences in elevation because of the shading, which is enhanced by the presence of each of the triangles which make up the TIN. Being able to recognize these slight differences in elevation proved to be useful considering the general lack of relief in Houston-area topography. Figure 3.8: DEM (left) and TIN (right) of the same area Of interest is the fact that the TIN of the modeling region occupies 6.4 megabytes of memory while the floating-point DEM of the same area resides in 14.5 megabytes. 25 3.2 Preparation of NEXRAD Data The Princeton Environmental Institute (the Institute) supplied all rainfall data for the 16-18 October 1994 storm, which began about 1530 hrs on the 16 th and ended around 2230 hrs on the 18 th . As received from the Institute, the data were zipped into three files, with each day of the storm covered by one file (e.g. 1016.zip, 1017.zip, 1018.zip). These zip files contain many different spreadsheet-readable text files, each of which held instantaneous rainfall information for a given time in a geographic map projection. To determine the time associated with a given text file, note that the names are set up in an mmdd_hhmmss format. For instance, 1016_152947 indicates that the data in that file are for a little before 1530 hrs on October 16 th . Though not separated by a constant interval, these rainfall data files were separated by approximately six-minute intervals. For example, the first three files in 1016.zip were 1016_152947 1016_153431 1016_154031 Note that the first and the second file are separated by four minutes and forty-four seconds while the second and third files were separated by exactly six minutes. Each of these files holds ten thousand points of rainfall intensities (mm/hr) at one-square-kilometer resolution in a latitude-longitude map projection. Furthermore, the ten thousand points formed a 100- x 100-cell grid with each cell assigned a 26 latitude, a longitude, and a rainfall intensity. Figure 3.9 shows how the rainfall data are oriented spatially. Figure 3.9: Area covered by NEXRAD rainfall data The preprocessing of the rainfall data contained two distinct aspects. The first aspect considered reprojection problems between the original geographic data and their SHG counterparts. The second required that the rainfall data be converted from a column-based format into a grid-based format. The end goal of the preprocessing of the NEXRAD data was a set of files which conformed to ArcView?s ASCII raster file format, which looks like ncols integer nrows integer xllcorner real 30.39? N, 96.29? W 29.50? N, 96.29? W 30.39? N, 95.27? W 29.50? N, 95.27? W 27 yllcorner real cellsize integer NODATA_value integer cell_value cell_value cell_value, etc where ncols = the number of columns in the grid nrows = the number of rows in the grid xllcorner = x-value of the position of the lower-left-hand corner of the grid yllcorner = y-value of the position of the lower-right-hand corner of the grid cellsize = the cell size of the grid NODATA_value = number assigned to a cell with no data in it The number of rows and the number of columns both equal 100 since the Institute indicated that this was the case. In addition, the Institute also specified that the cellsize equals 1000 meters since the rainfall data are at 1-km 2 resolution. Next, the NODATA_value was arbitrarily set to -99999. The xllcorner and the yllcorner values, however, required a bit more thought. The lower-left-hand corner values in the SHG map projection were obtained by taking the coordinates of the lower-left hand corner of the rainfall intensity grids in the geographic map projection (e.g. from 29.50? N, 96.29? W) and converting them into the corresponding SHG map projection coordinates. An Avenue program named defgage.ave that converts a point shape file from a geographic map projection to any Albers equal area map projection was used. This program, which may be found in the Appendix, calculated that the xllcorner value equals -28715 meters while the yllcorner value equals 713439. Using a cell size of 1000 meters, the upper-right-hand corner of the grid resides at (70985 m,813139 m). Unfortunately, visual inspection of a grid using these header values reveals that a cell size of 1000 meters is not correct. Consider Figure 3.10, which is in the SHG map projection. 28 Figure 3.10: Inadequacy of 1000 meter cell size in ASCII files 29 The top image in Figure 3.10 shows the lower-left-hand corner of a rainfall grid using a cell size of 1000 meters while the bottom image shows the upper-right- hand corner of the same grid. The dotted grid was made by taking the latitude/longitude points from the original rainfall data set and reprojecting them into the SHG map projection directly. Each dot represents a rainfall intensity from the original rainfall data sets. The square-based grid underneath the dotted grid was created using the header file described earlier. Inspection of the lower image in Figure 3.10 reveals that the right-most column in the salmon grid does not contain any intensity values (i.e. green dots). Notice how the right-most column, whose boundaries can be determined by looking at the difference in coloration in the upper right-hand-corner of the lower image, does not have any green dots in it. Each dot must have its own unique cell in the square-based grid. Keep in mind that the dot does not have to fall in the middle of its unique cell; it merely has to fall somewhere inside of it. Hence, initial guesses as to a solution figured that changing the cell size value in the header for the grid would solve this problem. After some trial and error, a cell size of 997 meters generated a grid whose upper-right-hand corner is shown in the following image. 30 Figure 3.11: Upper-right-corner of rainfall grid using cell size of 997 meters Clearly the problem has been solved. Since each cell in the far-right column has its own unique dot, the reprojection from a geographic map projection to the SHG map projection was now successful. This discrepancy between the grid resolution as defined by the Institute and the grid resolution used in the SHG map projection may be explained as follows. The difference resulted from the map projection conversion from the original geographic map projection to the SHG map projection. Where the geographic map projection takes into account the curvature of the earth, the SHG map projection, being an Albers map projection, is flat. Thus, the conversion from the geographic map projection to the SHG map projection warranted a change in cell size to ensure that each of the ten thousand ?rain gauges? in the NEXRAD rainfall grid was represented by a unique cell in that grid. With this problem overcome, the second piece of the preprocessing of the NEXRAD data could be tackled. 31 Completion of the second aspect of the preprocessing necessitated solving three primary problems. First, the ten thousand points in each data set, while originally oriented in a single column, needed to be reorganized in a 100- X 100-cell grid. Second, the irregular time steps between each of the rainfall files had to be converted into a single consistent time step. A time step of ten minutes was chosen since this value would be small enough to capture the details of the storm, yet large enough that routing problems do not arise from too small a time step. Finally, as HEC-HMS requires rainfall depths when using grid-based rainfall, conversion from intensities to depths was needed. A MATLAB computer code was written to solve these problems; these programs, which are described below, may be found in their entirety in the Appendix. Since these rainfall data could always be represented using matrices whether in column or grid form, MATLAB was chosen since a primary strength of MATLAB is its ability to manipulate matrices. The pseudocode of the programs is listed below. Program 1 ? convert.m ? Changes time step from irregular values to a constant ten-minute value. Then switches rainfall intensities to rainfall depths. ? Load original rainfall data ? Load times associated with each of the separate rainfall files ? Remove first two columns from rainfall data (latitude/longitude values) ? Determine ten-minute value closest to time of first data set of the day -- i.e. 21:42:39 yields 21:40:00 32 ? Reorganize data from irregular time steps into regular ten-minute time step o Include data sets until combined time-interval of data sets exceeds ten minutes o Convert intensities in each data set to depths of rain since start time of set o Sum the rainfall accumulated in the ten-minute time step o Store rainfall in remaining time for the next time step o Find data sets in the next ten-minute interval and keep looping until finished ? Store results in a temporary text file. Program 2 ? makefile.m ? Switches output of program 1 from column format into 100x100 grids. Assigns header required by ArcView ASCII grids to top of each new file. ? Load output from program 1 ? Take each column from this output and send it to a filer function: o Convert to a 100x100 grid. Start building grid at lower left- hand corner, building rows before columns. o Attach ASCII header to top of grid. This header is described earlier in this section. o Output results to file. Application of these two programs to the Institute rainfall intensity data sets resulted in 331 grids, each of which holds the rainfall depth values for one ten-minute interval in the storm. With the pre-processing of the NEXRAD data complete, work on the digital line graph data sets could begin. 3.3 Preparation of Digital Line Graph Data for the Stream Network A typical problem when working with DEM data, particularly when considering an area as flat as the Houston region, stems from a lack of resolution in 33 the DEM. In fact, this dilemma represented one of the major issues facing this project. With slopes around five feet per mile, the modeling area proved to be too flat to describe adequately the stream network. Initial attempts at determining the stream network based on the DEM alone resulted in streams which made no sense based on inspection of the area. Consider Figure 3.12: Figure 3.12: Demonstration of inadequacies in original DEM 34 The black lines represent the delineated stream network on the unaltered DEM. However, the exit points from each reservoir on the stream network, marked by the black arrows, clearly do not agree with the delineated exit points. To mitigate this problem, 1:100,000-scale digital line graph (DLG) data were used. Essentially, these data are digital representations of cartographic information. Data files of topographic and planimetric map features are derived from cartographic source materials using manual and automated digitizing methods. Intermediate or 100,000-scale DLG data are derived from U.S. Geological Survey (USGS) 1:100,000-scale 30- by 60-minute quadrangle maps. If these maps are not available, Bureau of Land Management planimetric maps at a scale of 1:100,000 are used. The 1:100,000-scale DLG data distributed by the USGS are in the DLG- Level 3 (DLG-3) format. A DLG-3 file contains a full range of attribute codes, has full topological structuring, and has passed certain quality-control checks. (source http://www.ce.utexas.edu/prof/maidment/gishyd97/library/websites/dlg1doc.htm) The DLG data source for this project was a CD-ROM from the USGS entitled ?1:100,000-Scale DLG Data Hydrography and Transportation?, which contains 1993 DLG information over Texas and Oklahoma including the 1:100,000-scale Beaumont and Houston maps. This CD-ROM provided three different types of raw DLG data? political boundary data, hydrography data, and transportation data. Of these three data sets, the hydrography set possessed the greatest value. While the political boundary data did little more than show county boundaries and while the transportation data displayed major transportation routes, the hydrography data provided an excellent source for delineating the stream network. Before discussing 35 the benefits of these data, it is important to note that some preprocessing was required before the data could be used in an ArcView format. Conveniently , an AML for Arc/Info was available for preprocessing the DLG data. This AML may be found at http://www.ce.utexas.edu/prof/maidment/ gishydro/docs/amls/dlg.htm. As this is a generalized piece of code, the exact code used to process the hydrography data for the Houston area may be found in the Appendix. In a nutshell, this code first converts the data from the DLG optional format (as found on the CD-ROM) into coverages which Arc/Info recognizes. From here, the code creates the topology of the data set using the ?build? command. Since each of the coverages created only represents a 15? map with a border around it, the separate sets had to be combined and the borders removed. Finally, the coverages had to be reprojected into the SHG map projection. Once this series of steps was completed, the data sets were usable for this project. For convenience, a piece of the completed DLG data is shown in Figure 3.13. Note the berms of the two reservoirs near the center of the image. 36 Figure 3.13: DLG hydrography features of the study area 37 While this data set proved to helpful, it suffers from the fact that while the features are spatially correct, there is no labeling system to identify what each feature is in the set. For instance, the two reservoirs in the system, located toward the upper- left-hand corner of the above image are no different from a data point of view from any of the other lines. While this posed a problem, it was eventually overcome. As the discussion in this section is confined to data preprocessing, solutions to this problem will be discussed later in this thesis. 3.4 Preparation of USGS Flow Data During the storm event of interest, the USGS maintained seven flow gauges and two reservoir storage gauges in the study area. The table below lists the names of these gauges as well as their positions in latitude and longitude coordinates. Gauge Name Latitude Longitude Cypress Creek at Katy-Hockley Road near Hockley, TX 29? 57? 00" N 95? 48? 29" W Cypress Creek at House-Hahl Road near Cypress, TX 29? 57? 32" N 95? 43? 03" W Barker Reservoir near Addicks, TX 29? 46? 11" N 95? 38? 49" W Bear Creek near Barker, TX 29? 49? 50" N 95? 41? 12" W Langham Creek at West Little York Rd near Addicks, TX 29? 52? 01" N 95? 38? 47" W Addicks Reservoir near Addicks, TX 29? 47? 28" N 95? 37? 24" W Buffalo Bayou at West Belt Dr at Houston, TX 29? 45? 43" N 95? 33? 27" W Buffalo Bayou at Piney Point, TX 29? 44? 48" N 95? 31? 24" W Buffalo Bayou near Katy, TX 29? 44? 35" N 95? 48? 24" W Table 3.1: Names and locations of pertinent USGS flow and storage gauges 38 At http://txwww.cr.usgs.gov/databases.html, which is within the USGS Texas web page, one may request via email archived data from the USGS. Thankfully, Dr. Fred Liscum honored this request and sent materials via mail to CRWR. The data he sent document the water elevation above a given datum at a given gauge as a function of time. The time interval for these tables is one hour. Data were requested for the time period covering midnight on 1 October 1994 to midnight on 15 November 1994. In USGS terminology, the initial tables with water surface elevations are called primary gauge heights. To translate these primary gauge heights into actual flow values, the USGS uses rating tables. Before these rating tables may be used, however, the gauge heights themselves need to be converted from primary gauge heights to adjusted gauge heights. The relationship between water surface elevation and flow at a particular site is not constant over time due to factors such as vegetation growth, the accumulation of sediment in the channel, and the scouring of the channel from erosion. To address this issue, the USGS periodically uses field observations linking water surface elevation and flow to create variable shift tables. These tables provide the adjustments required to shift the primary elevation values to their associated adjusted elevation values. As an example, consider the Katy gauge from the above table. From here, consider three arbitrary primary elevations which may have occurred on 20 October, such as 21.00 feet, 27.00 feet, and 33.00 feet. Now consider the variable shift table for this gauge, which is partially reproduced below. 39 DATE/TIME INPUT SHIFT INPUT SHIFT INPUT SHIFT --------------------------------------------------------------------------------- 09 10 2300 23.00 1.22 24.50 1.22 26.50 0.00 10 18 2300 23.00 1.22 24.50 1.22 26.50 0.00 --------------------------------------------------------------------------------- 10 18 2400 23.00 1.24 24.00 1.24 31.00 0.00 02 28 0400 23.00 1.24 24.00 1.24 31.00 0.00 --------------------------------------------------------------------------------- 02 28 0500 22.00 1.13 24.00 1.13 31.00 0.00 ? ? ? ? ? ? ? NB: The shifts during each time interval are fixed. The USGS uses field observations to update the shift values from time interval to time interval. Table 3.2: Partial variable shift table for Buffalo Bayou near Katy, TX Since the date for this discussion is 20 October, the time interval of interest on the variable shift table is 2400 hrs on 18 October 1994 to 0400 hrs on 28 February 1995. For primary elevation values less than 23.00 feet in this interval, the shift is 1.24 feet. Thus, the first arbitrary primary elevation of 21.00 feet would have an adjusted value of 22.24 feet. Similarly, the third value of 33.00 feet would have an adjusted value of 33.00 feet since primary elevations above 31.00 feet do not need to be adjusted. The middle arbitrary value, 27.00 feet, requires a somewhat more complex process for it to be adjusted. As the shift at a primary elevation of 24.00 feet is 1.24 feet, and the shift at a primary elevation of 31.00 feet is 0.00 feet, an interpolation procedure is used to convert from primary to adjusted elevation between 24 and 31 feet. In this case, 27.00 feet is three-sevenths of the way from 24.00 feet to 31.00 feet. Thus the shift is (1-3/7)*1.24 or 0.71 feet, and the adjusted value would 40 then be 27.71 feet. The conversion of all primary elevation data to adjusted elevation data occurred using an Excel spreadsheet. Once the adjusted gauge elevations were available for each of the gauges (note that for a reservoir the primary and adjusted elevations equal each other), USGS rating tables were employed to deduce gauge flows (or storage values) for each gauge. Rating tables contain, for each adjusted elevation value, a corresponding flow value (or storage value); simply take a given adjusted gauge elevation, look at the rating table, and find the associated flow (or storage) value. The process of converting all adjusted elevation values to flow (or storage) by hand would require an extensive amount of time. To expedite this procedure, a database containing all of the rating table data was created. From there, MATLAB programs were again employed. This simple program, findflow.m, may be found in the Appendix; its pseudocode appears below: ? Read adjusted elevation data for a given gauge ? Read rating table data for the same gauge ? Loop through all adjusted elevations o For each elevation, find the two elevations on the rating table on each side of the given elevation o Using interpolation and the information from the rating table data, convert the adjusted elevation to the proper flow value ? End loop Once this program was executed on each of the adjusted gauge values, flow and storage values were available. Values generated by these programs could be checked 41 for accuracy since the data from the USGS included hourly average flow values. Note that flow values had units of ft 3 /sec and storage values had units of acre-feet. An additional way to view this process of generating flows from elevations is graphically. First consider that two relationships are available for the USGS? elevation vs. time and flow vs. elevation. The top two charts in Figure 3.14 reveal this relationship for the Katy gauge from midnight on 17 October 1994 to midnight on 25 October 1994. The key is to realize that the following equality also applies: The flow vs. elevation relationship provides a single number by which to multiply every point on the elevation vs. time chart such that a new relationship, flow vs. time, is produced. As a result, because flow is proportional to elevation, the flow vs. time curve looks like the elevation vs. time curve. However, the primary difference is that since larger elevations are multiplied by larger flows, the flow vs. time curve is higher and narrower than its elevation vs. time counterpart. Time Flow Elevation Flow Time Elevation =* Eqn. 3-1 42 Figure 3.14: Graphical relationship among elevation vs. time (top), flow vs. elevation (middle), and flow vs. time (bottom) Elevation vs Time 25 27 29 31 33 35 37 0 50 100 150 200 250 Time (days) E l evation (ft) Flow vs Elevation 0 500 1000 1500 2000 2500 20 25 30 35 40 Elevation (ft) Flow (cfs) Flow vs Time 0 500 1000 1500 2000 2500 0 50 100 150 200 250 Time (days) F l o w (cfs) 43 3.5 Acquisition of Anderson Land Use/Land Cover data The TNRIS offers Anderson Land Use/Land Cover data for the region of interest through its FTP site (ftp://www.tnris.state.tx.us/pub/GIS/). Data for the region were originally collected by the USGS from 1:100,000 and 1:250,000 maps and then converted to Arc/Info format by the EPA in a geographic map projection. The TNRIS then took these data and converted them into the Texas State Mapping System Lambert map projection. Two land use/land cover data sets were downloaded from the TNRIS--Beaumont and Houston. To expedite downloading, the TNRIS stores data sets of this type in a zipped format on its web site, so the first step for preparing the data involved unzipping the files using PKUNZIP. Unzipping a given file, such as houston.zip, generated a file named houston.e00. Using the ESRI program IMPORT 71, these *.e00 files were converted into files which ArcView could readily recognize. Before use in ArcView, however, these data sets needed to be converted from the Texas State Mapping System Lambert map projection into the SHG map projection. This was accomplished within the Arc/Info environment. Once reprojected, the data were viewed in the ArcView system. Viewed with ArcView, the land use data sets were initially monochromatic and uninformative. To get around this, an ArcView legend file was created. The parameters within this legend follow. 44 Value Label Symbol Color 0-9 Unknown White 10-19 Urban Red 20-29 Agriculture Yellow 30-39 Range Land Green 40-49 Forest Dark Green 50-59 Water Blue 60-69 Wetlands Light Blue 70-79 Barren Grey Table 3.3: Legend for land use data Once this process was complete the Land Use/Land Cover data were easy to read and use. A cursory check of the data set indicated that the agriculture label dominated most of the area in the Addicks/Barker reservoir region, although small percentages of wetlands and forests were present. The following table summarizes the relative presence of each of the different land-use types in the modeling area. Value Label Percentage 0-9 Unknown 0.00 10-19 Urban 13.5 20-29 Agriculture 67.9 30-39 Range Land 0.00 40-49 Forest 17.1 50-59 Water 0.29 60-69 Wetlands 0.99 70-79 Barren 0.22 Table 3.4: Percentages of land-use types in modeling area While some parts of the region have been developed since these data were accumulated, visual inspection of aerial photography (see section 3.8) of the region reveals that most of the agricultural area has remained agricultural. 45 3.6 Acquisition of the STATSGO Data Set As with the Land Use/Land Cover data, the TNRIS has available STATSGO soils data for the entire state of Texas on its web site (ftp://tnris.state.tx.us/GIS/ land_use/statsgo/). After the data had been downloaded from the TNRIS site, they were moved into a UNIX environment since these files required the UNIX command gunzip for unzipping. The end result of using gunzip was a standard *.e00 file which could be prepared for use in ArcView using ArcView?s IMPORT 71 utility. IMPORT 71 produced a STATSGO soils data set in the Texas State Mapping System Lambert map projection for the entire state of Texas. After reprojecting the data within Arc/Info to the SHG map projection, the data were loaded into ArcView. To make the data more readable, the Legend Type in the Legend Editor was changed from single color to unique value, with muid (i.e. map unit identifier) being the field used to determine unique values. After preparing the data set for ArcView, the comp (i.e. soil component identifier) and muid tables, located in the INFO directory associated with the data set, were loaded to get an idea concerning what types of soils were present in the area. Conveniently, one particular muid--TX248--dominated the area upstream of the Addicks and Barker reservoirs. This muid contains ten separate components. The primary surface texture, making up 91% of the muid, is fine sandy loam. Additionally, 93% of this muid has a hydrologic soil group type of D, which means that it drains poorly. In a more general sense, consider the next table, which breaks down the entire modeling region by hydrologic soil group. 46 Hydrologic Group Name Percentage Present A 0.0479 B 0.144 C 10.7 D 89.2 Table 3.5: Breakdown of hydrologic groups in study area Note that this information plays a role later when deciding what sort of initial infiltration parameters should be used with the HMS model. 3.7 Preparation of the Curve Number Grid Traditionally, one combines information from land use/land cover data and soils data to calculate what is known as the curve number grid. This concept, established by the Soil Conservation Service, assigns a curve number, CN, such that 0 &1 ,I DQ DUHD KDV D &1 RI WKHQ DOO RI WKH UDLQZDWHU ZKLFK ODQGV LQ WKLV area will runoff; conversely, if CN equals 0, then the soil absorbs all water. Thus, one may use this grid to determine what amount of rainfall will turn into direct runoff in a storm. While the USACE requested that the Anderson Land Use/Land Cover data and the STATSGO soils data be acquired for this project, it did not specify that the curve number grid for this analysis had to be created from these two results. While that was the original plan for generating the curve number grid for this study, an alternate curve number grid was provided by Dr. Francisco Olivera at CRWR. Dr. Olivera had created a CD-ROM in May, 1998, entitled "U.S. Spatial Hydrology 47 Database" which contains a 250-meter curve number grid of the entire country with the following parameters for its map projection: Projection ALBERS Zunits NO Units METERS Spheroid CLARKE1866 Xshift 0.0000000000 Yshift 0.0000000000 Parameters 29 30 0.000 /* 1st standard parallel 45 30 0.000 /* 2nd standard parallel -96 0 0.000 /* central meridian 23 0 0.000 /* latitude of map projection?s origin 0.0 /* false easting (meters) 0.0 /* false northing (meters) Dr. Olivera suggested that this grid be used for this project since, in his opinion, it would have been unlikely that a more accurate grid could have been generated from the available land use/land cover and soils data. In addition, as indicated by this curve number grid, the curve number for most of the actual modeling area was 70. With this kind of consistency present in the curve number grid, it was deemed unnecessary that a new grid be made specifically for this project. 3.8 Acquisition of Digital Orthophoto Quarter Quadrangles for Study Area While digital orthophoto quarter quadrangles (DOQQs) exist for only parts of the state of Texas, the great majority of the data collection area and the entire modeling area for this study have been photographed. DOQQs, made available through the Texas Orthoimagery Program and the Texas Strategic Mapping Program, may currently be downloaded from the TNRIS web site (http://www.tnris.state.tx.us/ 48 DigitalData/DOQ_2.5m/2_5mdoqs.htm) in 2.5-meter resolution, although photo resolutions of one, ten, and 30 meters may be ordered. Earlier in this project, the TNRIS allowed users to download 30-meter DOQQs; these were used to save space since a DOQQ of this resolution occupies 175 Kb. This makes manipulation of these data much easier considering that a 2.5-meter DOQQ of the same area would occupy about 8,000 Kb while a one-meter DOQQ would be 149 Mb. Other background information about DOQQs include that these orthophotos are scanned aerial photos that have been corrected to removed distortions. Furthermore, they are derived from 1:40,000 National Aerial Photography Program (NAPP) images taken between 1994 and 1997. . The quadrangle names listed in Table 3.2 above, which contains the names of the individual 30-meter DEMs in the area, also apply to the DOQQ sets except that four DOQQs reside in a given USGS quadrangle sheet in the interest of keeping file size low. Of the forty-eight quads of interest, complete DOQQs exist for all but three--La Porte, League City, and Plum Grove. As these missing data sets reside in the corners of the data collection area, far from the modeling area itself, their loss has no negative impact on this study. Once downloaded, the DOQQs were converted from their original map projection, Universal Transverse Mercator (UTM), to the SHG map projection. This reprojection occurred in the Arc/Info environment. To convert color images from one map projection to another, the following set of commands is used in Arc/Info. 49 imagegrid original_DOQQ.tif ucode color project grid ucodec1 acodec1 txtoshg.txt project grid ucodec2 acodec2 txtoshg.txt project grid ucodec3 acodec3 txtoshg.txt gridimage acode.stk none new.tif tiff compression kill ucode all kill acodec1 all kill acodec2 all kill acodec3 all The first line converts the image file into a stack of three grids. These three grids represent the red, green, and blue content on a scale of zero to 255 for each pixel in the picture. The next three lines of code reproject each color grid from the UTM map projection to the SHG map projection. The projection file used follows. input projection UTM zone 15 units METERS datum NAD27 spheroid CLARKE1866 xshift 0.0000000000 yshift 0.0000000000 parameters output projection ALBERS datum NAD83 zunits NO units METERS spheroid GRS80 xshift 0.0000000000 yshift 0.0000000000 Parameters 29 30 0.000 /* 1st standard parallel 45 30 0.000 /* 2nd standard parallel -96 0 0.000 /* central meridian 23 0 0.000 /* latitude of map projection?s origin 0.00000 /* false easting (meters) 0.00000 /* false northing (meters) end 50 Finally, the gridimage command line converts the reprojected color stack of grids back into a usable image file. This process represents the technique used to reproject one image. However, in this case, multiple images needed to be reprojected. In addition, these individual images had to be merged into a single image. To achieve this single image required converting all of the individual images into their own color stacks, merging all of the color stacks into a single color stack using the merge command with Grid, and then converting this large color stack into a single digital orthophoto of the entire region. As an example of what a piece of this final image looked like, consider Figure 3.15, which is an aerial photo of the Addicks Reservoir. 51 Figure 3.15: Digital orthophoto of Addicks Reservoir 52 Chapter 4: Terrain Processing Using CRWR-PrePro Dr. Maidment and Dr. Olivera?s research team at CRWR has worked on the CRWR-PrePro code for several years. This code, written for the ArcView environment in the Avenue programming language, considers a fundamental problem with GIS in general. Figure 4.1: Linking GIS Data and modeling software Few mechanisms currently exist which fill the role of box 2 in the above diagram. CRWR-PrePro, however, breaks down GIS data and creates an output file which is the input for box 3, which in this case is HEC-HMS. By using CRWR- PrePro to analyze the terrain data sets described in the previous chapter, one can save considerable time and energy when creating a model to be used in HMS. The following sections in this chapter describe the use of CRWR-PrePro for creating a basic model for HMS. 4.1 Descriptions of Functions within CRWR-PrePro Dr. Olivera maintains a web site at CRWR, http://www.ce.utexas.edu/ prof/olivera/prepro/prepro.htm, which documents the most recent versions of CRWR- Box 1 GIS Data Box 2 ????? Box 3 Modeling Software 53 PrePro and includes directions on how to use CRWR-PrePro. The work done for this project incorporates the version prepro03.apr from this web page. To begin using CRWR-PrePro, an ArcView user loads a CRWR-PrePro project file. As long as the user has the spatial analyst extension within ArcView, CRWR-PrePro may be executed. To develop a HMS model, CRWR-PrePro proceeds through a series of steps. The image below shows these steps in order. Figure 4.2: CRWR-PrePro menu The first primary step in the list above, Fill Sinks, analyzes the DEM and looks for sinks, or pits. Sinks on a DEM are places where water gets trapped. In other words, each of the elevation points surrounding a pit is higher than the elevation 54 of the pit itself; therefore, pits trap water and hinder flow. It is necessary in an analysis of this type that if water falls on to any of the cells in the DEM, the water is able to flow to the edge of the DEM via some route. Thus, the Fill Sinks function eliminates all pits from the DEM. While it is necessary for all pits to be removed, the inaccuracies of the DEM data set will create some pits that are not really there. Fortunately, this is not a problem since pits of this nature are minor, and smoothing them out of the DEM will not significantly alter the flow patterns in the DEM. It should be noted, however, that in some cases, a pit in the DEM should not be removed, such as one which represents a salt lake with no outlets. It is up to the user to ensure that only superfluous pits are removed. The second function creates from the DEM a flow direction grid. When water lands on one of the grid points on the DEM, it may flow one of eight different directions to an adjacent grid cell depending upon which direction offers the steepest descent. This flow direction grid indicates which direction water will flow from each point on the elevation grid. Thus, from this flow direction grid, a stream network may be created. Next, CRWR-PrePro makes a flow accumulation grid from the flow direction grid. Each value on the flow accumulation grid represents the number of grid cells that the grid cell in question drains. This flow accumulation grid plays an important role in the Stream Definition calculation section of CRWR-PrePro. When the user selects the threshold of the stream definition, the user informs CRWR-PrePro how many cells must be upstream of a given cell for that cell to be considered part of the 55 stream network. Obviously, all cells belong to the stream network, but only the cells that drain significant (i.e. a threshold number) of cells should be considered a river or a stream or some other sort of drainage mechanism. For this project, the stream definition threshold chosen is 8000 cells. While this number is arbitrary, it appears to create a dense enough stream network while at the same time not overcrowding the GIS screen with superfluous stream lines. With the stream network established, CRWR-PrePro then moves toward developing watershed boundaries. The Stream Segmentation (Links) algorithm breaks up each stream into its corresponding links in which each stream reach between tributaries is labeled as a unique entity. For instance, a stream shaped like a Y would have three links, two upstream, one downstream. After identifying all of the links with the stream network, CRWR-PrePro assigns an outlet to the downstream cell of each of the links. CRWR-PrePro carries out these two steps of identifying links and assigning outlets since each outlet will eventually be the exit location of a sub-watershed. Figure 4.3 illustrates how CRWR-PrePro separates a stream network into links and outlets. The left image in the figure represents a stream network while the right image shows a collection of links. Each link is drawn with a different color and possesses an outlet at its downstream cell. 56 Figure 4.3: A stream network broken into links The Add Outlets option allows the user to enter additional outlet points within the stream network. This function makes CRWR-PrePro more flexible in that it allows the user to specify certain points within the stream network as sub-watershed outlet points even though these additional points are not at the downstream ends of outlets. This function has great value for this project in that the locations of the USGS gauges may be properly placed within the stream network. This was accomplished by creating a shapefile of the USGS gauge points in a geographic map projection since the USGS specifies the locations of all of its gauges in decimal degrees on its web page. Using this shape file as input, one may add outlets to a CRWR-PrePro model. Schematically, adding outlets breaks links into smaller links. To create the sub-watersheds themselves, one invokes the Sub-Watershed Delineation section of CRWR-PrePro. This set of calculations takes each of the outlets, both those which CRWR-PrePro established when creating outlets from links 57 and those which the user added, determines which cells are upstream of that outlet but downstream of the nearest upstream outlet, and delineates this set of cells as a sub- watershed. CRWR-PrePro was used successfully to establish the stream network and the shapes of the sub-watersheds within the DEM of interest. While this is an excellent start to creating a surface-water flow model of a given area, no parameters have been assigned to any of the watersheds. Since it is difficult to assign parameters to a group of cells in a raster environment within ArcView, CRWR-PrePro now converts the streams and sub-watersheds into vector data sets so that tables better describing the system may be created and assigned to the appropriate elements of the model. Hence, the next step of the model, Vectorize Streams and Watersheds, does just that. By converting the streams and sub-watersheds into a vector format, not only can CRWR- PrePro now more easily assign data to the different watersheds in the modeling area, but the user may also manipulate the GIS data more easily. For example, the user may now merge sub-watersheds quickly using the next tool in the CRWR-PrePro battery. Merging sub-watersheds is usually necessary since a few extremely tiny sub- watersheds, with an area equal to only a few cells on the DEM, are created. Since water may flow through these minute watersheds more quickly than the ten minutes specified as the analysis time interval, small sub-watersheds may skew model results. Hence, it is in the best interest of the model that they be eliminated. In this case, all sub-watersheds smaller than two km 2 were merged with the next sub-watershed upstream. 58 The next two sets of calculations within CRWR-PrePro, Soil Group Percentages and Curve Number Grid, are required only if the user wishes to create a curve number grid based on land use/land cover and soils data. Since this model incorporates a curve number grid from a separate source (see Chapter 3.7), there was no need to invoke these two functions. In the Calculate Attributes section, CRWR-PrePro determines the remaining parameters it requires before it is able to create a HMS input file. The types of parameters calculated in this section include sub-watershed slopes, curve numbers, and lag times. The end goal of this block of calculations is to determine both abstraction and lag time data. The heart of the CRWR-PrePro system resides in the next set of calculations. Within the HMS Schematic section, CRWR-PrePro creates the input basin model file for the HMS processor. Features of the HMS model which CRWR-PrePro creates include sub-watershed connectivity, sub-watershed areas, and average sub-watershed curve numbers. In addition, Muskingum routing parameters (K and X) are included in the model. While the user has to insert additional data for the Buffalo Bayou into the HMS system before executing the model, the majority of the model will be prepared by the CRWR-PrePro code. 4.2 Use of CRWR-PrePro to Create the HMS Basin Model Using the thirty-meter DEM as described in Chapter 3.1, an exploratory run was made to get a general idea of where any troublesome areas in the DEM were 59 located. This preliminary investigation revealed obviously incorrect stream networks in certain areas on the DEM, especially in the vicinity of the reservoirs; therefore, it was decided that burning the streams into the DEM would be necessary. The best available stream data for the region were the DLG hydrography data described in Chapter 3.3. The process of burning in streams is given in the following Figure 4.4 . 60 Figure 4.4: Burning a stream network into a DEM 61 Essentially, burning in streams involves taking the stream network and converting it into a grid in the same map projection as the DEM onto which the streams will be burned. Once the stream grid is dropped onto the DEM, all points on the DEM which are not part of the stream network are raised an arbitrary value; in this case five hundred feet was used. Burning in the streams accomplishes two things. First, it ensures that a more accurate stream network will be created from the DEM. Second, it guarantees that once water enters a part of the stream network, it will stay within the network and eventually flow off of the DEM. Before burning in the streams took place, a grid of the stream network itself had to be created. Though the DLG data are the best available data for burning in the streams, this data set has its limitations. The data set?s primary weakness is that none of its features is labeled. Since hydrography features include rivers, streams, canals, ditches, irrigation channels, and pipelines, it was not easy to determine which aspects of the hydrography data should be included in the actual stream network. Therefore, for the next attempt with CRWR-PrePro, all data from the hydrography data set were eliminated except for those which were considered natural streams. These streams stood out compared to the rest of the data in that they exhibited a meandering quality whereas non-natural flow paths possessed man-made qualities such as straightness and sharp angles. While this natural stream network did not include some stream paths which were part of the network, this network nonetheless provided a good starting point for CRWR-PrePro analysis. Figure 4.5 shows a representative portion of the natural network used in this part of the analysis. 62 Figure 4.5: Natural stream network used when burning in steams After filling in the sinks, the flow direction and accumulation calculations occurred without incident. When delineating the stream network, a stream threshold value of eight thousand cells was chosen via trial and error. By choosing eight thousand as the threshold, CRWR-PrePro created a stream network which gave good insight as to the shape of the stream network without being overly complex. From this point, the link and outlet grids were created. Incorporating the flow and reservoir storage gauges into the CRWR-PrePro system required first building a shapefile of the gauge locations. From an Avenue script, a table of the gauges was made which held the location of each gauge in the geographic map projection. This program then reprojected these points into the SHG map projection. This new reprojected shapefile was then used to add the gauges as outlets in the CRWR-PrePro model. 63 The next several steps in the CRWR-PrePro process delineated the watershed boundaries. After converting these results from raster format to vector format, and after clipping out the appropriate watersheds, the modeled region looked like this. Figure 4.6: Watershed delineation using natural stream network Unfortunately, a cursory inspection of the results at this time indicated that the stream network which was burned into the original DEM lacked enough detail to model the region of interest properly. Consider the following table. Gauge Name USGS gauge area (km 2 ) Calculated gauge area (km 2 ) % Error Buffalo Bayou at Piney Point 821 793 -3.4 Addicks Reservoir 352 342 -2.9 Barker Reservoir 332 371 11 Langham Creek 63.7 58.3 -8.5 Bear Creek 55.7 89.3 60 Buffalo Bayou near Katy, TX 164 196 20 Table 4.1 Area comparisons using natural stream network 64 While the overall size of the watershed compared favorably with the official USGS value, the Bear Creek and Katy gauge areas exhibited substantial error. In an attempt to improve upon these results, the stream network burned into the original DEM was scrapped in favor of one that took into account the man-made aspects of the hydrography data set. Any piece of the hydrography data set that looked like it belonged to the stream network was included as part of the stream network. As before, CRWR-PrePro analyzed this stream network to the point where a vector image of the watersheds could be viewed. The image below shows the watershed system delineated using this stream network as in input. Figure 4.7: Delineated network using a combination of natural streams and man-made channels 65 Again, calculated drainage areas represent the best parameter for determining whether the watersheds have been delineated properly. The following table reveals how well this stream network fared. Gauge Name USGS gauge area (km 2 ) Calculated gauge area (km 2 ) % Error Buffalo Bayou at Piney Point 821 816 -0.6 Addicks Reservoir 352 375 6.5 Barker Reservoir 332 330 -0.6 Langham Creek 63.7 90.7 42 Bear Creek 55.7 47.7 14 Buffalo Bayou near Katy, TX 164 155 -5.5 Table 4.2: Area comparisons using a combination of natural streams and man-made channels While this represents a clear improvement over the first attempt at modeling the system, the large error present for the Langham Creek gauge indicates that a better match might be achieved. To create yet another stream network for the burning process required some method of determining which pieces of the hydrography data set actually belonged in the stream network. Using the 3-D Analyst extension within ArcView, a TIN of the DEM was created. Via careful inspection of the TIN, a third stream network was created by looking at the hydrography data set draped over the TIN and determining whether various aspects of the hydrography data belonged to the network or not. A piece of the hydrography information belonged to the network if the TIN indicated that a flow channel was indeed associated with that given piece. To get a better idea of how this network was created, consider the following image. 66 Figure 4.8: DLG data draped over TIN to show which parts of the DLG data belong in the network This image shows a four-km 2 area within the modeling zone with the black line representing a piece of the DLG data and the background displaying a piece of the TIN. Inspection of this image reveals that the TIN shows, via valleys in the graphic, where the actual channel is located within the DLG data set. Hence, the nearly horizontal line, with no gully under it, should not be included in the stream 67 network. To determine which parts of the DLG data to include in this newest stream network, careful inspection of the entire modeling area using the TIN and the DLG data sets was performed. Only those segments of the DLG data with depressions beneath them were included in this network. When burning this stream network onto the DEM, a slightly different process was used which required two steps of burning instead of one. The stream network for this trial can be considered to be made up of two different networks--a natural network and a man-made network. While the natural network includes the tortuous aspects of the major drainage paths, the man-made network includes the smaller, artificial parts of the network such as irrigation ditches. Since the area of interest is so flat, the DEM elevations were raised 500 hundred feet above the man-made network and 1000 feet above the natural network. This way, once water flowed into the man-made network, it remained in the network until it flowed off of the DEM. Furthermore, once water flowed into the natural network, not only did it remain in the network, but it never reverted back to the man-made network. Using this network, the following watersheds were delineated. 68 Figure 4.9: Delineated system using stream network based on TIN and DLG data With this new network, CRWR-PrePro yet again delineated a set of watersheds. The table below shows how well the modeled drainage areas match up with the USGS drainage areas. Gauge Name USGS gauge area (km 2 ) Calculated gauge area (km 2 ) % Error Buffalo Bayou at Piney Point 821 792 -3.5 Addicks Reservoir 352 348 -1.1 Barker Reservoir 332 322 -3.0 Langham Creek 63.7 92.8 45 Bear Creek 55.7 42.5 -24 Buffalo Bayou near Katy, TX 164 162 -1.2 Table 4.3: Area comparisons using stream network derived from TIN 69 From these results, it appears that modeling the area drained by the Langham Creek and Bear Creek gauges may be more difficult than imagined. No matter which of the three stream networks provided the foundation for the CRWR-PrePro model, both of these gauges exhibited significant error between their actual and modeled drainage areas. However, if one considers the combined areas of these two gauges, the error equals about -12%. Nonetheless, these results indicate that for regions as flat as the area west of Houston, a thirty-meter DEM is at least partially inadequate, and more accurate terrain information, such as that from LIDAR, is needed. Regardless of the erroneous results concerning the Langham Creek and Bear Creek gauges, the stream network made from the DLG and TIN data was chosen as the best available stream network for the following reasons. Whereas the first attempt completely ignored the DLG data and the second attempt over relied on it, this attempt carefully considered each segment of the DLG data for inclusion in the stream network via direct investigation of the TIN and the DEM. Furthermore, since the Langham Creek and the Bear Creek gauges exist in the vicinity of each other, it is possible that any shortcomings resulting from the DEM would be localized to the area where these two gauges are located, especially since the errors associated with the other gauges all fall under four per cent. It is particularly interesting to note that the Langham Creek and Bear Creek gauges, with their incorrect drainage areas, both fall within the Addicks Reservoir basin, which has an overall drainage area error of about one per cent. 70 Having decided on this third network as the network of choice, it was possible to move forward within the CRWR-PrePro section of the modeling process. The next CRWR-PrePro step involved adjusting the sizes of the sub-watersheds within the region such that all basin areas were manageable. Using the Merge Sub-Watersheds function on the CRWR-PrePro menu allows the user to combine easily one sub- watershed with another. This came in handy since CRWR-PrePro created nine sub- watersheds with areas less than two km 2 . Sub-basins this small can lead to quirky results in the HMS because their lag times approach the ten-minute time step of the model. Therefore, each sub-watershed smaller than two km 2 was merged with the sub-watershed directly upstream from it. On the other end of the spectrum, CRWR-PrePro created some sub- watersheds which were exceptionally large. Since a higher number of sub-basins allows for better modeling, all sub-basins larger than thirty km 2 were split in half. To do this, seven additional outlets were placed in the same shapefile that held the gauge location information. CRWR-PrePro was then executed again from the Add Outlets menu choice so that a new set of sub-basins could be created in which no sub-basin exceeded the thirty-km 2 criterion. 71 Figure 4.10: Basin used by CRWR-PrePro to create the HMS basin model The final step before creating the HMS basin model required calculating a variety of parameters for the now-complete set of sub-watersheds using the Calculate Attributes option. These parameters addressed abstractions, lag times, sub-watershed velocities, and Muskingum routing parameters. For the abstractions calculation, CRWR-PrePro offers two choices, SCS or Initial/Constant. With a curve number grid readily available, the SCS method was chosen, and CRWR-PrePro generated average curve numbers for each of the sub-basins. Generally, the average curve number for a given basin was around 70 upstream of the reservoir exits and around 90 downstream of the reservoir exits. For calculating lag time, CRWR-PrePro offers two options-- 72 SCS and length/velocity. Again, with a curve number grid available, the SCS equation for generating lag times was used. ) *1900 )9)/1000((**100 (*6.0 5.0 7.08.0 S CNL t ? = Eqn. 4-1 where t = lag time in minutes L = length of the longest flow path in feet CN = average sub-basin curve number S = average sub-basin slope expressed as a percentage Note that the lag time = 0.6 * T c , or the time of concentration of a watershed. Calculation of the Muskingum routing parameters first required computing the sub-watershed velocities, which represent how quickly water flows through the stream channels. While each sub-watershed technically has its own velocity value, the same value was assigned to every sub-watershed. In several ways, the watershed as a whole is homogenous; slopes, curve numbers, and land usage, for example, do not vary significantly across much of the watershed. As a result, a single velocity seemed appropriate. To generate this velocity, the hydrographs of two gauges on Buffalo Bayou (Piney Point and West Belt Drive) which reside near each other were used. The period of record analyzed for this project was 1 October 1994 through 15 November 1994. By analyzing the flow data using the CORREL function in Excel, it was possible to estimate how long it took for water to travel between the gauges. 73 Correlation Between Piney Point and West Belt Drive Gauges 0.9 0.92 0.94 0.96 0.98 1 01234 Lag (hours) Corre la t i on Figure 4.11: Correlation to determine average watershed velocity By taking the flow data from each gauge and shifting them appropriately, the above graph is generated. The time when the correlation factor reaches its highest value corresponds to the time it takes water to travel from one gauge to the other. Since it is possible to measure this travel distance within ArcView, a velocity may easily be generated. In this case, the flows from the two gauges of interest best match up when they are shifted by 1.1 hours (3960 sec) compared to each other. Since the flow distance between the two gauges is 5500 m, the watershed velocity is 5500 m/3960 sec or 1.39 m/s. This value was used to calculate the Muskingum K parameter, which equals the length of the stream flow path in the sub-watershed over this velocity value. To enter sub-watershed values into CRWR-PrePro, one merely creates a text file with the grid code (i.e. sub-watershed id value), sub-watershed velocity, and sub-watershed Muskingum X value. For this model, X equaled 0.2 for all sub-watersheds since 0.2 is typically acceptable for natural streams. When 74 CRWR-PrePro is ready for this file, it prompts the user for this file?s location. To generate the number of Muskingum reaches in a stream, CRWR-PrePro uses the higher result from the following equations: 1)) 60** (**2(1 + ? = VT L XtruncateN Eqn. 4-2 1) 60***3 (2 + ? = VT L truncateN Eqn 4-3 where N1 = number of reaches using Eqn. 4-2 N2 = number of reaches using Eqn. 4-3 X = Muskingum X value = 0.2 L = length of stream flow path in meters ?T = model time step in minutes V = sub-watershed stream velocity = 1.39 m/s Completion of the Muskingum calculations marked the end of the parameter calculation section. With all parameters generated, it was now possible to create the basin model for the HMS. To do this, one invokes the HMS Schematic option on the CRWR-PrePro menu. This portion of the code uses the nodes in the stream network to identify junctions and sinks using a node classification scheme developed by Hellweger and Maidment (1997). By definition, junctions and sinks also serve as sub-basin outlets. A sub-basin and its associated outlet are linked together based on the proximity of the outlet point to the sub-basin polygon boundary. Finally, pertinent streamlines are identified as reach elements by the type of nodes at the beginning and end of each streamline. In addition, this section of the code also creates symbolic point and symbolic line attribute tables to store information related 75 to the connectivity among hydrologic features. All of this information is then written into a format which conforms to that of an HMS basin file (Olivera, 1998). 76 Chapter 5: NEXRAD Processing The preprocessing of the NEXRAD data, as described in Chapter 3.2, created text-based grids which ArcView could easily import for viewing. To get a better idea of the quantity of rain which fell over the model area during the 16-18 October 1994 event, consider the following image: Figure 5.1: Rainfall totals (inches) over the Buffalo Bayou This image indicates that rainfall values of four to nine inches fell over the rainfall area. Compared to other regions around Houston, the Addicks/Barker reservoirs avoided the full strength of this event. For instance, one does not have to travel far to the northwest of this image to encounter areas where the rainfall quantities consistently exceed twenty inches. Nonetheless, while the reservoirs may 77 not have received an incredibly large amount of rain, the storm still provides ample rainfall for this research. 5.1 Preparing the Rainfall Grids for the HMS The HMS system allows a user to include grid-based rainfall data in a model only if the data are in the Hydrologic Engineering Center Data Storage System (DSS) format. This format, created by the HEC in 1979, was designed primarily for water resource applications. The strength of this format resides in its ability to transfer data in blocks of continuous data. The DSS system considers each block of continuous data a single element, which makes searches of the data set easier and more efficient. Tom Evans of the HEC supplied a program called ai2dssGrid for converting rainfall data in ArcView?s ASCII-grid format into the DSS format. His program requires that the data be in one of two map projections, one of which is the SHG map projection. To convert from ASCII-grid to DSS, one simply types in a UNIX environment for each rainfall depth grid ai2dssGrid input dssfile pathname sdate stime edate etime gridtype where input = name of input ASCII grid containing cumulative rainfall data over one time step dssfile = name of output DSS file (e.g. rain.dss) pathname = apart/bpart/cpart/dpart/epart/fpart/ within DSS file = /bayou/buff/PRECIP///radar/ sdate = starting date of rainfall data in input ASCII grid stime = starting time of rainfall data in input ASCII grid edate = ending date of rainfall data in input ASCII grid etime = ending time of rainfall data in input ASCII grid gridtype = map projection of ASCII grid = SHG 78 As the rainfall in the October, 1994 event covered 331 ten-minute intervals, each of which has a grid as described in Chapter 3.2, a batch file with 331 calls of ai2dssGrid was used to create a DSS file named rain.dss. With this file completed, one simply needs to inform HMS of the file?s location for inclusion in a HMS model. 5.2 The ModClark Parameter File with the HMS The HMS requires a ModClark parameter file for use in conjunction with grid-based rainfall in the DSS format. This parameter file contains values for hydrologic properties of cells defined by the intersection of a set of watershed boundaries with cells in the SHG grid. The intersection of the watershed boundaries with the SHG grid necessitated that all ArcView data sets be reprojected in the SHG map projection. A discussion of the SHG grid follows. 5.2.1 The SHG Grid In an effort to standardize grid-based precipitation reporting practices, the HEC created the Standard Hydrologic Grid, an Albers map projection which covers the conterminous United States. The HEC officially recognizes resolutions of 10,000 m, 5,000 m, 2,000 m, 1,000 m, 500 m, 200 m, 100 m, 50 m, 20 m, and 10 m, although the resolution may technically be set to any value. Typically, the HEC prefers as a default a 2,000-m grid, although, in this project, a grid resolution of 997 m has been chosen to accommodate the available rainfall data (See Chapter 3.2). Each cell in the grid possesses a unique pair of integer indices (i,j) for identification. These indices represent the number of grid cells the southwest corner 79 of a given cell is from the origin of the grid, which is located at (23? N, 96? W). The formulae for determining the indices for a given cell follow: )( cellsize easting floori = Eqn. 5-1 )( cellsize northing floorj = Eqn. 5-2 where floor(x) = the largest integer less than or equal to x. As an example, consider that the HEC, located in Davis, CA, resides at (38? 35? N, 121? 45? W). This position translates to -2185019 m easting and 2063359 northing. Using the default grid resolution of 2,000 m, i = floor(-2185019/2000) = floor(-1092.5) = -1093 and j = floor(2063359/2000) = floor(1031.7) = 1031. 5.2.2 Use of GridParm to Create the ModClark Parameter File In addition to the DSS file, the HMS system also requires a ModClark parameter file which, for each sub-watershed, lists the i- and j-indices of the SHG cells within the sub-basin, the average travel distance in kilometers from each SHG cell within the sub-basin to the sub-basin outlet, and the area of each SHG cell which resides in the sub-basin in square kilometers. The following text shows a piece of the parameter file used in this model: SUBBASIN: 86 GRIDCELL: 25 761 0.13176 0.00454754 GRIDCELL: 26 761 0.25619 0.0250194 GRIDCELL: 24 762 5.25758 0.00559098 GRIDCELL: 22 765 7.74665 0.964091 GRIDCELL: 23 765 6.52171 0.994009 GRIDCELL: 24 765 5.4878 0.743663 80 GRIDCELL: 23 767 7.36279 0.165136 GRIDCELL: 24 767 7.03907 0.196063 GRIDCELL: 25 767 6.96628 0.00252854 END: SUBBASIN: 87 GRIDCELL: 28 760 4.31344 0.0163324 GRIDCELL: 29 760 3.70745 0.273056 GRIDCELL: 30 760 2.70715 0.552455 GRIDCELL: 33 766 6.40111 0.0675318 GRIDCELL: 31 767 8.31459 0.0024277 GRIDCELL: 32 767 8.04311 0.0761606 END: Notice that the largest value in the area column is 0.994009 km 2 , which is derived from the cell size of 997 meters. The HEC offers a block of AML scripts collectively called GridParm which allows a user of the HMS system to develop the ModClark parameter file using Arc/Info. A general outline of the procedure used within GridParm resides below: ? Collect DEM ? Remove pits from DEM ? Calculate flow direction, flow accumulation, and stream grids ? Delineate watersheds from terrain grids and basin outlet locations ? Combine watersheds and parameter data with grid cells to create ModClark file Since the HEC wrote these programs before Avenue for ArcView existed, GridParm generates its terrain grids from the DEM using Arc/Info AMLs. In this case, however, the first few steps may be skipped since GridParm may use the terrain grids created within the CRWR-PrePro environment. GridParm then delineates the watersheds using the terrain grids and a basin-outlet-locations file. This input file 81 lists all of the sub-basin outlet locations by latitude and longitude. For this process, GridParm was slightly modified to accept sub-basin locations as points in meters on the SHG since the sub-basin outlets locations were readily available on the SHG from the CRWR-PrePro environment. Having finished the watershed delineation, GridParm then created the ModClark parameter file. Again, a slight modification of GridParm occurred to accommodate the unorthodox 997-meter cell size of the rainfall grids. Except for the merging of the sub-watersheds with the SHG, almost all of the calculations performed by GridParm also take place within the CRWR-PrePro environment. Therefore, it seemed logical that writing a GridParm procedure for inclusion in CRWR-PrePro would be wise. Hence, the next chapter discusses the creation of a GridParm process for CRWR-PrePro using Avenue. 5.2.3 Creation of ModClark Parameter File in ArcView Using Avenue Several aspects of the GridParm code suggested that the creation of the ModClark parameter file within CRWR-PrePro would save time and effort. The most obvious of these is that both systems derive the majority of their information from the same data set?the DEM. Additionally, learning how to run GridParm in UNIX requires learning how to use an entirely new set of programs on a different platform. Furthermore, not all CRWR-PrePro users have access to a UNIX machine, so if they want to use grid-based rainfall, they will have to create the ModClark parameter file some other way. Finally, getting GridParm to generate a ModClark parameter file 82 compatible with a HMS basin model created in CRWR-PrePro requires altering some of the GridParm code. Therefore, because of all these reasons, inclusion of the GridParm procedure in CRWR-PrePro seemed natural. Because CRWR-PrePro creates a variety of data sets which do not exist in the GridParm procedure, the CRWR-PrePro version of GridParm attacks the ModClark procedure differently from the way GridParm does. Consider, for instance, how both methodologies determine the average flow length from a given SHG cell to the sub- basin outlet. In the GridParm technique, flow lengths are calculated by following the flow-direction grid to the sub-basin outlet from each point in the sub-basin. In CRWR-PrePro, however, there exists the FLDStoWO, or Flow Length DownStream to Watershed Outlet, grid. This grid contains the distance from each cell downstream to the sub-basin outlet. Therefore, to determine the average distance of all of the points in each SHG cell to the exit of the sub-basin, one needs to determine which cells are in a given SHG cell and then average these distances. Below resides the pseudocode for creating the ModClark parameter file with Avenue in the CRWR-PrePro environment. ? Ask for watershed shapefile theme ? Ask for sympoint theme, which may be used to link the name of a sub-basin in ArcView with the name of the same sub-basin in HMS. ? Ask for flow-length-downstream-to-watershed-outlet grid ? Ask for watershed grid ? Ask for SHG cell size ? Ask for name/location of output file 83 ? Take HEC id values from sympoint and assign them to watersheds in vector theme ? Determine extent of view ? Create and draw SHG shapefile over the watersheds ? Intersect SHG shapefile with watershed shapefile ? For each watershed o Note which parts of which SHG cells fall in that watershed o Convert these SHG cells into a grid o Using watershed grid and FLDStoWO grid, compute mean travel distance ? End ? Write all results to output parameter file To understand better this pseudocode, consider Figure 5.2. 84 Figure 5.2: Visual representation of ModClark code 85 The first image represents an arbitrarily delineated watershed area. The ModClark code sizes this watershed basin and then creates a grid of SHG cells which are a user-specified size. This grid is large enough to cover the area of interest. From there, the ModClark routine combines these two data sets to create a third set, which shows how the SHG cells themselves fall within each of the sub-watersheds. From here, the program may then calculate the necessary parameters for the ModClark parameter file. Once the program was completed, a button, , was placed on the ArcView toolbar which allows the user to execute the program by depressing this button. Clearly, this method provides an easier option for creation of the ModClark parameter file. In fact, besides being easier, the Avenue code outperformed its UNIX counterpart, completing the parameter file as much as 50% faster. 86 Chapter 6: Model Calibration Within HEC-HMS Once CRWR-PrePro had finished generating the basin model for the HMS, the creation of the Buffalo Bayou model moved out of the ArcView environment and solely into the HMS package. Getting the basin file from CRWR-PrePro into the HMS required using the ?Import? function from the ?Edit? menu. The image below shows what this basin model looked like upon import. Figure 6.1: Basin model as imported into HEC-HMS Before it could run the model, HMS required several additional sets of information. These included the location of the file created by the ModClark Avenue code, the initial loss of rainfall via absorption into the ground at the beginning of the event, the T c and R values for each sub-watershed for creating unit hydrographs, base 87 flow data, creation of the Precipitation model, establishment of the Control Specifications, and gauge discharge data for comparison purposes. The following section describes each of these requirements individually. 6.1 Additional Data for the HMS Model Some of the additional data for the model entailed merely a quick step or two in the HMS. For example, informing the HMS as to the location of the grid-cell file from the ModClark program meant going to the ?Basin Model Attributes? section under the ?File? menu in the basin editor and browsing to the proper file location. In addition, establishment of the Control Specifications took only a minute or two since that section merely needed to know the model start and end times. In this case, modeling began at 1520 hrs on 16 October right before the rainfall for this event began. A conclusion time of midnight on 25 October was chosen for two reasons. First, this allowed for the rainfall to run off the flat, Houston-area terrain. Next, a second rain storm, which this model does not consider, commenced in the early hours of 25 October. Additionally, this attribute section required a modeling time-step, which had already been established in PrePro as ten minutes. Finally, completing the Precipitation Model within the HMS, which uses the grid-based precipitation method, entailed specifying the location of the DSS file of gridded rainfall and listing the pathname parts in this file (A = bayou, B=buff, F = radar). Other aspects of the additional information, however, required either more time to enter, more time to consider what to enter, or both. Entering the gauge 88 discharge data, for instance, requires little more than entering the data by hand for each time step and letting the HMS know which element within the Basin model is associated with that particular gauge. Since the USGS flow data had already been processed, these data were readily available. The base flow in each sub-watershed was set to zero since this model considered a single, large event with flows that well exceeded any base flows by a substantial margin. Initial losses for runoff calculations were also set to zero since the area received enough rain in the day or two before the event to consider the soil saturated. Furthermore, as mentioned in Chapter 3.6, most of the soil in this region is in Hydrologic Group D, which means that it drains poorly. Finally, one must also consider the impervious cover of a given area when constructing a rainfall-runoff model such as this. Consider the following figure, which shows the Anderson land use-land cover for the modeled region. 89 Figure 6.2: Land use over the modeled region Except for the red block (urban), which occupies the eastern edge of the area of concern, agriculture (yellow) dominates the land use in the region. A tour of the area by the author confirmed that this is the case although some blocks of suburbia have developed in the area, particularly in the upper-central region of the model. Other small regions of forest (green) and wetlands (light blue), both with very low impervious cover, reside in the region. In fact, except for minor areas of urban development in areas outside of the eastern edge, little urbanization of the watershed exists upstream of the reservoirs. In addition, every sub-watershed upstream of the reservoirs possesses at minimum a majority of non-urban use. Therefore, for the practical purposes of this model, especially for parameter estimation, one may assume 90 that a level of about five per cent imperviousness is acceptable for the region upstream of the reservoirs. For the region downstream of the reservoirs, 85% imperviousness levels may be used. This selection is based on information from the United States Department of Agriculture TR-55 manual on urban hydrology, which states that commercial areas average 85% imperviousness. The HMS requires one final set of information to complete the model. The two primary parameters associated with the creation of a unit hydrograph using the ModClark method, time of concentration (T c ) and storage (R), had to be entered for each of the sub-basins. Calculating these parameters involved trying several different techniques before the most acceptable solution made itself apparent. Originally, data taken from a study of the Buffalo Bayou region completed in a 1977 USACE- Galveston study were used. This study listed for each of its sub-basins the area, T c , and R. T c and R values were calculated using many storms over several years. These data provided a starting point for determining modern T c and R values. Using the data from this 1977 study, a chart of T c vs. area shows a linear relationship between these two sub-basin parameters. 91 Tc vs Area in Buffalo Bayou Region 0.5 0.75 1 1.25 1.5 1.75 2 0 5 10 15 20 25 30 Area (km 2 ) T c (h r) Actual Data Best-fit Line Figure 6.3: T c vs. area for the modeling region A perusal of these data revealed that the R value of a given sub-basin approximately equaled twice the T c value for the same basin. In fact, the average T c value of all of the sub-basins equals 1.10 hrs while the R values average 2.14 hrs; therefore, an initial guess of the R values for a given sub-basin of two times the T c value seemed acceptable. Unfortunately, reservoir release data could be located for this study. Many efforts to obtain these data from the USACE-Galveston office resulted in nothing more than a computer program which calculated exit flows from the reservoirs if one knew how many of the release gates had been opened and to what heights these gates had been raised. Without these crucial data, it was impossible to accurately predict the flow moving downstream of the reservoirs. Since the version of HMS used in this 92 study can only work with uncontrolled reservoirs, this made modeling anything downstream of the reservoirs next to impossible. More in-depth studies in the future would require adequate reservoir release schedules to generate acceptable results downstream of the reservoirs. 6.2 Execution of the HMS Model Initial runs using these parameters were not encouraging. Consider the hydrograph below for the Buffalo Bayou gauge near Katy, TX, which was chosen since the modeled and actual drainage areas for this gauge compare the most favorably of all of the gauges considered. Figure 6.4: Hydrograph of Buffalo Bayou near Katy, TX gauge using original parameters 93 In this image, as with all hydrographs displayed in this thesis, the red line represents the observed runoff while the blue line represents the modeled flow. Clearly, a limitation exists in the parameters used for this model, especially since other gauges in the area exhibited similar disparities between the two curves. Inspection of this plot reveals that while the magnitudes of the two sets of flows differ substantially, the timing of the hydrograph peaks matches well. For instance, the calculated peak flow occurs at 2350 hrs on 18 October, which compares favorably with the observed time of 2300 hrs on 18 October. In addition, the observed volume equals 3.6 inches while the calculated volume equals 3.2 inches. This difference is well within the range of potential error for NEXRAD rainfall data. Based on these reasons, the inaccuracies in the above image probably stem from incorrect R values. Higher R values, which are related to longer storage times in a watershed, would attenuate the calculated hydrograph. In an attempt to construct a better match between observed and calculated flows, a series of equations was used to generate new T c and R values for the watershed. Using an equation from The Handbook of Hydrology (1993), T c was calculated according to 5.0 14.2 * * 3 2 S nL Tc = Eqn. 6-1 where T c = time of concentration in minutes L = length of the longest open channel in feet 94 n = Manning?s coefficient S = dimensional slope Since CRWR-PrePro calculated L and S earlier in the modeling process, Manning?s coefficient remained the sole unknown variable for this equation. Considering that straight streams have a coefficient of 0.03 and winding streams have a coefficient of 0.04, 0.035 appeared to be a valid estimation of the streams in the area. Calculation of the R values for each of the sub-basins employed an equation developed in Hydrology and Floodplain Analysis (Bedient and Huber, 1992) which had been developed for watersheds in the Houston area. 706.0 )(* S L CRTc =+ Eqn. 6-2 where T c = time of concentration in hours R = R in hours C = 7.25 if per cent development is under eighteen per cent L = length of longest open channel in miles S = slope in feet/mile Using these new equations to determine the T c and R values for the Buffalo Bayou at Katy, TX gauge resulted in the following hydrograph. 95 Figure 6.5: Hydrograph of Buffalo Bayou near Katy, TX using T c and R equations While this is an improvement, problems persist. As before, the peaks of the curves match up well with a calculated peak time of 0100 hrs on 19 October. Note that this is two hours later than the peak using the previous set of T c and R values. Furthermore, while the peaks better resemble each other in this scenario, the R values remain too low for a good match to be achieved. In an effort to match better the two curves, the R values for the Katy gauge were systematically doubled and then tripled. While this is a crude technique for generating a more accurate calculated hydrograph, the lack of an additional technique for determining R values necessitated this procedure. Furthermore, one could 96 postulate that varying the R values is an acceptable form of model calibration since its technique of generation is the among the defensible of all of the parameters. Using R values which equaled three times the R values generated by the Bedient R equation yielded the following hydrograph. Figure 6.6: Hydrograph of Buffalo Bayou near Katy, TX using 3*Bedient R values Using these parameters creates a good match between the two hydrographs. In addition to the peak times differing by two hours, the peak flows themselves vary by thirty cfs (2350 cfs observed vs. 2380 cfs calculated). A downfall to these new results, however, is the volumetric difference between the two curves with the 97 observed curve still containing 3.6 inches of runoff while the calculated curve moves 3.1 inches of runoff. This is most likely attributed to errors in the NEXRAD data the abstraction losses were set to zero based on assumed antecedent conditions. One final note about this hydrograph is that it is disconcerting that an arbitrary technique was used to adjust the R values so that the peaks would better agree. At this point, considering the match achieved at the Katy gauge, the Langham Creek and the Bear Creek gauges should each be acknowledged. Using parameters generated similarly to those created for the Katy gauge (i.e. multiplying R by a factor of three), HMS calculated the following hydrographs for the Bear Creek and Langham Creek gauges respectively. 98 Figure 6.7: Hydrograph of the Bear Creek gauge using 3*Bedient R values Without a doubt, a fundamental error in the model prohibits a good match for this gauge. While Figure 6.7 shows the results when R is multiplied by three, other multipliers were tried as well, and each one resulted in hydrographs which did not agree very well with each other. In addition, the effective rainfall for this gauge equals 8.7 inches of runoff, while the calculated value totaled 3.0 inches for an error of -66%. This differs substantially from the -24% error between the observed (55.7 km 2 ) and calculated (42.5 km 2 ) areas. Along the same vein, consider the Langham Creek hydrograph. 99 Figure 6.8: Hydrograph of the Langham Creek gauge using 3*Bedient R values Here again significant disparities occur between the observed and calculated values for this plot. For instance, the peak flow, while similar in magnitude for each hydrograph, varies by over a day. In addition, the shapes themselves do not agree well. Interestingly, however, the calculated runoff volume of 3.8 inches agrees favorably with the observed 3.9 inches. This could be attributed to levels of development that have occurred in this watershed according to the aerial photography of the region, though the likely culprit is the error between the modeled and official drainage areas. In addition, as with the Bear gauge, other R multipliers were tried. 100 This indicates that while multiplying the R values for the Katy gauge by three provided a good match, it appears that that multiplier is unique to that watershed. In conclusion, consider the following table, which summarizes the results at the three flow gauges addressed. Gauge Obs Peak Time Calc Peak Time Obs Runoff (in) Calc Runoff (in) Obs Peak Flow (cfs) Calc Peak Flow (cfs) Obs Drainag e Area (km 2 ) Calc Drainag e Area (km 2 ) Katy 18 th @ 23:00 19 th @ 01:00 3.6 3.1 2350 2380 164 162 Bear 20 th @ 01:00 19 th @ 02:00 8.7 3.0 760 540 55.7 42.5 Langham 18 th @ 05:00 19 th @ 02:00 3.9 3.8 1540 1510 63.7 92.8 Table 6.1: Summary of data at flow gauges in the HMS model While the model successfully duplicated the observed results at the Katy gauge, significant problems exist at the other two gauges. In response to the problems present at the Bear and Langham gauges, some potential reasons for the results are presented in the next chapter. 101 Chapter 7: Conclusions and Future Work The application of CRWR-PrePro in conjunction with the HMS to create a rainfall-runoff model of the Addicks/Barker Reservoir area provided a unique experience and generated valuable insights, particularly into modeling an area as flat as the terrain around western Houston. The rest of this chapter addresses conclusions of the results generated from the CRWR-PrePro/HMS tandem; in addition, some comments on potential future work are presented. 7.1 Modeling Conclusions The peculiarities of the results of the rainfall-runoff model require some additional discussion. On one hand, the Katy gauge exhibits acceptable matches in hydrograph peak magnitudes, peak times, general curvature, total runoff volume, and drainage area after an arbitrary adjustment of the R values. On the other, the Bear and Langham gauges possess glaring discrepancies in many of these parameters. Clearly, the differences between the actual and calculated drainage areas for the problematic gauges play a primary role in producing such conflicting results. Apparently, the delineation of the watersheds themselves must have played a part in fomenting the errors. However, a point to ponder is why the HMS had the most success modeling the Katy gauge. A possible explanation follows. The aforementioned 1977 report performed on the Addicks and Barker Reservoir region noted that on occasion flow from Cypress Creek overflows into the northern portion of the Addicks Reservoir watershed. As part of this observation, the report provided 102 a series of rating tables which, for subdivided portions along the Addicks/Cypress watersheds, indicated how much flow Cypress Creek could handle before overflow occurred. In fact, consider the following map, which provides a close up of the region pertinent to this discussion. Figure 7.1: Distribution of rainfall depths in inches The 1977 report mentioned that Cypress Creek at Katy-Hockley could sustain flows of at least 9500 cfs before overflowing. Moreover, the report also determined that no water would flow into the Addicks Reservoir at the Cypress Creek at House- Hahl gauge location. In fact, the minimum flow which would cause Cypress Creek to flood into the Addicks reservoir at any location was 4500 cfs. According to the flow data supplied by the USGS for this project, the maximum flow at either location on Cypress Creek was 2400 cfs. Interestingly enough, however, oral communication with James Doan at the HEC indicated that eyewitnesses did report overflow from the Cypress Creek watershed to the Addicks watershed during this particular rain event. 103 In addition, oral communication with Fred Liscum of the USGS, Houston Office revealed that the area in this particular region of the study zone is so flat that the direction of the wind can sometimes cause water to flow between the two watersheds. The Katy sub-basin, however, resides far enough away from the Cypress Creek such that no crossover flow occurs. Some other differences exist between the drainage area of the Katy gauge and the drainage areas of the Bear and Langham gauges. The average slope as calculated in CRWR-PrePro equals 0.00091 for the entire modeling region, 0.00082 for the Katy gauge, 0.00056 for the Langham gauge, and 0.00078 for the Bear gauge. In addition, the error between actual and calculated drainage areas for the entire region was - 3.5%, 1.2% for the Katy gauge, -24% for the Bear gauge, and 45% for the Langham gauge. Thus, in this case, a sub-basin?s average slope directly relates to how accurately the area of a delineated watershed will be. One should realize, however, that while the Katy and Bear regions had similar overall slopes, significant differences in their errors appeared. Overall, while this model does not establish a definitive cutoff for modeling with thirty-meter DEMs, one should be leery of using DEMs at this resolution if the slope is under 0.0008. If there exists a critical slope that sends a warning flag when dealing with thirty-meter DEMs then perhaps another critical parameter may be established as well. In this case, one could consider the average stream length per a given area in a potential modeling area. Using the ModClark parameter file, one may perform such calculations and determine if there indeed exists a cutoff value. Clearly, the less the 104 average stream length per area, the straighter the stream. Calculations were performed to determine this value for the entire watershed and for each of the three gauges considered earlier. The following table holds the results from these calculations. Location of Interest Location Area (km 2 ) Average Stream Length (km/km 2 ) All of modeling area 883 5.86 Bear gauge 42.5 5.55 Langham gauge 92.8 5.63 Katy gauge 162 4.68 Table 7.1: Average stream length per unit area As the table indicates, the area which owned the best modeling success, namely the Katy region, exhibits a substantially lower stream length per square kilometer. Apparently, an additional reason for the success in the Katy area resides within this parameter. While this chart does not give a definitive cutoff value, one should be leery of stream lengths per area larger than five km/km 2 when using thirty- meter DEMs. Finally, one should take into account that finer-resolution DEMs should generate higher stream lengths per area since a stream would have a greater opportunity to meander in such a grid, so caution should be used when considering this parameter for other grid resolutions. Finally a brief note about the NEXRAD rainfall data. While the NEXRAD data more than likely contributed somewhat to the differences between observed and calculated results, significant discussion was not spent on why this may be the case 105 since it appears that the resolution issues in the DEM data dominated the error propagation. The poor resolution of these data, in retrospect, indicates that they were not suitable for use in this study. However, in the case of the NEXRAD data, there is no reason to believe that were not applicable to this project, especially considering that the data had been quality controlled with local rainfall information at the Princeton Environmental Institute before they were made available. 7.2 Recommendations for Future Work Essentially, the combination of low slopes within a thirty-meter DEM worked together to make modeling the Buffalo Bayou region difficult. Future efforts to develop more reliable models could employ several different tools. First and most obvious, terrain models of higher resolution would prove invaluable. While Texas has recently completed thirty-meter DEMs of the entire state, which many thought would be accurate to model even exceptionally flat regions such as the southeastern Texas coast, finer resolution terrain maps have begun to appear. The Houston Advanced Research Center (HARC), located in The Woodlands, TX, has generated DEMs on the order of ten-foot resolution using airborne light detection and ranging (LIDAR) technology. By measuring the time it takes for a pulse of light to travel from an airplane to the ground and back, HARC has developed LIDAR-based maps of all of Harris County. While certain issues persist in the use of this technology, such as filtering out noise and management of such large quantities of data, these data will undoubtedly assist greatly in future modeling attempts. 106 Future improvements could also be achieved via the use of high-resolution DOQQs. These aerial photographs could provide excellent insight into the actual stream network in a given region. Used in conjunction with terrain data, one could generate a high-quality stream network by burning in streams based on information gleaned from DOQQs data sets. Furthermore, aerial photography provide more modern land use information than the Anderson data sets do. A more general suggestion for improving the rainfall-runoff models involves better parameter data sets. For example, although the curve number varied only slightly over most of the region considered for this study, those data were nonetheless derived from 500-meter grids. Completion of a national SSURGO data set would allow the generation of higher-quality curve number grids. In addition, better techniques for determining abstractions would assist in modeling. Convenient techniques for developing Green & Ampt parameters would prove useful. Fortunately, for this model of a single event, a rainstorm in the days before established level III antecedent moisture conditions. At best, this was a fortuitous situation. Finally, if one employs the ModClark method, a more accurate technique for calculating T c and R values needs to be developed. In fact, perhaps the ModClark parameter file could include T c and R values for each of the SHG cells within a given sub-basin. As far as parameter determination is concerned, this project probably would have benefited most greatly from better equations for calculating T c and R. 107 APPENDIX A COMPUTER CODE OF INTEREST This Appendix contains the various blocks of computer used for this thesis. Avenue program defgage -- This program converts a set of points from a geographic map projection to an Albers map projection. User sets parameters of Albers map projection. ?-------------------------------------------------------- ------ ?-------------------------------------------------------- ------ ? Name: swbp.definegages ? Written by Seann Reed ? Headline: ? Self: ? Returns: ? Description: Create a point shape file from ? locations specified ? in a table. ? ? Topics: ? Search Keys: ? Requires: ? History: ? ?-------------------------------------------------------- --?------------------------------------------------------ ---- theProject=av.GetProject theView=av.GetActiveDoc theDocs=theProject.GetDocs tabList=List.Make for each d in theDocs if (d.Is(Table)) then tabList.Add(d.getname) end end ?--- IDENTIFY INPUT TABLE intablename=msgbox.choiceasstring(tabList,"Choose table with Lat/Lon Values","Table") if (intablename=nil) then 108 exit end intable=theproject.finddoc(intablename) invtab=intable.getvtab infields=invtab.getfields ?--- IDENTIFY INPUT FIELDS latfield=msgbox.choiceasstring(inFields,"Choose the latitude field.","Latitude") if (latfield=nil) then exit end lonfield=msgbox.choiceasstring(inFields,"Choose the longitude field.","Longitude") if (lonfield=nil) then exit end idfield=msgbox.choiceasstring(inFields,"Choose id field.","ID") if (idfield=nil) then exit end ?--- ?define output map projection if desired ?--- projpoints=true labels={"central meridian","lower standard parallel","upper standard parallel","reference latitude","false easting","false northing","spheroid"} defaults={"- 100","27.4167","34.9167","31.1667","1000000","1000000","# spheroid_grs80"} inlist=msgbox.multiinput("Input Albers parameters. Cancel to get ouput in geographic coordinates.","Projection Parameters",labels,defaults) if (inlist.count=0) then projpoints=false else albprj=albers.make(rect.makenull) albprj.setcentralmeridian(inlist.get(0).asnumber) albprj.setlowerstandardparallel(inlist.get(1).asnumber) albprj.setupperstandardparallel(inlist.get(2).asnumber) albprj.setreferencelatitude(inlist.get(3).asnumber) albprj.setfalseeasting(inlist.get(4).asnumber) 109 albprj.setfalsenorthing(inlist.get(5).asnumber) albprj.setspheroid(inlist.get(6).asenum) end ?--- READ AND PROCESS DATA OutFileName=FileDialog.Put("Albstat".asfilename,"*.shp"," Output Shape File" ) if(OutFileName=Nil)then exit end OutFileName.SetExtension("shp") OutFTab=FTab.MakeNew(OutFileName,point) outTheme=Ftheme.make(outftab) ?CREATE FIELDS FOR THE NEW POINT TABLE outFields=List.Make outfields.Add(Field.Make("Albstat#",#field_long,9,0)) outFields.Add(Field.Make("Latitude",#field_decimal,8,6)) outFields.Add(Field.Make("Longitude",#field_decimal,8,6)) outFieldsc=outFields.DeepClone outftab.addfields(outfieldsc) theView.addtheme(outTheme) if(outFtab.CanEdit)then outFtab.SetEditable(true) else msgbox.info("Can?t edit the output theme.","Error") exit end ?IDENTIFY FIELDS FOR WRITING shpField=outFtab.findfield("shape") oidfield=outftab.findfield("Albstat#") olatfield=outftab.findfield("Latitude") olonfield=outftab.findfield("Longitude") for each rec in invtab id=invtab.returnvalue(idfield,rec) lat=invtab.returnvalue(latfield,rec) if ((id>0) and (lat>0)) then ?lat=invtab.returnvalue(latfield,rec).asstring ?lon=invtab.returnvalue(lonfield,rec).asstring lat=invtab.returnvalue(latfield,rec) lon=invtab.returnvalue(lonfield,rec) ?latdeg = lat.left(2).asnumber + (lat.middle(2,2).asnumber/60) + (lat.right(2).asnumber/3600) 110 ?londeg = lon.left(2).asnumber + (lon.middle(2,2).asnumber/60) + (lon.right(2).asnumber/3600) newrec=outFtab.AddRecord ?pt=point.make(londeg*(-1),latdeg) pt=point.make(lon,lat) if (projpoints) then ptp=pt.returnprojected(albprj) outFtab.Setvalue(shpField,newrec,ptp) else outFtab.Setvalue(shpField,newrec,pt) end outftab.setvalue(oidfield,newrec,id) outftab.setvalue(olatfield,newrec,lat) outftab.setvalue(olonfield,newrec,lon) end end outftab.seteditable(false) Avenue program txdot.modclark -- This program creates the ModClark parameter file for use with grid-based rainfall in the HEC- HMS ?******************************************************** ************************** ?Name: modclark ? ?Description: Avenue version of the Hydrologic Engineering ? Center?s GridParm program, which creates the ? ModClark parameter file to be used in HEC- HMS. ? ?Author: Seth Ahrens and Ximing Cai ? ?History: Modified for CRWRPre-Pro by Brian Adams ? Modified by Francisco Olivera for CRWR-PrePro (prepro04.apr) on 5/17/99 111 ?******************************************************** ************************** ? ?get the view and working directory ? ? theview=av.GetActiveDoc thedisplay = theview.getdisplay theDir=av.GetProject.GetWorkDir ? ask user for watershed shapefile theme; make sure choice is plausible iptheme = MsgBox.List(theview.GetThemes,"Select ?Watershed? polygon theme","ModClark") theftab = iptheme.getftab theshapef = theftab.findfield("shape") theshape = theftab.returnvalue(theshapef,0) if (theshape.getclass.getclassname = "polygon") then ipfound = true end if (not ipfound) then msgbox.error("Watershed theme has to be a polygon theme.", "ModClark") exit end ? ask user to choose appropriate sym point shapefile ? this theme contains the hec id numbers which will be assigned to the ? watershed theme selected above. also make sure user chooses valid file. pttheme=MsgBox.List(theview.GetThemes,"Select ?SymPoint? point theme","ModClark") ptftab = pttheme.getftab pthecf = ptftab.FindField("hecid") if (pthecf = nil) then 112 MsgBox.Warning("SymPoint theme does not have?hecid? field -- Exiting Program","ModClark") exit end ? get flow downstream to outlet grid flowdist=MsgBox.List(theview.GetThemes,"Select ?Flow Length to Watershed Outlet? grid theme","ModClark") flowgrid = flowdist.GetGrid ? get watershed grid wtshdthm=MsgBox.List(theview.GetThemes,"Select ?Watershed? grid theme","ModClark") watgrid = wtshdthm.GetGrid aFN=av.GetProject.GetWorkDir.MakeTmp("wshgrd", "") ? rename data set watgrid.Rename(aFN) ?ask user for the cell size. cellsize=Msgbox.input("Enter cell size in meters for the precipitation mesh.", "ModClark", "2000").AsNumber ? ask user for name of output file cellfname=Msgbox.input("Enter the path and nameof the grid-parameter output file", "ModClark-- OUTPUT FILE", thedir.asstring).AsFileName cellfile= linefile.make(cellfname, #FILE_PERM_WRITE) ? take original watershed shapefile theme. find all fields and ? add hec id values from sym point shapefile if not done already. wshed=iptheme wshdftab=wshed.GetFtab 113 wshpf=wshdftab.Findfield("Shape") widf=wshdftab.Findfield("Id") wgridcodef=wshdftab.Findfield("Gridcode") wlngflowf=wshdftab.Findfield("Lngflwpth") wslopef=wshdftab.Findfield("Slope") wshvelf=wshdftab.Findfield("Wshvel") wlagtimef=wshdftab.Findfield("Lagtime") wshdftab.seteditable(true) hecidf = wshdftab.FindField("hecid") if (hecidf = nil) then hecidf = field.make("hecid",#FIELD_DECIMAL,16,0) wshdftab.addfields({hecidf}) hecidf = wshdftab.FindField("hecid") end ?----take hec id values from sym point shapefile. add these ?----values to the watershed ftable. ptgridf = ptftab.FindField("gridcode") for each rec in ptftab ptid = ptftab.ReturnValue(pthecf,rec) ptgrid = ptftab.ReturnValue(ptgridf,rec) for each x in wshdftab wgridcode = wshdftab.ReturnValue(wgridcodef,x) if (wgridcode = ptgrid) then wshdftab.setvalue(hecidf,x,ptid) end end end wshdftab.seteditable(false) ? find total area of watersheds. tarea=0 for each w in wshdftab 114 wshp=wshdftab.ReturnValue(wshpf, w) tarea=tarea+wshp.ReturnArea end ? ? set area tolerance, a ratio of the total area. if a watershed ? is smaller than this area, then the script ignores it. atol=0.000001 ?------------create precipitation mesh shapefile.-------- ------ ? first calculate the range for the mesh based on the extent of the watershed polygon theme arect=wshed.ReturnExtent theheight=arect.GetHeight thewidth=arect.GetWidth xcenter=arect.ReturnCenter.Getx ycenter=arect.ReturnCenter.Gety xminf=xcenter-(thewidth/2.0) yminf=ycenter-(theheight/2.0) xmaxf=xcenter+(thewidth/2.0) ymaxf=ycenter+(theheight/2.0) ?calculate grid boundaries in even cells cmin=(xminf/cellsize).Truncate rmin=(yminf/cellsize).Truncate xmin=cellsize*cmin - cellsize ymin=cellsize*rmin - cellsize xmax=cellsize* ( ( xmaxf/cellsize).Truncate ) + cellsize ymax=cellsize* ( ( ymaxf/cellsize).Truncate ) + cellsize ? calculate grid size in rows ncol=(xmax-xmin)/cellsize nrow=(ymax-ymin)/cellsize 115 ? generate a polygon shapefile of the shg cells. av.showmsg("generate a polygon shape of the cells") polycell=FN.merge(thedir.asstring,"shg") pcftab=FTab.MakeNew(polycell, polygon) pcshpf=pcftab.FindField("shape") pcftab.Seteditable(true) idfld = field.make("cell_id", #FIELD_DECIMAL, 16, 0) areafld=field.make("area", #FIELD_DECIMAL, 36, 6) xfld=field.make("shg_x", #FIELD_DECIMAL, 24, 0) yfld=field.make("shg_y", #FIELD_DECIMAL, 24, 0) pcftab.addfields( {idfld, areafld, xfld, yfld} ) id=0 for each r in 0 .. (nrow-1) keepgoing = av.setstatus((r/nrow)*100) row=ymin+ (r*cellsize) col=xmin numrow=rmin+r numcol=cmin for each c in 0 .. (ncol-1) pt1=point.make(col, row) pt2=point.make(col,(row+cellsize) ) pt3=point.make((col+cellsize),(row+cellsize) ) pt4=point.make((col+cellsize), row) theshape=polygon.make( { {pt1, pt2, pt3, pt4} } ) therec = pcftab.addrecord id=id+1 pcftab.setvalue(pcshpf, therec, theshape) pcftab.setvalue(idfld, therec, id) pcftab.setvalue(areafld, therec, theshape.ReturnArea) pcftab.setvalue(xfld, therec, numcol) pcftab.setvalue(yfld, therec, numrow) numcol=numcol+1 col=col+cellsize end 116 end pcthm=Ftheme.Make(pcftab) theview.AddTheme(pcthm) pcthm.Setvisible(true) pcftab.seteditable(false) ? ? create a shapefile of shg cells dropped over watershed theme. ? overlayfn=FN.merge(thedir.asstring,"wshshg") wgftab=FTab.MakeNew(overlayfn, polygon) wgshpf=wgftab.FindField("shape") wgftab.Seteditable(true) wgshpf=wgftab.Findfield("Shape") wgidf=field.make("wshd_id", #FIELD_DECIMAL, 20, 0) cgidf=field.make("cellid", #FIELD_DECIMAL, 20, 0) wggridf=field.make("gridcode", #FIELD_DECIMAL, 20, 0) wglngflowf=field.make("lngflwpth", #FIELD_DECIMAL,20, 4) wgslopef=field.make("slope", #FIELD_DECIMAL, 16, 4) wgshvelf=field.make("wshvel", #FIELD_DECIMAL, 36, 4) wglagtf=field.make("Lagtime", #FIELD_DECIMAL, 36, 4) wgareaf=field.make("area", #FIELD_DECIMAL, 36, 6) wshgxf=field.make("shgx",#FIELD_DECIMAL,10,0) wshgyf=field.make("shgy",#FIELD_DECIMAL,10,0) wgftab.addfields({ wgidf, cgidf, wggridf, wglngflowf, wgslopef, wgshvelf, wglagtf, wgareaf,wshgxf,wshgyf} ) ? actually drop shg cells over watersheds in this block. 117 av.showmsg("Overlay the cell shape on the watershed shape") rec=0 for each c in pcftab keepgoing = av.setstatus((rec/pcftab.GetNumRecords)*100) cshp=pcftab.ReturnValue(pcshpf, c) cid=pcftab.ReturnValue(idfld, c) xval=pcftab.ReturnValue(xfld,c) yval=pcftab.ReturnValue(yfld,c) for each p in wshdftab wshp=wshdftab.ReturnValue(wshpf, p) wid=wshdftab.ReturnValue(widf, p) wgridcode=wshdftab.ReturnValue(wgridcodef, p) wlngflw=wshdftab.ReturnValue(wlngflowf, p) wslope=wshdftab.ReturnValue(wslopef, p) wshvel=wshdftab.ReturnValue(wshvelf, p) wlagtime=wshdftab.ReturnValue(wlagtimef, p) wgshp=cshp.ReturnIntersection(wshp) if(wgshp.ReturnArea> (atol*tarea) ) then therec = wgftab.addrecord wgftab.SetValue(wgshpf,therec,wgshp) wgftab.SetValue(wgidf, therec, wid) wgftab.SetValue(cgidf, therec, cid) wgftab.SetValue(wggridf, therec, wgridcode) wgftab.SetValue(wglngflowf, therec, wlngflw) wgftab.SetValue(wgslopef, therec, wslope) wgftab.SetValue(wgshvelf, therec, wshvel) wgftab.SetValue(wglagtf, therec, wlagtime) wgftab.SetValue(wgareaf, therec, wgshp.ReturnArea) wgftab.SetValue(wshgxf,therec,xval) wgftab.SetValue(wshgyf,therec,yval) end end rec=rec+1 end 118 wgthm=Ftheme.Make(wgftab) theview.AddTheme(wgthm) wgthm.Setvisible(true) wgftab.seteditable(false) av.showmsg("write parameter output file") ? for each watershed in the original watershed shapefile theme: ? - from overlay, take portions of shg cells which fall in that watershed. ? - these extracted shg cells will be in their own shapefile. ? - convert this shapefile to a grid ? - for each cell id, use a zonalstat to determine mean travel distance. ? - use zonalstatstable to link each cell id with its proper travel distance ? - add these data to ftable associated with shg cells overlapped onto watershed shapefile. ? - output data for the given watershed to the modclark file for hms. ? num = 0 wgbitmap=wgftab.getselection tempVtabFN=FN.merge(thedir.asstring,"tmpstat.dbf") for each rec in wshdftab ?pull out the shg cells associated with one watershed. ?use bitmap to select cells and then convert directly to a grid ? this stores the shg cellid/gridnumber/shape of each ? shg cell in the watershed of interest shedgrid = wshdftab.ReturnValue(wgridcodef,rec) 119 wgbitmap.clearall wgftab.updateselection for each x in wgftab gridnum = wgftab.ReturnValue(wggridf,x) if (gridnum = shedgrid) then wgbitmap.set(x) wgftab.updateselection end end ? convert celltemp to a grid using the ftable cellidgrid = grid.makefromftab(wgftab,prj.makenull,wgftab.findfield("c ellid"),nil) afn=FN.merge(thedir.asstring,"subwsh.dbf") theftab = wgftab.export(afn,shape,TRUE) theftab.seteditable(true) aFN=av.GetProject.GetWorkDir.MakeTmp("cellid", "") ? rename data set cellidgrid.Rename(aFN) outtheme=gtheme.make(cellidgrid) outtheme.setname("cellid") ? for each cellid in the watershed, compute mean travel distance zonefield = theftab.findfield("cellid") distgrid = flowgrid.zonalstats(#GRID_STATYPE_MEAN,theftab,Prj.MakeNu ll,zonefield,False) trunkdistgrid=distgrid.int trunkdistgrid=trunkdistgrid.combine({cellidgrid,watgrid}) 120 outtheme=gtheme.make(trunkdistgrid) tempvtab=trunkdistgrid.getvtab table.make(tempvtab) fieldlist=tempvtab.getfields ? add mean cell travel distances to theftab cellvaluef = fieldlist.get(3) wcellf = theftab.findfield("cellid") wcellgf = theftab.findfield("gridcode") thedistf = fieldlist.get(2) dfield = theftab.findfield("cell dist") theftab.join(wcellf,tempvtab,cellvaluef) outtheme=ftheme.make(theftab.clone) ? theview.addtheme(outtheme) ? add a field to the overlay theme for storing the mean cell travel dist. keepgoing = av.setstatus(( (num+2)/wshdftab.GetNumRecords)*100) ? hec id from original whecidf = wshdftab.Findfield("hecid") whecid=wshdftab.ReturnValue(whecidf, rec) ? output data to modclark file starting here outstring="SUBBASIN: "+whecid.AsString cellfile.writeelt(outstring) fieldlist=theftab.getfields shgxf = fieldlist.get(9) shgyf = fieldlist.get(10) distf = fieldlist.get(13) areaf = fieldlist.get(8) 121 for each g in theftab dist = theftab.returnvalue(distf,g) if (dist.isnull) then continue end shg_x = theftab.returnvalue(shgxf,g) shg_y = theftab.returnvalue(shgyf,g) area = theftab.returnvalue(areaf,g) str1="GRIDCELL: " str2=shg_x.asstring str3=shg_y.asstring str4=(dist/1000).asstring str5=(area/1000000).asstring outstring=str1+ str2+" "+str3+" "+str4 +" "+str5 cellfile.writeelt(outstring) end outstring="END:" cellfile.writeelt(outstring) num=num+1 theftab.unjoinall end wgftab.seteditable(false) cellfile.close ? ? end of the program ? 122 MATLAB program adjuster.m -- This program converts adjusted gauge elevations to flow values using data available in USGS rating tables. % load all height data for the gauges % heights has 8 colums load heights -ascii % load all data from the usgs rating tables % each rating has two colums % first column is for elevation % second column is for flow load ratings1 -ascii % lenheight is number of flow values to be calculated lenheight = length(heights(:,1)); for x = 1:lenheight value = heights(x,1); notfound = 1; pos = 1; while notfound > 0 if value > ratings(pos,1) pos = pos + 1; else notfound = -1; end end interp = (value - ratings(pos-1,1) ) / 0.01; flow(x) = ratings(pos-1,2) + (ratings(pos,2) - ratings(pos-1,2))*interp; end flow = flow?; 123 MATLAB program convert.m ? This program is the exact code used to analyze the first batch of data from the first day of the storm. To analyze other pieces of the storm, the same general code was used, but the data sets were different. clear all; % load data load zero.dat; load aout001.dat; % aout files represent original data sets. load aout002.dat; % each aout file is at a different time. load aout003.dat; load aout004.dat; load aout005.dat; load aout006.dat; load aout007.dat; load aout008.dat; load aout009.dat; load aout010.dat; load aout011.dat; load aout012.dat; load aout013.dat; load aout014.dat; load aout015.dat; load aout016.dat; load aout017.dat; load aout018.dat; load aout019.dat; load aout020.dat; load aout021.dat; load aout022.dat; load aout023.dat; load aout024.dat; load aout025.dat; load aout026.dat; load aout027.dat; load aout028.dat; 124 load aout029.dat; load times.txt; % load file which contains times associated with % all of the aout files. times = times(1:30) %only need first thirty times since only first %data sets are considered for this run. times = times?; % grab intensity values from each file aout001=aout001(:,3); aout002=aout002(:,3); aout003=aout003(:,3); aout004=aout004(:,3); aout005=aout005(:,3); aout006=aout006(:,3); aout007=aout007(:,3); aout008=aout008(:,3); aout009=aout009(:,3); aout010=aout010(:,3); aout011=aout011(:,3); aout012=aout012(:,3); aout013=aout013(:,3); aout014=aout014(:,3); aout015=aout015(:,3); aout016=aout016(:,3); aout017=aout017(:,3); aout018=aout018(:,3); aout019=aout019(:,3); aout020=aout020(:,3); aout021=aout021(:,3); aout022=aout022(:,3); aout023=aout023(:,3); aout024=aout024(:,3); aout025=aout025(:,3); aout026=aout026(:,3); aout027=aout027(:,3); aout028=aout028(:,3); aout029=aout029(:,3); 125 % make grid of all rainfall grids % keep time data separate % the zero data set was entered because inspection of the data % revealed that at some point during the storm, no rain was % falling. hence, to properly block the data into constant % ten minute intervals, these zero data sets were required. newgrid = [aout001,aout002,aout003,aout004,aout005,aout006,aout007]; newgrid = [newgrid,aout008,aout009,aout010]; newgrid = [newgrid,aout011,aout012,aout013,aout014,aout015,aout016]; newgrid = [newgrid,aout017,aout018,aout019,aout020]; newgrid = [newgrid,aout021,aout022,aout023,aout024,aout025,zero]; newgrid = [newgrid,aout026,aout027,aout028,aout028]; % take time array and convert into % decimal time in minutes from the % start of the storm % first determine nearest ten-minte value % before storm began % i.e. 21:42:39 yields 1340 since 9.40 pm % is 1340 min past midnight rawtime = times(1); hours = floor(rawtime/10000); mins = floor((rawtime - hours*10000) / 100); secs = rawtime - hours*10000 - mins*100; inittime = hours*60 + floor(mins/10)*10; outputt = [?inittime for this run is ? num2str(hours) ?.? num2str(floor(mins/10)*10)] 126 % calculate first value of decimal time in minutes % to keep the next loop simpler. dectime(1) = mins + secs/60 - floor(mins/10)*10; %need to make sure time calcs are not messed %up if storm carries over to other days daynum = 1; % this loop converts times into minutes past time of first data set. for i = 2:length(times) rawtime = times(i); if ( times(i) < times(i-1) ) daynum = daynum + 1; end hours = floor(rawtime/10000); mins = floor((rawtime - hours*10000) / 100); secs = rawtime - hours*10000 - mins * 100; dectime(i) = (daynum-1)*24*60 + hours*60 + mins + secs/60 - inittime; end % now regroup rainfall data into ten-minute grids colout = 1; sizenew = size(newgrid); highcolin = sizenew(2); colin = 1; colout = 1; resid = 0; for colin = 1:length(dectime) if (colin == 1) rainval = 1; totime = dectime(colin); resid = totime; 127 colin = colin + 1; rainval = 2; else totime = dectime(colin) - 10*(colout-1); rainval = rainval + 1; colin = colin + 1; end % NOTE: USING RAINVAL CONCEPT FAILS IF TOTIME > 20 if (totime > 10) % not while since totime < 20 rainval = rainval -1; if (colout == 1) % don?t want to call dectime(0) totime = totime - 10; outtime(colout) = colout * 10; rainfall = newgrid(:,(colin-1)- rainval+1)/60*dectime((colin-1)- rainval+1); for x = 2:(rainval-1) timeinter = dectime((colin-1)- rainval+x) ? dectime((colin-1)-rainval+x-1); rainfall = rainfall + newgrid(:,(colin-1)- rainval+x)/60*timeinter; end resid = dectime(colin-1) - colout*10; timeinter = dectime((colin-1)) - dectime((colin-2)) ? resid; rainfall = rainfall + newgrid(:,(colin ? 1))/60*timeinter; outgrid(:,colout) = rainfall; colout = colout + 1; else totime = totime - 10; outtime(1,colout) = colout * 10; rainfall = resid/60*newgrid(:,(colin-1)- rainval); for x = 1:rainval 128 timeinter = dectime((colin-1)- rainval+x) ? dectime((colin-1)-rainval+x-1); rainfall = rainfall + newgrid(:,(colin-1)- rainval+x)/60*timeinter; end resid = dectime(colin-1) - colout * 10; rainfall = rainfall - resid/60*newgrid(:,(colin-1)); outgrid(:,colout) = rainfall; colout = colout + 1; end rainval = 1; end end if (rainval == 1) outtime(1,colout) = colout * 10; rainfall = resid/60*newgrid(:,(colin-1)); outgrid(:,colout) = rainfall; else outtime(1,colout) = colout * 10; rainfall = resid/60*newgrid(:,(colin-2)); timeint = dectime(colin-1) - dectime(colin-2); rainfall = rainfall + timeint/60*newgrid(:,(colin- 1)); outgrid(:,colout) = rainfall; end finalout = [outtime;outgrid]; newval = size(finalout); maxi = newval(2); % write output to a temporary text file. % got this info from 'help fprintf' fid = fopen('biggrid.txt','w'); midstr = blanks(1); for i = 1:maxi midstr = str2mat(midstr,'%'); midstr = str2mat(midstr,'8'); 129 midstr = str2mat(midstr,?.?); midstr = str2mat(midstr,?3?); midstr = str2mat(midstr,?f?); midstr = str2mat(midstr,blanks(1)); end midstr = str2mat(midstr,?\?); midstr = str2mat(midstr,?n?); begstr = str2mat(?f?); begstr = str2mat(begstr,?p?); begstr = str2mat(begstr,?r?); begstr = str2mat(begstr,?i?); begstr = str2mat(begstr,?n?); begstr = str2mat(begstr,?t?); begstr = str2mat(begstr,?f?); begstr = str2mat(begstr,?(?); begstr = str2mat(begstr,?f?); begstr = str2mat(begstr,?i?); begstr = str2mat(begstr,?d?); begstr = str2mat(begstr,?,?); endstr = str2mat(?,?); endstr = str2mat(endstr,?f?); endstr = str2mat(endstr,?i?); endstr = str2mat(endstr,?n?); endstr = str2mat(endstr,?a?); endstr = str2mat(endstr,?l?); endstr = str2mat(endstr,?o?); endstr = str2mat(endstr,?u?); endstr = str2mat(endstr,?t?); endstr = str2mat(endstr,setstr(39)); endstr = str2mat(endstr,?)?); endstr = str2mat(endstr,?;?); outstr = str2mat(begstr,setstr(39),midstr,setstr(39),endstr)? eval(outstr) MATLAB program findflow.m -- This program converts adjusted gauge elevations to flow values using data available in USGS rating tables. 130 % load all height data for the gauges % heights has 8 colums load heights -ascii % load all data from the usgs rating tables % each rating has two colums % first column is for elevation % second column is for flow load ratings1 -ascii % lenheight is number of flow values to be calculated lenheight = length(heights(:,1)); for x = 1:lenheight value = heights(x,1); notfound = 1; pos = 1; while notfound > 0 if value > ratings(pos,1) pos = pos + 1; else notfound = -1; end end interp = (value - ratings(pos-1,1) ) / 0.01; flow(x) = ratings(pos-1,2) + (ratings(pos,2) - ratings(pos-1,2))*interp; end flow = flow?; AML for developing DLG data ? This program was used to process raw hydrography data of the Houston area. 131 dlgarc optional ho3hyf01 anghyd1 dlgarc optional ho3hyf02 anghyd2 dlgarc optional ho3hyf03 anghyd3 dlgarc optional ho3hyf04 anghyd4 dlgarc optional ho3hyf05 anghyd5 dlgarc optional ho3hyf06 anghyd6 dlgarc optional ho3hyf07 anghyd7 dlgarc optional ho3hyf08 anghyd8 build anghyd1 line build anghyd2 line build anghyd3 line build anghyd4 line build anghyd5 line build anghyd6 line build anghyd7 line build anghyd8 line reselect anghyd1 angwatr1 line # line res rpoly# > 1 N Y res lpoly# > 1 N N reselect anghyd2 angwatr2 line # line res rpoly# > 1 N Y res lpoly# > 1 N N reselect anghyd3 angwatr3 line # line res rpoly# > 1 N Y 132 res lpoly# > 1 N N reselect anghyd4 angwatr4 line # line res rpoly# > 1 N Y res lpoly# > 1 N N reselect anghyd5 angwatr5 line # line res rpoly# > 1 N Y res lpoly# > 1 N N reselect anghyd6 angwatr6 line # line res rpoly# > 1 N Y res lpoly# > 1 N N reselect anghyd7 angwatr7 line # line res rpoly# > 1 N Y res lpoly# > 1 N 133 N reselect anghyd8 angwatr8 line # line res rpoly# > 1 N Y res lpoly# > 1 N N append hydromap angwatr1 angwatr2 angwatr3 angwatr4 angwatr5 angwatr6 angwatr7 angwatr8 project cover hydromap areahyd utmalb.prj build areahyd line clip areahyd border hydarcs line &return 134 APPENDIX B SUPPLEMENTAL DATA T c and R values from the USACE-Galveston hydrologic study of the Addicks and Barker Reservoir regions. Addicks Reservoir Sub-watershed Area (mi 2 ) T c (hr) R (hr) 10.19 1.67 3.25 3.19 1.00 1.94 3.22 1.04 2.02 2.80 1.02 1.98 3.10 1.00 1.94 2.43 0.93 1.81 2.14 0.93 1.81 2.46 0.93 1.81 2.42 0.98 1.90 1.39 0.81 1.57 4.20 1.06 2.05 2.34 0.96 1.86 2.81 0.96 1.86 2.61 0.96 1.86 2.02 0.93 1.81 1.99 0.91 1.77 2.15 0.93 1.81 2.60 0.98 1.90 2.79 0.98 1.90 2.52 0.98 1.90 1.87 0.96 1.86 1.47 0.89 1.72 3.16 1.26 2.44 1.69 0.89 1.72 1.15 0.78 1.52 1.31 0.65 1.26 5.48 1.36 2.63 2.97 1.00 1.94 1.98 0.39 1.81 1.87 0.86 1.67 1.21 0.78 1.52 3.24 1.08 2.09 2.21 0.91 1.77 135 1.54 0.81 1.57 0.72 0.72 1.40 0.98 1.00 1.94 3.91 0.98 1.90 3.25 0.98 1.90 2.43 0.96 1.86 2.54 0.91 1.77 1.38 0.81 1.57 1.49 0.84 1.62 2.00 1.18 2.29 0.93 0.72 1.40 0.67 0.72 1.40 21.89 6.05 11.75 3.27 1.00 1.94 Barker Reservoir Sub-watershed Area (mi 2 ) T c (hr) R (hr) 19.03 1.77 3.43 2.34 0.91 1.77 1.95 0.89 1.72 0.94 0.75 1.46 3.49 0.96 1.86 3.69 1.26 2.44 1.54 0.84 1.62 1.50 0.84 1.62 1.37 0.75 1.46 6.94 1.46 2.83 3.30 1.04 2.02 4.24 1.08 2.09 2.80 0.98 1.90 4.06 0.84 1.62 4.53 1.15 2.23 2.91 1.02 1.98 2.64 1.06 2.05 1.71 0.86 1.67 3.19 0.75 1.46 5.48 1.26 2.44 2.52 0.96 1.86 2.34 0.93 1.81 1.42 0.81 1.57 1.23 0.81 1.57 136 0.47 0.75 1.46 11.14 1.48 2.81 2.03 0.91 1.77 1.53 0.78 1.52 1.14 0.78 1.52 1.68 1.08 2.09 0.50 0.75 1.46 25.20 7.40 14.40 2.10 1.10 2.13 1.24 0.81 1.57 0.36 0.75 1.46 0.39 0.78 1.52 137 REFERENCES Bedient, P. B. and W. C. Huber, Hydrology and Floodplain Analysis, 2nd Ed., Addison-Wesley Publishing Co., New York, NY, 1992. Bernard Johnson Inc., "Addicks Reservoir Watershed Management Study.". Prepared for the Harris Country Flood Control District, April,1992. Chow, V.T., Maidment, D.R., and Mays, L.M., Applied Hydrology, McGraw-Hill, Inc., New York, NY, 1988. Hellweger, F. and Maidment, D.R., "HEC-PREPRO: A GIS Preprocessor for Lumped Parameter Hydrologic Modeling Programs," CRWR Online Report 97-8. The University of Texas at Austin, Department of Civil Engineering. Austin, TX, 1997. http://www.crwr.utexas.edu/online.html Hill, J.M., "Accuracy and Application Assessment of the Airborne LIDAR Topographic Mapping System (ALTMS)," Houston Advanced Research Center, Woodlands, TX, October, 1996. HEC, "GridParm: Procedures for Deriving Grid Cell Parameters for the ModClark Rainfall-Runoff Model," U.S. Army Corps of Engineers Hydrologic Engineering Center, Davis,CA, November,1996. HEC, "ModClark Model Development for the Muskingum River Basin, OH," U.S. Army Corps of Engineers Hydrologic Engineering Center, Davis,CA, October, 1996. HEC, "A Pilot Application of Weather Radar-Based Runoff Forecasting, Salt River Basin, MO," U.S. Army Corps of Engineers Hydrologic Engineering Center, Davis,CA, May, 1996. Maidment, D.R.,Editor, Handbook of Hydrology, McGraw-Hill, New York, NY, 1993. Olivera, F.,"HEC-PrePro v. 2.0: An ArcView Pre-Processor for HEC?s Hydrologic Modeling System," presented at the 1998 ESRI User?s Conference, San Diego, CA., 1998. http://www.ce.utexas.edu/prof/olivera/header.htm 138 Peters, J.C., and Easton, D.J. "Runoff Simulation Using Radar Rainfall Data." Water Resources Bulletin of the American Water Resources Association, 32 (1996), 753-760. Reed, S.R. and Maidment, D.R., "A GIS Procedure for Merging NEXRAD Precipitation Data and Digital Elevation Models to Determine Rainfall-Runoff Modeling Parameters," CRWR Online Report 95-3. The University of Texas at Austin, Department of Civil Engineering. Austin, TX, 1995. http://www.crwr.utexas.edu/online.html Hellweger, F. and Maidment, D.R., "HEC-PREPRO: A GIS Preprocessor for Lumped Parameter Hydrologic Modeling Programs," CRWR Online Report 97-8. The University of Texas at Austin, Department of Civil Engineering. Austin, TX, 1997. http://www.crwr.utexas.edu/online.html U.S. Army Enginnering District, Galveston, "Hydrology of the Buffalo Bayou and Tributaries," U.S. Army Corps of Engineers, Galveston, TX, August, 1977. 139 VITA Seth Robert Ahrens was born in Harrisburg, Pennsylvania on 15 November 1973, the son of Sara Ann Jefferson Ahrens and Christian Voss Ahrens. After completing his work at Camp Hill High School, Camp Hill, Pennsylvania, in 1992, he entered the University of Virginia in Charlotteville, Virginia. In January, 1995, he transferred to Rice University in Houston, Texas where he received a Bachelor of Science in Chemical Engineering in May, 1997. In August, 1997, he entered the Graduate School at the University of Texas at Austin. Permanent Address: 21 S. 26 th St. Camp Hill, PA 17011-4612 This thesis was typed by the author.