i CRWR Online Report 04-08 Geospatial Description of River Channels in Three Dimensions by Venkatesh Mohan Merwade, PhD. Graduate Research Assistant And David R. Maidment, PhD. Principal Investigator August 2004 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.shtml ii Copyright by Venkatesh Mohan Merwade 2004 iii Acknowledgements I would like to thank Dr. David Maidment, my advisor, for his support, guidance and new visions. I would also like to thank the Texas Water Development Board for funding and supporting my PhD research. Thanks go Barney Austin and Tim Osting from the Texas Water Development Board for their ideas and interest in my work. iv Table of Contents List of Tables.......................................................................................................vii List of Figures.....................................................................................................viii List of Figures.....................................................................................................viii Chapter 1 Introduction...........................................................................................1 1.1 Introduction..........................................................................................1 1.2 Background and Motivation.................................................................2 1.3 Objectives ............................................................................................9 1.4 Hypotheses And The Scope of Work .................................................11 1.5 Contributions From This Research and its Significance.....................13 1.5.1 Anisotropic considerations in the spatial interpolation of river channel bathymetry ...................................................................14 1.5.2 Description of channel bathymetry using cross-sections and profile-lines ...............................................................................14 1.5.3 Analytical model for description of river channels....................14 1.6 Dissertation Outline ...........................................................................16 Chapter 2 Literature Review................................................................................17 2.1 Definition of GIS and River Channel Terms ......................................17 2.2 Instream Flow Studies........................................................................22 2.3 Hydrodynamic Models.......................................................................26 2.3.1 One-dimensional Hydrodynamic Modeling ..............................26 2.3.2 Two-dimensional Hydrodynamic Modeling..............................30 2.3.3 Three-dimensional Hydrodynamic Modeling............................35 2.4 GIS and Hydrodynamic modeling......................................................35 2.5 GIS And River Channel Morphology.................................................39 2.6 Channel Bathymetry Over Large Spatial Domain ..............................42 v 2.7 Meandering of River Channels...........................................................45 2.8 Hydraulic Geometry of River Channels .............................................47 2.9 Summary............................................................................................50 Chapter 3 Data Collection....................................................................................54 3.1 Introduction........................................................................................54 3.2 Brazos River ......................................................................................55 3.3 Guadalupe River ................................................................................57 3.4 Sulphur River.....................................................................................58 3.5 Data Collection ..................................................................................59 3.5.1 Multibeam Echo-Sounder Data .................................................60 3.6 Data Description ................................................................................61 Chapter 4 Development of a Geospatial Structure for River Channels ................65 4.1 Introduction........................................................................................65 4.2 Three-dimensional Representation of Polylines .................................65 4.3 Curvilinear Orthogonal Coordinate System .......................................68 4.4 FishNet...............................................................................................71 4.5 Methodology......................................................................................72 4.5.1 Identifying the Thalweg Location .............................................74 4.5.2 Procedure to Assign (s,n,z) Coordinates....................................81 4.5.3 Data Transformation .................................................................86 4.5.4 Spatial Interpolation of Bathymetry ..........................................88 4.5.5 Creating a FishNet...................................................................122 4.5.6 Transforming the FishNet .......................................................124 4.6 Summary..........................................................................................127 Chapter 5 An Analytical Model for Extrapolation of River Channel Bathymetry ...............................................................................................128 5.1 Introduction......................................................................................128 5.2 Conceptual Model............................................................................129 vi 5.3 Methodology for Developing RCMM..............................................131 5.3.1 Data Transformation and Normalization .................................133 5.3.2 Relationship Between Channel Planform and Thalweg Location ..................................................................................136 5.3.3 Relationship Between Thalweg Location and Cross-sectional Form 143 5.3.4 Rescaling the normalized cross-sections .................................161 5.3.5 Creating profile-lines using cross-sections ..............................167 5.4 Application of the RCMM to Site 1 of the Brazos River..................168 5.5 Verification of RCMM.....................................................................174 5.6 Application of RCMM to Guadalupe River and Sulphur River in Texas................................................................................................180 5.6.1 Hydraulic Geometry Relationships..........................................180 5.6.2 Thalweg Prediction .................................................................182 5.6.3 Description of cross-sections...................................................189 5.7 Regional Application of RCMM......................................................192 5.8 Modeling of Prediction Errors..........................................................194 5.9 Summary..........................................................................................198 Chapter 6 Conclusions and recommendations ...................................................200 6.1 Three-dimensional Description of River Channels Using Cross- sections and Profile-lines .................................................................200 6.2 Analytical Model for Extrapolation of River Channel Bathymetry ..202 6.3 Recommendations For Future Work ................................................205 Appendix A .......................................................................................................207 References .........................................................................................................225 Vita ..................................................................................................................237 vii List of Tables Table 3.1 Raw bathymetry data. The elevation is recorded with respect to the mean sea level ..........................................................................62 Table 4.1: Ranking of spatial interpolation methods for the sample dataset on along Site 1 of Brazos River............................................120 Table 4.2: Ranking of spatial interpolation methods for the entire dataset along Site 1 of Brazos River.............................................................121 Table 5.1: Summary of different analytical forms for fitting channel cross-section.....................................................................................156 Table 5.2: Summary of regression analysis for fitting a combination of beta distributions to the cross-section data .......................................157 Table 5.3: Beta parameters for different thalweg locations .......................160 Table 5.4. Hydraulic geometry parameters for the Brazos River. A is the drainage area; Q 2 is 2-year return period flow; a, c, k, and b, f, m, are the coefficients and power terms of equations (5.20), (5.21), and (5.22), respectively ....................................................................165 Table 5.5: Parameters of radius of curvature ? thalweg location relationship for Brazos and Guadalupe River...................................186 Table 5.6: Semivariogram parameters for prediction errors at Site 1 and Site 2 along Brazos River.................................................................196 Table 5.7: Semivariogram parameters for prediction errors Guadalupe and Sulphur Rivers...........................................................................197 viii List of Figures Figure 1.1 Instream flow study segments in Texas .......................................5 Figure 1.2 Bathymetry data representation using points, raster, and TIN......7 Figure 1.3 Bathymetry data as points, cross-sections and a finite-element mesh.....................................................................................................8 Figure 1.4 A wireframe model representation of sphere ...............................9 Figure 1.5 Three-dimensional description of a river channel using cross- sections and profile-lines....................................................................11 Figure 1.6 Three-dimensional description of a river channel using cross- sections, profile-lines and hydro-volume ...........................................15 Figure 2.1 Spatial representation of simple features....................................18 Figure 2.2 (a) Data representation using raster; (b) Data representation using TIN ...........................................................................................19 Figure 2.3 (a) Channel planform definition; (b) Cross-sections and profile-lines........................................................................................21 Figure 2.4 Basic components of instream flow studies ...............................23 Figure 2.5 River system schematic in HEC-RAS........................................28 Figure 2.6 A typical HEC-RAS cross-section .............................................29 Figure 2.7 Finite element mesh for two-dimensional hydrodynamic modeling ............................................................................................32 Figure 2.8 Relative importance of user defined variables in RMA2 performance (From Donnell et. al., 2001) ..........................................34 Figure 2.9 Representation of river channel in Cartesian coordinate system ................................................................................................40 Figure 3.1 Location of Brazos, Guadalupe and Sulphur watersheds in Texas..................................................................................................55 Figure 3.2 Location of Site 1 along the Brazos River..................................56 Figure 3.3 Location of Site 2 along the Brazos River..................................57 Figure 3.4 Location of study area along the Guadalupe River ....................58 Figure 3.5 Location of study area along the Sulphur River.........................59 Figure 3.6 Bathymetry data collection using depth sounder and GPS.........60 ix Figure 3.7 a) Bathymetry data as point features in GIS; b) Attribute table linked to the data................................................................................62 Figure 3.8 Data density description at Site 1 and Site 2 along the Brazos River ..................................................................................................63 Figure 3.9 a) Real cross-section with a steep slope; b) Cross-section measured using a single beam depth sounder.....................................64 Figure 4.1 PolylineZM representation in ArcGIS. The m coordinate is a measure value along the line and the z coordinate stores the elevation.............................................................................................66 Figure 4.2 Relative and absolute measures in ArcGIS. ...............................67 Figure 4.3 Orthogonal curvilinear (s,n) coordinate system for river channels .............................................................................................69 Figure 4.4 Surface representation in ArcGIS using FishNet........................71 Figure 4.5 a) Site 1 on the Brazos River; b) Small section of Site 1 used for illustrations...................................................................................74 Figure 4.6 Flow chart for identifying thalweg.............................................76 Figure 4.7 Procedure for the thalweg identification in GIS. (a) Input data; (b) Raster surface created using bathymetry points; (c) Normal lines created for the initial line at 15m interval; (d) initial line is moved to the deepest points along all the normals; (e) Final result with initial line changed to thalweg; (f) Vertical profile of a normal showing the location of initial line and the deepest point (thalweg)............................................................................................78 Figure 4.8 Properties of normal line extracted from the channel bed (raster surface) ...................................................................................79 Figure 4.9 Braided river channel separated into two active channels A and B by an island..............................................................................80 x Figure 4.10 Assigning (s,n) coordinates to a bathymetry point P. (a) Point P is a bathymetry point in the proximity of segments AB and BC of the centerline; (b) The association of P with AB or BC is decided based on angles made by P with endpoints of AB and BC; (c) P is closest to BC because it makes acute angles with endpoints of BC; (d) P makes acute angles with both BC and CD. PG is smaller than PF, and therefore P is closest to BC; (e) P cannot be uniquely associated with either AB or BC because it makes one obtuse angle with both segments; (f) Straight line segments are converted to Bezier curves to maintain tangent continuity at each vertex of the centerline.......................................................................82 Figure 4.11 Polyline to piecewise Bezier curve conversion ........................85 Figure 4.12 Polyline with very sharp angles ...............................................86 Figure 4.13 Attribute table of bathymetry points with (s,n) coordinates .....87 Figure 4.14 Data transformation from (x,y,z) to (s,n,z). (a) Bathymetry data in (x,y) coordinates; (b) Bathymetry data in (s,n) coordinates ....87 Figure 4.15 Inverse distance weighting.......................................................89 Figure 4.16 River channel geometry ...........................................................91 Figure 4.17 Elliptical Inverse Distance Weighting......................................92 Figure 4.18 Spline Function........................................................................94 Figure 4.19 Ordinary kriging ......................................................................98 Figure 4.20 Semivariogram plot................................................................100 Figure 4.21 Semivariogram envelope for anisotropic kriging ...................104 Figure 4.22 Different search neighborhoods for anisotropic kriging.........105 Figure 4.23 Bathymetry points for creating training and test dataset ........106 Figure 4.24 Training and test datasets for spatial interpolation .................107 Figure 4.25 Test to find an optimum power for inverse distance weighting .........................................................................................108 Figure 4.26 Test to find an optimum search neighborhood for inverse distance weighting............................................................................108 Figure 4.27 Test to find an optimum search neighborhood for elliptical inverse distance weighting ...............................................................110 Figure 4.28 Test to find an optimum anisotropy ratio for elliptical inverse distance weighting ...............................................................111 xi Figure 4.29 Test to find an optimum value of ? for regularized spline......112 Figure 4.30 Test to find an optimum number of points for regularized spline with ? = 0...............................................................................113 Figure 4.31 Test to find an optimum value of ? for tension spline............114 Figure 4.32 Test to find an optimum number of points for tension spline.114 Figure 4.33 Effect of tension parameter on spline interpolation................115 Figure 4.34 A spherical semivariogram model with range, sill, and nugget equal to 55m, 2.4, and 0.094, respectively............................116 Figure 4.35 Test to find an optimum number of points for ordinary kriging..............................................................................................116 Figure 4.36 Semivariogram envelope for anisotropic kriging ...................117 Figure 4.37 Test to find an optimum number of points for anisotropic kriging..............................................................................................118 Figure 4.38 Summary of RMSE results for all spatial interpolation methods............................................................................................119 Figure 4.39 FishNet tool interface.............................................................123 Figure 4.40 FishNet generated from a raster surface.................................124 Figure 4.41 Transformation of a point from (s,n) coordinate system to (x,y) coordinate system. (a) Point P as one of the vertices of a line in (s,n) coordinates; (b) Mapping of P in (x,y) coordinates with respect to segment BC of the centerline. ..........................................125 Figure 4.42 Hydraulic FishNet..................................................................126 Figure 5.1 Conceptual model for RCMM. (a) Channel planform and thalweg; (b) cross-sectional forms for different thalweg locations...129 Figure 5.2 Flowchart for the development of RCMM...............................132 Figure 5.3 A typical cross-section with (s,n,z) coordinates. P(n i ,z i ) is any point in a cross-section; n L and n R are the n coordinates of left and right bank, respectively; Z b and Z d are the bank elevation and thalweg elevation, respectively ........................................................134 Figure 5.4 Data normalization for a cross-section. (a) Bathymetry data in original coordinates; (b) Bathymetry data in normalized coordinates.......................................................................................136 Figure 5.5 Locations of segments used for computing the radius of curvature ..........................................................................................138 xii Figure 5.6 Radius of curvature calculations to quantify the meandering shape of the channel.........................................................................139 Figure 5.7 Relationships between thalweg location and radius of curvature ..........................................................................................140 Figure 5.8 Map with radius of curvature signs at meandering bends.........142 Figure 5.9 Bathymetry data showing a cross-section ................................144 Figure 5.10 Power function fitted to the left hand side of the cross- section shown in Figure 5.9..............................................................144 Figure 5.11 Power function fitted to the right hand side of the cross- section shown in Figure 5.9..............................................................145 Figure 5.12 Quadratic function fitted to the channel cross-section............146 Figure 5.13 Third-degree polynomial fitted to the channel cross-section..147 Figure 5.14 Interpolating and smoothing splines ......................................148 Figure 5.15 Smoothing spline fitted to the cross-section...........................149 Figure 5.16 Beta distribution.....................................................................150 Figure 5.17 Beta distribution fitted to the cross-section ............................151 Figure 5.18 Combination of two beta distributions ...................................152 Figure 5.19 Combination of Beta distributions fitted to the channel cross-section.....................................................................................153 Figure 5.20 Gamma density functions for different (a,b) pairs .................154 Figure 5.21 Gamma distribution fitted to the cross-section.......................155 Figure 5.22 Combination of gamma distributions fitted to the channel cross-section.....................................................................................156 Figure 5.23 Beta parameters for different thalweg locations.....................159 Figure 5.24 Hydraulic geometry relationship between (a) average width; (b) average depth; (c) average velocity and flow of the Brazos River at Richmond ...........................................................................163 Figure 5.25 Gaging stations locations along the Brazos River ..................164 Figure 5.26 Relationship between 2-year return period flow (Q 2 ) and drainage area (A) for the Brazos River.............................................166 Figure 5.27 Generation of profile-lines using cross-sections.....................168 xiii Figure 5.28 Thalweg prediction using the radius of curvature-thalweg location relationship (equations 5.5 and 5.6). Circled areas show locations where the thalweg prediction is imprecise. Sections A-A and B-B are plotted in Figures 5.30 and 31 ......................................170 Figure 5.29 Thalweg prediction using the non-dimensional radius of curvature-thalweg location relationship (equations 5.7 and 5.8) ......171 Figure 5.30 Cross-section described by RCMM at section A-A shown in Figure 5.28.......................................................................................172 Figure 5.31 Cross-section described by RCMM at section B-B shown in Figure 5.28.......................................................................................173 Figure 5.32 Three-dimensional description of a river channel in the form of cross-sections and profile lines ....................................................174 Figure 5.33 Observed data and model output at Site 2 ..............................176 Figure 5.34 Observed data and model output at Site 1 ..............................177 Figure 5.35 Histogram of prediction errors at Site 1. Group 1 covers 84 percent of the data with low prediction errors, and Group 2 cover data with high prediction errors........................................................178 Figure 5.36 Comparison of obverted elevation and model output for group 1 at Site 1 ...............................................................................178 Figure 5.37 Spatial distribution of prediction errors at Site 1 of Brazos River. (a) Lightly and darkly colored points cover the area with prediction errors in Group 1 and Group 2, respectively. (b) Measured thalweg and the thalweg modeled by RCMM..................179 Figure 5.38 Location of study area along the Guadalupe River ................181 Figure 5.39 Location of study area along the Sulphur River .....................182 Figure 5.40 Thalweg prediction along the Guadalupe River using RCMM.............................................................................................183 Figure 5.41 Thalweg prediction along the Guadalupe River using equations 5.7 and 5.8........................................................................184 Figure 5.42 Thalweg prediction along the Sulphur River using equations 5.7 and 5.8........................................................................................185 Figure 5.43 Radius of curvature versus thalweg location for the Guadalupe River ..............................................................................187 Figure 5.44 Radius of curvature versus thalweg location for the Sulphur River ................................................................................................187 xiv Figure 5.45 Thalweg prediction along the Guadalupe River after incorporating the sinuosity criteria. Circled areas show locations where the thalweg prediction is imprecise........................................188 Figure 5.46 Thalweg prediction along the Sulphur River. Circled areas show locations where the thalweg prediction is imprecise ...............189 Figure 5.47 Spatial description of prediction errors for the Guadalupe River. Group 1 includes points with low prediction errors, and Group 2 includes points with high prediction errors.........................191 Figure 5.48 Spatial distribution of prediction errors for the Sulphur River. Group 1 includes points with low prediction errors, and Group 2 includes points with high prediction errors.........................192 Figure 5.49 Lower reach of the Brazos River............................................193 Figure 5.50 River channel description by RCMM at a location near Richmond gaging station. The cross-sections are approximately 500m apart. ......................................................................................194 Figure 5.51 Semivariogram of prediction errors (across flow, n- direction) for Site 2 on Brazos River................................................195 Figure 5.52 Semivariogram of prediction errors (along flow, s-direction) for Site 2 on Brazos River ................................................................196 1 Chapter 1 Introduction 1.1 INTRODUCTION Geographic Information Systems (GIS) are widely used in hydrology for data storage, representation, and analysis. Standardized datasets are available from national agencies such as the United States Geological Survey and the Environmental Protection Agency to represent the surface terrain and other hydrographic features like watershed boundaries, water bodies, and river channel networks for doing hydrologic studies. The surface terrain and the hydrographic features are represented using standardized formats. For example, surface terrain is generally represented using a Digital Elevation Model (DEM). Similarly, the hydrographic features such as gaging stations, river channel networks, and watersheds are represented using points, lines, and polygons, respectively. Unlike the land surface terrain and hydrography, however, there is no standardized way of representing the three-dimensional structure of river channels. The research presented in this dissertation provides a way of representing the three-dimensional structure of river channels in a geospatial environment. This dissertation has two main components that deal with river channel bathymetry. The first component (described in chapter 4) deals with developing a procedure to represent the three-dimensional structure of river channels in the form of a mesh of cross-sections and profile-lines. The cross-sections transverse the flow, and the profile-lines run along the flow. A network of cross-sections and profile-lines thus provides a three-dimensional mesh for representing channel 2 bathymetry. The second component (described in chapter 5) deals with developing a model for describing the river channel bathymetry over large spatial domains by using the field data collected for small-scale (few kilometers long) study reaches. The model, called River Channel Morphology Model (RCMM), is based on deriving relationships among different channel characteristics such as the channel planform (shape of river in plan-view), the thalweg location, and the cross-sections. In addition, River Channel Morphology Model characterizes the channel planform and cross-sections by analytical functions. Therefore, by knowing only the planform of the channel which is available for large spatial domains, the cross-sections and profile-lines can be described in three- dimensions. 1.2 BACKGROUND AND MOTIVATION The term channel bathymetry, as used in this research, refers to the elevation of the riverbed with respect to the mean sea level. Similar to land surface terrain data that are necessary for many hydrologic studies, channel bathymetry data are necessary for conducting hydrodynamic studies of river channels. River channels are complex systems, and hydrodynamic models are used to describe the distribution, direction, and magnitude of velocities and depths for different flow conditions. The processes caused by flowing water within a river channel are mostly fluvial, such as sediment transport, erosion, deposition, and change in channel planform. The flowing water within a river channel also provides habitat for fish species. Both the fluvial processes and the habitat conditions for fish species can be related to the hydrodynamic conditions within a 3 river channel. For example, the velocity distribution of flowing water within a river channel may suggest which areas are susceptible to erosion. A change in flow condition alters the hydrodynamics, which in turn alters the fluvial processes and habitat conditions. In this way, the results from hydrodynamic models are used to make decisions related to fluvial processes and fish habitat conditions. Instream flow is defined as the minimum flow necessary to maintain an ecologically sound environment in river channels (Wentzel, 2001). An ecologically sound environment includes the diversity and productivity of different types of fish species. In addition to ecological aspects, instream flow also restores economic and aesthetic values of rivers that are useful for recreation and navigation. In 2001, the 77 th session of the Texas Legislature passed Senate Bill 2, which established the Instream Flow Program for Texas Rivers (Texas Instream Flow Studies, 2002). Senate Bill 2 directed the Texas Water Development Board (TWDB), the Texas Parks and Wildlife Department (TPWD), and the Texas Commission on Environmental Quality (TCEQ) to perform engineering and scientific studies to determine minimum flow conditions for maintaining an ecologically sound environment for Texas Rivers. The three important components of instream flow studies are: hydrologic and hydraulic evaluation by TWDB, biological evaluation by TPWD, and water quality evaluation by TCEQ. In compliance with the Senate Bill 2, priority study areas were identified based on potential water development projects, water rights permitting issues, and other factors such as location of existing structures. The priority study areas 4 include the lower Guadalupe River, lower Brazos River, lower San Antonio River, middle Trinity River, lower Sabine River, and middle Brazos River (Figure 1.1). At the time Senate Bill 2 was passed, there were cooperative instream flow studies along the lower Guadalupe River and the lower Brazos River to assess the instream flow requirements for fish and to assess the downstream impact of proposed Allens Creek reservoir, respectively. Marvin Nichols Wright Patman Upper Sulphur Red Bols d?Arc Upper Sabine Prairie Creek Toledo Bend Sam Rayburn Eastex Neches B. A. Steinhagen Livingston Bedias Middle Trinity Middle Brazos Little River Lower Colorado Allens Creek Lower Brazos per Guadalupe Canyon Lower Guadalupe Lower San Antonio San Antonio Austin Dallas Fort Worth Lower Sabine Priority Study 2 n Up d Tier Special Study Existing Reservoirs Proposed Reservoirs Figure 1.1 Instream flow study segments in Texas As a part of the cooperative instream flow studies along the lower Guadalupe River and lower Brazos River, the TWDB has collected flow data, stage data, and bathymetry data. These data are required to perform 5 6 hydrodynamic simulations. Among flow, stage, and bathymetry data, bathymetry data play a major role in the hydrodynamic modeling process (Donnell et. al., 2001). The bathymetry data are collected using a depth sounder. In this method, a depth sounder combined with a Global Positioning System (GPS) is mounted on a boat, and the bathymetry is recorded as a set of scattered (x,y,z) points as the boat moves over the river bed. Even though the depth sounding technique provides a more detailed description of river bathymetry compared to traditional cross- sectional surveys, there are two issues associated with it. First, the riverbed, which is a surface, is usually not represented using points because they do not provide a continuous surface. The bathymetry points can be interpolated to create a raster grid or a Triangulated Irregular Network (TIN) (Figure 1.2). A raster stores the data using a matrix of square cells, and a TIN stores the data as an integrated set of nodes and edges. However, storing the data using a raster grid or a TIN is not desirable because they demand more disk space and computing power. 03015 Meters 4 Bathymetry data as points Bathymetry data as a raster grid Bathymetry data as TIN Figure 1.2 Bathymetry data representation using points, raster, and TIN Second, even though the data can be stored in different forms as shown in Figure 1.2, hydrodynamic models require data as lines. For example, a one-dimensional model such as HEC-RAS requires input bathymetry data in the form of cross- sections. RMA-2, on the other hand, is a two-dimensional model and requires bathymetry data in the form of finite element mesh (Figure 1.3). Transferring the data among different users and different models in different formats results in loss of information related to channel bathymetry. Therefore, storing the data using a suitable format is an issue with channel bathymetry, and the need to develop an efficient and standardized dataset for three-dimensional representation of river bathymetry is a primary motivation for this research. 7 Bathymetry Points Cross-sections Finite element mesh Figure 1.3 Bathymetry data as points, cross-sections and a finite-element mesh When a structure such as a reservoir is built across a river channel, the area of influence upstream and downstream of the structure has to be analyzed in order to make decisions related to the instream flow conditions. This approach requires bathymetry data for the entire river segment in order to conduct hydrodynamic simulations. Collecting data for the entire river segment requires an enormous amount of effort by trained personnel, and is practically not feasible with currently available technology and the funding associated with instream flow projects. Therefore, instead of performing the instream flow studies for the whole river segment, a short representative reach is selected for analysis. The instream flow decisions for the whole river segment are then made based on the analysis for the representative reach. Such regional decisions based on local studies of short representative reaches are debatable. Another motivation for this research is the inability to collect or acquire detailed bathymetry data over a large spatial domain. If bathymetry data are available over large spatial domains, instream 8 flow studies can either be carried out over these large areas, or decisions based on local scale studies can be verified. 1.3 OBJECTIVES As mentioned earlier, there are two issues associated with channel bathymetry. The first issue is about representing the channel bathymetry in an efficient manner where the data are available as (x,y,z) points, and the second issue is about the unavailability of channel bathymetry over large areas. The ability to store the vector data as three-dimensional features (details explained in section 4.2) enables representation of continuous surfaces using wireframe models (Yang et. al., 1994). Figure 1.4 shows a simple wireframe model for a sphere. Figure 1.4 A wireframe model representation of sphere A traditional wireframe model stores data in the form of vertices and edges for any given object. A wireframe model can also be created by using a network of three-dimensional lines to represent surfaces. Therefore, an efficient way of storing and representing channel bathymetry is by using three-dimensional lines, cross-sections and profile-lines. 9 10 Channel planform refers to the two-dimensional (x,y) plan-view of a river channel in space. Each river cross-section used to represent a channel bed has a typical shape, and this shape can be related to the channel planform. For example, in a meandering river channel, the cross-sections at meandering bends are asymmetric (thalweg is nearer to one bank than the other) compared to the cross- sections along a straight reach (thalweg is near the channel center). If the shape of the cross-sections and the channel planform are interrelated using analytical functions, it may be possible to describe the channel cross-sections by knowing the channel planform. Channel planform data are readily available from digital orthophoto quadrangles (DOQ) or from the national hydrography dataset (NHD). Therefore, by knowing only the channel planform, cross-sections can be defined over large spatial domains. This will overcome the second issue of unavailability of channel bathymetry over large areas. Based on the above discussion, this research has two main objectives: 1. Given the channel bathymetry as (x,y,z) points, develop a procedure to store the data in the form of cross-section and profile-lines. A network of cross- sections and profile-lines forms a three-dimensional mesh to represent the channel bathymetry (Figure 1.5). Figure 1.5 Three-dimensional description of a river channel using cross-sections and profile-lines 2. Develop a model in which the channel cross-sections and channel planform are represented using analytical functions. Such a model will allow inter- relating cross-sections with channel planform. Therefore, by knowing only the channel planform, cross-sections can be described. The objectives listed above are accomplished by using two datasets collected by TWDB along the Brazos River in Texas. The analytical model developed to accomplish the second objective is calibrated by using the data along the Brazos River in Texas, and is verified by applying it to the same river (an area upstream of the calibration site), the Guadalupe River, and the Sulphur River in Texas. 1.4 HYPOTHESES AND THE SCOPE OF WORK Based on the objectives listed earlier, this research is trying to answer two questions: 1. Given the river channel bathymetry as (x,y,z) points, how to create a standardized representation in the form of cross-sections and profile-lines? 11 12 2. Given the standardized representation of channel bathymetry, how can this information be used to predict the three-dimensional structure of river channels in areas with no bathymetric data? The above two questions can be categorized into descriptive and predictive types. The first question, which falls into the descriptive category, involves using the available data and techniques, with or without modifications, to come up with an answer. Therefore, it does not require any hypotheses. The second question, which falls into the predictive category, involves using the data at one location to make predictions at another location. Therefore, in addition to using the available data and techniques, it requires some hypotheses. River channels evolve through the process of erosion and deposition, which change both the channel planform (plainimetric shape of the river channel) and the cross-sectional shape (Leopold et. al., 1964; Knighton, 1998). Erosion at one bank and the deposition at the other cause the river channels to meander, which in turn change the shape of river cross-sections. Therefore, the change in channel planform and the cross-sectional shape are interrelated. In addition to change in the shape, the size (cross-sectional area) of cross-sections also change to accommodate the varying flow conditions in the channel. The following two hypotheses are made to answer the second question: i. The shape of river cross-sections is related to the channel planform. ii. The size of river cross-sections is related to the flow in the river channel. Besides channel planform and flow, the river geometry is also influenced by the geological characteristics such as substrate material, floodplain deposits, and 13 several other factors such as climate, land use types, etc (Knighton, 1998). In order to answer the second question completely, the proposed approach should also include geology, climate, and all other factors in its hypotheses. This, however, increases the scope of work beyond a single dissertation research. Therefore, this dissertation is focused only on studying the role of channel planform and flow in describing the river cross-sections. In addition, the research presented in this dissertation is applicable only to the main stem of meandering river channels. This excludes braided river channels, urban streams or constructed channels, and any other channels that are not explicitly meandering. This research also does not include branched river systems or the influence of tributaries on the main stem. The research presented in this dissertation is framed within the context of ArcGIS 8.x, a GIS software package from the Environmental Systems Research Institute (ESRI). ArcGIS stores spatial and temporal data in a relational database. In addition, ArcGIS also offers customization capabilities using Visual Basic and ArcObjects. ArcObjects is the development platform for ArcGIS and provides an infrastructure to build customized applications based on existing components (Waltuch et al., 2001). Visual Basic and ArcObjects are used to develop customized computer routines in this research. 1.5 CONTRIBUTIONS FROM THIS RESEARCH AND ITS SIGNIFICANCE There are three main contributions from this research, and they are briefly discussed below. 14 1.5.1 Anisotropic considerations in the spatial interpolation of river channel bathymetry This research demonstrates the importance of anisotropy in the channel bathymetry data in the application of spatial interpolation schemes. A modification to inverse distance weighting, a spatial interpolation scheme described in chapter 4, is suggested that takes into account the anisotropy in the channel bathymetry. Channel bathymetry plays a major role in hydrodynamic simulations. Therefore, creating an accurate bathymetric surface using spatial interpolation schemes contributes towards better hydrodynamic modeling results. 1.5.2 Description of channel bathymetry using cross-sections and profile- lines A procedure is developed that uses the channel bathymetry data collected as (x,y,z) points to create the channel representation in the form of cross-sections and profile-lines. The channel bathymetry in the form of cross-sections and profile-lines can be used to populate the river channel component of the Arc Hydro data model. Storing the channel bathymetry data as linear features in Arc Hydro data model enables the data to be shared and used by the community in a standardized form. In addition, cross-sections and profile-lines can serve as input to one, two, or three-dimensional hydrodynamic models. 1.5.3 Analytical model for description of river channels An analytical model is developed in this research to create a three- dimensional description of river channels using cross-sections and profile-lines. The analytical model, also called the River Channel Morphology Model, is based on the channel planform and its geometry. Therefore, by knowing only the planform, the channel can be described in three-dimensions. River Channel Morphology Model is a useful tool for creating the three- dimensional description of river channels where there are no bathymetry data. The overall contribution from this research is the new way of geospatially representing river channels. This research describes river channels using cross-sections and profile-lines (Figure 1.6), and such a description is oriented with respect to the direction of flow in the river. At macro level, the area between the two cross- sections can be treated as a control volume (Hydro-Volume), and this concept of Hydro-Volume opens doors for developing hydrodynamic models and water quality models (mass balance models) for river channels using GIS. Figure 1.6 Three-dimensional description of a river channel using cross-sections, profile-lines and hydro-volume 15 16 1.6 DISSERTATION OUTLINE This dissertation is organized into six chapters. The first chapter provides background information about the instream flow studies in Texas, the motivation for this research, the objectives of this research, and the contributions from this research. The second chapter provides a literature review to understand the current state of knowledge and to identify the gaps that this research aims to fill. The third chapter provides an overview of the study area, the data collection procedure for river channel bathymetry using a depth sounder, and a brief description of the data. The fourth chapter describes the procedure for developing a standardized three-dimensional description of river bathymetry using cross- sections and profile-lines. The fifth chapter describes the River Channel Morphology Model for creating a three-dimensional description of river channel bathymetry by inter-relating the cross-sectional form to the river channel planform. The final chapter, chapter six, discusses the conclusions and recommendations for future work based on this research. 17 Chapter 2 Literature Review A literature review is presented to provide a background relating to issues with the bathymetry of river channels. This chapter is divided into five sections. In the first section (2.1), basic terms used in this research are defined. The second section (2.2) illustrates the importance of hydrodynamic modeling in decision making processes related to river channels by using instream flow studies as an example. Section 2.3 discusses different types of hydrodynamic models and their data requirements. Section 2.4 presents an overview of issues specific to GIS and hydrodynamic modeling studies. Section 2.5 briefly discusses some issues associated with using Cartesian coordinate system for analyzing river channels. Section 2.6 provides a background to illustrate the need for representation of river channel bathymetry over large spatial domains. Sections 2.7 and 2.8 present a discussion on meandering of river channels and hydraulic geometry. Finally, section 2.9 presents a discussion to summarize the research gaps and formulate the goals for this research. This chapter uses some GIS and river channel morphology terms that may be new to readers who are not familiar with GIS or river morphology research or both. Therefore, before the literature review, some of the terms that are used in GIS and river morphology research are first defined in section 2.1. 2.1 DEFINITION OF GIS AND RIVER CHANNEL TERMS Some terms as defined in Zeiler (1999) and Knighton (1998) that are used in this research are presented in this section. ? Feature: A feature is a spatial entity such as a point, line or polygon. A single point is called a point feature, a single line is called a line feature, and a single polygon is called a polygon feature. A set of connected lines is called a polyline feature (Figure 2.1). Point Features Line Feature Polyline Feature Polygon Feature Vertices End points Figure 2.1 Spatial representation of simple features ? Feature Class: A feature class is a collection of features with the same type of geometry. For example, a point feature class can only store point features. A feature class is a table with each row representing a single feature and each column representing a field. The values stored in the fields are the attributes of each feature. A feature class also stores the geometry of features as attributes in the shape field. A feature cannot be stored outside a feature class. ? Feature Dataset: A feature dataset is a collection of feature classes that share the same coordinate system. A feature class can also exist by itself outside of feature datasets. ? Object class: An object class, which does not store any geometry information, is a table that stores only descriptive information. As such, a row in an object class does not represent any geographic feature. An object class 18 can be associated with features in a feature class, but it cannot store or represent features by itself. For example, time-series records of discharge measurements stored in a table is an object class. ? Geodatabase: A geodatabase is the top-level unit of geographic data. It is a collection of feature datasets, feature classes, and relationship classes. A relationship class is a table that stores relationships between objects and features. For example, a gaging station (point feature) can be related to an object class (time series table) through a relationship class. ? Vector: Vector data represent features as points, lines, and polygon features (Figure 2.1). ? Raster: A raster is a regularly spaced matrix of cells with each cell having associated attribute information (Figure 2.2). For example, a digital elevation model (DEM), which is a raster dataset, has an elevation value associated with each cell. 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 19 1 2 3 4 5 a) Raster b) Triangulated Irregular Network (TIN) Cell Cell Size Nodes Edge Face Figure 2.2 (a) Data representation using raster; (b) Data representation using TIN 20 ? TIN (Triangulated Irregular Network): A TIN is an integrated set of points (nodes) with elevations and triangles with edges. A TIN is therefore made of points (nodes), lines (edges), and polygons (faces) (Figure 2.2). ? Channel planform: Channel planform refers to the two-dimensional (x,y) plan-view of a river channel in space (Figure 2.3a). ? Channel cross-section: A river channel cross-section is a ground profile transverse to the flow direction (Figure 2.3b). ? Channel profile-line: A river channel profile-line refers to the ground profile parallel to flow direction (Figure 2.3b). ? Thalweg: A thalweg is a polyline connecting the points of lowest bed elevation in a river channel (Figure 2.3b). ? Meander length: The distance along the channel at the meandering bend is called meander length (M L in Figure 2.3a) ? Valley length: The straight line distance between any two points along a river channel is called valley length (Figure 2.3a) ? Sinuosity: The meander length divided by the straight-line valley length is the sinuosity of the meandering bend. Valley length ML x y A Flow direction a) b) Cross-section Profile lines Thalweg Three-dimensional view at A Figure 2.3 (a) Channel planform definition; (b) Cross-sections and profile-lines 21 22 2.2 INSTREAM FLOW STUDIES Human activities alter the flow regime in a river channel. Consumptive water use from a river channel for irrigation, industrial, and municipal purposes decrease the amount of water available for organisms living in a river channel. Even non-consumptive water use like hydroelectric power generation can significantly alter the flow patterns within a river channel, thus affecting the fish community. In order to minimize human impacts due to consumptive and non- consumptive use of river water, a minimum amount of water has to be reserved for maintaining healthy ecological conditions within a stream channel. The amount of water reserved for this purpose is referred as instream flow (Wentzel, 2001). Since the mid 1970s, the instream flow technique has matured from setting fixed minimum flows to incremental methods in which aquatic habitats are quantified as a function of stream discharge (Stalnaker et. al., 1995). According to Wentzel (2001), a typical instream flow study has four basic components: 1) habitat descriptions and biological sampling, 2) hydrodynamic modeling, 3) habitat model, and 4) decision making (Figure 2.4). The first component, habitat descriptions and biological sampling, includes identifying different types of fish species, and the conditions under which they live. These conditions include substrate type, water depth and velocity, temperature, etc. The hydrodynamic modeling component includes modeling of spatial and temporal variations in depth and velocity for different flow-rates. A habitat model combines the results of hydrodynamic model and biological sampling to determine the spatial extent of the habitat for different fish species under varying flow conditions. Hydrodynamic Modeling Habitat Descriptions/ Biological Sampling Habitat Model Instream Flow Decision Making Figure 2.4 Basic components of instream flow studies Instream flow studies vary based on their criteria for describing habitats, assumptions and techniques for hydrodynamic and habitat modeling, and the decision making process. For example, the habitat model proposed by Austin and Wentzel (2001) for instream flow studies in Texas is coupled with GIS to provide a visual output for each analysis to make the decision making process easier. All the components shown in Figure 2.4 are equally important for instream flow studies. However, only the hydrodynamic modeling component, which is the motivating factor behind this research, is discussed in the following sections. Based on flow assessment methods, Karim et. al. (1995) group instream flow methods into three major categories: 1) historic flow methods, 2) hydraulic rating methods, and 3) habitat rating methods. Historic flow methods use recorded or estimated flow statistics at a single point along a river channel. These methods 23 24 assume a strong relationship between the natural flow regime and the biological health of a river. One of the goals of these methods is to recommend a minimum flow that is within the historic flow range because it is assumed that aquatic life has survived such flows in the past. The minimum recommended flow could be average monthly flow, one in five year low flow, flow equaled or exceeded 96% of time, etc. (Jowett, 1997). The historic flow methods are the easiest and least sophisticated of all the methods. The Tennant (1976) method, also known as the Montana method, is the most widely known historic flow method. The Tennant method correlates habitat quality with various percentages of mean annual flow (MAF). The flow conditions in the Tennant method range from 10 percent MAF for severely degraded habitat conditions to 200 percent MAF for flushing flows. Hydraulic rating methods relate hydraulic parameters of channel cross- sections, such as width, depth, velocity and wetted perimeter, with discharge. These methods assume a relationship between certain hydraulic parameters and the available fish habitat. The most commonly used hydraulic parameter is wetted perimeter. Variations in hydraulic parameters are established by measurements at different flows, prediction from cross-section surveys and stage-discharge relationships, or calculation of water surface profiles for different flows (Richardson, 1986). Habitat rating methods are an extension of the hydraulic rating methods. These methods involve relating the hydraulic conditions associated with different flows to specific habitat types. Unlike historic rating methods and hydraulic rating methods, habitat rating methods involve hydrodynamic modeling of river systems 25 to predict water depth and velocity throughout a river reach. First, habitat suitability criteria are developed for different hydraulic conditions. Results from hydraulic models are then compared with these criteria to determine the area of suitable habitat for the target fish species. Since the habitat rating methods are physically based, these methods are considered more reliable. Physical habitat simulation (Milhous et. al., 1989) also called PHABSIM is a widely used habitat rating instream flow method in the United States. Since the advent of PHABSIM, which uses a one-dimensional hydrodynamic model, several habitat rating methods employing two-dimensional hydrodynamic methods have been developed. Leclerc et. al. (1991), Ghanem et. al. (1996), Waddle et. al. (1997), and Parasiewicz and Dunbar (2001) are some of the examples that use two-dimensional hydrodynamic models with habitat rating methods. In general, the instream flow studies have become more sophisticated over time due to increased understanding of the biological processes, the hydrodynamics of river systems, and the availability of faster computers and advanced measurement techniques. While the shift from a one-dimensional habitat model to a two-dimensional model was being made, several studies comparing these two methods were published. The comparison between a one dimensional PHABSIM and a two- dimensional habitat model involves using the same biological model thereby relating the change in the model performance only to the river hydrodynamics. Ghanem et al. (1996), for example, used the bathymetry data collected for a one- dimensional model to run both one- and two-dimensional models. While others, 26 Leclerc and Lafleur (1997), for example, used the bathymetry data collected for a two-dimensional model to compare the two approaches. The conclusion, however, was the same in all the studies: the two-dimensional habitat model performed better than the one-dimensional model. Advantages provided by two-dimensional models include fewer velocity measurements, easier calibration, and flexibility with respect to modeling under different scenarios. Besides advantages, two-dimensional models were also reported to be very sensitive to channel bathymetry (Bovee, 1996; Tarbet and Hardy, 1996). The accuracy of any two-dimensional model is, therefore, highly dependent on the accuracy of channel bathymetry and on the characteristics of finite element grids used in the model. 2.3 HYDRODYNAMIC MODELS Hydrodynamic models are used to describe the distribution, direction, and magnitude of velocities and depths for different flow conditions. The terms hydraulic and hydrodynamic modeling are used interchangeably within the engineering community. In this dissertation, the term hydrodynamic modeling is used throughout. Based on the equations used and the computational domain, hydrodynamic modeling can be categorized into three main groups: 1) one- dimensional, 2) two-dimensional, and 3) three-dimensional. Each of these groups is briefly discussed in the following sections. 2.3.1 One-dimensional Hydrodynamic Modeling One-dimensional models are the simplest option available for modeling the flow conditions within a river channel. They are best suited for modeling the 27 flow conditions within a river channel that is described by a set of stream cross- sections. The simplest forms of one-dimensional models solve one-dimensional energy equations to compute water surface elevations at each cross-section for steady gradually varied flow conditions. These models can also accommodate momentum equations for rapidly varied flow conditions, such as a hydraulic jump. Slightly more sophisticated one-dimensional models simulate unsteady- state flow conditions in river channels and solve cross-sectional averaged Saint- Venant equations to route the flow hydrograph and compute water surface elevation at each cross-section. For example, HEC-RAS (USACE, 2002b), a commonly used one-dimensional hydrodynamic model, has the capability to perform both steady and unsteady state simulations. A typical one-dimensional model, HEC-RAS in this case, requires geometric and flow data. Basic geometric data consist of the river system schematic and the cross-section data. A river system schematic (Figure 2.5) defines how various reaches are connected and the direction of flow among the reaches. Junction 1 Junction 2 Junction 3 Reach 1 Reach 2 Reach 4 Tributary Tributary Tributary Cross-section Cross-section Cross-section Mcs M s Reach 3 Figure 2.5 River system schematic in HEC-RAS The point of intersection of two reaches is called a stream junction. The cross- section data include ground profiles transverse to the flow direction. Figure 2.6 shows a typical cross-section for HEC-RAS. Each reach on the river system schematic has a unique identifier. Each cross-section on the river system schematic is associated with a reach identifier, and a station identifier (M s ), which specifies distance along the reach. The cross-reach distance along the cross- section is specified by station numbers (M cs ). Looking downstream, the left end of a cross-section has the lowest station number and the right end has the highest station number. 28 M s M cs Figure 2.6 A typical HEC-RAS cross-section The flow data for HEC-RAS consists of flow regime, discharge information, initial conditions and boundary conditions. The flow regime is specified as subcritical, supercritical, or mixed. Discharge information includes at least one flow value along every reach on the river system schematic. The initial and boundary conditions are specified in terms of initial water surface elevations at the upstream and downstream end, flow hydrograph, or discharge-rating curve. Besides geometry and flow data, additional information such as substrate type to calculate Manning?s n is also required. Besides HEC-RAS, MIKE 11 (DHI, 2000) from the Danish Hydraulic Institute and FLDWAV (Fread and Lewis, 1988) from the National Weather Service are also used for one-dimensional modeling. 29 30 One-dimensional models have traditionally been used for flood plain mapping, and are still commonly used for this purpose. The model set-up is easy and the computations are fast. However, they are limited by their ability to model two-dimensional characteristics such as channel meander, velocity currents, flooding in flat urban environments, lateral distribution of flows, etc. Some limitations of one-dimensional models are overcome by using two-dimensional models. 2.3.2 Two-dimensional Hydrodynamic Modeling One-dimensional models work in the longitudinal direction, and two- dimensional hydrodynamic models go one step further by adding the lateral component of river hydrodynamics to the system. Two-dimensional hydrodynamic models solve depth-averaged mass and momentum equations to compute water surface elevations and velocities. In other words, at each point, three items are computed: water depth, and velocities in two directions (x and y). Most two-dimensional models operate under hydrostatic assumption, which means the accelerations in the vertical direction are negligible. The spatial or computational domain of the river system is divided into a set of elements, and for each element, the velocity vectors are assumed to point in the same direction over the entire depth of the water column at any instant of time. The computational domain of two-dimensional hydrodynamic models can be represented by using either (x,y) coordinates or channel oriented curvilinear orthogonal coordinates. In the curvilinear orthogonal coordinate system, the data are plotted with respect to the flow direction. The hydrodynamic models, however, take the data in (x,y) 31 coordinates, and the computations are performed in curvilinear orthogonal coordinates (Hodges and Imberger, 2001). The details of curvilinear orthogonal coordinate system are presented in chapter 4. The two-dimensional models use either the finite-difference method or the finite-element method to solve the energy and momentum equations. The choice of solution method depends on several factors such as the computational speed, model set-up, and shape of computational domain. However, in most instances, finite element methods are usually preferred over finite difference methods due to their ability to model complex geometries in an efficient manner (Jennings, 2003; Rathburn and Wohl, 2003). Commonly used two-dimensional hydrodynamic models, such as RMA2 (Donnell et. al., 2001) and FESWMS (Froehlich, 1992) are based on finite element technique. MIKE 21 (DHI, 2001) from the Danish Hydraulic Institute is an example of a finite-difference model. Typical data requirements for two-dimensional models consist of geometry data, boundary conditions, and calibration data. The geometry data mainly consist of channel bathymetry in the form of finite element mesh specified by (x,y) coordinates at each mesh node (Figure 2.7). The geometry data has to be detailed in order to represent the spatial variations in the channel bed. Generally, the geometry data collected on site are in the form of irregularly arranged points, but the hydrodynamic models require the bathymetry in the form of finite element mesh with the bed elevation specified at each node (Figure 2.7). Fortunately, most models have the capability to interpolate the bathymetry data to get the elevation at each mesh node. The representation of finite element mesh using (x,y) coordinates avoids extra attributes such as station numbers and reach identifiers as used by one-dimensional models. The collection of geometry data for two-dimensional models, however, demands more resources compared to one-dimensional models. Finite element nodes Finite element mesh Bathymetry points Figure 2.7 Finite element mesh for two-dimensional hydrodynamic modeling The most commonly used boundary conditions for two-dimensional models are downstream water surface elevation and upstream flow-rate. Calibration data for two-dimensional models mainly include point measurements of depth and velocity. These measurements are compared with model estimates, and the parameters are adjusted to match the measured values. Additional inputs 32 33 to two-dimensional models consist of eddy viscosity values and roughness coefficients for various substrate types. Two-dimensional models have both advantages and disadvantages compared to one-dimensional models. Advantages include two-dimensional distribution of flow and velocity, and the ability to simulate complicated flow patterns. The main disadvantage of two-dimensional models is that they are not as powerful as one-dimensional models for simulating flow across control structures such as weirs, pumps, tidal gates, etc. This disadvantage of two-dimensional models can be overcome by combining them with one-dimensional models to use the best of both models (Runge and Oslen, 2003). The crucial factor while using two-dimensional models is an accurate description of channel bathymetry. The RMA2 reference manual (Donnell et. al., 2001) suggests that geometry and study design are the most significant factors in the model application (Figure 2.8). Geometry and study design mainly include the description of mesh, which in turn depends on the channel bathymetry. Geometry and Study Design 60% Boundary Conditions 20% Roughness 10% Viscosity 6% 4% Other Figure 2.8 Relative importance of user defined variables in RMA2 performance (From Donnell et. al., 2001) Attempts to address the issue of channel geometry with two-dimensional models have been made both by scientists and engineers. Some of these attempts include refining the mesh resolution to account for small-scale topographic features, such as boulders and woody debris (Crowder and Diplas, 2000; Gilvear et. al., 1999). All these studies conclude that higher resolution produces better model results. Although mesh resolution is important, there are some basic issues that need to be addressed with the channel bathymetry. According to French and Clifford (2000), these issues are related to the quantity and quality of the terrain data and the tools that are used to interpolate these data onto computational meshes. Even though the channel bathymetry plays a significant role in the model output, the procedures by which the terrain data are used to interpolate the bathymetry onto mesh nodes are poorly documented. There are, however, exceptions such as Carter and Shankar (1997) who suggest that kriging algorithms 34 35 work better for interpolating river channel bathymetry. In general, collecting the terrain data, interpolating the data, and developing a suitable finite element mesh are some of the crucial points related to two-dimensional hydrodynamic modeling. 2.3.3 Three-dimensional Hydrodynamic Modeling Three-dimensional models are similar to two-dimensional models, except that the governing equations are not depth-averaged. Besides the two horizontal dimensions (x and y), the vertical dimension (z) is also modeled by introducing a number of layers in the water column. Advantages of three-dimensional models include their ability to model the vertical distribution of flow and velocity, which include stratification, diffusion and dispersion processes. Three-dimensional models, however, require additional computational time. Given the increase in computational speed of computers over past few years, computational time may not be an issue with three-dimensional models in the future. The data requirements for three-dimensional models are similar to two- dimensional models with additional measurements of velocity fields in vertical direction (Fissel et. al., 2002). The issues related to bathymetry data with two- dimensional models also hold true with three-dimensional models. 2.4 GIS AND HYDRODYNAMIC MODELING GIS provides a platform to create, store, analyze, and visualize spatial data in vector and raster form. The geometric data required by all hydrodynamic models can be created and stored inside GIS. One of the basic issues with the GIS data is how to use it to run hydrodynamic simulations. Djokic et. al., (1994) 36 developed a GIS interface to export the GIS data as input to HEC-RAS (then HEC-2), and then import the output results from HEC-RAS to GIS for visualization. This effort was carried forward by Ackerman et. al., (1999), which led to the development of HEC-GeoRAS (USACE, 2002a). HEC-GeoRAS is a GIS extension that allows GIS users to create and export HEC-RAS geometry files such as channel centerline and cross-sections using a digital terrain model. The output from HEC-RAS is also available in GIS format for visualization. Development of HEC-GeoRAS was a significant step forward because it allowed users to visualize and animate results to study flood events for various scenarios (Azagra et. al., 1999; Snead, 2000). It also led to addressing some of the issues associated with floodplain mapping by using GIS capabilities. These issues mainly include the quality of terrain data (Werner, 2001) and the integration of HEC-RAS geometry data with the surface terrain (Tate et. al, 2002). Hydrodynamic modeling packages such as Surface Water Modeling System (SMS) and the MIKE series have GIS modules to create a finite element mesh for two- and three-dimensional models. These modeling packages, however, do not have the ability to process the raw data and deal with issues such as integrating two different terrains to describe the study area. Due to the inability of most models to process raw data, GIS is now widely used to pre-process raw data for two- and three-dimensional models. A typical pre-processing step in GIS to prepare the input data for a two- or three-dimensional model include one of the following: identifying errors in the measured bathymetry and correcting them; creating a digital terrain model based on point measurements; extracting only a 37 small portion of the available terrain dataset; integrating two terrains to describe the study area; and creating stream centerline or introducing additional features such as banklines, levees, etc. Since GIS can handle many kinds of pre-processing tasks, it has become a common data analysis tool for many hydrodynamic models. This is true with hydrologic models as well (Biftu and Gan, 2001). In addition, GIS has been used to integrate hydrologic and hydrodynamic models for studies involving flood forecasting (Anderson, 2000). In addition to supporting different types of data such as raster and vector, GIS also has capabilities to convert the data from one form to another. For example, a raster dataset such as a digital elevation model that has elevation values for every cell can be converted to a vector dataset (points). The resulting vector data will store the corresponding elevation values as attributes. Conversely, points with elevation attributes can be interpolated using spatial interpolation techniques to create a raster grid or a TIN. Similarly, contour lines can be included while creating a raster or a TIN. In addition to points and lines, polygons can be used. In the case of raster datasets, polygons can be used only to define a spatial boundary, while TIN polygons can define three-dimensional objects such as buildings in addition to spatial boundaries. The ability of GIS to handle several kinds of data and the ease with which the data can be processed and used for different studies certainly leads to one question. Is it possible to develop a data model that can store the data inside GIS in a standard format? The data model in question should also have some tools specifically designed to meet the demands of the modeling community and have 38 the ability to communicate with different models. The Center for Research in Water Resources at the University of Texas at Austin and the Environmental Systems Research Institute (ESRI) established a GIS and Water Resources Consortium to answer the above question. A team effort by this consortium led to the development of Arc Hydro (Maidment, 2002). Arc Hydro is a geospatial and temporal data model for water resources that operates within the ArcGIS environment, which can be used for integrating hydrologic and hydrodynamic simulation models. Arc Hydro divides water resources data into five components: network, drainage, channel, hydrography and time series. The channel component of Arc Hydro includes three-dimensional description of river channels in the form of cross-sections and profile lines. The Arc Hydro data model also has a toolset, consisting of four main functions to perform basic watershed analysis. These functions are: 1) Terrain Processing, which deals with basic processing of digital elevation models (DEM) such as filling sinks, burning streams, creating flow direction gird, flow accumulation grid, etc; 2) Watershed Processing, which deals with delineating watersheds and sub- watersheds based on points; 3) Attribute Tools, which assign key attributes to Arc Hydro feature classes; and 4) Network Tools, which deal with geometric networks. Network tools are used to connect lines or points with watersheds and to assign flow direction to geometric network features. 39 The current version of Arc Hydro toolset does not include tools for time series and river channels. In the last two years, the time series component has matured in two parts. The first part deals with retrieving online time series data from USGS and other sites, while the second part deals with sharing or transferring time series data among hydrologic/hydraulic models. The time series tools are currently independent of the main Arc Hydro toolset. The new version of Arc Hydro is likely to have the time series component or the time tools may stay independent as they are now. The river channel component of Arc Hydro is equally important to hydrology compared to other components, but besides HEC- GeoRAS, which can use Arc Hydro channel data as input, Arc Hydro does not have any generic toolset for river channels. The generic geo-processing or editing tools available in GIS can be used to pre-process any kind of data required for river channel studies. These tools, however, are not case specific, and may involve considerable manual editing to get prepare a dataset for a specific application. Therefore, developing a toolset for river channels is necessary to make the GIS capabilities more efficient for river channel applications. Development of a procedure that will eventually lead to populating the river channel component of Arc Hydro with cross-sections and profile-lines is the first step towards a generic toolset for river channels in GIS. 2.5 GIS AND RIVER CHANNEL MORPHOLOGY GIS provides a spatial framework for representing, analyzing, and visualizing earth and its resources. GIS treats earth and its resources as spatial objects, and their geographic locations are described by using Cartesian coordinates (x,y). Although most geographic features can be handled in Cartesian coordinates, some natural features (example, river channels), however, can be better handled in a different coordinate system. One of the shortcomings of using Cartesian coordinate system for handling river channels can be demonstrated by a simple example of calculating distance between two points along a river channel. In Figure 2.9, distance between points P and Q in Cartesian coordinates is d pq , but the actual distance between P and Q is the flow-length, L pq , which is greater than d pq . Using d pq for calculating travel-time, slope, etc along river channels may produce unreliable results. P Q d Y pq L pq X Figure 2.9 Representation of river channel in Cartesian coordinate system As mentioned earlier, bathymetry data for river channels are collected in the form of (x,y,z) points, and these points are then spatially interpolated to create a continuous surface. This is true even for data collected in the form of cross- 40 41 sections. Another issue associated with using Cartesian coordinate system for analyzing river channels is the spatial interpolation of channel bathymetry for creating continuous surfaces (Burroughes and George, 2001; Goff and Nordfjord, 2004). The interpolation techniques used for spatial interpolation do not take into account the flow-oriented morphology of river channels, and therefore produce unsatisfactory results. The issue of Cartesian coordinate system for handling river channels is well known, and therefore, the concept of flow-oriented curvilinear orthogonal coordinate system (s,n) is sometimes used for analyzing river channels (Fukoka and Sayre,1973; Holley and Jirka, 1986). In (s,n) coordinate system (described in detail in section 4.3), s is the distance along and n is the distance across the river channel. The (s,n) coordinate system is used in various river channel studies. Nelson and Dungan (1989), and Johannesson and Parker (1989) used the (s,n) coordinate system to formulate a predictive model for the evolution of meandering rivers. Goff and Nordfjord (2004) used (s,n) coordinate system for spatial interpolation of river channel morphology to conclude that the spatial interpolation in transformed coordinates provide realistic results compared to results in Cartesian coordinates. With regard to mathematical modeling, (s,n) coordinate system offers computational efficiency by mapping a sinuous river form into a rectilinear computational space. For example, Ye and McCorquodale (1997), Sinha et. al. (1998), and Chau and Jiang (2001) used (s,n) coordinates for improving the computational efficiency of numerical modeling in the context of river channels and estuary studies. 42 Since the geospatial framework in GIS uses Cartesian coordinate system, several tools available for spatial analysis in GIS are inadequate to handle the complex morphology of river channels. Therefore, using a coordinate system that is referenced with respect to the flow direction in a river channel is necessary for developing a geospatial structure for river channels in GIS. 2.6 CHANNEL BATHYMETRY OVER LARGE SPATIAL DOMAIN Besides a standardized dataset for storing river channel bathymetry, another problem associated with channel bathymetry is the availability of data over large areas (several hundred miles). The typical cost associated with collecting channel bathymetry using a single-beam depth sounder and GPS (details given in chapter 3) is approximately $1,750 per river mile. The channel bathymetry can also be collected using a multi-beam depth sounder, side scan sonar and GPS (details given in chapter 3), but the cost associated with this technique ($2,500 per river mile) is even higher than the single-beam depth sounder. The minimum cost associated with single-beam depth sounder and multi-beam depth sounder is $6,000 and $15,000, respectively. In addition, the collection of channel bathymetry data requires a considerable amount of time expended by trained personnel, and is therefore limited to short channel reaches. For any kind of study that involves bathymetry data collection, a short representative reach (less than 10 miles) is selected, and the results from this short reach are then applied for the entire study area. For example, the instream flow studies in Texas that use such an approach have study areas that cover several hundred miles of main-stem river segments (Figure 1.1). Similarly, studies 43 involving contaminant transport and soil erosion modeling, which also use local scale studies to make decisions at regional scales, are open to argument (Van Winkle et. al., 1998; Addiscott and Mirza, 1998; Renschler and Harbor; 2002). If additional bathymetry data are available, then these data can be used to verify the decisions that are made based upon the studies on the representative reaches. Therefore, the ability to acquire bathymetry data over large spatial domain is a key to successful regional scale studies. The availability of remotely sensed satellite data, such as LIDAR, has proven to be a very useful for regional scale studies, but their inability to penetrate deep and turbid water limits their use in river channel studies. On one hand, the collection of detailed bathymetry data over large spatial domain using conventional techniques is limited by available resources. On the other hand, techniques such as LIDAR that are used to map the regional topography cannot measure the channel bathymetry due to their inability to penetrate deep and turbid water. One of the options (besides having bigger budgets) for obtaining bathymetry over a large spatial domain is to extrapolate available bathymetry data on short reaches to create a regional description of channel bathymetry. As mentioned earlier, the spatial interpolation schemes do not account for the complex morphology of river channels thus giving unsatisfactory results. Therefore, interpolating available bathymetry to fill the gaps is an issue by itself, and the idea of extrapolation makes the problem more complex. Satisfactory interpolation of bathymetry, however, can be achieved by considering the spatial arrangement of data and channel geometry (Carter and Shankar, 1997; 44 Burroughes and George, 2001; Goff and Nordfjord, 2004). The extrapolation of bathymetry data, on the other hand, requires more than just the spatial arrangement and the channel geometry. River channels form by the interaction of flowing water with the earth surface, and they develop and adjust over time through the combined process of erosion and deposition. Erosion of the riverbed and deposition of the sediments are in turn influenced by the channel planform and the discharge in the channel (Knighton, 1998). Therefore, various factors such as the geomorphology of the river channels, channel planform, geology, discharge, etc. play a major role in shaping the channel bathymetry. Geomorphology involves the study of evolution and adjustment of river channels over time. The traditional approach to the evolution of river channels involves comparison of the present morphology with that recorded by previous measured sources. The change in morphology is then related to field indicators such as vegetation and urbanization. Studies involving cross-sections and topographic maps are localized in extent due to the unavailability of extensive data at a watershed scale (Gregory et. al., 1992). With the availability of aerial photographs, image processing tools and GIS, morphological comparisons are now carried out using aerial photographs and satellite-based remotely sensed data (Werritty and Leys, 1999; Winterbottom, 2000). In addition, satellite remote sensing offers the possibility of studying spatial and temporal variations in geomorphology from the fine scale of mapping pebbles to the regional scale of a single catchment to the global scale of the world?s topography (Yang et. al., 1999; Mertes, 2002; Stein et. al., 2002). 45 The GIS data used for geomorphologic studies, such as aerial photographs and remotely sensed data, are also shared by hydrologists for several applications. For example, remotely sensed data are used in distributed hydrologic modeling where the model is distributed based on the numbers of cells (Schultz, 1993; Biftu and Gan, 2001). As a result of sharing of common datasets and computing tools, hydrology and geomorphology are on a converging path. Besides the traditional approach of using cross-section surveys, the use of GIS has emerged in investigating the relationships between the channel morphology and watershed characteristics (Miller et. al., 1996; Harman et.al, 1999). The morphological variables include channel width, radius of curvature, meander wavelength and mean depth, while the watershed characteristics include drainage area, maximum flow length, stream order, and relief. The ability of GIS to easily compute the above listed variables makes it popular in watershed studies. 2.7 MEANDERING OF RIVER CHANNELS Irrespective of scale or boundary material, most river channels have an inherent tendency to meander. The initiation of meander requires bank retreat, which alternates from one side of the channel to other thus developing a series of bends. According to Callander (1978) periodic deformation of channel bed and the development of sinuous thalweg are necessary precursors to erosion of channel banks. However, the literature available on the physics of channel meandering suggests different theories for the development of meanders. Most theories fall into two broad categories: meandering as a result of interaction between the flow and a mobile channel bed; and meandering due to inherent properties of turbulent flow (Knighton, 1998). Development of alternate bars deflects the flow against opposite banks thus initiating sinuous pattern in a river channel (Callander, 1978). Studies by Nelson and Smith (1989) suggest that such development occurs best in wide, shallow channels with relatively coarse bed material. Experiments by Schumm et. al., (1987) concluded that irrespective of bed material (bedrock or alluvial), the following conditions are necessary for development of a stable meandering pattern: a) the availability of bedload ample enough to develop a point bar; and b) a mechanism which can deposit and stabilize point bars, such as cohesive bank material, bank stabilizing vegetation, etc. Similar studies carried out by Parker and Johannesson (1989), and Seminara and Tubino (1989) support these conclusions. With regard to meandering as a result of turbulent flow, the periodically reversing helicoidal flow influence the pattern of erosion and deposition thus initiating meandering in river channels. According to Einstein and Shen (1964), and Shen and Komura (1968), the helicoidal flow can result in the formation of meandering thalweg and alternating bars, which will eventually lead to channel meandering. Yalin (1992) suggest that turbulent structure of the flow in a river channel is auto-correlated, and similar flow perturbations are approximately separated by equal distance ( ?2 times width), leading to alternate regions of high- speed and low-speed flow and to alternate regions of deposition and erosion. 46 47 An important outcome or an element of river meandering is the flow pattern at the meandering bends. According to Knighton (1998), the main features of that flow pattern are: ? superelevation of the water surface against the concave (outer) bank; (b) a transverse current directed towards the outer bank at the surface and towards the inner bank at the bed to give a secondary circulation additional to the main downstream flow; and (c) a maximum-velocity current which moves from near the inner bank at the bend entrance to near the outer bank at the bend exit, crossing the channel through the zone of greatest curvature?. These flow pattern at the meandering bends also induce cross-stream variations in the boundary shear stress and velocity fields (Dietrich, 1987), which in turn cause and propagate the asymmetries in cross-sectional shapes. Channel meandering, therefore, is a complex phenomenon that involves interaction among flow fields, banks, and bed material. It can be initiated by either ample bank load or helicoidal flow pattern or bank erosion. Due to point bars that initiate the meandering of river channels and the helicoidal flow pattern, the cross-sections at meandering bends undergo deposition at the inner banks and erosion at the outer banks. This process of erosion and deposition in conjunction with variations in boundary shear stresses and velocity fields at the meandering bends creates asymmetric cross-sectional shapes at meandering bends. 2.8 HYDRAULIC GEOMETRY OF RIVER CHANNELS The interaction between the mobile water and the earth surface gives rise to several channel patterns of channel planform. The magnitude of discharge and the type of flow regime influence the channel geometry, which includes both the 48 channel planform and the channel cross-sectional shape. This interdependence between the flow and river geometry is the underlying principle of the hydraulic geometry approach (Leopold and Maddock, 1953). The hydraulic geometry approach assumes that discharge is the dominant independent variable, and the dependent variables (average channel width, average depth, average velocity) are related to it in the form of simple power functions as shown below: W = aQ b (2.1) d = cQ f (2.2) v = kQ m (2.3) where W is average width, d is average depth, , v is average velocity, Q is discharge, and a, c, k, b, f and m are numerical coefficients. Discharge (Q) is the product of mean velocity (v) and cross-sectional area of flow (A). Thus, A = w x d (2.4) Q = w x d x v (2.5) which can also be written as Q = aQ b x cQ f x kQ m (2.6) Thus, b + f + m = 1 (2.7) and, a x c x k =1 (2.8) Hydraulic geometry relationships can be developed at a single cross- section (at-a-station hydraulic geometry) or can be developed across several cross-sections along a river (downstream hydraulic geometry). Since the introduction of the concept, hydraulic geometry relationships have been developed on several rivers over years (Schumm, 1960; Park 1977; Hey 1978: Knighton 1985; Miller et. al., 1996; Harman et. al., 1999). Most studies related to hydraulic geometry relationships are focused on finding and refining the exponents of equations 2.1 to 2.3. Hydraulic geometry relationships are useful because they can be used to predict the dimensions of channel morphology by using the flow data. Rosgen (1994) used hydraulic geometry as one of the criteria for classifying stream channels into different types by studying the effect of cross- sectional shape (width/depth ratio), slope, sinuosity, and meander geometry on hydraulic geometry relationships. Besides the prediction of channel morphology and its dimensions, hydraulic geometry relationships can also be used to understand the dynamics of river systems. For example, Leopold et. al. (1964) developed a m/f ratio using the exponent of hydraulic geometry relationships. In this ratio, m is the rate of increase of velocity with discharge, and f is the rate of increase of depth with discharge. A higher m/f ratio indicates increase in sediment load with discharge and vice versa. Another set of empirical relationships that are known for predicting channel dimensions are the regime equations (Lacey, 1930; Blench 1971). Regime equations in their original form can be given as: R v f 2 73.0 = (2.9) 3/1 47.0 ? ? ? ? ? ? ? ? = f Q R (2.10) 49 50 3/16/5 ? =A 2/1 66.2 QP = 3 S = 26.1 fQ (2.1) (2.12) /56/1 00053.0 fQ ? (2.13) where f is Lacey silt factor, v is velocity in feet/sec, R is hydraulic radius in feet, A is cross-sectional area in square feet, P is wetted perimeter in feet, S is slope, and Q is discharge in cubic feet per second. Regime equations are based on the hypothesis that channels adjust their dimensions (width, depth, slope) until they are in equilibrium with incoming discharge and sediment load. The term ?regime? in the regime equations is understood to mean that the silting and scouring are balanced over time. The main difference between hydraulic geometry and regime equations is the Lacey silt factor, which accounts for sediment load in the regime equations. Although developed mainly for irrigation canals in India, the regime equations have been widely used in natural channels with slight modifications (Simons and Albertson, 1960; Chang, 1980; Savenije, 2003; Eaton et. al., 2004). 2.9 SUMMARY Projects such as instream flow studies require coupling of hydrodynamic models with biological data. Due to their importance in the overall study design, hydrodynamic models are gaining importance in decision making processes related to river channels. The importance of hydrodynamic modeling has also led to the use of two- and three-dimensional models instead of one-dimensional models. In addition, GIS is used to couple the results from hydrodynamic models with other data to support the decision making process. For example, in the case 51 of instream flow studies, GIS is used to couple hydrodynamic modeling results with biological data to support the decision making process. The popularity of hydrodynamic models (especially two- and three- dimensional) in studies related to river channel has led to several investigations involving the testing of model sensitivity to model parameters. Results from these studies show that quality of channel bathymetry data plays a major role in the performance of hydrodynamic models. Raw bathymetry data that are used for hydrodynamic modeling studies are usually pre-processed using GIS. However, GIS currently does not have tools that can take into account the complex morphology of river channels while analyzing the data. In addition, there is no standardized way of representing the three-dimensional description of river channel bathymetry in GIS. Therefore, a procedure that will take into account the orientation of channel morphology to create a standardized three-dimensional description of river channels in GIS needs to be developed. All the hydrodynamic models use the channel bathymetry in the form of polyline features. For example, cross-sections that are input to one-dimensional models are polylines. Similarly, a finite element mesh used for two- and three- dimensional modeling is also a network of polylines. Therefore, the idea of creating a standardized description of river channels in the form of polyline features is worth exploring. A standardized dataset for channel bathymetry is important for storing the data and sharing the data among different users and models without putting too much effort into pre-processing. Arc Hydro already has a channel component to store the data in vector form using cross-sections and 52 profile lines. A procedure needs to be developed that will lead to the extraction of these features from the existing data and store them in Arc Hydro. A standardized description of river bathymetry using polyline features may resolve one issue associated with using channel bathymetry in GIS for hydrodynamic modeling. However, the main issue related to channel bathymetry is the availability of the data. Collection of channel bathymetry data is usually limited to short river reaches due to the considerable amount of time and resources involved. However, the results of the studies from these short reaches are sometimes used to make decisions over large areas. For example, a typical instream flow study is carried out over a reach of not more than 10 km, but the result from these studies are applied to river segments that may be a hundreds kilometers long. It is difficult to verify such decisions over large spatial domains, and therefore, such decisions are open to argument. Considering the limitations associated with the collection of bathymetry data over large spatial domains, a procedure needs to be developed to describe the channel bathymetry in areas with no bathymetric data. This will prompt additional hydrodynamic studies thus providing a better base for making the decisions. Describing the channel bathymetry in areas with no bathymetric data is a complex problem. However, if the information that is available and the factors that shape the channel geometry are taken into consideration, it is possible to formulate a procedure to create the channel description in areas with no bathymetry data. The process of meandering interrelates channel planform with cross-sectional asymmetries. The concept of hydraulic geometry provides a basis 53 for predicting river channel morphology by knowing the discharge in the channel. Therefore, factors such as channel planform and discharge play a vital role in adjusting the channel geometry, and they can be used to create the description of river channels. In addition, the description of channel bathymetry in the form of polyline features can be used to study the trends in the channel bathymetry, and then fit analytical forms to these trends. These analytical forms in combination with information about channel planform and discharge can be used to describe the channel bathymetry in areas with no bathymetry data. 54 Chapter 3 Data Collection 3.1 INTRODUCTION This research mainly uses the channel bathymetry data collected by the Texas Water Development Board (TWDB) along the lower Brazos River in Texas. Additional datasets that are used include the channel bathymetry along the lower Guadalupe River and the Sulphur River in Texas, also collected by TWDB (Figure 3.1). The bathymetry data are supplemented with the Digital Orthophoto Quadrangles (DOQ) and the channel boundaries for all the study areas. This chapter gives a brief description about the study areas, the data collection methodology, and a brief overview of the data. 0 400200 Kilometers Guadalupe Basin Brazos Basin Sulphur Basin Guadalupe River at Rain water Ranch Site 1 Brazos on river Site 2 Brazos on river Sulphur River Site Figure 3.1 Location of Brazos, Guadalupe and Sulphur watersheds in Texas 3.2 BRAZOS RIVER The Brazos River (Figure 3.1) is the largest river basin in the State of Texas with a drainage area in excess of 45,000 square miles. The Brazos River begins in eastern New Mexico near the Texas border, and flows southeast through the center of Texas, discharging into the Gulf of Mexico near Freeport, Texas. The study area (Figure 3.1) is located downstream of the proposed Allens Creek Reservoir (-96 o 5?56?, 29 o 40?59?), which is NE of Wallis County, TX on FM 1093, and south of Simoton County, TX. The data are available for two reaches 55 along the Brazos River. The first reach, which is just downstream of the proposed Allens Creek Reservoir, is about 7.5 kilometers long (Figure 3.2), and the average width of this reach is about 100 meters. The Brazos River, in the vicinity of the first study reach, is deeply incised, and the meandering nature of the river has created numerous oxbow lakes. Brazos River Brazos Basin 040200 Kilometers 021 Kilometers Allens Creek Reservoir (proposed) Study area Figure 3.2 Location of Site 1 along the Brazos River The second reach, which is downstream from Hempstead, Texas, is about 50 kilometers long (Figure 3.3), and the average width along this reach is about 50 meters. Hereafter, the first study reach is referred as Site 1, and the second study reach is referred as Site2. Site 1 and Site 2 are 84 river-miles and 4 river- miles upstream of the Gulf of Mexico, respectively. 56 02010 Kilometers Brazos River Brazos Basin 4 0 400200 Kilometers Allens Creek Reservoir Site Site 1 Site 2 Figure 3.3 Location of Site 2 along the Brazos River 3.3 GUADALUPE RIVER The Guadalupe River originates in Western Kerr County in Texas, and has a drainage area of 6,700 square miles. The study area (Figure 3.4), which is located near Seguin, is 1.2 km long and about 35m wide. 57 Guadalupe River Guadalupe Basin 01050 Kilometers 020100 Meters Seguin Figure 3.4 Location of study area along the Guadalupe River 3.4 SULPHUR RIVER The Sulphur River begins in northeast Texas and eventually flows into the Red River in Arkansas. The study area (Figure 3.5), which is located about 50 km east of the City of Texarkana, is about 1.4 km long and 30m wide. 58 Sulphur River Sulphur Basin 01050 Kilometers 0250125 Meters Texarcana Figure 3.5 Location of study area along the Sulphur River 3.5 DATA COLLECTION The bathymetry data are collected using a boat-mounted acoustic single beam depth sounder which is coupled to a GPS unit (Figure 3.6). The GPS unit is equipped with an antenna that is capable of receiving Differential GPS (DGPS) corrections. The depth sounder sends out a sound pulse in the water, and the pulse gets reflected back to its source as an echo after hitting the channel bottom. The time interval between the initiation of a sound pulse and the echo returned from the channel bed is used to determine the depth of the channel bed from the source (depth sounder). The depth sounder is accurate to approximately ?0.1cm, and depths can be measured in water as shallow as 0.3m. The DGPS, on the other hand, provides the latitude and longitude with ?1m horizontal accuracy. With this setup, which provides latitude, longitude, and water depth every second, the boat can be moved throughout the river reach to collect bathymetry data at any 59 resolution. The measurements are stored on a computer, which is connected to the depth sounder and the GPS. Satellites GPS antenna Depth sounder Boat Computer Reference station Water surface Channel bed Figure 3.6 Bathymetry data collection using depth sounder and GPS Since the depth sounder provides depth relative to the water surface, the output from the depth sounder is subtracted from the water surface elevation to get the channel bathymetry. Therefore, the whole set-up outputs a set of points with location (x,y) and elevation (z) attributes. The bathymetry data at all the study areas are collected using the same set-up. 3.5.1 Multibeam Echo-Sounder Data In addition to single beam depth sounder data, the TWDB has collected the bathymetry for a small section (about 5 km) at Site 2 of Brazos River using a multibeam echo-sounder. Multibeam eco-sounder works on the same principle as the single beam sounder, but instead of capping one location, the multibeam echo- sounder caps several locations at one ping. Effectively, the job of a single beam 60 61 echo-sounder is performed at several different locations on the channel bottom at once (Geyer, 1992). The multibeam echo-sounder gives a very high resolution of bathymetry data (one point every 0.5 meter), but these data are not used in this research for two main reasons. First, the multibeam data were collected only for a small section (about 5 km) of Site2. Second, working with such a high-resolution data (1.5 million data points for a river section 100m wide and 5km long) as a whole demands considerable computing power and time in ArcGIS. For example, it takes more than a minute just to display the data on the computer screen in ArcGIS. Computing power, however, may not be a problem with the newer version (9.x) of ArcGIS. 3.6 DATA DESCRIPTION The raw bathymetry data at all the study areas are in the form of a text file as shown in Table 3.1. These data can be displayed as points inside GIS. Once the data are displayed inside GIS, they can be exported and stored in a geodatabase as a feature class. Figure 3.7 shows a snapshot of the bathymetry data stored in a feature class and the corresponding attribute table. Even though the data are now transformed to GIS features, the raw data shown in Table 3.1 are still preserved in the attribute table as shown in Figure 3.7. The ?ELEV? field in the attribute table stores the elevation. The data can be accessed, selected, and exported by querying the values in the ?ELEV? field. 62 Date Time Latitude Longitude Elevation 512 512 512 512 512 512 512 512 2001 21:22:08 29.658416 -96.033064 10.416 2001 21:22:09 29.658361 -96.033087 10.436 2001 21:22:10 29.658319 -96.033094 10.344 2001 21:22:11 29.658271 -96.033107 10.238 2001 21:22:12 29.658223 -96.033119 10.190 2001 21:22:13 29.658176 -96.033129 10.202 2001 21:22:14 29.658129 -96.033138 10.208 2001 21:22:15 29.658082 -96.033145 10.128 Table 3.1 Raw bathymetry data. The elevation is recorded with respect to the mean sea level 0 30 15 Meters BathymetryPoints 4 a) b) Figure 3.7 a) Bathymetry data as point features in GIS; b) Attribute table linked to the data The quality of bathymetry data, collected as points, can be judged by the number of points measured per square meter (data density). The bathymetry data density at Site 1 along the Brazos River, Guadalupe River and Sulphur River are denser (@5 points/100m 2 ) compared to Site 2 (1 point/100m 2 ) along the Brazos River. A snapshot of data at Site 1 and Site 2 along the Brazos River shows the difference in the data density as depicted in Figure 3.8. The data for Site 1 was collected for detailed two-dimensional hydrodynamic studies, which require finer resolution bathymetry. The data for Site 2 was collected for studying the increase in the intrusion of salt water from the Gulf of Mexico in the lower reaches of Brazos River. The model for studying salt-water intrusion uses a coarse mesh, and does not require a fine resolution bathymetry. 03015 Meters BathymetryPoints 4 Site 1 Site 2 Figure 3.8 Data density description at Site 1 and Site 2 along the Brazos River In addition to data density, there is one more point regarding the quality of the data. A single beam depth sounder cannot measure steep banks, which are common at Site 2. This scenario is illustrated in Figure 3.9. 63 -10 -8 -6 -4 -2 0 0204060 Dista nce a cross cha nne l (m) E l ev at i o n ( m ) Cross- section -10 -8 -6 -4 -2 0 204060 Dista nce a cross cha nne l (m) E l ev at i o n ( m ) Cross- section a) A B b) Figure 3.9 a) Real cross-section with a steep slope; b) Cross-section measured using a single beam depth sounder Figure 3.9(a) shows a cross-section with a steep bank (AB). The depth sounder, however, cannot measure this bank due to the steep slope thus recording the data as shown in Figure 3.9(b). As a result, the dataset for Site 2 at Brazos River has several locations where it does not describe the real shape of the cross-sections. This limitation has to be taken into consideration while using the data. 64 65 Chapter 4 Development of a Geospatial Structure for River Channels 4.1 INTRODUCTION This chapter describes a procedure for developing a standardized three- dimensional dataset of river channel bathymetry. The result is a vector dataset in the form of cross-sections and profile-lines. The main procedure is discussed in the methodology section (section 4.5) and the earlier sections (section 4.2 to 4.4) cover some of the basic concepts that are used in developing the procedure. 4.2 THREE-DIMENSIONAL REPRESENTATION OF POLYLINES To create a geospatial structure for representing river channels using polylines it is necessary to understand how polylines are stored in ArcGIS. In ArcGIS, a polyline is a set of line segments that may or may not be connected (Zeiler, 1999). The polyline discussed in this research refers to a set of connected line segments as shown in Figure 4.1. A standard GIS polyline stores the (x,y) coordinates of its vertices. In ArcGIS, besides these two basic coordinates, a polyline can also store two extra coordinates: m and z (Figure 4.2). Figure 4.1 PolylineZM representation in ArcGIS. The m coordinate is a measure value along the line and the z coordinate stores the elevation. When an ArcGIS polyline stores (x,y,m), it is called a Polyline M feature, and when it stores all four coordinates (x,y,z,m), it is called a Polyline ZM. The three- dimensional polyline that is used in this research is of type Polyline ZM. There are two types of measure coordinates that can be assigned to any polyline: relative or absolute. In the case of relative measure, as used in the National Hydrography Dataset, measure values are assigned at the two end points, and the intermediate values are interpolated between these values based on the length of the polyline. For example, consider a polyline that has a length of 15m (Figure 4.2). The starting and the ending points of this polyline are assigned a measure of zero and 100, respectively. In the case of relative measure, a point on the polyline that is located at a distance of 3m from the starting point has a 66 measure of 20 ? ? = 20100 ? ? ? ? * 15 3 . Similarly, a point that is located at a distance of 6m from the starting point will have measure of 40, and so on. In the case of absolute measure, the polyline has measure values that are assigned with reference to a fixed point along the line. The fixed point that is used as a reference to assign absolute measures is usually the starting or the end point of the line. For example, if the starting point of the same 15m line, discussed in previous paragraph, is used as a reference point, then the end point will have a measure of 15m (Figure 4.2). Original Length (m) 0 3 6 9 12 15 Relative measures 0 20 40 60 80 100 Absolute Measures (m) 0 3 6 9 12 15 Polyline MZ Figure 4.2 Relative and absolute measures in ArcGIS. A point that is from the starting point will have a measure of 3m, a point that is 6m awa is used as a measure of have a measure of 12 do not hav These meas measures are 67 3m away y will have a measure of 6m, and so on. However, if the end point reference point (with a value of 0m), the starting point will have a 15m. In this case, a point that is 3m away from the starting point will m (15 minus 3). In addition, unlike relative measures which e measurement units, absolute measures have measurement units. urement units are the same as that of the line to which absolute assigned. Absolute measures can therefore be assigned in meters or 68 feet or miles or kilometers as necessary. Absolute measures in meters are used in this research. By default, the measures are assigned from upstream end to the downstream end (in the digitized direction), but this can be reversed, if desired. In this research, the measures are assigned from the upstream end to downstream end. 4.3 CURVILINEAR ORTHOGONAL COORDINATE SYSTEM The curvilinear orthogonal (s,n) coordinate system assigns coordinates to the points in a river channel with reference to its centerline such that all the points in the river channel are oriented in the direction of flow (Fukoka and Sayre,1973; Holley and Jirka, 1986). Engineers and geomorphologists often use channel-fitted (s,n) coordinate system to model river channel processes using mathematical models. The (s,n) coordinate system offers computational efficiency by mapping a sinuous river form into a rectilinear computational space thus simplifying the discretization of the model equations (Ye and McCorquodale, 1997; Sinha et. al., 1998; Chau and Jiang, 2001; Wadzuk and Hodges, 2001). Because of its significance in mathematical modeling of river channels, the curvilinear orthogonal coordinate system (s,n) is also used in this research for developing the three-dimensional structure of river channels in a geospatial environment. In the curvilinear orthogonal coordinate system (Figure 4.3), s is the distance along the centerline and n is the perpendicular distance from the centerline. n1 n 2 P Q Centerlin e Bankline s ( s = 0, n =0) P ( s 1 , n 1 ) Q ( s 2 , n 2 ) s 1 s 2 Figure 4.3 Orthogonal curvilinear (s,n) coordinate system for river channels The centerline, which runs in the direction of flow has the s coordinate equal to zero at the beginning (upstream end) of the channel, and is equal to the length of the centerline at the downstream end of the channel. The s coordinate for any point in a channel is always positive. Looking downstream, all the points lying to the left hand side of the centerline have negative n-coordinates, and all the points lying to the right hand side of the centerline have positive n-coordinates. For example, in Figure 4.3, point P has a negative n-coordinate, and point Q has a positive n-coordinate. There are no standard sign conventions for (s,n) coordinate system. Therefore, the sign conventions used in this research may not conform to those used in other studies. In addition, the definition of centerline as a reference for assigning (s,n) coordinates is arbitrary. The thalweg, which is the traditional way of defining a centerline (Chang and Hill, 1976; Lisle, 1987; Oetter et. al., 2004) can be used in the (s,n) coordinate system. Thalweg is flow-independent and does not change for a given bathymetry data. Thalweg, however, has one shortcoming: it is not 69 70 smooth, and usually has abrupt changes in directions along the flow. This shortcoming can be overcome by filtering (removing unwanted segments) the thalweg to get a smoother line. A simple user-defined geometric centerline running through the middle of the channel can also be used as a reference for assigning (s,n) coordinates(Goff and Nordfjord, 2004). However, to define a reasonable geometric centerline, the knowledge of flow (for locating banks) is necessary. Some studies even use one of the banks as a reference to avoid negative n-coordinates (Nelson and Dungan, 1989; Johannesson and Parker, 1989). A thalweg or a centerline or a bank can introduce singularities (intersecting of cross-sections) in a very sinuous channel, and sometimes just any smooth line that will serve the purposes of assigning (s,n) coordinates is used (Wadzuk and Hodges, 2001). To develop a consistent procedure for assigning (s,n) coordinates, a unique reference similar to the origin in the Cartesian coordinate system is required. Since (s,n) coordinates are measured with respect to the centerline, it is necessary to define a unique centerline. Thalweg is used as the centerline in this research. Using the thalweg as a centerline does not provide a unique reference because the channel bed changes over time thus changing the thalweg location. Nevertheless, the thalweg is still a better solution because it is flow independent, and does not change for a given dataset. The curvilinear orthogonal coordinate system can have its origin (s=0, n=0) either at the upstream end of the channel or at the downstream end of the channel. In this research the upstream end is taken as the origin, and the sign conventions discussed in previous paragraph are used for assigning the coordinates. As mentioned earlier, the (s,n) coordinates can also be applied by using a different origin (downstream end) and different sign conventions, if desired. 4.4 FISHNET The term FishNet as used in GIS refers to a map of a three dimensional surface represented as a set of three-dimensional lines (PolylineZM) that are parallel to the x- and y-axes. A FishNet is similar to a wireframe model that is used in computational geometry. The surface FishNet tool available in ArcGIS can be used to create a FishNet whose elevation values are obtained from a raster or a TIN surface (Figure 4.4). High : 12 Low : 4.6 Fishnet Raster grid FishNet FishNet in 3D Raster Figure 4.4 Surface representation in ArcGIS using FishNet It is well known to GIS users that surfaces from raster and TIN are very slow to display or render. Therefore, surfaces such as raster or TINs are sometimes converted and stored as FishNet for faster rendering and visualization. In 71 72 addition, a FishNet also provides easy navigation guidelines for visualizing the three-dimensional surface. 4.5 METHODOLOGY The desired result is a three-dimensional vector dataset for storing river channel bathymetry in the form of cross-sections and profile-lines. The surface FishNet tool available in ArcGIS provides a network of three-dimensional lines, but these lines are parallel to the Cartesian (x,y) coordinate axes and are not oriented in the direction of flow (Figure 4.4). Therefore, the lines parallel to the x- axis are not cross-sections and the lines parallel to the y-axis are not profile-lines. However, the concept of (s,n) coordinate system can be used to convert a regular FishNet into a network of cross-sections and profile lines. When the data are mapped in the (s,n) coordinate system, the flow direction in the river channel is parallel to the s-axis and transverse to the n-axis. So, if the FishNet is generated after the data are mapped in to the (s,n) coordinates, the resulting FishNet will have lines parallel to the s- and n-axes, which are actually parallel and transverse to the flow. Therefore, a FishNet in (s,n) coordinates gives a network of cross-sections (lines parallel to n-axis) and profile-lines (lines parallel to s-axis). If the FishNet generated in (s,n) coordinates is transferred back to the Cartesian coordinates, the orientation of the FishNet lines with respect to the flow direction is still preserved giving a network of cross- sections and profile-lines. The basic methodology for creating a three-dimensional vector dataset thus involves four steps: mapping the bathymetry data (x,y,z points) into the 73 (s,n,z) coordinate system, interpolating the data to create a surface, creating a FishNet, and then transferring the FishNet to Cartesian coordinate system. However, GIS does not have some of the functions necessary to execute the methodology in the above listed four easy steps. For example, GIS does not have tools to assign (s,n,z) coordinates to (x,y,z) bathymetry points. The basic methodology is slightly modified and is executed using following steps: 1. Develop a procedure to identify the thalweg to serve as a reference for assigning (s,n,z) coordinates. 2. Develop a procedure to assign (s,n,z) coordinates and map the data using (s,n,z) coordinates. 3. Interpolate the z values to create a raster grid in (s,n,z) coordinates. 4. Create a FishNet using the raster grid in (s,n,z) coordinates. 5. Transfer the FishNet to Cartesian coordinates to get a network of cross- sections and profile-lines. This methodology is similar to Goff and Nordfjord (2004), and is developed using Site 1 on the Brazos River (Figure 3.2). For the purposes of illustration, the graphics presented in this chapter cover only a small section of the Site 1 (Figure 4.5). b) a) Figure 4.5 a) Site 1 on the Brazos River; b) Small section of Site 1 used for illustrations 4.5.1 Identifying the Thalweg Location A procedure is developed to identify the thalweg location based on three- dimensional discrete bathymetry points. The input data for the procedure are (x,y,z) bathymetry points, an arbitrary centerline (polyline feature) defined by the user, and the boundary of the channel (polygon feature). The initial arbitrary centerline can be a blue line from a digital map hydrography, a centerline interpreted from an aerial image or a centerline digitized over the bathymetry points. However, as explained later, care must be taken to make sure that the individual segments of the arbitrary centerline are at least 100m in length. In this research, the arbitrary centerline digitized over the bathymetry points is used. Similarly, the channel boundary (left and right bank lines) can be digitized from an aerial image. In this research, the channel boundary is defined by using the digital orthophoto quadrangle (DOQ) for Site 1 on the Brazos River. The initial 74 75 boundary digitized using the DOQ, however, did not cover all the bathymetry points because the DOQ is more than five years old and the river changes its path over time. Therefore, the initial boundary was modified slightly to cover all the measured bathymetry points. The thalweg is a line running through the deepest part of a riverbed, and the riverbed is a continuous surface. Therefore, to locate the thalweg using discrete points, it is important to first interpolate the discrete points to create a surface for the channel bed. An accurate bathymetric surface is necessary to define a unique thalweg. As discussed later, creating an accurate bathymetric surface from discrete points using spatial interpolation methods is a challenge by itself. However, it is possible to create a thalweg that is accurate enough to serve the purposes of the curvilinear orthogonal coordinate system. In addition, if the same spatial interpolation method is used, the thalweg identified by the procedure can be assumed to be unique. A flow chart of the operations needed to identify the thalweg is shown in Figure 4.6 and the corresponding steps are illustrated in Figure 4.7. As mentioned earlier, the procedure requires three inputs: the bathymetry points, an arbitrary centerline, and the boundary polygon (banklines). First, the z values of all the bathymetry points are interpolated inside the channel boundary to create a raster grid using the inverse distance weighting method in ArcGIS (Figure 4.7b). (1) Bathymetry points, arbitrary centerline, and channel boundary (3) Is the centerline dense enough? (4) Are there sharp angles in the centerline Filter s out unwanted egments (7) Is it the last vertex of the line? Densify the line (5) Create normal at each vertex (6) Find the deepest point for the normal and move the original vertex to this deepest point (8) Join all the deepest points to get the thalweg (2) Create a raster surface from the points using inverse distance weighting yes No No Yes yes No Figure 4.6 Flow chart for identifying thalweg After a raster grid is created for the channel bed, picking a cell with the smallest z value in each row, and joining these cells will give a thalweg. This process, however, will not produce a smooth thalweg (no abrupt disjoints) desired for the (s,n) coordinate system. Therefore, instead of going through each row in the raster, the procedure goes through the vertices of the input arbitrary line, and 76 77 relocates them to follow the deepest part of the channel bed. Again, if the input arbitrary polyline has very short segments (<25m), the resulting thalweg may not be smooth. On the other hand, if the input arbitrary line has very long segments (200m), the resulting thalweg may be imprecise. By trial and error, a segment length of 100m is found adequate to provide a reasonable number of vertices for locating the thalweg for the study area. The length of individual segments in the input arbitrary line, however, may vary for other channels depending on the width and sinuosity of the channel. If the initial centerline has segments that are too long, then the line is modified such that each segment is approximately 100m long. The initial centerline is also checked for angles between the segments. Any angle between 0-135 degrees or between 225-360 degrees is considered a sharp angle and if two adjoining segments make a sharp angle with each other, then these segments are removed and modified to achieve a smoother initial centerline line. To identify the deepest part along the channel bed, normal lines (lines perpendicular to the centerline) are created at all the vertices of the centerline (Figure 4.7c), which are cross-sections covering the entire width of the channel. Each normal line is a three-dimensional polyline (PolylineZM) which has a vertex for each underlying cell of the raster grid (the channel bed), and the elevations, z, for each vertex are extracted from the underlying raster grid (Figure 4.8). 78 Input Steps 6, 7, 8 Steps 3,4,5 Output Step 2 a) d) e) c) b) 4 6 8 10 12 0 20 40 60 80 Initial location Thalweg Transverse distance (m) El eva t i o n ( m ) Normal lines Initial line Thalweg f) Figure 4.7 Procedure for the thalweg identification in GIS. (a) Input data; (b) Raster surface created using bathymetry points; (c) Normal lines created for the initial line at 15m interval; (d) initial line is moved to the deepest points along all the normals; (e) Final result with initial line changed to thalweg; (f) Vertical profile of a normal showing the location of initial line and the deepest point (thalweg) Therefore, the initial centerline now has one normal line at each of its vertices. To identify the thalweg, and relocate the initial centerline, all the vertices of each normal line are sorted in the ascending order of their corresponding z values. Then the vertex that has the minimum z value is picked as the deepest point along that normal. The deepest points along all the normal lines are similarly found and then joined to form a three-dimensional polyline, which is the thalweg (Figure 4.7d). For purposes of illustration on a small data set, the normal lines shown in Figure 4.7 are created at 15m intervals instead of 100m that is actually used to locate the thalweg on Site 1. Arbitrary line Normal line Properties of normal line Vertex with lowest elevation Figure 4.8 Properties of normal line extracted from the channel bed (raster surface) Depending on the number of segments in the input centerline, the resulting channel thalweg location may vary slightly. This problem can be overcome by providing an input centerline whose segments are longer than 100m and then splitting the line into segments of standardized length, such as 100m, thus giving the same result for different initial centerlines. The procedure to identify the thalweg is developed for meandering channels. The procedure may not work for braided channels (Figure 4.9) unless the channel is separated into parts (A and B in Figure 4.9), and is applied to each part individually. The application of this procedure to braided channels is, however, not tested in this research. 79 Bank lines Thalweg B A Figure 4.9 Braided river channel separated into two active channels A and B by an island After the thalweg is identified, the next step is to use it as a reference for assigning (s,n,z) coordinates to the (x,y,z) bathymetry points. A procedure to assign (s,n,z) coordinates in ArcGIS is described in the next section. 80 81 4.5.2 Procedure to Assign (s,n,z) Coordinates A procedure is developed to assign (s,n,z) coordinates to (x,y,z) bathymetry points using thalweg as the reference line. Hereafter, the centerline refers to the thalweg identified in the previous section (section 4.5.1). First, the centerline is converted to a piecewise Bezier curve to have a common tangent for all the adjacent segments in the centerline. A Bezier curve is a curve determined by a set of control points (vertices of the polyline) in which each point of the curve is calculated from a polynomial function that uses the coordinates of the control points as parameters (Bezier, 1993). Second, the centerline is assigned absolute measures using the upstream end as the reference point. Therefore, each vertex of the centerline is assigned a measure value equal to the flow length from the upstream end of the centerline. The s-coordinate at each vertex is equal to the measure value at that point. The z-coordinate remains unchanged in both (x,y,z) to (s,n,z), therefore, the procedure involves only the (s,n) coordinates. The procedure to assign (s,n) coordinates is described with reference to Figure 4.10. 82 a) d) A B C P P A ? ? B C ? ? C G B ? P F D ? P C A ? ? B ? ? b) e) C P G B ? ? P ? ? B C A c) f) Figure 4.10 Assigning (s,n) coordinates to a bathymetry point P. (a) Point P is a bathymetry point in the proximity of segments AB and BC of the centerline; (b) The association of P with AB or BC is decided based on angles made by P with endpoints of AB and BC; (c) P is closest to BC because it makes acute angles with endpoints of BC; (d) P makes acute angles with both BC and CD. PG is smaller than PF, and therefore P is closest to BC; (e) P cannot be uniquely associated with either AB or BC because it makes one obtuse angle with both segments; (f) Straight line segments are converted to Bezier curves to maintain tangent continuity at each vertex of the centerline 83 The concepts presented in this section remain the same regardless of whether the centerline is a Bezier curve or not. However, to illustrate the importance of Bezier curve, the steps involved in the computation of (s,n) coordinates using a regular centerline with straight segments are presented. Consider a hypothetical bathymetry point, P, in a river channel (Figure 4.10a). The task of assigning (s,n) coordinates to P involves finding the perpendicular distance of P from the centerline (n-coordinate) and the flow distance of P from the upstream end (s- coordinate). Since the centerline is a polyline (a set of segments), it is easier to consider individual segments of the centerline while calculating the (s,n) coordinates. The problem then reduces to the point P and the segment of the centerline that is closest to P. The point P in Figure 4.10a lies in the vicinity of two segments, AB and BC. If the lines connecting P to both ends of a segment form acute angles with the segment, then the point P must be closer to that segment than to any other segment in the polyline. To determine whether P is closer to AB or BC, the angles between P and the lines joining vertices A, B, C are computed. Angle ? (Figure 4.10b) is an acute angle while ? is an obtuse angle, therefore P is not closest to segment AB. Angles ? and angle ? in Figure 4.10b are both acute angles, so P is closest to BC. The perpendicular distance from P to the closest segment of the polyline (segment BC) is PG (Figure 4.10c). This is calculated as: PG = PB Sin ? (4.1) The distance BG (in Figure 4.10c) is BG = PB Cos ? (4.2) 84 The distance PG is the n-coordinate of P, while the s-coordinate is the sum of BG and the s coordinate at vertex B. In some cases, such as the one shown in Figure 4.10d, it is possible that P forms acute angles with more than one segment. In such cases, the segment that gives the shortest perpendicular distance to the point P is selected. For example, in Figure 4.10d, P forms acute angles with both BC and CD. However, the perpendicular distance from P to BC, PG, is smaller than the perpendicular distance from P to CD, PF. Therefore, point P is referenced to BC and not to CD. A problem arises when a bathymetry point lies near the intersection of two segments. In this case, it is difficult to associate that point with either segment because it does not form acute angles with any of them and cannot have unique (s,n) coordinates. This problem is illustrated in Figure 4.10e. In this figure, P does not form acute angles with both AB and BC, and therefore cannot be uniquely associated with either of the segments. This issue arises because the point B has two normals: one perpendicular to AB and the other perpendicular to BC. Any point that lies between these two normals cannot be associated to any of these segments. This case can be accounted for by converting the polyline into a smooth curve, such as a piecewise Bezier curve. Therefore, the centerline is converted to a piecewise Bezier curve before assigning the (s,n) coordinates. The ?smooth? method in ArcObjects is used to approximate the polyline by a Bezier curve. When the polyline is converted into a piecewise Bezier curve, each individual segment of the polyline is converted into a Bezier curve such that at each vertex, the adjoining Bezier curves satisfy the tangent continuity condition (at the vertex, the curves have common tangent line). Figure 4.11 shows a polyline with regular segments, Bezier curve segments, and common tangents at each vertex. After a polyline is converted to a piecewise Bezier curve, the transition is smooth at the intersection of segments, providing a unique normal at any point on the curve, as shown in Figure 4.10f. The s-coordinate of P is equal to the s-coordinate of point B, and the n-coordinate of point P is the normal distance PB defined from the Bezier curve. Polyline Bezier curves Tangent Figure 4.11 Polyline to piecewise Bezier curve conversion As with the thalweg procedure, this procedure for assigning (s,n) coordinates is not applicable to braided river channels (Figure 4.9). In addition, although many possible scenarios are incorporated in assigning (s,n) coordinates, this procedure may still not work when the centerline/thalweg has very sharp 85 angles (Figure 4.12). In such a case, if a point lies in the middle of two sharp segments, even transforming the polyline to a Bezier curve may create problems because the point cannot be uniquely assigned to either segment. However, such an extreme condition was not encountered in this research, and may not be common for simple meandering channels. If such a condition arises, some manual editing may be necessary to remove the sharp angles of the centerline. P Figure 4.12 Polyline with very sharp angles 4.5.3 Data Transformation The procedure described in the previous section is used to assign (s,n) coordinates to the bathymetry data. The (s,n) coordinates are stored in the attribute table of the bathymetry points along with the original attributes as shown in Figure 4.13. 86 New attributes ? (s,n) coordinates Figure 4.13 Attribute table of bathymetry points with (s,n) coordinates The (s,n) coordinates of each point are used to transform the bathymetry data from Cartesian coordinates (Figure 4.14a) into the curvilinear orthogonal coordinate system (Figure 4.14b), where the bathymetry points are oriented along the flow direction. a) b) s o n - + n y s o +n -n s n x Figure 4.14 Data transformation from (x,y,z) to (s,n,z). (a) Bathymetry data in (x,y) coordinates; (b) Bathymetry data in (s,n) coordinates 87 88 The actual process of using ArcGIS to plot the data in the (s,n) coordinate system using the (s,n) coordinates from the attribute table is described in Appendix A. In addition to bathymetry points, the boundary of the channel is also transformed from the Cartesian coordinates to (s,n) coordinates. Similar to a polyline, a polygon is also a set of connected segments in which the starting point of the first segment is the same as the ending point of the last segment (closed polyline). Therefore, in the case of the boundary polygon, the vertices of the polygon are actually transformed from (x,y) coordinates to (s,n) coordinates. 4.5.4 Spatial Interpolation of Bathymetry To create a FishNet, the data should first be converted to a continuous surface (raster). Since the quality of the resulting FishNet is dependent on the interpolated surface, it is important to make sure that the interpolated surface is a close approximation of the real channel bed. In ArcGIS, the spatial interpolation techniques available with the spatial analyst extension are used to create a continuous surface (a raster grid) using discrete points. This section discusses the methods that are available with the spatial analyst extension. Specifically, inverse distance weighting, splines, and ordinary kriging are discussed. In addition, to account for anisotropy in the river channel bathymetry, a method called Elliptical Inverse Distance Weighting is developed in this research. As the name suggests, it is a modified version of inverse distance weighting method. Inverse Distance Weighting Inverse distance weighting (McCoy and Johnston, 2001) assumes that the observations (elevation, rainfall measurement, temperature measurement, etc.) that are close to one another are more alike than those that are further apart. To estimate a value for a new location (z 0 ), which lies in the vicinity of observed values, a search neighborhood around the new location is defined and a weighted average is taken of the observed values within the defined neighborhood (Figure 4.15). The search neighborhood is a circle, and the weights that are used for averaging are a decreasing function of distance. The simplest weighting function is the inverse power (1/d p ), hence the name inverse distance weighting. z0 z1 z2 z4 z3 d2 d3 d4 d3 Figure 4.15 Inverse distance weighting In Figure 4.15, if z 0 is the interpolating point (point with no observation), z 1 , z 2 , z 3 , z 4 are the observed values at a distance of d 1 , d 2 , d 3 , d 4 , respectively from the interpolating point, the estimate using inverse distance weighting is calculated as 0 ?z ? ? = = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = n i p i n i p i i d d z z 1 1 0 1 ? , with n = 4 in this case. (4.3) 89 90 Inverse distance weighting has two parameters: the search neighborhood, and the power (p) used for weighting function. The search neighborhood, which is a circle, can be defined either by specifying the radius of the circle or by specifying the number of observations (n) to be included in the averaging. In the latter case, the search neighborhood will expand until it covers the specified number of observations. The second parameter, power, is changed to increase or decrease the influence of neighboring observations. If the power is zero, the inverse distance weighting collapses to a normal averaging method, and all the observations will have equal influence for estimating a value for the interpolating point. As the power increases, the influence of the farther observations on the interpolating point decreases. If the power is high, only surrounding observations closer to the interpolating point will have influence on the estimated value. In most applications, a power of two is used for calculating weights and accordingly, the procedure is called inverse square distance weighting. Elliptical Inverse Distance Weighting River channels have anisotropic terrain (Figure 4.16) due to which the bed-slope along the flow (?x in Figure 4.16) is smaller than the bed-slope across the flow (?y in Figure 4.16). A A B B x z y Plan view of a river channel Section A-A Section B-B ?x = z/x ?y = z/y z Figure 4.16 River channel geometry As a result, the variations in the bathymetry along the flow are smaller than that are across the flow. The circular search neighborhood used in inverse distance weighting does not take into account this anisotropy. Therefore, the search area has to be modified from a circle to an ellipse. The major axis of the ellipse is parallel to the flow and the minor axis is perpendicular to the flow (Figure 4.17) allowing more observations that lie along the flow to be included in the interpolation. Figure 4.17 shows the same points as shown in Figure 4.15 with an elliptical search neighborhood. 91 z0 z1 z2 z5 z3 d2 d1 d5 d3 d4 z4 a b Figure 4.17 Elliptical Inverse Distance Weighting Besides having an elliptical search neighborhood, it is also necessary to have a scheme for assigning elliptical weights. For example, let a and b be the major and minor axes of the ellipse (search neighborhood) in Figure 4.17. If all the measured bathymetry points to be included in the interpolation are less than or equal to b units away from the interpolating point then there is no difference between regular inverse distance weighting and elliptical inverse distance weighting. However, if some of the measured bathymetry points to be included in the interpolation are more than b units farther from the interpolating point, then these points are assigned weights that are smaller than the weights assigned to the points that are less than b units away from the interpolating point. Therefore, to increase the weights of the points that lie along the flow direction and are farther than b units away from the interpolating point, equation 4.3 is modified as below: ? ? = = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = n i p i n i p i i d d z z 1 1 0 1 ? if d i < b (4.4) 92 ? ? = = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = n i p i i p i d d 1 1 1 ?? n i z z 0 ? if d i > b (4.5) f i i e d =?d (4.6) where In equation 4.6, e f is any positive number (? 1) to account for elliptical distance. Anisotropy in the data is defined as the ratio of variance in the bathymetry along n-axis to the variance in the bathymetry along s-axis. )( )( s n r zVar zVar a = (4.7) The bathymetry data used in this analysis (Site 1 at Brazos River) have an anisotropy of 3, but instead of using a single fixed value, a r can also be used as an variable (Tomczak, 1998). The ratio of major axis to minor axis of the elliptical search neighborhood is kept equal to the anisotropy ratio. Elliptical inverse distance weighting, therefore, has four parameters: the search neighborhood (defined by number of points, n), the power (p e ), elliptical distance factor (e f ), and the anisotropy ratio (a r ). Splines The term spline arises by an analogy with a draftsman?s aid of that name, a thin metal or wooden strip, which is bent elastically so as to pass through certain points on constraints. The curve taken up by a physical spline is the one which minimizes its internal strain energy. Interpolating the data using splines is also based on similar assumptions. The interpolation function should pass through the data points, and at the same time should be as smooth as possible. A spline curve 93 is a parametric, composite (piecewise) curve formed with polynomial functions satisfying specified continuity conditions at the intersection of two adjoining pieces (Dierckx, 1993). A polynomial function has the form f(x) = a m-1 x m-1 + a m-2 x m-2 +?+ a 1 x + a 0 , where x is the independent variable, m is a non-negative integer (degree of the polynomial function), and a m-1 , a m-2 ,?, a 0 are the coefficients of the polynomial function. The most commonly used continuity is C 1 continuity, which means the slope of the curve (the first derivative of the function) is continuous. To achieve C 1 continuity at the ends or knots, the polynomial functions must be of at least degree three (m = 3). Therefore, the splines that are commonly used are cubic splines, and hereafter the term spline refers to a cubic spline. To achieve higher order continuity (C 2 , C 3 , ..), higher degree polynomial functions must be used. Figure 4.18 shows a spline curve defined for a set of five control points. The set of data points that are used to define the curve are called control points, and the points that lie on the curve are called knots. 0 1 2 3 4 5 6 7 8 9 0123456 0 1 2 3 4 5 6 7 8 9 0123456 Knots Control points Spline z x Interpolating Spline Approximating Spline Figure 4.18 Spline Function 94 In the case of interpolating spline that passes strictly through all the data points, the control points are the same as knots. In the case of approximating or smoothing splines where the curve does not necessarily pass through all the data points, knots are different than control points (Figure 4.18). Even with interpolating splines, the condition that the spline should strictly pass through the data points is sometimes relaxed to achieve a smoother curve. This section deals with the discussion of interpolating splines only. A discussion about smoothing splines is presented in section 5.3.3 of chapter 5. Fitting a 2D spline surface to a set of discrete input data points is similar to fitting a thin metal sheet that is constrained not to move at the data points (Mitas and Mitasova, 1988; 1999). The idea is to define a function f(x,y) which passes through the input data points and minimizes the bending energy function given by the following equation : ( )dxdyffffI yyxyxx ?? ? ++= 2 222 )( 2 2 .. f (f cbyaxyxT ++=),( (4.8) where ? is the two-dimensional Euclidian space (area of interest), and is the square of the second-order derivative of the function given by equation 4.9. ),(),(), 1 yxTrrRtyx n j jj += ? = (4.9) (4.10) In equations 4.9 and 4.10, n is the number of input data points, R(r,r j ) represents a basis function, and T(x,y) represents the local trend function. The coefficients t j , a, b, and c are determined by solving a set of linear equations for all the data points. The basis function R(r,r j ) is designed to obtain minimum curvature, and can be 95 expressed in different forms. When R(r, r j ) is given by equation 4.11, equations 4.8 and 4.9 lead to a thin plate spline function. 96 )log(*),( 2 jj ddrrR j = ) (4.1) ()( 2 j 2 jj yyxxd ?+?=where is the distance from any point (x,y) to the j th point (x j ,y j ). Thin plate spline function works well with gradually varying surfaces such as a flat sloping terrain with no spikes. However, with data that are not gradually varying, the plate stiffness causes the function to overshoot in regions where the data points have large gradients. This is undesirable especially with channel bathymetry, which has lots of noise associated with the data. The stiffness of a thin plate spline can be relaxed by adding the first order derivatives into the energy equation (equation 4.12) and the resulting function is called a thin plate spline with tension. ( ) ( ){ }dxdyffffffI yyxyxxyx ?? ? ++++= 2 222222 2)( ? T( (4.12) The trend function T(x,y) for tension spline is given by equation 4.13, and the basis function is defined by equation 4.14. (4.13) ayx =), () ? ? ? ? ? ? ++? ? ? ? ? ? ?= ? ? ?? j j j dKc d rrR 0 2 2 ln 2 1 ),( (4.14) Where, parameter ? defines the weight of tension, and it can tune the surface from a stiff plate (? = 0, low tension) into an elastic membrane (high tension), K 0 is the modified Bessel function of zero order, and c is a constant equal to 0.577215. Another version of thin plate spline is the regularized spline, which improves the analytical properties of thin plate spline by adding third and higher order derivatives into the energy function (equation 4.15). ( ) ( ){ } ?? ? ++++= 2 ....2)( 22222 dxdyfffffI xxxyyxyxx ? (4.15) The trend function for regularized spline is same as the thin plate spline (equation 4.10), and the basis function is given by equation 4.16. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ++? ? ? ? ? ? + ? ? ? ? ? ? ?+? ? ? ? ? ? = ?? ? ?? 2 ln1 2 ln 42 1 ),( 0 2 2 jjjj j d c d Kc dd rrR (4.16) The parameter ? defines the weight of the third derivatives in the energy function with higher values providing a smoother surface. The regularized spline in ArcGIS, however, includes only the third order derivative. ArcGIS has the option of using either the thin plate spline with tension or regularized spline while using the spline function for spatial interpolation. Besides ? and ?, each method has the number of points (n) as the second parameter for interpolation. Ordinary Kriging Oridinary kriging, which is a geostatistical approach to interpolation, is similar to inverse distance weighting in that it assigns weights to the observed values to estimate a value for an arbitrary location usually over a grid of locations (Isaaks and Srivastava, 1989). However, unlike inverse distance weighting, the weights are assigned based not only on the distance between the prediction location and the observation points, but also on the spatial correlation among the observation points. The spatial correlation is characterized through the use of the 97 semi-variogram model, which provides a measure of variance as a function of distance between the observation points. The semi-variance (? ij ) between two points, located at (x i ,y i ) and (x j ,y j ), with observations z i and z j , respectively is given as ( ) 2 )( 2 ji ij zzE ? =? ()() (4.17) and the distance between these two points is given as 22 jijiij yyxxd ???= (4.18) The kriging procedure to predict a value (z 0 ) at any location using measured values at known locations is described below. Let (x 1 ,y 1 ), (x 2 ,y 2 ), ?, (x n , y n ) be the locations with known values z 1 , z 2 ,?, z n , respectively. The problem is to predict z 0 at (x 0 ,y 0 ). Similar to inverse distance weighting, the kriging procedure also interpolates the data based on the number of observations in a circular neighborhood (Figure 4.19). zo z1 z3 z2 d01 d03 d12 d02 d23 d13 Figure 4.19 Ordinary kriging The following steps are involved in kriging: 98 99 1. Calculate the distance (d ij ) between each pair of known locations using equation (4.18). 2. Calculate the semivariance (? ij ) for each pair using equation (4.17). 3. Plot the semivariogram. A semivariogram is a plot of semivariance versus the distance (? ij versus d ij ) for each pair with semivariance as the ordinate (Figure 4.20). Generally, with real datasets that have thousands of observation points, a procedure called binning is employed. As the number of the observation points increase, the number of pairs of locations also increases rapidly making the computations unmanageable. In the binning procedure, all the pairs that lie in a specified range of distance are grouped together (in a bin) and the average values of d ij and ? ij for each bin are used for computation. 4. Fit a function to the semivariogram plot. This process is called semivariogram modeling. Semivariogram modeling is similar to fitting a least-squares line in regression analysis. Any type of function that can serve as a reasonable model can be used to fit the semivariogram. For example, a spherical type that rises at first and then levels off for larger distances beyond a certain range is commonly used (Figure 4.20). Three terms are important in semivariogram modeling: the range, sill, and nugget. The distance from the origin to a point where the semivariogram first flattens out is called as the range. 0 0.5 1 1.5 2 2.5 0 50 100 150 200 250 300 d ? Range Sill Nugget Semi-variogram model Empirical Semi-variogram Figure 4.20 Semivariogram plot All the locations that are separated by distances shorter than the range are spatially correlated, whereas the locations further apart are not (spatially independent). On a semivariogram plot, the range is measured along the horizontal (x) axis. When a semivariogram model attains the range, the corresponding value on the y-axis is called as the sill. The value along the y- axis at a point where the semivariogram intercepts with the y-axis is called as the nugget. At zero separation distance (d ij = 0), the semi-variance is zero. However, measurements tend to vary at infinitesimal separation distances. The nugget, also referred to as nugget effect, can be attributed to measurement errors or to spatial variations at micro-scales smaller than the sampling distances. The sill minus the nugget is referred to as the partial sill. The semivariance for each pair is again calculated, but this time, the fitted semivariogram model is used for calculations. For example, for a spherical model, the expression is 100 ? ? ? ? ? ? ? ? ? ? ? ? ?+= 0 22 ) rr s ?? ?? ?? ?? 3 13 ( hh h? sh for 0 ? h ? ? r, and (4.19a) ??? += 0 0 ?z ? = = N i ii zz 1 0 ? ? 0 ? 1 1 = ? = N i i ? () for h > ? r (4.19b) where, ? 0 is the nugget, ? s is the sill value, h is the distance between the two locations (d ij ), and ? r is the range of the model. Using expression (4.19), the semivariance (? ij ) for each pair is calculated. 5. Calculate the weights for interpolation. As mentioned earlier, kriging assigns weights to the observed values to estimate/predict a value for a given location (prediction location). Therefore, the prediction, , for z 0 at any location (x 0 ,y 0 ) is calculated as (4.20) Where z i is the observed value at a location (x i , y i ) and ? i is the weight associated with that location. The weights are calculated such that the difference between the true value (z 0 ) and the predicted value ( z ) is as small as possible subject to . That is, minimize the statistical expectation of the following expression, ( ) 2 00 2 ?zzE ?=? ? ? ? ? ? ? ? ? ? ? ?= ? 2 0 2 n ii zzE ?? 1 1 = = N i i ? (4.21) ? ? ? ? ?? =1i (4.22) Subject to: ? (4.23) 101 102 where, ? 2 is the variance. The condition that the sum of the weights, ? i , is equal to one makes the predictor unbiased. Hence, ordinary kriging is called the best linear unbiased estimator. Equations 4.22 and 4.23 are a constrained minimization problem which can be converted to an unconstrained problem using a Lagrange multiplier (?). Therefore, the problem transforms to minimizing () ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= ?? == 1, 1 2 1 0 N i i n i iii zzEL ????? ()[] (4.24) The minimization problem can then be solved by partially differentiating L with respect to ? i and ?, and equating the partial derivatives to zero. 0 , = ? ? i i L ? ?? ()[] (4.25) 0 , = ? ? ? ?? i L ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 . . . . . . * 01...1 1 .... .... .... 1... 10 101 1 111 ? ? ? ? ? ?? ?? nnnn n (4.26) The solution to this problem in the matrix form (Jenson et. al., 1997) can be written as ? ? = g (4.27) 103 1? ?=? 0 ? ? ? ? ? ? ? ? ? = ? = 1 1 N i i ? The weights are then determine by solving g 6. After the weights are determined, z can be calculated using equation 4.20. The kriging procedure without imposing the constraints on the weights is called as simple kriging. In this research, ordinary kriging is used. Ordinary kriging that can also take into account the directional influences in the data (anisotropy) is called anisotropic ordinary kriging (Eriksson and Siska, 2000). Unlike simple ordinary kriging, anisotropic kriging uses a two-dimensional semivariogram model (semivariogram envelope) to calculate kriging weights (Figure 4.21). In other words, simple ordinary kriging uses only one model for fitting an empirical semivariogram while anisotropic kriging uses a series of models for different search directions. For example, if a bathymetry point to be used in interpolation is located at an angle of 90 degrees from the interpolating point, a semivariogram corresponding to that direction is used for calculating the weight, and if the interpolating point makes a different angle then a different semivariogram model from the semivariogram envelope is used. Semi-variogram envelope Semi-variogram model for 90 degrees (m) (m 2 ) Figure 4.21 Semivariogram envelope for anisotropic kriging The number of semivariogram models in a semivariogram envelope depends on the variability in the data. Figure 4.21 has 15 semivariogram models for directions ranging from zero degree to 90 degrees. The lowermost model corresponds to 90 degrees and the uppermost model corresponds to zero degrees. Like elliptical inverse distance weighting, anisotropic kriging also uses an elliptical search neighborhood. The length of minor axis of the elliptical search neighborhood is equal to the range of the uppermost semivariogram model in the semivariogram envelope, and the length of the major axis is equal to the range of the lowermost semivariogram model. The geostatistical analyst in ArcGIS (Johnston et. al, 2001) also has the capability to divide the elliptical search neighborhood into different sectors and specify the minimum number of points to be included for interpolation from each sector (Figure 4.22). 104 a) Ellipse with one sector b) Ellipse with four sectors c) Ellipse with four sectors d) Ellipse with eight sectors a d c b a dc ba g f d e cb a h Figure 4.22 Different search neighborhoods for anisotropic kriging In this research, an ellipse with four diagonal sectors (Figure 4.22c) is used for anisotropic kriging. This means that the search neighborhood (ellipse) is defined by the user to include a minimum number of points in each a, b, c and d sectors of Figure 4.22c while performing the interpolation. Approach for comparing spatial interpolation methods 105 i z? To compare different interpolation methods, the original dataset is split into two sub-datasets: the training dataset and the test dataset. The training dataset is used to create the surface, and the performance of the interpolation method is evaluated by comparing the interpolated values with the measured values in the test dataset using the root mean square error (RMSE). If there are n data points in the test dataset with observed values z i , i = 1,2,?, n, and are the corresponding values estimated by the interpolation method, RMSE is calculated as () ? = ?= n i ii zz n RMSE 1 2 ? 1 (4.28) There are two ways to create a test dataset. One way is to select random points from the original dataset, and another way is to select the points manually. Random selection works well with regularly spaced data such as LIDAR in which the removal of one point creates a gap at that location. However, in the case of channel bathymetry data which are not irregularly spaced (Figure 4.23), random selection does not provide the best training dataset. Figure 4.23 Bathymetry points for creating training and test dataset In the case of channel bathymetry that is collected using a single beam depth sounder, removal of points where the data are clustered together does not make much difference. Instead the focus should be on those areas where the data are not collected. For example, in Figure 4.23, removing few random points from the green circle does not significantly help to study how the interpolation is working in the area with red circle, which is the main area of concern with regard to spatial interpolation. In this research, the data are manually separated to create the training and test datasets. All the data points that lie along the direction of flow are selected to create areas similar to the red circle shown in Figure 4.23. In addition, the bathymetry at Site 1 along the Brazos River is fairly uniform; all the data lie along 106 one big meandering bend. Therefore, only a representative area (Figure 4.5b) is selected for comparison purposes, and the results obtained from this small area are then verified for the entire dataset (Figure 4.5a). Figure 4.24 shows the original dataset of the representative area, the training dataset, and the test dataset in the (s,n,z) coordinate system. Representative Dataset Training Dataset Test Dataset Figure 4.24 Training and test datasets for spatial interpolation Comparison of Spatial Interpolation Methods The interpolation methods discussed earlier, namely, inverse distance weighting, elliptical inverse distance weighting, kriging (ordinary and anisotropic), and spline (regularized and tension) are compared using the RMSE values. In addition, for each method, the influence of different parameters is tested to arrive at optimum values. The goal of this exercise is to come up with an interpolation method that gives the least RMSE so that it can be used to create a surface (channel bed) from the bathymetry points. Each method is discussed individually and the results are summarized in the end to arrive at a conclusion. Inverse Distance Weighting Inverse distance weighting has two parameters: number of points (n) and the power (p). First, an analysis is done to verify the optimum power for a fixed search neighborhood, which is specified by using eight points. As shown in 107 Figure 4.25, a power of two gave the least RMSE which is not surprising given that this is the default value, and is widely used in inverse distance weighting. 0.7 0.7 0.7 0.7 0.768 R M SE ( m ) 60 62 64 66 024681012 Power (p) Figure 4.25 Test to find an optimum power for inverse distance weighting Second, the search neighborhood which is specified by the number of points is tested. Figure 4.26 shows the results. 0. 0. 0.70 0.75 0.80 R M SE ( m ) 60 65 20406080 Number of points (n) Figure 4.26 Test to find an optimum search neighborhood for inverse distance weighting According to the results in Figure 4.26, the model performs better for more data points (bigger search neighborhood). The RMSE for the worst-case scenario (small search neighborhood with five points) is 0.8m, and for the best- 108 case scenario (search neighborhood with fifty points) is 0.65m thus giving an overall improvement of 0.15m in model results. After a gradual decreasing trend, the RMSE stays fairly constant for more than 50 points. The results make sense given the nature of the training dataset (Figure 4.24). When the interpolating points are located close to the boundary, the search neighborhood is more influenced by the data along the boundary, giving unrealistic results. As the search neighborhood expands other data points that do not lie along the boundary are also included giving better results. It is important to notice the influence of search neighborhood on the interpolation scheme. Generally, the inverse distance interpolation scheme is applied by using default parameters. The default value for the number of points is 12 in ArcGIS which gives a RMSE of 0.75m (0.10m higher than the best result of 0.65m). Elliptical Inverse Distance Weighting Elliptical inverse distance weighting has four parameters: number of points (n), power (p e ), elliptical distance factor (e f ), and anisotropy ratio (a r ). The power is not tested again as the weighting function ? ? ? ? ? ? p d 1 for elliptical inverse distance weighting method is the same as that of the inverse distance weighting method. First, the search neighborhood is tested with a r = 3 (anisotropy ratio in the data) and e f = 1 (simplest case). The results are shown in Figure 4.27. 109 0. 0. 0. 0. 0.53 R M SE ( m ) 48 49 50 0.51 52 0 20406080 Number of points (n) Figure 4.27 Test to find an optimum search neighborhood for elliptical inverse distance weighting As expected, the results are significantly better than the inverse distance weighting scheme. The improved results suggest that the anisotropy in the data should be taken into consideration while using the spatial interpolation schemes. The result for the worst-case scenario (0.52m) is better than even the best-case scenario of the regular inverse distance weighting approach (0.65m). The RMSE for the best-case scenario (20 points) is 0.49m, which is only (0.52-0.49) 0.03m lower than the worst-case scenario. Unlike regular inverse distance weighting, the RMSE for elliptical inverse distance weighting increases for more than twenty points. The elliptical search neighborhood expands more in the direction of flow. Therefore, sometimes points that are too far away from the interpolating point are included in the interpolation causing higher RMSE. To overcome this negative influence of points that are farther away from the interpolating point, the length of the major axis of the elliptical search neighborhood is restricted to be less than 100m (50m on each side of the interpolating point in the direction of flow). As reported in kriging 110 technique later (Figure 4.34), the bathymetry points are spatially correlated for a distance of up to 50m in the direction flow as observed with the semivariogram models. Therefore, a length of 100m is used to restrict the major axis of the elliptical search neighborhood. Even if the search neighborhood is defined to include 50 points, the procedure includes only those points that are enclosed by an ellipse with a major axis of 100m. With n = 20, and e f = 1, the sensitivity of anisotropy ratio is tested, and the results are shown in Figure 4.28. An anisotropy ratio of 7 gives the lowest RMSE. 0. 0. 0. 0.51 0.52 R M SE ( m ) 48 49 50 024681012 Anisotropy Ratio (a r ) Figure 4.28 Test to find an optimum anisotropy ratio for elliptical inverse distance weighting Finally, with n = 20 and a r = 7, the sensitivity of e f is tested. In order to give higher weights to the points that are farther than the minor axis of the elliptical search neighborhood (Figure 4.17), generally e f is taken equal to anisotropy ratio (Tomczak, 1998). The RMSE for e f = 7 is 0.5015m which is higher than the RMSE (0.4832m) with e f = 1 (Figure 4.28). A possible explanation for such a behavior may be that a value of e f = a r or higher brings distant points too close to the interpolating point thus overshadowing the influence of local bathymetry 111 112 bd i =? points. Values of e f lower than a r produced lower RMSE, and by trial and error it is found that e f = d i /b (d i is the distance between the interpolating point and observed bathymetry points, and b is the length of minor axis of the elliptical search neighborhood) gives the least RMSE (0.4743m). The expression e f = d i /b changes equation 4.6 to . Regularized Spline Regularized spline also has two parameters: the number of points (n) and ?. The parameter ? defines weight of the third derivatives in the energy function (equations 4.15 and 4.16) with higher values providing a smoother surface. For a fixed number of points (12), the influence of ? on the regularized spline is tested (Figure 4.29). The RMSE is the lowest for ? = 0, it then increases with ?, and finally starts to fall down for higher values. The lowest RMSE for ? = 0 is 3.33m, which is very high compared to inverse distance weighting and elliptical inverse distance weighting. 3 4 5 6 RM S E ( m ) 0 0.2 0.4 0.6 0.8 1 Weight (?) Figure 4.29 Test to find an optimum value of ? for regularized spline For ? = 0, the RMSE does reduce significantly as the number of points included in the interpolation increase (Figure 4.30). The lowest RMSE for regularized spline is found to be 1.74m with ? = 0 and n = 50 points. Overall, the RMSE results for the regularized spline are too high compared to other methods. 0 2 4 6 RM S E ( m ) 102030405060 Number of points (n) Figure 4.30 Test to find an optimum number of points for regularized spline with ? = 0 Tension Spline Tension spline also has two parameters: the number of points (n) and ?. In this case, however, the parameter ? defines the weight of tension. The change in tension tunes the surface from a stiff plate (low tension) into an elastic membrane (high tension). For a fixed number of points (12), ? in tension spline does seem to work well compared to ? of regularized spline (Figure 4.31). However, ? has to be significantly high (>25) for the tension spline to work well. The lowest RMSE obtained for the tension spline is 0.65m with ? = 50. 113 0.6 65 0.7 0 50 100 150 200 250 Weight (?) 0. 0.75 RM S E ( m ) Figure 4.31 Test to find an optimum value of ? for tension spline The number of points reduces the RMSE from 0.67m for five points to 0.54m for 50 points (Figure 4.32). The tension spline therefore performs better than inverse distance weighting and regularized spline. 0.50 0.55 0.60 0 10 2030405060 Number of points (n) 0.65 0.70 RM S E ( m ) Figure 4.32 Test to find an optimum number of points for tension spline For the tension spline, ? = 0, which is not shown in Figure 4.31, gave the highest RMSE of 5.72. A value of zero for ? makes the tension spline to behave like a regular thin plate spline and the steep gradients create spikes in the interpolated surface. These spikes are removed by introduction of ? to create a smooth surface (Figure 4.33). 114 Tension spline ? = 0 Tension spline ? = 50 Figure 4.33 Effect of tension parameter on spline interpolation Ordinary Kriging The most important aspect of kriging techniques is the description of semivariogram model to estimate the weights used for interpolation. The geostatistical analyst extension in ArcGIS has an interface to model the semivariogram. A spherical semivariogram model is used, and the optimized parameters of the model (range, sill, nugget) as estimated by the interface seemed to fit the data well. An example semivariogram modeled by ArcGIS is shown in Figure 4.34. 115 (m) (m 2 ) Figure 4.34 A spherical semivariogram model with range, sill, and nugget equal to 55m, 2.4, and 0.094, respectively Since the geostatistical analyst in ArcGIS optimizes the semivariogram model, the only parameter that remains to be tested with the ordinary kriging is the search neighborhood, which is specified by the number of points. The RMSE produced by ordinary kriging for different number of points is shown in Figure 4.35. 0. 0. 0. 0.51 0.53 0.55 RM S E ( m ) 45 47 49 0 10 2030405060 Number of points (n) Figure 4.35 Test to find an optimum number of points for ordinary kriging The lowest RMSE produced by ordinary kriging is 0.48m with 40 points. Similar to elliptical inverse distance weighting, there is not much difference 116 among the RMSE results produced by ordinary kriging for different number of points. The difference between the RMSE produced by the best-case scenario and the worst-case scenario is only (0.53-0.48) 0.05m. Anisotropic Kriging As mentioned earlier, anisotropic kriging creates a semivariogram envelope (Figure 4.36), which is a series of semivariogram models for different directions. Depending on the angle any bathymetry point makes with the interpolating point, a semivariogram model corresponding to the same direction from the semivariogram envelope is used for calculating the kriging predictions. Like elliptical inverse distance weighting scheme, the anisotropic kriging uses elliptical search neighborhood for interpolation. The major and minor axes of the elliptical search neighborhood are adjusted depending upon the anisotropy in the data. Semi-variogram envelope Semi-variogram model for 90 degrees (m) (m 2 ) Figure 4.36 Semivariogram envelope for anisotropic kriging 117 Similar to ordinary kriging, the geostatistical analyst in ArcGIS optimizes the semivariogram envelope and the anisotropy ratio. The only parameter that remains to be tested is the search neighborhood, which is specified by the number of points. The RMSE produced by anisotropic kriging for different search neighborhoods defined in terms of the number of points is shown in Figure 4.37. Anisotropic kriging outperformed any other interpolation method by producing the lowest RMSE even by using only five number points. The lowest RMSE produced by anisotropic kriging is 0.36m for 15 points. 0. 0. 0.40 0.42 RM S E (m ) 36 38 0 102030405060 Number of Points (n) Figure 4.37 Test to find an optimum number of points for anisotropic kriging Summary of Results from Spatial Interpolation Methods Figure 4.38 summarizes the results from different spatial interpolation schemes in one chart. The RMSE results for regularized spline are too high (>1.5m) and are not shown in the figure. 118 0.3 0.4 0.5 0.6 0.7 0.8 5 10 15204050 Number of Points (n) 0.9 RM S E (m ) IDW EIDW T. Spline Kriging An. Kriging Figure 4.38 Summary of RMSE results for all spatial interpolation methods The anisotropic kriging performed the best among all the methods examined in this research. Detrending of Data Detrending refers to removal of trends from the data. River channel bathymetry has trends in directions along the flow and across the flow (Figure 4.16). It is reported in the literature (Carter and Shankar, 1997) that trends in the data introduce bias, and these trends must be removed while interpolating the data. Anisotropic consideration does take into account these trends, but methods that do not take anisotropy into consideration are affected by the trends. Therefore, a final step in the comparison of spatial interpolation techniques would be to detrend the data, and then compare the RMSE results. Following steps are involved: 119 120 1. Define the trends in the bathymetry data using analytical functions. Chapter five deals with fitting analytical functions to bathymetry data in detail. For the purposes of spatial interpolation, the geostatistical analyst in ArcGIS is used to fit third-degree polynomial functions to the trends in the bathymetry data. 2. Subtract the trend from the measured bathymetry to get residual errors. 3. Perform spatial interpolation on residual errors with optimum parameters, and then add the trend to the interpolated surface to get the final result. Table 4.1 shows the comparison between the results obtained from different spatial interpolation methods before and after detrending. Interpolation Method RMSE Before Detrending (m) RMSE After Detrending (m) Rank Based on RMSE Inverse Distance Weighting 0.65 0.76 5 Elliptical Inverse Distance Weighting 0.47 0.47 2 Regularized Spline 1.74 3.41 6 Tension Spline 0.54 0.82 4 Ordinary Kriging 0.48 0.51 3 Anisotropic Kriging 0.36 0.38 1 Table 4.1: Ranking of spatial interpolation methods for the sample dataset on along Site 1 of Brazos River As seen in Table 4.1, the RMSE increased for all the methods except for elliptical inverse distance weighting. These results contradict the conclusions of Carter and Shankar (1997) who suggest detrending of data before applying the interpolation scheme. Regardless, anisotropic kriging outperformed all other methods with a lowest RMSE of 0.38m. Based on the results before detrending, the interpolation 121 schemes can be ranked with anisotropic kriging being the first (lowest RMSE), and regularized spline being the last (highest RMSE). Spatial Interpolation of the Entire Dataset So far, all the results are obtained for a representative dataset at Site 1 (Figure 4.5b). In this section, the results obtained by applying the spatial interpolation methods to the entire dataset of Site 1 on Brazos River (Figure 4.5a) are presented. The entire dataset was also split into training dataset and test dataset in the same way as described earlier (Figure 4.24). About 25 percent of the data points that lie along the flow are removed to form the test dataset. Only the optimum parameters obtained on the representative dataset are used to verify if similar results can be obtained for the entire dataset. In other words, it is hypothesized that the dataset used in earlier analysis is representative of the entire dataset. The training dataset for the entire area is interpolated without detrending, and RMSE is calculated for each procedure using the test dataset. The RMSE results are summarized in Table 4.2. Spatial Interpolation Method RMSE (m) Rank Inverse Distance Weighting 0.53 5 Elliptical Inverse Distance Weighting 0.32 2 Regularized Spline 0.59 6 Tension Spline 0.45 4 Ordinary Kriging 0.44 3 Anisotropic Kriging 0.31 1 Table 4.2: Ranking of spatial interpolation methods for the entire dataset along Site 1 of Brazos River 122 For the entire dataset, anisotropic kriging performed the best with a lowest RMSE of 0.31m. In addition, based on the RMSE, all the methods have the same ranks as established using the representative dataset. This justifies the hypothesis that the representative dataset used for the initial analysis is indeed a representative sample. Based on the spatial interpolation analysis, the anisotropic kriging performed better than any other techniques examined in this research. The elliptical inverse distance scheme, which also takes into account the anisotropy in the bathymetry data performed second next to anisotropic kriging. Compared to anisotropic kriging, the RMSE produced by elliptical inverse distance weighting scheme is only one centimeter higher for the entire dataset. Therefore, elliptical inverse distance weighting is a promising method given the complexity of the kriging procedure. However, in this research, anisotropic kriging is used to create an interpolated surface for generating the FishNet discussed in the next section. 4.5.5 Creating a FishNet The FishNet described in section 4.4 can be created from the interpolated surface by using the FishNet tool (Figure 4.39) available in ArcGIS. Figure 4.39 FishNet tool interface The FishNet tool has two main inputs: the surface from which the FishNet has to be generated, and the spacing between the lines in the FishNet. The interpolated surface created in the previous section is used as an input to the FishNet tool. The spacing between the lines depends on the amount of details required to be incorporated into the FishNet. A FishNet with smaller spacing captures more details compared to a FishNet with larger spacing. Figure 4.40 shows a FishNet with 15m spacing overlaid on a raster surface. 123 n s Figure 4.40 FishNet generated from a raster surface. The FishNet created in the Cartesian coordinate system has the lines parallel to the x- and y-axes (Figure 4.4). Similarly, the FishNet created in the (s,n) coordinate system has the lines parallel to the s- and n-axes. Since s- and n- axes are oriented in the direction of flow and across the flow, respectively, the FishNet in (s,n) coordinates gives a network of cross-sections and profile-lines. 4.5.6 Transforming the FishNet To get a three-dimensional description of river bathymetry using cross- sections and profile-lines in the Cartesian coordinate system, the FishNet created in (s,n) coordinates is transformed back to the (x,y) coordinate system. Each cross-section and profile-line in the FishNet is transformed individually. The coordinate transformation procedure described in section 4.5.2 is developed for points. Therefore, to transform a line from (s,n) coordinates to (x,y) coordinates, all its vertices are first converted into points. The points are transformed back to (x,y) coordinates, and are then joined to form a line. If a point P(s p ,n p ) has to be transformed to (x,y) coordinates, three items are required with reference to the centerline in (x,y) coordinates (Figure 4.41): the segment of the centerline closest 124 to point P, angle made by P with the closest segment, and distance of P from the starting point of the closest segment. To locate the segment closest to P, s P is evaluated with respect to the measure values at all the vertices of the centerline. For example, with reference to Figure 4.41b, m B and m C are the measures at vertices B and C, respectively. So, if m B ? s P < m C , then P is closest to BC, and the angle made by P with BC (?) can be computed as follows, ? ? ? ? ? ? ? ? ? = ? BP p ms n 1 tan? (4.29) 125 C ? ? B G P x Thalweg C B A n y Vertices P(sp,np) s a) b) Figure 4.41 Transformation of a point from (s,n) coordinate system to (x,y) coordinate system. (a) Point P as one of the vertices of a line in (s,n) coordinates; (b) Mapping of P in (x,y) coordinates with respect to segment BC of the centerline. The distance of P from the starting point of segment BC is computed by using by using equation 4.29, and point P can then be mapped into (x,y) coordinates by using BP and ? with reference to BC. )sin(? = p n BP (4.30) After all the vertices of a single line in (s,n) coordinates are transformed to (x,y) coordinates, they are joined to form a line. This process is repeated until all the lines in (s,n) coordinates are transformed to (x,y) coordinates. The relative orientation of lines in the FishNet is preserved when the data are transformed from (s,n) coordinates to (x,y) coordinates. Therefore, in the (x,y) coordinates, the lines perpendicular to the centerline are cross-sections while the lines parallel to the centerline are profile-lines (Figure 4.42). 2D Hydraulic FishNet 3D Hydraulic FishNet Figure 4.42 Hydraulic FishNet In addition, when the cross-sections and profile-lines intersect in (x,y) coordinates, they are mutually perpendicular. This is a useful result that can serve as input to hydraulic models in two different ways. For two or three-dimensional models, the result can be used as it is to define a finite element mesh. For one- 126 127 dimensional models the cross-sections can be separated, and then used as an input to define the one-dimensional channel geometry. 4.6 SUMMARY This chapter answers the first question (Given the river channel bathymetry as (x,y,z) points, how to create a standardized representation in the form of cross-sections and profile-lines?) posed in section 1.4. The goal of developing a procedure for storing the river channel bathymetry in a standardized vector dataset is achieved by creating a three-dimensional channel-oriented mesh comprising of cross-sections and profile-lines. The procedure uses the concept of (s,n) coordinate system to transform a regular ArcGIS FishNet into a channel- oriented network of cross-sections and profile-lines. The issue associated with calculations of (s,n) coordinates for the bathymetry points lying adjacent to the intersection of polyline segments is resolved by replacing the polyline segments with Bezier curves. Since the quality of FishNet is dictated by the quality of underlying raster grid, interpolation techniques available with the spatial analyst extension of ArcGIS are tested for their ability to predict a reasonable channel bed from discrete bathymetry points. It is found that the anisotropy in the river bathymetry data plays a major role in the interpolation process. Anisotropic kriging, which accounts for anisotropy in the data performed better than any other techniques examined in this research. A modified inverse distance scheme developed in this research, elliptical inverse distance scheme, which also takes into account the anisotropy in the bathymetry data performed second next to anisotropic kriging. Chapter 5 An Analytical Model for Extrapolation of River Channel Bathymetry 5.1 INTRODUCTION Chapter 4 describes a procedure to create a three-dimensional description of a river channel. The procedure uses a set of bathymetry points as an input, and the end result is a network of cross-sections and profile-lines. The FishNet description of channel bathymetry can be considered as an analytical model described by the function: 128 ),(),(?),( yxyxzyxz ?+= ),(? yxz ),(? yxz (5.1) Where z(x,y) is the measured bathymetry, is the mean surface for the channel bed (major trend in the bathymetry, deterministic component) and ?(x,y) are the departures from the mean (stochastic component). If the three-dimensional features of the FishNet (cross-sections and profile-lines) are described by using analytical forms then it may be possible to use the analytical form of the channel to estimate the channel bathymetry at locations where there are no bathymetric data. This chapter deals with developing an analytical model for the channel bathymetry. This analytical model hereafter will be referred to as River Channel Morphology Model (RCMM). Both and ?(x,y) in equation 5.1 are equally important to provide a meaningful description of channel bathymetry. However, due to the complex nature of the problem, only the first term (deterministic component) is addressed in this research. In addition, as mentioned in chapter 1, the three-dimensional structure of river channels is influenced by channel planform, flow, geology, land use, climate, etc. RCMM, however, takes only the channel planform and flow into consideration. A conceptual model for the development of RCMM is discussed in section 5.2. Section 5.3 describes the development of RCMM, following the application of RCMM in sections 5.4, 5.5, and 5.6. Section 5.7 briefly discusses the modeling of ?(x,y) term in equation 5.1 in the context of stochastic analysis. 5.2 CONCEPTUAL MODEL RCMM relates channel characteristics, namely the thalweg location and the cross-sectional shape, with the channel planform. Looking downstream, the thalweg location refers to the distance of the thalweg from the left bank of the channel. RCMM is based on a conceptual model that is explained with reference to Figure 5.1. C C C C B BB B A A A A a) b) Figure 5.1 Conceptual model for RCMM. (a) Channel planform and thalweg; (b) cross-sectional forms for different thalweg locations 129 130 As mentioned in section 2.7, in meandering river channels, the cross-sections at meandering bends undergo deposition at the inner banks and erosion at the outer banks. This process of erosion and deposition in conjunction with variations in boundary shear stresses and velocity fields at the meandering bends creates asymmetric cross-sectional shapes by moving the thalweg towards the outer banks. Therefore, the development of cross-sectional asymmetry in river channels can be related to the change in the channel planform (Leopold et. al., 1964; Knighton, 1998). RCMM assumes the cross-sectional asymmetries are related to channel planform, and is based on the following characteristics about meandering channels: ? The thalweg has a typical pattern that follows the channel planform (Figure 5.1a). For example, at the meandering bend, the thalweg is close to the outer bank, while the thalweg is approximately in the center when the channel is straight (not meandering). ? Depending on the thalweg location, the cross-sections have asymmetric or symmetric forms (Figure 5.1b). For example, when the thalweg is close to the bank, the cross-sections have asymmetric form. If the thalweg is approximately at the center of the channel, the cross-sections have a symmetric form. In summary, the direction of channel meander relates to the thalweg location, which in turn helps to define the cross-sectional asymmetry. Therefore, if the knowledge about the channel planform is available, it may be possible to predict the cross-sectional form. This conceptual model is simple, and there are 131 exceptions to it. For example, as shown later in the results, the thalweg does not always lie in the center when the channel is straight. However, the idea is to build a model based on simple concepts, and study the results so that the model can be modified in future to incorporate additional information. RCMM is developed and calibrated by using the data for Site 2 of the Brazos River, and is then verified by applying it to Site 1 of the Brazos River (Figure 3.1). The data at Site 1 covers a short reach that includes only one meandering bend whereas the data at Site 2 covers 50 km of Brazos River with several meandering bends. Therefore, the data at Site 2 is preferred to develop RCMM. As mentioned in chapter 1, RCMM is applicable only to meandering channels with an alluvial bed, and it does not take into account the effect of tributaries. The methodology to develop RCMM is described in the next section. 5.3 METHODOLOGY FOR DEVELOPING RCMM Based on the conceptual model discussed in previous section, RCMM relates the shape of the channel planform with the cross-sectional form. Therefore, by knowing only the shape of the channel planform, the three- dimensional structure of the channel can be described. The flowchart for the development of RCMM is shown in Figure 5.2. The detailed procedure to develop RCMM is discussed step-by-step in the following sections. (Input) Bathymetry Points and channel boundary (1) Transfrom the data from (x,y,z) to (s, n, z) (2) Normalize the data so that the width and depth of the channel are equal to one (3) Locate the thalweg (4) Calculate the radius of curvature for short river segments using the left bank (5) op relationship between the thalweg location and radius of curvature Devel (6) analytical form to describe normalized cross-sections nd relate their parameters to thalweg locations Use an a (7) evelop hydraulic geometry relationships to get width and depth of channel from the flow data D (8) Use the width and depth computed in (7) to rescale the normalized cross-sections (9) in the cross-sections with profile-lines to create a three-dimensional description of the river channel Jo Data are in original dimensions Data are in normalized dimensions Figure 5.2 Flowchart for the development of RCMM 132 133 5.3.1 Data Transformation and Normalization This section describes steps 1 and 2 of the flowchart shown in Figure 5.2. Application of RCMM involves using the data at one site for calibration and the data at other locations for verification. Therefore, the application of RCMM involves dealing with different locations and different scales. To make this process generic, the data are handled in a system that is independent of location and scale. After processing, the data can be transformed back to their original coordinates. The first transformation to the data involves converting the data from Cartesian coordinates (x,y,z) to curvilinear orthogonal coordinates (s,n,z) as described in section 4.5.3. Data transformation from (x,y,z) to (s,n,z) coordinate system makes the data independent of location because in the (s,n,z) coordinate system, all river channels are referenced with respect to the flow direction, irrespective of their location and planform in the Cartesian coordinate system (Figure 4.14). In this research, data normalization refers to making the data independent of scale. This is accomplished by converting the data to a non-dimensional form, in which the width and the depth of the channel are unity. For example, Figure 5.3 shows a typical cross-section with regular (s,n,z) coordinates. In Figure 5.3, looking downstream, n L is the n coordinate at the left bank with a negative sign, n R is the n coordinate at the right bank with a positive sign, and Z b is the bank elevation with respect to mean sea level. The s coordinate is a measure along the channel, and it is the same at any location for a given cross-section. After transforming the data to a normalized domain, the n coordinate for any point Comment: Applying what? The verb applying usually has an object with it. becomes n*, which is zero at the left bank and is equal to one at the right bank. Similarly, the z coordinate for any point becomes z*, which is zero at the banks and is equal to one at the thalweg. d = Zb - Zd 0 nL nR Zb (-) (+) P(ni, zi) Zd W = abs(nL) + nR Figure 5.3 A typical cross-section with (s,n,z) coordinates. P(n i ,z i ) is any point in a cross-section; n L and n R are the n coordinates of left and right bank, respectively; Z b and Z d are the bank elevation and thalweg elevation, respectively The data are transformed using equations 5.2 and 5.3 as described below. With reference to Figure 5.3, for any bathymetry point P(n i ,z i ), the non- dimensional coordinates are: n i * = (n i ? n L )/W (5.2) z i * = (Z b ? z i )/d (5.3) Where W = width of the channel d = maximum depth (thalweg) Similarly, if t is the thalweg location (distance from the left bank to the thalweg), then t* (thalweg location in normalized domain) is calculated as t* = |n L /W| (5.4) 134 135 The conversion of bathymetry data to a normalized form can be carried out in two different ways: global and local. Here the term global refers to the entire channel, and local refers to a short reach within the channel. In the case of global conversion, the widest section and the deepest point in the river channel are used for W and d in equations 5.2 and 5.3, respectively. The global conversion process is undesirable because the deepest point in the entire channel may be several times deeper than the rest of the bathymetry data. This makes the entire channel in non-dimensional form relatively flat compared to the deepest point. Therefore, local conversion is used. In the case of local conversion, the channel is divided into a number of short reaches, and the data are converted for each reach individually. Finding a reach length for local conversion depends upon the density of the bathymetry data. An optimum reach length that captures enough data points to define a cross-section adequately must be determined. The main criteria used to find an appropriate length for a reach is to have enough data points to be able to locate the thalweg. Reaches ranging from 25m to 300m were analyzed, and a 200m long reach (@ 200 points) was found satisfactory for Site 2. For Site 1, which has a higher density of bathymetry data, a 50m long reach (@ 200 points) was found satisfactory for normalization. Among all bathymetry points for each reach, the points with minimum n-coordinate (n L ) and maximum n-coordinate (n R ) are used to calculate the width (W = abs[n L ] + n R ). Similarly the points with minimum z-coordinate (Z d ) and maximum z-coordinate (Z b ) were used to calculate the depth (d = Z b -Z d ). Figure 5.4 shows the data for a reach in the original form and in the non-dimensional form. The conversion only changes the scale while preserving the original shape of the cross-section. 136 0 0.25 0.5 0.75 1 0 Normalized width, n* E l eva t i o n , z * 0.25 0.5 0.75 1 Bathymetry points 17 . 4 19 . 9 22.4 -20 5 30 55 80 Width, n (m) E l e v a ti o n , z (m ) Bathymetry points a) b) Figure 5.4 Data normalization for a cross-section. (a) Bathymetry data in original coordinates; (b) Bathymetry data in normalized coordinates The entire bathymetry data at Site 2 (50 km long) is transformed and normalized to make it independent of location and scale. The transformation does not necessarily mean mapping the data in the (s,n*,z*) coordinate system. It just means that all the bathymetry points are assigned (s,n*,z*), and these coordinates are stored as their attributes for use in further analysis. 5.3.2 Relationship Between Channel Planform and Thalweg Location This section describes steps 3, 4, and 5 of the flowchart shown in Figure 5.2. First, the thalweg for Site 2 is identified by using the methodology discussed in section 4.5.1. Looking downstream, the thalweg location refers to the distance of the thalweg along a cross-section from the left bank of the channel. The thalweg location can also be calculated using the right bank. In this research, however, the left bank is used because it has the least n coordinate in the (s,n,z) 137 coordinate system. In the normalized domain, where the width of the channel is unity, the thalweg location (t*) is always between zero and one. The thalweg location dictates the asymmetries in the cross-sectional forms (Figure 5.1), and these asymmetries have been related to flow by Knighton (1981; 1982; 1984). The main goal of Knighton?s work was to develop asymmetry indices to study the adjustment of channel geometry over time, and these indices were based on extensive field measurements. In this research, the radius of curvature of the left bank, which is easy to compute, is used to quantify the cross-section asymmetry. Radius of curvature can be computed by using either banks or the thalweg, but the radius of curvature of the thalweg is usually different than the radius of curvature of the channel. Between left and right bank, left bank is preferred for the same reason mentioned earlier (least n-coordinate). The radius of curvature, which is used as an indicator for the channel planform, is then related to the thalweg location (t*). With reference to the conceptual model in Figure 5.1, if the radius of curvature of a particular reach is small (meandering bend), the thalweg is close to the bank with an asymmetric cross-section. If the radius of curvature is large, the thalweg is approximately at the center of the channel with a nearly symmetric cross-section. A GIS procedure is developed to compute the radius of curvature, and relate it with the thalweg location. The boundary of the channel (polygon) is manually split into two banks, with only the left bank (looking downstream) used in computations. The left bank is divided into a number of segments such that the length of each segment is approximately equal to the meander wavelength, which is approximately 10-14 times the width of the channel (Knighton, 1998). The average width of the channel at Site 2 is about 60m, and the average length of the segments is about 650m. The bank is divided manually to make sure that each segment does actually represent a bend and the quality of the data at each segment is good (Figure 3.6). The locations of segments along Site 2 (26 segments along a 50 km long reach) that are used for computing the radius of curvature, and a typical segment used for computations are shown in Figure 5.5. Reach Location [ 063 Kilometers Channel Boundary 4 [ [ [ [ [ [ [ [ [ [ [ [ [ [ [ [ [ [ [ [ [ [ [ 012060 Meters Bathymetry Points Thalweg Boundary Left Bank Figure 5.5 Locations of segments used for computing the radius of curvature 138 For each segment, two items are computed: 1) the radius of curvature and 2) the distance from the midpoint of the segment to the thalweg. For each segment, the circle of curvature is computed using the two end points (i and k in Figure 5.6) and the mid point (point j). r2 r1 C1 C2 Left boundary p i j k C2 r2 Figure 5.6 Radius of curvature calculations to quantify the meandering shape of the channel The point of intersection of the perpendicular bisectors of sub-segments ij and jk defines the center of the circle of curvature (c 2 ), while the radius of this circle is the radius of curvature (r 2 ) for that particular segment. Looking downstream, if the center of curvature is to the right hand side of the segment (circle c 1 ), the radius of curvature is assigned a positive value, and the radius of curvature for circle c 2 is assigned a negative value. Therefore, a positive radius of curvature 139 means the channel is meandering to the left, and a negative radius of curvature means the channel is meandering to the right. Thus, the value of the radius of curvature not only quantifies the meandering channel planform, but also indicates whether the channel is meandering to the right or left. For each segment, the radius of curvature and the thalweg location (t*) are computed. Figure 5.7 shows the results that relate the radius of curvature and the thalweg location. 0.00 0.25 0.50 0.75 -15000 -10000 -5000 0 5000 10000 15000 Radius of Curvature (m) 1.00 Th a l w e g l o c a t i on, t * Model Observed Data R 2 = 0.8238 R 2 = 0.8717 Figure 5.7 Relationships between thalweg location and radius of curvature Figure 5.7 has following relationships: t* = -0.076*ln(r) + 1.21 0 < r < 12500m (5.5a) t* = 0.5 r ? 12500m (5.5b) t* = 0.087*ln(abs(r)) ? 0.32 -12500m < r < 0 (5.6a) t* = 0.5 r ? -12500m (5.6b) Where, t* and r are thalweg location and radius of curvature, respectively. 140 141 Equations 5.5and 5.6 relate a dimensional term (radius of curvature) to a non-dimensional thalweg location (t*). It can be assumed that wider channels have higher radius of curvature and compared to narrower channels. The relationship between the width and the radius of curvature, however, is debatable (Hey, 1976). In addition to equations 5.5 and 5.6, two more equations (5.7 and 5.8) are developed with non-dimensional radius of curvature (r* = r/W). t* = -0.096*ln(r*) + 0.922 0 < r* < 150 (5.7a) t* = 0.5 r* ? 150 (5.7b) t* = 0.092*ln(abs(r*)) ? 0.07 -150 < r* < 0 (5.8a) t* = 0.5 r* ? -150 (5.8b) Equations 5.7 and 5.8 will be used only to test the dependence of radius of curvature on channel width. Figure 5.8 shows how the radius of curvature changes between positive and negative values in the middle part of Site 2. 031.5 Kilometers Boundary (+) (+) (+) (+) (+) (+) (+) (-) (-) (-) (-) (-) (-) Location with positive radius of curvature Location with negative radius of curvature (-) (+): (-): Figure 5.8 Map with radius of curvature signs at meandering bends For the negative values of radius of curvatures (looking downstream and the channel meandering to the right), the thalweg location predicted by equation 5.5 (or 5.7) will be between 0.5 and 1.0. For the positive values of the radius of curvature (looking downstream and the channel meandering to the left), the thalweg location predicted by equation 5.6 (or 5.8) will be between zero and 0.5. Thus, equations 5.5 and 5.6 can be used to locate the thalweg by knowing the radius of curvature of the channel planform. For high values of r (which means a straight channel), t* is assigned a value of 0.5 to avoid the overlapping of domains defined by equations 5.5 and 5.6. Based on the calculations for the Brazos River, an upper bound of r = 12500m is used for equations 5.5 and 5.6. The upper bound on the radius of curvature may, however, vary for other river channels. 142 143 5.3.3 Relationship Between Thalweg Location and Cross-sectional Form To develop a relationship between the channel planform and the thalweg location, the radius of curvature is used as an indicator of the channel planform. Likewise, the cross-sectional shape has to be quantified with an analytical form. The concept of fitting an analytical form to cross-sections is not new. Attempts have been made to fit analytical forms to glacial valley cross-sections. For example, James (1996) compared power functions and second-order polynomials for glacial valley cross-sections in three Sierra Nevada valleys. In this research, several analytical forms including power functions, polynomials, splines and probability density functions (Gamma and Beta), were considered for fitting the channel cross-sections. For initial analysis, the analytical forms are evaluated mainly based on their ability to fit the bathymetry data (maximum R 2 ), number of parameters, and the ease with which they can be used or applied. After a suitable analytical form for describing the cross-sections is found, the significance of each parameter is tested by doing a non-linear regression analysis using SAS (Cody and Smith, 1997). A brief description of all the functions that are considered for fitting channel cross-sections is given below by using one sample cross-section as an example. Power Functions Power functions have the form f(x) = ax m . Power functions are simple and easy to compute with only two parameters: a and m. In the case of channel bathymetry, f(x) is z*, and x is the n* coordinate. Power functions work well with symmetric cross-sections such as glacial valley cross-sections used by James (1996), but are not very useful to represent the shape of river channel cross- sections. In addition, power functions can describe only one side of the channel cross-section from the thalweg. Therefore, the cross-section has to be split at the thalweg, and each side of the cross-section has to be fit with an individual power function. Figure 5.9 shows bathymetry data for a cross-section. 0.0 0.2 0.4 0.6 0.8 1. 0 0.0 0.2 0.4 0.6 0.8 1.0 n* z* Figure 5.9 Bathymetry data showing a cross-section The data are split at the thalweg and the power functions for each side are shown in Figures 5.10 and 5.11. The functions are fit using the least squares regression technique. z* = 0.9349(n*) 0.8625 R 2 = 0.9558 0 0.2 0.4 0.6 0.8 1 0.0 0.2 0.4 0.6 0.8 1.0 n* z* Power Function Bathymetry Points Figure 5.10 Power function fitted to the left hand side of the cross-section shown in Figure 5.9 144 z* = 0.1965(n*) -6.3469 R 2 = 0.6633 .0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1 0 n* z* Power Function Bathymetry Points Figure 5.11 Power function fitted to the right hand side of the cross-section shown in Figure 5.9 Even though the R 2 for the left part is greater than 0.9, the power function does not seem to describe the channel shape very well (Figure 5.10). In addition, there is no control over the parameters, and they can take any values between ?? +? to . Fitting two different functions, and then joining them so as to represent a continuous surface adds more complexity to an already complex problem. Polynomial Functions A polynomial function has the form f(x) = a m-1 x m-1 + a m-2 x m-2 +?+ a 1 x + a 0 , where x is the independent variable, m is a non-negative integer, and a m-1 , a m-2 ,?, a 0 are the coefficients of the polynomial function. The non-negative integer, m, is also called the degree of the polynomial function. Depending on the degree, polynomial functions have different names. A polynomial function of zero degree (m=0) has only one constant term (a 0 ). If a 0 =0, the zero degree polynomial is called a zero polynomial. If a 0 ? 0, the polynomial function is called a constant. A polynomial function of degree one has the form f(x) = ax + b, and is called a linear function. A polynomial function of 145 degree two has the form f(x) = ax 2 +bx+c, and is called a quadratic function. A third degree polynomial function has the form f(x) = ax 3 +bx 2 +cx+d. In general, when a polynomial function is plotted on a graph, the number of bends in the graph is equal to the degree of the polynomial minus one. Therefore, a linear function (zero number of bends) is not suitable to describe a cross-section. A quadratic function does produce one bend, but it does not provide enough flexibility to fit the cross-sectional shape (Figure 5.12). z* = -2.1407(n*) 2 + 2.6012(n*) - 0.1553 R 2 = 0.6163 0 2 4 6 8 0 -0.2 0. 0. 0. 0. 0. 1. 0.0 0.2 0.4 0.6 0.8 1.0 n* z* Quadratic Function Bathymetry Points Figure 5.12 Quadratic function fitted to the channel cross-section A third-degree polynomial fitted to the bathymetry data is shown in Figure 5.13. 146 z* = -6.93(n*) 3 + 8.64(n*) 2 - 1.83(n*) +0.21 R 2 = 0.8447 0.0 0.2 0.4 0.6 0.8 1.0 n* 0.0 0.2 0.4 0.6 0.8 1. 0 z* Third-degree Polynomial Bathymetry Points Extra Bend Figure 5.13 Third-degree polynomial fitted to the channel cross-section A third-degree polynomial does fit well to the bathymetry data, but it produces one extra bend (red circle in Figure 5.13). The problem of the extra bend can be tackled by joining two polynomials, but this will introduce at least three more terms in the new analytical form. In addition, the problem of continuity at the joints needs to be addressed. Instead, spline functions which fit piecewise polynomial functions maintaining the continuity at joints are considered. Smoothing Spline The discussion of spline functions in section 4.5.4 is limited to interpolating splines where the function passes through the observed data points. The spline function discussed in this section is an approximating or smoothing spline where the spline function does not necessarily pass through the observed data points (Figure 5.14). 147 0 1 2 3 4 5 6 7 8 0123456 9 0 1 2 3 4 5 6 7 8 9 0123456 Knots Control points Spline z x Interpolating Spline Approximating Spline Figure 5.14 Interpolating and smoothing splines Let (x 1 , z 1 ), (x 2 , z 2 ), ?, (x n , z n ) be the n number of pairs of observations. Let z(x) be a B-spline function that is fitted to the observed data. A B-spline uses basis functions (weights having values from zero to one) to describe the polynomial functions. It is desired to fit a smoothing spline through (x 1 ,z 1 ), (x 2 , z 2 ), ?, (x n , z n ). The problem is to find a spline z(x) for which ()dxxz n x x ? = 1 2 )(? )=? (5.9) is minimal, subject to the condition that ( Szzw n i i ii ?? ? =1 2 )?( (5.10) where ? is a smoothing norm (energy function), w i are the weights, and S is a specified non-negative number, also called as smoothing factor. The goal is to find the optimum number of knots and the parameters of z(x). The solution for equations 5.9 and 5.10 is a compromise between the following two approximation objectives: ? to obtain an approximating function that is as smooth as possible (? is small) 148 ? to obtain an approximating function that fits the data as close as possible (? is small) The smoothing factor, S, controls the extent to which these objectives are satisfied. The CURFIT algorithm available in the public domain (www.netlib.org) fits a smoothing spline using equations 5.9 and 5.10. A cubic B-spline with ten knots fitted to the cross-section is shown in Figure 5.15. 0.0 0. 0. 0. 0. 1. 0.0 0.2 0.4 0.6 0.8 1.0 n* z* 2 4 6 8 0 Smoothing Spline R 2 = 0.8905 Bathymetry Points Figure 5.15 Smoothing spline fitted to the cross-section The smoothing spline fits well to the bathymetry data (R 2 = 0.8905). However, the number of parameters (10) is relatively higher compared to power and polynomial functions. Regardless, spline is the best fit so far compared to power and polynomial functions. Beta Distribution A random variable x is said to have a standard beta distribution with parameters ? and ? if the probability density function (pdf) of X is given by ,0,0,10,)1( ),( 1 ),|( 11 >>< ?, and is symmetric when ? = ? (Figure 5.16). Figure 5.16 Beta distribution 150 The beta distribution, however, has two shortcomings. First, as shown in Figure 5.16, it is relatively flat at one of the tails. A flat tail is undesirable because it indicates zero cross-sectional area towards the end of the river channel cross- sections. Second, in the normalized domain, z* < 1 for channel cross-sections. This condition (z* < 1) is violated with a single beta distribution. In other words, it is difficult to fit a pdf to a channel cross-section while preserving its unit area property (integral of a pdf is equal to one) given by: ? ? ?? =1)( dxxf (5.14) Therefore, the pdf has to be multiplied by a factor (0><< ? ? ? ? ? ? ?= ?? kxkxx B xf ?? ?? ?? ?? (5.15) Equation 5.15 is used to fit a beta distribution to the cross-section as shown in Figure 5.17. 0.0 0.2 0. 0. 0. 1. 0.0 0.2 0.4 0.6 0.8 1.0 n* z* R 2 = 0.8313 4 6 8 0 Figure 5.17 Beta distribution fitted to the cross-section 151 The parameters of the beta distribution are determined by using the Newton- Rhapson optimization technique available with the SOLVER in MS Excel. The objective function is the squared difference between the observed values and the model predictions. The goal is to estimate the parameters of beta distribution that minimize the objective function. Although a beta distribution multiplied by k fits well to the bathymetry data, the problem of flat tail (red circle in Figure 5.17) still remains. In addition, the beta fit (R 2 = 0.8313) is not as good as the spline function fit (R 2 = 0.8905). The shortcomings of a single beta distribution are overcome by combining two beta distributions (Figure 5.18). Since the main reason behind using two beta distributions is to deal with the flat tail of a single beta distribution, a symmetric beta distribution (? = ?) is used to add only one extra parameter to equation 5.15. A symmetric beta (? = ?) is added to an asymmetric beta (? ? ?), and multiplied by k to maintain z*< 1 and overcome the flat tail problem. Figure 5.18 Combination of two beta distributions 152 The combination of two beta distributions is designated as a f(z), and is calculated as: f(z) = {f(x|? 1 ,? 1 ) + f(x|? 2 ,? 2 )}*k, where ? 1 ? ? 1 , ? 2 = ? 2 , 0>?= ?? baxexbaxf bxa )(? ab a (5.17) where ?(a) is a standard gamma function given by equation 5.13. Figure 5.20 shows a series of gamma probability density functions for several (a,b) pairs. The parameter a controls the peakedness of the distribution, and the parameter b controls the spread of the distribution. Accordingly, parameters a and b are called shape and scale parameters, respectively. 0.0 02468 x f(x | a , b ) 0.2 0.4 0.6 0.8 1.0 a = 1, b = 1 a = 2, b = 0.5 a = 2, b = 1 a = 2, b = 2 Figure 5.20 Gamma density functions for different (a,b) pairs To overcome the problem of unit area property of a pdf (equation 5.14), the gamma distribution is also multiplied by k (0>? ? ? ? ? ? ? ? ? ? = ?? kbaxkex ab baxf bxa a (5.18) The gamma distribution fitted to the cross-section using equation 5.18 is shown in Figure 5.21. 154 0. 0. 0. 0. 0. 1. n* z* 0 2 4 6 8 0 0.0 0.2 0.4 0.6 0.8 1.0 R 2 = 0.6729 Figure 5.21 Gamma distribution fitted to the cross-section Similar to beta distribution, the parameters of gamma distribution are also determined by using the Newton-Rhapson optimization technique available with the SOLVER in MS Excel. Like a single beta distribution, a single gamma distribution also has a flat tail problem (red circle in Figure 5.21). In addition, the fit (R 2 = 0.6729) is not as good as a single beta distribution (R 2 = 0.8313). Therefore, a combination of two gamma distributions is used for fitting the channel cross-section. In addition, a = b for one of the distributions to minimize the number of parameters. The combination of two gamma distributions is designated as a f(z), and is calculated as: f(z) = {f(x|a ? ?b ? ) + f(x|a 2 ?b 2 )}*k, where a ?? ? ? b ?? , a 2? = ? b 2 , and k < 1. (5.19a) Which for fitting a cross-section looks like {}kbanfbanfz *,|*(),|*(*? 2211 += (5.19b) Like smoothing spline and beta distributions, a combination of gamma distributions also fits well to the bathymetry data (Figure 5.22). The fit, however, 155 is better in the case of smoothing spline and beta distributions based on the R 2 value. 0. 0. 0. 0. 0. 1. 0.0 0.2 0.4 0.6 0.8 1.0 n* z* 0 2 4 6 8 0 Gamma Distribution R 2 = 0.8717 Bathymetry Points Figure 5.22 Combination of gamma distributions fitted to the channel cross- section Summary of Analytical Forms Table 5.1 gives a summary of different analytical forms used in this research to fit a normalized channel cross-section. Analytical Form Number of Parameters R 2 Power Function* 2 0.8096 Quadratic Polynomial 3 0.6163 Third Degree Polynomial 4 0.8447 Smoothing Spline 10 0.8905 Single Beta Distribution 3 0.8313 Two Beta Distributions 4 0.8974 Single Gamma Distribution 3 0.6729 Two Gamma Distributions 4 0.8717 *: R 2 averaged for both left and right part. Table 5.1: Summary of different analytical forms for fitting channel cross-section 156 157 From Table 5.1, it is clear that among different analytical forms considered, a combination of two beta distributions outperforms all other functions in terms of better fit (highest R 2 ). In addition, a combination of two beta distributions has only four parameters which are easy to estimate. Finally, the overall fit and the significance of each parameter for a combination of two beta distributions is verified by performing a non-linear regression analysis, and the results are shown in Table 5.2. Overall R 2 = 0.89 P value = < 0.0001 Mean Square Error (MSE) = 0.043 Estimate Std. Error t ratio ? 1 6.63 0.219 30.27 ? 1 2.65 0.071 37.32 ? 1 / ? 1 1.63 0.073 22.33 K 0.22 0.005 44 Table 5.2: Summary of regression analysis for fitting a combination of beta distributions to the cross-section data From Table 5.2, it is clear that a combination of two beta distributions have statistically significant predictive capability (P<0.0001), and all the variables are statistically significant (low standard error). Therefore, a combination of two beta distributions is used for fitting the channel cross-sections so that the cross- sectional form can be related to the thalweg location. Relating Channel Cross-sections to Thalweg locations This section describes step 6 of the flowchart shown in Figure 5.2. Once the channel cross-sections are defined using a combination of two beta distributions, the next step is to relate the parameters of beta distributions to 158 different thalweg locations. To relate beta distributions (channel cross-sections) to different thalweg locations, the bathymetry points are selected manually for nine different thalweg locations and beta parameters (? 1 , ? 1 , ? 2 , ? 2 and k) estimated by performing non-linear regression analysis. The nine different thalweg locations are 0.1 to 0.9 in increments of 0.1 units, and they are selected from the same 26 numbers of reaches that are used for the radius of curvature computations (Figure 5.5). So, for example, if the thalweg is located at about 0.3 units from the left bank, then the bathymetry points covering one reach length (200m) are selected manually for one computation. This gives a beta cross-section for a thalweg location that is 0.3 units from the left bank. This process is repeated for all nine thalweg locations. In some of the reaches, the effect of macroscopic features such as point bars, riffle, etc. is more pronounced compared to others, and such reaches are not used for estimating beta parameters. For the same reason, if two or more reaches have the same thalweg location, then the reach that gives a best fit (least sum of squares of residuals) is used for estimating the beta parameters. However, for each thalweg location, the parameters are slightly modified to generalize the shape of the cross-section at those locations. The generalized parameters are useful for computing the beta distributions at other river channels. The beta parameters estimated for different thalweg locations are shown in Figure 5.23 and Table 5.3. 0 2 4 6 8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t* 10 ? 1 , ? 1 , ? 2 , ? 2 , k ?1 ?1 ?2, ?2 k Figure 5.23 Beta parameters for different thalweg locations Figure 5.23 can be used to describe a cross-section using a combination of two beta distributions for any thalweg location. Although the value of k does not change considerably for different thalweg locations, k is still kept as a variable rather than assuming it as a constant. This may prove useful while applying the model to other river channels. 159 160 Overall Model Results Individual Parameter Results t* R 2 MSE p-value Parameter Estimate Std. Error t ratio ? 1 1.33 0.038 35 ? 1 6.61 0.345 19 ?2/?2 1.65 0.091 18 0.1 0.89 0.008 < 0.0001 k 0.19 0.050 4 ?1 2.49 0.156 16 ? 1 6.16 0.481 13 ?2/?2 1.31 0.062 21 0.2 0.78 0.019 < 0.0001 k 0.23 0.005 46 ? 1 1.73 0.073 24 ?1 3.04 0.140 22 ? 2 /? 2 1.47 0.048 31 0.3 0.88 0.008 < 0.0001 k 0.28 0.003 93 ? 1 3.34 0.276 12 ?1 3.46 0.186 19 ?2/?2 1.15 0.044 26 0.4 0.75 0.028 < 0.0001 k 0.29 0.005 58 ?1 4.83 0.350 14 ?1 4.06 0.261 16 ?2/?2 1.08 0.037 29 0.5 0.84 0.012 < 0.0001 k 0.26 0.004 65 ?1 3.30 0.195 17 ?1 1.94 0.085 23 ?2/?2 1.53 0.049 31 0.6 0.88 0.008 < 0.0001 k 0.31 0.005 62 ?1 6.39 0.509 13 ?1 3.00 0.183 16 ?2/?2 1.03 0.035 29 0.7 0.77 0.013 < 0.0001 k 0.25 0.005 50 ?1 5.01 1.018 5 ?1 1.39 0.212 7 ?2/?2 2.13 0.451 5 0.8 0.87 0.008 < 0.0001 k 0.22 0.007 31 ?1 7.91 0.358 22 ?1 1.99 0.043 46 ?2/?2 4.02 0.211 19 0.9 0.79 0.016 < 0.0001 k 0.21 0.050 4 Table 5.3: Beta parameters for different thalweg locations 161 5.3.4 Rescaling the normalized cross-sections This section describes steps 7 and 8 of the flowchart shown in Figure 5.2. The cross-sections defined using a combination of beta distributions are in (s,n*,z*) coordinate system with unit width and unit depth. To use these cross- sections for describing real river channels, first, they have to be rescaled to fit field dimensions, and second, transferred back to the Cartesian (x,y,z) coordinate system. River channel cross-sectional form adjusts over time through the process of erosion and deposition to accommodate the varying flow. Since the discharge increases downstream, the width and depth of the channel should similarly vary (Knighton, 1998). This is the underlying principle of hydraulic geometry relationships. It has been a common practice to describe the spatial variability of channel-width and mean channel depth using the flow data (Moody and Troutman, 2002; Harman et. al., 1999; Miller et.al., 1996). The same approach is used here to rescale the normalized data by developing hydraulic geometry relationships. As discussed in section 2.8, hydraulic geometry relates the independent variable (flow) to dependent variables (average channel-width, average channel depth, and average velocity) through simple power forms as shown below: W = aQ b (5.20) d = cQ f (5.21) v = kQ m (5.2) where W, d, v, and Q are respectively average width, average depth, average velocity and discharge. Although the data are initially normalized by using the 162 thalweg depth, as shown later in the results, the average depth predicted by the downstream hydraulic geometry relationships is found to be adequate for rescaling the data. The United States Geological Survey (USGS) is the primary source of data on all the rivers in the United States. The data include time-series of stream levels, steamflow, and cross-section measurements for more than 850,000 gaging stations. To develop the hydraulic geometry relationships for the Brazos River, flow data and cross-section measurements are downloaded from the USGS website (http://nwis.waterdata.usgs.gov/usa/nwis/measurements). The USGS performs cross-section measurements, which include channel-width, mean cross- sectional area and mean velocity, to update the rating curves for the gaging stations. The criteria used in the selection of gaging stations takes into account the data quality requirement for flow data studies, and this introduces some bias in the cross-sections surveyed at the gaging stations. For example, the cross-sections at gaging sites are less susceptible to erosion and deposition. However, hydraulic geometry relationships have been developed using gaging station data (Harman, et. al., 1999; Dodov and Foufoula, 2003). In addition, as will be shown later, these data are found adequate for this research. A typical relationship between average depth (d), average width (W), average velocity (v) and the flow (Q) for a gaging station at Richmond (# 08114000) is shown in Figure 5.24. w = 95.654Q 0.1206 R 2 = 0.8164 100 100 1000 10000 100000 Flow, Q(cfs) a = 95.654 b = 0.1206 1000 A ver ag e Wi d t h , w ( f eet ) a) d = 1.4895Q 0.2537 2 1 10 100 A v er ag e D e p t h , d ( f eet ) c = 1.4895 f = 0.2537 R = 0.8672 100 1000 10000 100000 Flow, Q (cfs) b) y = 0.0094x 0.5894 R 2 = 0.9576 0.1 1 10 Ve lo c i t y , v (fp s ) 100 1000 10000 100000 Flow, Q (cfs) k = 0.0094 m = 0.5894 c) Figure 5.24 Hydraulic geometry relationship between (a) average width; (b) average depth; (c) average velocity and flow of the Brazos River at Richmond 163 Hydraulic geometry relationships are developed for ten gaging stations (Figure 5.25) using the measurement data from the USGS, with the results summarized in Table 5.4. 0 200100 Kilometers 05025 Kilometers 08116650 08114000 Site 1 Site 2 Gaging Stations Brazos River Brazos Watershed 08116650 08114000 08111500 08108700 08098290 08096500 08093100 08091000 08090800 08089000 Figure 5.25 Gaging stations locations along the Brazos River Drainage areas (A) in Table 5.4 are obtained from the USGS. Q 2 (two-year return period flow) for each station is determined from the flow duration curves, also obtained using the USGS flow data. The use of Q 2 is explained later. 164 165 Gauging station number A (km 2 ) Q2 (m 3 /s) a b c f k m a*c*k b+f+m 08089000 36894.38 404.93 17.822 0.327 0.169 0.376 0.403 0.274 1.21 0.98 08090800 40587.70 539.44 24.632 0.321 0.252 0.350 0.204 0.291 1.27 0.96 08091000 42092.48 625.80 21.369 0.361 0.169 0.375 0.364 0.236 1.32 0.97 08093100 45785.81 549.35 54.609 0.196 0.140 0.451 0.130 0.353 1.00 1.00 08096500 51781.63 668.28 20.746 0.333 0.134 0.502 0.257 0.234 0.71 1.07 08098290 54053.05 835.35 189.41 0.075 0.036 0.567 0.136 0.364 0.94 1.01 08108700 76360.62 1172.32 97.624 0.114 0.488 0.354 0.025 0.513 1.17 0.98 08111500 88872.85 1472.48 120.33 0.103 0.373 0.394 0.023 0.492 1.04 0.99 08114000 92050.76 1580.08 95.654 0.121 1.490 0.254 0.009 0.589 1.34 0.96 08116650 92651.64 1461.15 60.153 0.157 0.065 0.582 0.241 0.272 0.95 1.01 Table 5.4. Hydraulic geometry parameters for the Brazos River. A is the drainage area; Q 2 is 2-year return period flow; a, c, k, and b, f, m, are the coefficients and power terms of equations (5.20), (5.21), and (5.22), respectively The parameters in Table 5.4 can be used to predict W, d, and v for any given flow at the gaging stations listed in the table. For any point along the river that does not lie exactly at any of the gaging station locations, the values obtained by using the parameters presented in Table 5.4 can be linearly interpolated with respect to the flow or upstream watershed area to predict W, d and v. The results (channel-width and depth) can then be used to rescale the normalized cross- sections. Since the rescaling depends on the flow, different flow conditions will produce different cross-sections. Bankfull discharge, which is defined as the discharge at which the channel is completely full, is commonly used in geomorphology to compute hydraulic geometry relationships. Bankfull discharge, however, does not have a consistent return period, and is reported to have return periods ranging from 1 to 2 years (Williams, 1978; Leopold et. al., 1964). Also, a bankfull channel cannot always be defined in the field (Knighton, 1998). To develop a consistent procedure, a 2-year return period flow is used as a standard to develop hydraulic geometry relationships and rescale the normalized cross- sections. In addition, the 2-year return period flow (Q 2 ) is related to upstream watershed area (A) as shown in Figure 5.26 to get the following relationship: Q 2 = 0.0004A 1.3284 (5.23) R 2 = 0.9703 10 0 10 0 0 30000 50000 70000 90000 110000 Drainage area, A (km 2 ) 10 0 0 0 D i sch ar g e , Q 2 (m 3 /s ) Figure 5.26 Relationship between 2-year return period flow (Q 2 ) and drainage area (A) for the Brazos River. In equation 5.23 Q 2 is in m 3 /s and A is in km 2 . Equation 5.23 can be used to express equations 5.20, 5.21, and 5.22 in terms of upstream watershed area, and the hydraulic geometry parameters (d, W, and v) can be predicted at any location along the Brazos River using the upstream watershed area. However, it should be noted that hydraulic geometry relationships work only for natural channels. For unnatural river channels (constructed channels), hydraulic geometry relationships do not apply. Under such circumstances, the width many not increase with the increase in drainage area downstream. Therefore, additional information such as 166 167 digital orthophoto quadrangles (DOQ) should be used to calculate the width. In addition, hydraulic geometry relationship gives just the depth of water; the elevation (z) should be calculated by subtracting the water depth from digital elevation models (DEM) or other survey data, if available. Thus, the hydraulic geometry relationships in conjunction with additional information such as DOQ and DEM can be used to rescale the normalized cross-sections. In the case of Site 2 of the Brazos River (Figure 5.25), there are no gaging stations located along the entire reach. The nearest gaging station (# 08114000) is located 40 km upstream. The DOQs of the area surrounding Site 2 are used to find the width of the channel, and the depth obtained from the hydraulic geometry relationship at gaging station # 08114000 is slightly increased to account for the increased upstream watershed area. After rescaling, the data can be transformed back to their original (x,y,z) coordinates by using the theory described in section 4.5.6. 5.3.5 Creating profile-lines using cross-sections This section describes step 9 of the flowchart shown in Figure 5.2. To create a complete three-dimensional description of river channels, the cross- sections must be joined by profile-lines. Profile-lines can be generated using two different approaches: 1) depth-based approach and 2) area-based approach. The depth-based approach divides the cross-section into regions with equal depths. If a cross-section that is 10m deep (from banks to the thalweg) has to be divided into three regions, the depth-based approach will introduce the profile lines at d/3 on each side (Figure 5.27). On the other hand, the area-based approach will divide the cross-section into regions with equal areas. After the profile-lines are generated, they are converted into Bezier curves for creating a smooth representation of the three-dimensional channel form. a) Depth-based profile-lines b) Area-based profile-lines d d 1 =d/3 d2=d/3 d3=d/3 A A3=A/3 A2=A/3 A1=A/3 Figure 5.27 Generation of profile-lines using cross-sections 5.4 APPLICATION OF THE RCMM TO SITE 1 OF THE BRAZOS RIVER The application of the RCMM involves starting with a blue line (centerline) obtained from the National Hydrography Dataset, and then converting this blue line into a three-dimensional description of a river channel through the following steps: 1. Using equations 5.20 and 5.21, predict the average width (W) and depth (d) for a 2-year return period flow. However, to verify RCMM using the data at Site 1 which was measured during a flow of 9700 cfs, the values of d and W are predicted for 9700 cfs instead of a 2-year flow. The 2-year return period flow is a standard flow that can be used at locations where there are no data. 2. Offset the blue line obtained from the medium resolution (1:100,000) NHD by a distance of W/2 on each side to establish the channel boundary. The medium resolution NHD does not always follow the actual river channel. 168 169 Therefore, the channel boundary is modified to accommodate additional information from aerial photographs. 3. Looking downstream, use the left bank to calculate the radius of curvature, and then locate the thalweg using the radius of curvature-thalweg relationship (equations 5.5 and 5.6). Since equations 5.5 and 5.6 are developed using the left bank as a reference, left bank is used for locating the thalweg. 4. From the thalweg locations, generate beta cross-sections. The parameters in Figure 5.23 can be used to generate beta distributions for any thalweg location. 5. Rescale the cross-sections described by beta distributions using W and d obtained from hydraulic geometry relationships (Step 1). Get the bank elevation from the digital elevation model (DEM) and subtract the rescaled bathymetry to get actual elevations. 6. Create profile lines from the cross-sections using the procedure described in section 5.3.5. Using the procedure listed in step 1, the values for d and W corresponding to a flow of 9700 cfs are predicted to be 4.25m and 94m, respectively. To locate the thalweg, the left bank of the channel is first divided into a number of segments (approximately 10*W). For each segment, the radius of curvature is computed, and then the thalweg location is identified by using equations 5.5 and 5.6. The thalweg identified by RCMM for Site 1 on the Brazos River is shown in Figure 5.28. 060300 Meters Legend Modeled Thalweg Measured Thalweg BoundaryPolygon 4 B B A A Figure 5.28 Thalweg prediction using the radius of curvature-thalweg location relationship (equations 5.5 and 5.6). Circled areas show locations where the thalweg prediction is imprecise. Sections A-A and B-B are plotted in Figures 5.30 and 31 170 0 600300 Meters Legend Measured Thalweg Modeled Thalweg BoundaryPolygon 4 Figure 5.29 Thalweg prediction using the non-dimensional radius of curvature- thalweg location relationship (equations 5.7 and 5.8) From Figures 5.28 and 5.29, it can be said that thalweg predicted by both the dimensional and non-dimensional form of relationships is approximately the same. This is not surprising given the relationships were developed using the data at the same river (Site 2 of Brazos). In general, RCMM creates a good description of thalweg location for most of the reach. However, there are locations (circled on Figure 5.28), where the model result is not in agreement with the observed data. The main reason for such deviation is the model?s tendency to locate the thalweg close to the center of the channel when the channel is straight. The model needs to be modified to incorporate additional information when the reach is straight. 171 After the thalweg is located, the beta parameters from Figure 5.23 are used to describe the cross-sections using a combination of two beta distributions (equation 5.16). The resulting cross-sections are then rescaled to fit the observed data. Figure 5.30 shows a cross-section described by RCMM at location A-A shown in Figure 5.28. 16 18 20 22 -20 0 20 40 60 80 n-coordinate (m) 24 z - c o or di n a t e ( m ) Observed Model Figure 5.30 Cross-section described by RCMM at section A-A shown in Figure 5.28 The cross-section described in Figure 5.30 fits well to the observed data at that particular location. However, this is not true at all the locations. For example, at cross-section B-B shown in Figure 5.28, the cross-section described by RCMM does not match completely with the observed data (Figure 5.31). 172 16 18 20 22 24 - 55 -35 -15 5 25 Normalized Width E l e v a ti o n (m ) Observed Model Figure 5.31 Cross-section described by RCMM at section B-B shown in Figure 5.28 The disagreement between the model and the data at cross-section B-B is a result of an abrupt change in the bathymetry just upstream. Cross-section B-B is located just downstream of a big dip in the channel bed, and the bathymetry at this location is not as smooth compared to the rest of the channel. The thalweg is in an abrupt transition from the right bank (looking downstream) towards the center of the channel. RCMM does not take into account such abrupt changes in bathymetry. RCMM only provides a mean surface for the channel bed, and it is obvious to expect some deviations of the observed data from the mean surface. Figure 5.31 provides one example where RCMM is unable to describe the cross- section precisely. Finally, after the cross-sections are generated, profile lines are generated to get a three-dimensional mesh as shown in Figure 5.32. The profile- lines are generated using the depth-based approach to divide the cross-sections into ten regions. 173 Figure 5.32 Three-dimensional description of a river channel in the form of cross- sections and profile lines 5.5 VERIFICATION OF RCMM RCMM is developed by using the data at Site 2 on the Brazos River, and is then applied to Site 1 on the same river (Figures 3.2 and 3.3). The performance of the model is studied by comparing the prediction errors on Site 2 and Site 1 along the Brazos River. The prediction errors are the deviations obtained by subtracting the observed data from the model results. The prediction errors are analyzed by using two measures: root mean square error (RMSE) and coefficient of determination (R 2 ). The network of cross-sections and profile-lines obtained from RCMM (Figure 5.32) can be interpolated to create a TIN (triangulated irregular network) surface. The prediction errors are then calculated by subtracting the observed bathymetry data from the TIN surface. If z 1 , z 2 , ?, z n are the observed values, and , , ?, z are the predicted values from the RCMM model, the RMSE of prediction is given by equation 4.28. The RMSE gives an average error of prediction, and is scale dependent. For example, for a very shallow river channel which is less than one meter deep, a RMSE of 0.5 is very high. On the other hand, for a deep channel (d >10 m), a 1 ?z 2 ?z n ? 174 RMSE of 0.5 is relatively low. Therefore, to make the RMSE scale independent, the standardized RMSE (RMSE*) is used (Colla et. al., 1999). The RMSE* is computed by dividing the RMSE by the standard deviation observed values as shown below: )( * i z RMSE RMSE ? = (5.24) is the standard deviation of observed values. where )( i z? It is desirable to have least RMSE* (close to zero), but since the RCMM model provides a mean surface for the channel bathymetry, it is unrealistic to expect a very low (less than 0.5m) RMSE. For comparison purposes, the RMSE* of the calibration site (Brazos River at Site 2) is taken as a benchmark, and the results from other site are compared with this value. The coefficient of determination (R 2 ) gives the proportion of the variations in the observed data that is described by the model prediction (Devore, 1995). Therefore, the higher the value of R 2 (close to one), the more successful is the model in explaining the variations associated with the observed data. The coefficient of determination is computed as: SST SSE R ?=1 2 () ? = ?= n i i zzi 1 2 ? (5.25) where SSE (5.26) is the sum of square of prediction errors, and 175 2 1 2 1 ? ? ? ? = ? = nn i i zSST 1 ? ? ? ? ? ? =i i z n (5.27) is the sum of squared deviations about the sample mean of the observed values. The RMSE* values at Site 1 and Site 2 along the Brazos River are 0.87 and 0.56, respectively. RCMM is developed using the data at Site 2, and is then verified by applying it to the data at Site 1. Generally, the calibration results are better than the verification results. Therefore, a higher RMSE* at Site 2 compared to Site1 is not surprising. Figure 5.33 shows the comparison of observed data and the model output at Site 2. The corresponding R 2 value is 0.68. Considering that RCMM models only the mean surface for the channel bed, an R 2 value of 0.68 seems reasonable. -15 -12 -9 -6 -3 0 M o de l O u t p ut ( m ) Bathymetry Points -15 -12 -9 -6 -3 0 Observed Elevation (m) R 2 = 0.68 Line of perfect fit Figure 5.33 Observed data and model output at Site 2 Figure 5.34 shows the comparison of observed data and the model output at Site 1 with a corresponding R 2 value of 0.25. 176 6 8 10 12 14 M o d e l O u tp u t (m ) 6 8 10 12 14 Observed Elevation (m) R 2 = 0.25 Bathymetry Points Line of perfect match Figure 5.34 Observed data and model output at Site 1 An R 2 value of 0.25 at Site 1 suggests that RCMM is not doing well in describing the channel bathymetry at Site 1. This issue is further investigated by looking at the prediction errors at Site 1 in detail. Figure 5.35 shows the histogram of prediction errors at Site 1. The histogram is skewed to the left with a skewness of ?1.22. In addition, most (more than 80 percent) of the prediction errors lie in the range of (-1.5, 1.5). The prediction errors can be split into two groups. The first group (Group 1) which accounts for 84 percent of the data includes those points that have prediction errors within a range of (-1.5,1.5), and the second group (Group 2) contains the rest of the data. The R 2 value for group 1 is 0.59 (Figure 5.36), which is reasonable considering an R 2 value of 0.68 for Site 2. Group 2 has a negative effect on the R 2 value of the entire dataset. To find out the reasons for the negative behavior of the data in Group 2, the spatial distribution of the errors is studied. 177 0 200 400 600 800 1000 -1 0. 3 -9. 6 -8. 9 -8. 3 -7. 6 -6. 9 -6. 2 -5. 5 -4. 8 -4. 2 -3. 5 -2. 8 -2. 1 -1. 4 -0. 8 -0. 1 0. 6 1. 3 2. 0 2. 7 3. 3 Prediction Error (m) 1200 N u m b e r of da t a point s Group 2 Group 1 Group 2 Figure 5.35 Histogram of prediction errors at Site 1. Group 1 covers 84 percent of the data with low prediction errors, and Group 2 cover data with high prediction errors 6 8 10 12 14 M o d e l O u tp u t(m ) R 2 = 0.59 68101214 Observed Elevation(m) Bathymetry Points Figure 5.36 Comparison of obverted elevation and model output for group 1 at Site 1 Figure 5.37a shows the spatial distribution of prediction errors. 178 179 00.525 Kilometers etry Points al 0. Bathym Residu -8.0 - -1.5 ! -1.5 - 1.5 1.5 - 4 !!!!!!!! !!!!!!!!! !!!!!!!!! !!!!!!!! !! !!!!!! !!! !!!!!!!! !!! !!!!! !!! ! !!! !!!! !!!!!!!!!!!!! !!!! !!! !! ! !! !! !!!!! ! !!!!!!!!!! !!!!!! !!!! !!!!!!!!!!!!!!!! !!!!!!! !!!!! ! ! !!!! !!!!! !!!!! !!! !!! !!!! !!!! !!! ! ! !!!!!!!!!! !!!!! !!!!! !!!!!!! !!!!! !!!! !!!!! !!! !!!!!!!!!!!!!!!! !!!!!!!! !!!! !!!! !!!! !!!!! !!!! !!!!!!! !!!!!!! !!!!! !!! !! !! !!!! !!! !!!! !!!!! !!!!! !!!!! !!! !!! !! !!!! ! !!!!! !!!!!!!! !!!!!! !!!!!!!!!!!! !!! ! !!! !! !!!!! !!!!! !!!!!! !!!! !!!!!! !!!!!!! !!!! !!!!!! !!!!! !! ! !!!! !!!!! !!!! !!!! !!! !!! !!! !!! !!!! !! !! ! !!!!!! !!! !! !!!!! !!!!!!! !! ! !!!! !!!!!!! !!!!!!! !!!!! ! !!! !!! !!!! !!!!!!! !! !!!! !!!!!!! !!! !! ! !!!! !!! !!! !!! !!! ! !!!!!!!!! !!!! ! !!!!!! !!! !!!! !!!!! !!! !!!! !!! !!!! !! !!!! !!!!!! !!!!!!!! !!!! !!! !!!!! !!! !!! !!! !!! !! !!! !!! !! !!!!!! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!! !!!! !!!! !! ! ! !! !!!! !!!! !!! !!!!! !! !!!!! !!!!! !!!!!! ! !!! !!! ! !!!! !!! !!! !!!! !! !!!!!!!!! !!!!!!! !!!! !!!! !!! !! !!! !! !! !!! !!! !!! !!! !!! !!! !!! !!!! !!!! !!!!! !!!! !!!! !!! !! !!! !! !!!! !!! !!! !!! !!! !!! !!!!!!! !!!! !!!! !!!! !!! !! ! !!! ! !!! !! !!!! !!!!! !! !!!! ! !!! !!! !!!! !! !!!! !!!! !!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!! !!!! ! !!!! !!!! !!!! !!!! !!!! !!! !!! !!! !!! !! ! !!!!!!! !!!!! !!!!!! !!!!!!!!!! !!!!! !!!! ! !!!!! !! !!!!!!! !!!!!!!! !!! !!! !!!!!!!!! !!!!!!!!!!!!!!!! !!!! !!!!!! !!!! !! !!!!!!!! !! !!! ! !!! !!! !!!!! !!!!!! !!! !!! ! !!!!!!!!!!!!!!! !!!!!! !!!!!! !! !!!! !!! !!! !!!! !!! !!! !!! ! !!!!! !!!!! !!!!! !!!!! !!!!! !!! !!! !! !!!!!!!!!!! !!!! !!! !!!! !!!!!!!! !!! !! !!! !!! !!! !!!!!!!!! !!! ! !! ! ! !!! !!!!!!!!!!!!!!!!!!!!! !!!! !!!!!!! !!!! !!!! !!!!!! ! ! !!! !!! !!! !! !!!! ! !!! !! !!!! !! !!! !!! !!! !! !!! !!!!!!!!!!!! !!!!! !!!! !! !!!! !!! !!! !!! !!!!! !! !!!!! !!!! !!!! !!!!! !!!!! !! !!!! !!!!! !!! !!!! !!!!! !!! !!!! !! !!! !! !!!!! !!!! !!! !!!!! !!!!!!!!!!! !!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!! !!!! ! !!!!!!!!!!!!! !!!! !!!! !!!!!!!! !!!! !!!! !!!!!! !!!! !!!! !!!! !!!! !!! !!! !!! !!!! !!!!! ! !!! !! ! !!!! !!! !!! !!! ! !!! !!!! !!! !!! !!!! !!!! !!! !!!! ! !!!!!!! !!!!!!!! !!!!! !!!!!!! !! !!!!! !! ! !!!! !! !!!! !!!! !! !! !!!!! !!! !!!!!!! !! !!!!!! !! !!!!! !!! !!! !!! !!!!! !!!!! !!!!!! !!!!!! !!!!!! !!!!! !!!! !!! !!!! !!!!!!! ! !! !! !!! !!!! !!! !!!!!! !!!! !!! !!! !!!! !!!!! !!! !!!!!! !!! !!!! !!!!!! !!!!!! !!!!! !!!!!! !!!! !!!!! !!!! !!!!!! !!!!! !!!!!! !!!!!! !!!!!! !!!!! !!!!! !!!!!! !!!!!!! !!!!!! !!!! !!!!!!! !!!!!!!!! !!!!!!!!! !!!!!! !!!!!!!!!!!!!! !!!!!!!!! !!!!!!!!!!! !!!!!!!!!!!!!! !!!! !!!!! !!!!!!! !!!!!!! !!!!! !!!! !!! !!!!! !!!!! !!!! !!!! !!! !!!!!! !!! !!!!! !! !!!!! !!!! !!! !! !!!! !!!! !!!!! !!! !!!!!! !! ! !!! !! !!! !!! !!! !!!!!!! !!!!! !!! !! !! !! ! !!!!! !!!! !!! !!! !!!! !!!!! !!! !!!! !!! !!!! !!! !! !!! !!!! !!! !!!! !!! !!!! ! !!! ! !!!!! ! !!!!! !!!! !! !! !!!!! !!! !!!! !! ! !!!! !!!!! !!!!! !!!!! !!!! !!!!! !!!!! !!!! !!!! !!!!! !!! !! !! !!! !! !! ! !!!!! !! !!!!!!!!! !!!!!!!! !!!!!! !! !!!! !!!!!! !!!!!!!!! !!!! !!! !!!! !!!!! !!!! !!!! !!!! !!! !!! !!! !!!!! !!!! !!!! !! !!!! !!! !!!!!!! !!!!!! !!!! !!! !!!! !!!! !!!! !!!! !!! !!!! !!!!!! !!!!!! !!!!!! !!!!!! ! !!!!! !!!!! !!!!!!! !!!!! !!!!! !!!!!! !!!!!! !!!!!!!!!!!!!!!!!!!!!!! !!!!!! !!!!!!! !!!!! !!!! ! !!!! !! ! !!!! !!!! ! !! !!!! !!!!! ! !!!!! !!! !! !!!! !!! !!!!!!!!!! !!!!! !!!! !!! !!!! !!!! !!! !!! !!! !!!! !!!!! !!!!! !!!! !! !! !! ! !!!! !!!! !!!!! !!! !! !!! !!! ! !!!! !! ! !!!! !!!! !!! !!!! !! !!!! !! !!!!!!! !!!! !!!!! !!!!!! !!! !!!! ! !!!!!!!!! !!! ! !!!!!! !! !!!!! !!!!! !! !!!!!! !!!!!!! !!! !!!!!!!! !!!! !! !! !!! !! !!! !!!!! !!!!! !!!! !!!!!! !!!!! !!!!! !!!!!! !!!! !!! !!!! ! !! !!!!!!!!!!!! !! !! !!!!!! !!!!! !!!!!!!!! !!! !!! !! !!! !!!!!! !!! !!!! !!!!! !!!! !! ! !!! !!!! ! !!! !!!!!! ! ! !!! !!! !!! !! !!! !!!!!!! !! !!!!! !!!!!!!!!!!!!!! !!! !!!! !! !!! !!!!!! !!!!!! !!!!!! !!!!!!!!!!!!!!! !!! !!! !!!!!!! !!!!!!! ! !! !!!!!! !! !!!! !!!!! !!!!!!! !!! !!! !! !! !!!!!! !!! !!!!! !!!!!!! !!!!!! !!! !!! !!! !!!!!!! ! !! ! ! !! !! !!!!! !! !!!!!!!!! ! !!! ! !!!!!!!! !!! !!!!! !! !! !!!! !!!! !!!!!!!!!!!! !! !! !!!! !!!! !!!! !!! !!!!!!!!!!!!!!!! !! !!! !! ! !! !! !!!! ! !!! !!! !! !! !! !!! !!!! !! !! !!!! !!! !!! ! !!!! !!! ! ! !!!!!! ! !!!!! !!! !! ! !!! !!! !!! !!! !!!!!!!!!! !!!!! ! !! !!!! !!! !!!!! !!!!! !!!!! !!! !!!! !!!!!! !!! !!!! !!!! !!!! !!!! !!!!!!! !!!!!!!!!! !!!! !!!!! !!!!! !!! !!! !!! !!!!!!!!!!!!!!! !!!!!!!!!!!! !!!!!! !!!!!! !!!!! !! !!!!!!!!!!! !!!!!!! !!! !!!! !!!! !!!!! !!!! !!! !!!!!!!!!!!! !!!!!!! !!!!!!!!!!!!!! !!!!!! !!!!! !!!!! !! !!!!! !!!!!!! !!!!!!!!!!!!!!! !!!!!! !!!!! !!! !!!! !!!!! !!!!!!!!! !!!!!!!! !!! !!!! !! !! ! !! !!!! !!! !!!!!!!! ! !!! !!!!!!!!!!!!!! !! !!! !!!! !!! !!!!!! ! !!!!!! !!!!! !!!!!!! !!! !!!!! !!!!!!!! !!!!!! !!! !!!!! !!!! !! !!!! !!! !! !!! ! !! !!!! ! !! !!!! !!! ! !! !!! ! !!!!! !!! !!! !! !!! !! !! ! ! ! !!! !! !! !!!!! !! !!! !! ! ! !! !!!! ! !!! !!!!! ! !!! !!!! !!! !!!! !! !!! !!! !!!!! !!!! !! !!!! !!!! !! !!!!!!!! !!!!! !!!!! !!! !!!!! !!!!!!!!! !!!!!!!!! !!!!!!! !!!!! !!!! !!!!!!!! !! !!!! !!!!! !!!! !!!!! !!! !!!! !!!! !!! !!! !! !!!!! !! !!!! !!!! !!! !!! !!! !!! !! ! !! ! !! !!! !! !! !!! !!! !!! !!! !!! ! ! !! ! ! ! ! ! ! ! !! ! !! !!!! !!! !! !!!!!!!! !! !!! !!!!! !!!! !!! !!! !!!! ! !!!!!!!!!!!!!!!!!!!!!!! !!!!!!!! !!!!!!!!! !! !! !! !!!!! !!!! !!! ! !! !! ! ! !!!! !!!! !! !! !!!!! ! ! !!! ! !! !!!! !!! !!!! !!! !! !!!! ! !!! !! !!!!! !! !!!! !! !!! !!! !!!! !!!! !!!! !!!! !!!!! !!!!! !!!! !!! !! !! !!! !!! !! !!! !! ! !! !!! !! !! !! !! !!! !!! !!! !!! !!! !!! !!! !!!! !!! !!! !!! !! !!! !!! !!! !! !!! !!! !!! !!! !!! !!! !!! !!! ! !!! !! !!! !! !!! !!! !!! ! !! !! !!! ! !!! !!! !! !!!! !! !! !!! !! !! !! !!!!! !!! !!!! !! !!! !!!!! !!!! !!!!! !!!! !!! !!!!! !!!! !!!! !!!!! !!!!!! !!!!! !!!!!!!! !!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! !!!!!! !!! !!! !!! !!! ! !!!! !!!!!!!!!!!! !!!!! !!!!! !!!! !!!! !! !!!! !!! !! !!! ! !!! !!! !!! !! !!! !! !! !!! !! !! !!! !! !!! !! !!! !!! !!! !!! !! ! !! !!! !!!! !!!! !!! !!!! !!! !!! ! ! !!! ! !!!! !!! !!! !!!! ! !! ! !!!! !! !!!!!!! !! !!! !!! !! !!!! !!!!! !!! !!!!!!! ! !!! !! !!!! !!!!!!!!! !! !!! !!! !! ! !!!!!!!!! !!!!!!! ! !! ! ! ! ! !!! ! ! !!!!!!!!!!! !!!!!! !! ! !!! ! ! ! ! ! !!!!! ! ! ! !!! !! !! !! ! ! !!! ! ! ! ! !! !!!! ! !! !! ! !!! !!!!! !!!!!!!! !!!!!!!!!!!!! !!! !! !!! !!!!! !!!!! !!!!!! !!!!!!!!! !!!!!!! !!!!! ! !!! !!!!!!! !!!!!! !!!!!!!!!!! ! !!!!! !!!!!!!!!!!!!!!! !!!!! !!!! !!!! !!!!! !!!! !! !!!!!!!!!!!!!!! !!! ! !!!!!! !! !!! !! ! !!! !!!!!! ! ! !!!! ! !!!! !!! !!!!! !!!! ! !!!! !!!!! !!!!! !! !!! !!! !!!!! !!!!! !!!!! !!!! !!!! !!!! !!! !!! !!! !! !!! !!!! !!! !!!! ! !! !!! ! ! !!!! !!!!! !!!! !! !! !! ! ! !!!!!! !!!!! ! !!! !!!! !!!! !!!! !!!!! !!!!! ! ! !!!!! !! !!!! !!!!!! !!!!!!!! !!!! !!!! !!!! !!!!!!! !!!!!!! !!!! !!!!!! !!!! !!!!!! !!!!!! !!!!! !!!!! !!! ! !! !!!! !!! !!!! !!!!!! !!!!!! !!!!!! !!! !!!! !!!! !!!!! !!!!!!!!!!!!! !!!!!!!!!!!!! !!!!! !!!!! !!!! !!!!!!!!! !!!!!!! !!!!!!!!!! !!!!! !!! !!!!!! !!!! !!!!! !!! !! !!!!! !!!!!!!! !!!!!!! !!!! !!!! !!!! !!!! !!!!! !!!!!! !!!!! !!!! !!!!!!!!!!!!!!!!!!!!!! !! !! !!!! !!!! !!!! !!! !! !!! ! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!!!!! ! ! !!! !!!! !!! !!! !!! ! !!! !!! !!! !! !!! !! !!!! !! !!!! !!! !! !!!! !!! ! !!! !!! !!! !!! ! !!! !! !!! !!! !!! !!! !! !!! !! !! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!!!!!!! !!! !!! !!! !!! !!! !!! !!! ! !!! !!! !!!! !! !!!! !!! !! !!!! !!!!! !! !!! !!!!! !!! ! !!! !! !!! !!! !! !! !!! !! ! !!! !!!!!!!! ! !!!! ! !!! !! !!! !!! !!! !!! ! !! !!! !!! ! !!!!!!!!! !!!! !!! !!!!!!!! !!!!!!!! !!!! !!! !!!! !!!!!! !!! !!!!! !!!!!! !!! !!!! !!!!!!!!!!!!!!!!! !!!!!!! !!! !!!! !!!! !!! !!! !!!! !!! !!! ! !!!!!! !!!!!!!!!!!!!!!! !!!!!!!!!!! !!! !!! !! !!! !!! !!! !!! !!!!! !!!!! !!!! !!!! !! !!!! !!!!!!! !!!!! !!!!! !!!!! !!! !!!!! !! !!!! !!! !!!!!!!!!!!!!!!!! !!!!!!!!!! !!!! !!!! !!!!! !!! !!! !!! !!! !!!! !!! !!!!!!!!!!!!!!!!!!!!!!!!!!! !!!! !!!!! !!! !!!! !!!! !!!! !!!!!! !!!!!!!!!!!!!!!!! !!!! !! !!!!!! !!! !!!!! !!! ! !!! !!! !!! !!!! !! !!! !!! ! !!! !!!! !!! !!!! !!!! !!!! !!! !! !!! !!!! !!!! !!!! !!!! !!!!!! !!! !!!! !!!!!!!!!!! !!! !! !! !!!! !!!! !! !!! !! !!! !!! !!!! !! !!! !! !! !! !! !! !! !!! !!!!!!! !!! !! !!! !! !!!! !!! !! !!! !!! !!!! !!! !! !!!! !!! !!! !! !!! !!! !!! ! !!!!! !!! !! !! !!! !!! !!!! !! !!!! !!!! !!! !!! !!! !!!! !!! !!! !!! !!!!! ! !!! !!! !!! !!!! !!!! !!! !!!! !! !! !! !!!! !!! !! !!! !!! !!!! !!!! !!! !!! !!!! ! !!! !!! !!! !!! !!!! !!!! !!!! !!!! !!! ! !!! ! !!!! !! !!!! !!!!! !! !!!!!!!!!!!! !!!!!!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!!!!!!!!!!!!!!!!!!! !!!!!!!!!!! !!! !!! !!!! !!!! !!!! !!!! !!! !!! !!! !!!!!!!!!!!!! !!!!!! !!!!!!! !!!!! !!!!!!!!!! !!!! !!! !!!! !!! !!! !!! !!!!! !!! !!! !!! !!! !! !! !!! !!!! !!!!! !!! !!! !!! !!! !!!! !!! !!! !!!! !!!! !!!!!!!!!!! !!!!!!!!! !!!!!!!!!!!!!!!!! !!! !!!! !!! !!!! !!! !!! !!!!! ! !!! !!!!!!!!!!!!!!!! !!!!!!!!! !!!!! ! !! !!! !!! !!! !! !! !! !!!! !! !!!! !!!! !!!! !!! !!! !!!! !!!! !!!!!! ! !!! !!!!!!! !!!!!!! !!!!!!! !!!!!! !!!!!! !!!!! !!!!!!!!! !!!!!!!!!!! !!!!!!!!!!!!! !!!!!!!!!!!! !!!!!!!!!!!!! !!!! !! ! ! !!!! !!!!!!! !!!! !!!!! !!! !!!! !!! !!! !! !!! !! !!! !!! !!!! !!! !!! ! ! !!!!!! !!!! ! !! !!!!!! !!! !!!! ! !!!! !!! !!!!! !!!!! !!!!!! !!!! !!!! !!!!!! !!!! !!!!! !!!! !!!!!!!! !!!!!!!!!!! !!!!!!!! !!!!!!!!!!!!!!!! !! !!!!!!!!!!!!! !!!!!!!!!!!!!!!!! !!!!!! !!!!!!! !!!!! !!!!!! !!!!!!! !!!!! !!!!!! !!!! !! !!! !!!! !!!! !!!!!! !!!!!! !! !!!!! !!! ! ! !! ! ! !! !!!!! !!!!!!! !!!! !!!!! !!!! !!!!! !!!! !! !!!!!! !!!!!!! !!! !!!!!!! !!!!!!! !!!!!!! !!!!!!! !!!!!!!!!!! !!!!!!!!! !!!!!!!! !! !!!!! !!!!! !!! !! !!! !!!!!!!!!!!!!!!! ! !! !! !! !! !!! !! !!! ! !!!! !!!! !!! !!! !!!!!! !!!! !! !!!! !!! !!!! !!!! ! !!!! !!! !!!! !!! !!! !!!!! !!! !!!! !!! !!! !! !!! !!! !!!! !!! !!!! !!! !!!! !!!! !!! !!! !!!! ! !! !!!! ! !!! ! !!!! !!! !! !! !!! !!! !!!! !!! !! !! !!! !!!! !!!! !! !!!!! !!! !!! !!! !! !!! !!!! ! !! !!!!! !!! !!!! !!! !!!! !!!! !!!! !! !!! !!!! !!! !! !!! !!!! !!!!! !!!!! !!! !! !!! !!!! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!! !!!!!! !!! !!!!!!!!!! !!!! !!!! !!!!! !!!!!!!!!!!!!!!! !!!!! !!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!! !!!! !!!!!!!!!!!!!!! !!!!!!!!!!!!!! !!!!!!!!!!!!! ! ! !!!!!!! !!!!!!!!!!!! !!!!!!!! !!!!!!! !!!! ! !!! !!!! !! !!! !!! ! !!! !! !! !!!! !!!!! !!!!!! !!! !!!!! !! !!! !!!! !!! !! !! !! !!!!!!!!!!! !!!! !!!!! !!! !! !! ! !! !!!!!! !!!!! !! !! !! !! !!! ! !!! !! ! !!! !!!!!! !!! !!! !! !!!! !!! !!! !!!!! !!! !!! !!!! !!! !!!!!! !! !! !!! !! !! !! !!! !!!! !!!!!! !!!! !!!!!!!! !!!!!! !!!! !!!! !! ! !!! !!!!!!!!! !!!!!!!! !!!!! !!!!!! !!!!!!! ! !!! !!!! !!! !!!! !!!!! !!!!!! !!!!! !!!!!! !!!! !!!!!! !!!! !! ! !!! !! ! !!!! !!!!!! !!!!!! !!!!! !! ! !!!!!! !! !!! ! !!!! !!! !!!! !! !! !! !!!! !!!!!! !!!!! !! !!! !!!! !!!!!!!! !!!!! !!! !! !!! !! !!!!!!!!!!!!!!!!!! !!!! !!!! ! !!!!!!!!!!!!!!! !! !!! !!!! !!! ! !!!! !!!!! !!!!!!!! !!! ! !! !! !!!!! !!! !! !!! !!!!!! !! !!!! !!! ! ! !!!!! !!!! !!! !!! ! !!! !! !!!! !!! !!!! !!!! !!! !!! !!! !!!!! !!! !!! !!! !!! !!! !!! !!! !!!! !! !! !! !!!!!!!!! !!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!! !!!! !!!!! !!!!!! !!! !!!!!!!!! ! !! !!! !! !! !!!!!!!! !!!!!!!!!! !!!!!!!!! !!!!! ! ! ! !!!!!!! !! !!!!!!!!!!!!!!!! !!! !!!!! !!! !!! !!!!!! !! !!!!!!!!!!! !! !!!!!!!!!!!!!!!!!!!!! !!!!! !!!! ! ! !!!!!!!! !!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! !!!! !! !!!! !!!! !!!! !! !! ! !!!! !!! !!!! ! !!! !!!!!! !!!!!!!!!! !! !!!! !!! !!! !!!!!!!!!!!!! !!!!!!!! !!!! !!! !!! !!!! !!! !!!!!! !!!! !!!!! !!!! !!!! !!! !!! !!! !! !!! !!!!!!!!!!!!!!!! !! !!! !! !!!!!!!!!!! !!!! !!!! !!!! !! !!! !!! !!! !!!!!!!!!!!!! ! !! !!!!!!!! !!! !!! ! !!!!!! !!!!!!! !!!! !!!! !!!!! !! !!! !!!!! !!!! !!!!! !!! !!!!!! !!! !!!! !!! !!! !!!!!! !!! !!! !!!!! !!!! !!!!! !!! !! !!!! !!! ! ! !!!!!! !!!!!! !!!!!! !!!!!! !!! !! !!!!!! !!! !!!!! !! !!!! !! !!! !!! ! ! !! !!! !!!! !!! !!! !!! !!!!!! !!! ! ! !!!!! !!! !! !!!!! !!!!!!! !!! !!!!!! !!!!! !!!!!! !!!!!! !! !! !! !!!!! !!!!!!! !! !!!!!! !! !!!!!!!! !! !!! !! !! ! !! !! !!!! ! !! !!!!!! !!!!!!! !!!! !!! ! !!! !!!! !! !!!!!!!!!! !!!!!!!!! !! !!!!!! !!!!!!! !! !!!!! !!!!!!! !!!!!!!! !!!!!!! ! !!! !!!!! !! !!!! !!!!!!!! !!!!! ! !!!! !!!! !!!!!!!!! ! !!! !!! !! !! !!!!!!!!!!! !!!!!!!! !!!!! !!!! !!!!!! !!!!!! ! !!!! !!! !!!!! ! !!!! !!!! !! !!! !!!! !! !!! !! !! !!!! !!!! !!! !!! !! !! !! !!! ! !!!! !!! !! !!! !!! !!!! !!! !! !!!! !!! !!! !! ! !!! ! ! !!! !!!! !! !!!!! !!!! !! !! !!! !!!! ! !!! !!!! !! !!! !! !!!! !!! ! !!! !! !!! !!! !! !! !!! !!! !!! ! !!! ! !!! !!! !!! ! !! !!!! !!! !! ! ! ! ! !!!! ! !!! ! ! !!!!!! !!! !!! ! !!!!!! !!!!! !! !!!!! !! ! ! ! !!!!!! !!!!!!! !!!!!!!! !!! !!!!! ! !! ! !! ! ! !! !! !!! !! ! !! ! !!! !!!! !!!!!! !! ! ! !!!! !! ! !!! !! !! !!!! !!! !! ! ! !!!!!!!!!!! !!!!!!!!! ! !!!!! !!!! ! !! !!! !! ! !!!!!!! !! !!!!! !! !!!!!!! !!!!!!!!!!!!!!!!!!!!!!!! !! !!!!!! !!!!! !! ! !! !! !!!! !!!!! !! ! !! !! ! !! !! !!!!!! !!! !! ! !! !!!! ! !!!! !! !! !! !! !!!!! ! !!! ! ! ! !!!! !! !! !! ! !! ! ! ! ! ! !!!!!! !! !!!!!! !!! ! !!!! !! ! !!! ! !!! !! ! ! !!! !!!!! !! !!!! !!! !!! !!! !!! !! !! ! !!! !!! !!!! !! !!! ! ! !!! ! !!! ! !!!! !!! !!!!! !! !!!!! !!! !!!! !! !! !!!!! !!!! !!!!!! !! !! !!!!!! ! !!! ! !!!! !!!! !! !!!! !!! !! ! !!!! !! !! !! ! !!! !!! !! ! !! !! !! !!! !! ! !! !!! !! !! ! !!! !!!!!! ! ! !!! !!! ! ! ! !!!! !!!! ! !!!! !!! !! !!! !!!! !! !!!!! ! !!!!!! !!!!!! ! !!! ! ! ! ! !!!!! !!!! !!!!!!!!!!!!!! !!!!!!! !!! !!! !! !! !!! !!! ! !!!!!!! !!!! !!!!!! ! ! !!!! !!! !!! !!!!! !! ! !!! !!! !! !!! !!! !! !!! !!!!! !! !! !!! !!! !! !! !! !!! !! !!!!!!! !!! !! !!!!! !!!!!! ! !! !!! !!! !!!! !!!!! !!! ! !! !!!!!!!!!!!! !! !!! !! !!!!!!!!!!! !!! !! !!!!!! !! !!!!!!!!! !!! !!! !! !!! ! !! !!!!! !! ! ! !!!! ! !!!!!!! !! !!!!! !! !! !!!!! !!! !! !!! !! !!!! !!! !!!! !!!! ! !!!!!!!! !! !!!!!! !!!! ! ! !!!!!! !! !! !! !!! !!!! !!! !!!! !! !!!! ! !!!! ! !!! !! ! !!!!!! !!!! !!! !!!! !! !!!!!!!! ! !!!!!!!!!!!!! !! !!!!!!! !! !!!!!! !!!! !!!! !!!!! !! !! ! !! !! ! !! !!!! !! ! ! ! !! !!! !!! !! ! ! !!! !!!! !! ! ! ! ! ! ! !!!! !!! !!!!!!!!!! !! !! !!!!! !!!! !!!! !!!!!!! !!! ! !! ! !! !!! !!! !! !!! !!! !! !!!! !!!!!! !! !!! !!! !! ! !!! ! !!!! !! !!! !! ! ! !! ! ! !!!!! !!! !!! !!!!! !!!!! ! !!!!! !!!! ! !!!! ! ! ! ! !!! ! !!!!!!!! !!!! !! !! !!! !! ! !!! !! !! !! !! !!! !!!!!!! !!! !!!!!!! !!!!! !! ! !!! !!! !!!! !! !! !!!!! !!! ! !!!! !!!! !! !! !!! ! !!! !!! !!! !!!!! !!! !! !!!! !! !!! !!!! !!!! ! !!!!! !! ! !!!!!!!!!! !! !!!! !!!!! !! !! !!!!!!!!!! !!! !! !!! !! !!! !!! !!! !!! ! !!!!! !! !! !!!!!!!!!! !!!!!!!!!!!! !!!! !!! !! !!!! !!! !! !! !! !!! !!!! !! ! !! !! !! !!!! !! !!!! !!!!! !!!!! !!!!! !!!! !! ! ! ! ! !! ! !!! ! ! !!! !!!! !! ! !! ! ! !!! !!! !! !! !! !! !!! !! ! !! !! !! !!! !!!! !! ! ! !!! !! ! !!!! !! !!!!! !!!! !!! !! !! !!! !!!! !!!!!! !!!!!!!! !!!!!!! !!!!!!!!!! !!!!!!!!!!!!! !!!!!! !!!!! !!!! !!!!!!!!! !!! !!!!!!! !!!!! !!!! !! !!! !! !! !! !! !! !! !! !! ! !! !! !! !!! !! !! !! !! !! !!! !! !!!! !!!! !!!!!! ! ! !! ! !! !! ! !! !! !! !!!!! !!!!!!!!!!!! !!!! !!!! !!!! !!! !!!! !!!!! !!!!!!!!!! !!!!!!!!!! !!!!!!!!!! !!!!!!!!! !!!! !! !! !! !!!! !!!! !! !! !! !!!!!! ! !! ! ! ! ! !!!! ! ! ! !!!!!!!!!! !!!! !! !!! !!! !!! !!! ! !!! ! !! ! !!! !!! !!!!!! ! !! !!! !! ! ! !!! ! !! ! ! !!! ! !!!! !! !!! ! !! !! !! !!!!!!!!! ! !! !!! !!! !! ! !!! !!!!! !! !!!!! ! !! !!!!!! ! ! ! !! ! !!!!! !! !!! !!! !! !!!!!!! !! !!! !!!! ! !! ! ! ! ! !!! ! !! ! !!! !!!! !!! !! !!!!!!!!!!!!! !!!!!!!!!! !!! !!!! ! ! !! !! !! ! ! !!! ! !!!!! !!! !! ! !! !! ! !! !! !!! !!!!! ! !!! !!!!! !! !!!!! !! !!!!!! !!! !! !!!! !!!!!!!!!!!!!!!!!!!! !!!!! !! !! !! !!!!!! !!!! !!! !! !!! !!! !! !!! !! !!! !! !! !! !!!! !! !! !!!! !!! !! !! !! !!!!! !!! !! !!!!!!! !! !! !!! !!!! !! !! ! ! !! !! !!!!!!! !! !! ! !! !!!! !!!!!!!! !!!! !! !!! !!! ! !!! !!!!! !!!! ! !!! !!! !!! !! ! !!!!! !!!!!! ! !!!!!!!!!! ! ! ! !!! !! !! !!!! !! !!!!! !! !! !!! !!! ! !! ! !! ! ! !!! !! ! !! ! ! ! ! !! ! ! !! !!! !!! !! !! !!!! !!! !!!!! !!!! !!! !!! ! !!!! !!!!!! !!!!! !!!! !!! !! !!! ! ! ! !!!! ! !!!!!! ! !!!! !!!! !!!!! !! !!! !!! ! !! !! !!! !!! !!! !!!!! !! !!!!!!!! ! !!!! !!!!!!!! ! !! !!!! !! !!! !!!! ! ! ! !!! ! ! !! !!!! !!!! !! !! !! ! !!! !!! !!! !!!!!!!!!!!!!!! ! !!! !! !! ! !!!!!! !!! !!!!! !! ! !! !!!! !!! !!!!!!! !!!!!!!! !! ! !!! !!!!!!!!! ! !! !! !!! !!!! !!! ! !!!! !!! ! !!! !! !!!! ! !!!!!! !!!!!! !!!!!!! !!!!! !!!! !! !!!!! !!! !! ! !!!!!!!! !! !!! !!! !!! !!!! !!! ! !! !!! !! !!! !! ! ! ! !!!! ! !! !! !!!!!!!!! !! !!!!!!!! !!!!!! !! !! ! !!!!!! ! !!! !!!!!!! !!! ! !! !!! !! 00.50.25 Kilometers Model Thalweg Measured Thalweg Boundary a) b) Figure 5.37 Spatial distribution of prediction errors at Site 1 of Brazos River. (a) Lightly and darkly colored points cover the area with prediction errors in Group 1 and Group 2, respectively. (b) Measured thalweg and the thalweg modeled by RCMM 180 From Figure 5.37a and b, there appears to be a correlation between the prediction errors and the thalweg discrepancies. Most of the prediction errors in Group 2 lie in those areas where the thalweg described by RCMM does not match closely with the actual thalweg measured in the field. This makes sense because the parameters of beta distributions that are used to describe the river cross-sections are related to the thalweg location. If the thalweg is not located precisely, RCMM is unable to create a reasonable description of cross-sections using the beta distributions. Therefore, an accurate prediction of thalweg location is a key in producing a good description of channel bathymetry using RCMM. 5.6 APPLICATION OF RCMM TO GUADALUPE RIVER AND SULPHUR RIVER IN TEXAS So far, the bathymetry data at Site 2 and Site 1 along the Brazos River are used in developing and testing RCMM, respectively. In this section, the application of the RCMM to two more rivers in Texas, Guadalupe and Sulphur, is discussed. This task is carried out to test the applicability of RCMM on rivers other than the Brazos River, which is used to develop the model. The location of the study areas on the Guadalupe and Sulphur River with respect to the Brazos River are shown in Figure 3.1. 5.6.1 Hydraulic Geometry Relationships As shown in Figure 5.38, the study area along the Guadalupe River which is located between the USGS gaging stations at New Braunfels (# 08168500) and Gonzales (# 08173900) is 1.2 km long with an average width of about 35m. Hydraulic geometry relationships are developed for these gaging stations and the parameters are linearly interpolated based on the upstream watershed area to get the hydraulic geometry relationships at the study area location. The average width and the depth at the study area are 33m and 2.3m, respectively. The width is slightly modified to incorporate additional information by using the DOQ of the study area. 0 5025 Kilometers E E E E E E Cuero Sattler Gonzales Victoria Kerrville New_Braunfels E E Hunt Comfort 020100 Meters Figure 5.38 Location of study area along the Guadalupe River The study area along the Sulphur River is 1.4 km long with an average width of about 30m. The Sulphur River has only one USGS gaging station in operation which is located near Talco, TX (# 07343200), and is about 100 km upstream of the study area (Figure 5.39). Due to the absence of more than one gaging station along the Sulphur River, the values of width and depth obtained 181 from the hydraulic geometry relationships at the gaging station are slightly modified to match the observed average values at the study area. The average values obtained by using the hydraulic geometry relationships at the Brazos and Guadalupe River do match with the observed data. Therefore, in the absence of information, using the observed data to slightly modify the hydraulic geometry results seems appropriate for the study area along the Sulphur River. Accordingly, the average width and the depth for the study area along the Sulphur River are 30m and 4.5m, respectively. 0 10050 Kilometers E Talco 0250125 Meters Figure 5.39 Location of study area along the Sulphur River 5.6.2 Thalweg Prediction First, the relationship between the dimensional radius of curvature and the thalweg location (equations 5.5 and 5.6) are used without any change to predict the thalweg for the reach along the Guadalupe River, and Figure 5.40 shows the 182 result. From Figure 5.40, it is clear that the relationship derived for the Brazos River does not work on the Guadalupe River. 0 10050 Meters Modeled Thalweg Measured Thalweg Boundary Figure 5.40 Thalweg prediction along the Guadalupe River using RCMM The Guadalupe River is narrower than the Brazos River, and the discrepancy between the measured thalweg and the thalweg predicted by RCMM may be due to the dimensional form of radius of curvature versus thalweg relationship. To verify this, the thalweg at the Guadalupe River is also predicted by using the non-dimensional equations 5.7 and 5.8. Figure 5.41 shows the result. 183 0 10050 Meters Modeled Thalweg Measured Thalweg Boundary Figure 5.41 Thalweg prediction along the Guadalupe River using equations 5.7 and 5.8 Although the prediction of thalweg by equations 5.7 and 5.8 is improved at some locations (Figure 5.41) compared to prediction by equations 5.5 and 5.6 (Figure 5.40), the overall result is still unsatisfactory. Similarly, for the Sulphur river, the thalweg prediction by using equations 5.7 and 5.8 is found to be unsatisfactory (Figure 5.42). Although the non-dimensional form of radius of curvature versus thalweg relationship provide improvement in the thalweg prediction compared to dimensional form, it is still unable to capture the details necessary for describing the relationships between the channel planform and thalweg location. 184 0 10050 Meters Measured Thalweg Modeled Thalweg Boundary Figure 5.42 Thalweg prediction along the Sulphur River using equations 5.7 and 5.8 According to the conceptual model used to develop RCMM (Figure 5.1), the thalweg location depends only on the meandering shape of the channel. This, however, seems to be untrue from results obtained by application of the model to Guadalupe and Sulphur rivers. The simple assumptions made in the conceptual model, however, need to be modified to make the RCMM generic. The first addition to RCMM is the incorporation of sinuosity. In Guadalupe and Sulphur Rivers, the measured thalweg lies in the middle of the channel even at meandering bends, which is in contrast to RCMM?s conceptual model. If Site 2 of the Brazos, Guadalupe, and Sulphur Rivers are split into segments that have a length of about 10 to 14 times the channel width (Meander wavelength = 10-14 times channel width), the average sinuosity values for the corresponding river segments are 1.2, 1.5, and 1.8, respectively. A 185 186 decrease in the channel slope causes an increase in the sinuosity (Knighton, 1998). Sinuosity is the ratio of channel length to the straight-line valley length (Figure 2.2). Therefore, as the slope (elevation difference/channel length) decreases, the channel length increases, and the sinuosity increases. A decrease in slope lowers the flow velocity thus reducing the eroding effects along the bank of the channel. So, one theory can be formulated by using the channel sinuosity. The increase in channel sinuosity (decrease in slope) reduces the eroding effects of flow along the banks thus inhibiting the movement of the thalweg towards the bank, as observed in the Guadalupe and Sulphur Rivers (Figure 5.40 ? 5.42). To fit the thalweg on the Guadalupe River, the parameters in equations 5.5 and 5.6 are manually modified. Equations 5.5 and 5.6 can be expressed as: t * = a*ln (r) + b 0 < r < r max (5.28) t * = c*ln(abs(r)) + d -r max < r < 0 (5.29) Table 5.3 shows the values of a, b, c, d, and r max for Brazos and Guadalupe Rivers. Name Sinuosity A b c D r max (m) Brazos 1.2 -0.076 1.21 0.087 -0.3210 12500 Guadalupe 1.5 -0.070 1 0.07 0.0061 1500 Table 5.5: Parameters of radius of curvature ? thalweg location relationship for Brazos and Guadalupe River The corresponding plot of equations 5.28 and 5.29 along with observed thalweg locations for the Guadalupe River are shown in Figures 5.43. In addition, the parameters a, b, c, and d for the Guadalupe River also fits the data for the Sulphur River (Figure 5.44). 0.00 0.25 0.50 0.75 -3000 -2000 -1000 0 1000 2000 3000 Radius of Curvature (m) 1. 0 0 Tha l w e g Lo c a t i on , t * Model Observed Data Figure 5.43 Radius of curvature versus thalweg location for the Guadalupe River 0.00 0.25 0.50 -3000 -2000 -1000 0 1000 2000 3000 Radius of Curvature (m) 0.75 1. 0 0 Tha l w e g Lo c a t i on , t * Model Observed Data Figure 5.44 Radius of curvature versus thalweg location for the Sulphur River Although the thalweg prediction along the Guadalupe River by using equations 5.28 and 5.29 is significantly better (Figure 5.45) compared to Figure 5.40, there are still some areas (circled in Figure 5.45) where the thalweg predicted by the model does not match with the observed thalweg. 187 0 50 100 Meters Boundary Measured Thalweg Modeled Thalweg Figure 5.45 Thalweg prediction along the Guadalupe River after incorporating the sinuosity criteria. Circled areas show locations where the thalweg prediction is imprecise The thalweg in a natural channel is a result of a combination of complex erosion and deposition processes that are not accounted for RCMM. The inability of RCMM to locate the thalweg can also be attributed to other factors such as geology, climate, etc that are not included RCMM. Figure 5.46 shows the thalweg prediction for the Sulphur River. 188 01050 Meters Measured Thalweg Modeled Thalweg Boundary Figure 5.46 Thalweg prediction along the Sulphur River. Circled areas show locations where the thalweg prediction is imprecise As mentioned earlier, the average sinuosity for Brazos is lower than the sinuosity of Guadalupe and Sulphur Rivers. The parameters in Table 5.5 can be related to sinuosity. For example, parameters a and d increase with sinuosity, and parameters b and c decrease with sinuosity. 5.6.3 Description of cross-sections The beta distributions describe the general shape of the cross-sections based on the thalweg location, and the general cross-sectional shape is assumed to be not very different for different channels. Therefore, the same parameters that are used for describing the channel bathymetry along the Brazos River are used 189 190 for the study reaches along the Guadalupe and Sulphur River. The RMSE* for the Guadalupe River is 1.16, which is relatively higher compared to the Brazos river. Likewise, the R 2 value is also very low (-0.34). The R 2 value is significantly affected by the data that are not well described by RCMM. For example, in the case of the Brazos River, the R 2 value improved significantly after the prediction errors are separated into two different groups. The first group contained the data with low prediction errors, and the other group contained the rest of the data. Similarly, for the Guadalupe River, the prediction errors are separated into two different groups. Group 1 contains the data with low prediction errors (-1.0,1.0), and Group 2 contains the rest of the data. The R 2 values for Group 1 and Group 2 are 0.43 and ?0.2, respectively. Similar to earlier observations with the Brazos River, most of the prediction errors in Group 2 lie in the areas where the thalweg prediction by RCMM is imprecise. Figure 5.47 shows the reach along the Guadalupe River with the data in Group 2 colored in dark. In Figure 5.47, the circled area has a big dip in the bathymetry, which the model is unable to describe precisely. Measured Thalweg Modeled Thalweg Group 1 Group 2 01050 Meters Figure 5.47 Spatial description of prediction errors for the Guadalupe River. Group 1 includes points with low prediction errors, and Group 2 includes points with high prediction errors The RMSE* for the Sulphur River is 0.69, and the R 2 value is 0.53. The RMSE* and R 2 values for the overall dataset along the Sulphur River are very good given the previous results along the Brazos and Guadalupe River. This is due to the fairly regular bathymetry of the Sulphur River without any dips. Even with this high R 2 value, the data are split into group 1 (with low prediction errors) and group 2 (with rest of the data) to see the spatial distribution of prediction errors. Similar to earlier observations along the Brazos and Guadalupe River, high prediction errors (circled in Figure 5.48) lie in the areas where the thalweg location described by RCMM is imprecise. 191 192 Measured Thalweg Modeled Thalweg Group 1 Group 2 Figure 5.48 Spatial distribution of prediction errors for the Sulphur River. Group 1 includes points with low prediction errors, and Group 2 includes points with high prediction errors 5.7 REGIONAL APPLICATION OF RCMM Sections 5.4 and 5.6 discuss the application of RCMM to small reaches along the Brazos, Guadalupe, and Sulphur Rivers. In this section, the application of RCMM to a 300 miles segment along the Brazos River is presented (Figure 5.49). Since there are no bathymetry data along this reach, the accuracy of the results is not verified. However, the data at the gaging stations, which are measured, are used to generate the three-dimensional river channel bathymetry. The idea is to demonstrate the applicability of RCMM in generating three- dimensional river channel bathymetry at regional scale. The channel description for the 300 miles segment along the Brazos River is generated to study the frequency of flooding of oxbow lakes along the area of interest, and such studies do not require a very accurate description of channel bathymetry. Studies like this (flooding of oxbow lakes) or other studies involving modeling of some water quality constituents do not require a very accurate river channel bathymetry, and an approximate description of channel form is sufficient. Therefore, even an approximate representation of river channel bathymetry produced by RCMM is useful. 0 10050 Kilometers 08116650 08114000 08111500 08108700 020100 Kilometers Gaging Stations Brazos River Brazos Watershed 08116650 08114000 08111500 08108700 08098290 08096500 08093100 08091000 08090800 08089000 Figure 5.49 Lower reach of the Brazos River 193 Figure 5.50 shows a view of the channel bathymetry that is produced by RCMM at a location near Richmond gaging station (#08114000). 0105 Kilometers 08116650 08114000 Figure 5.50 River channel description by RCMM at a location near Richmond gaging station. The cross-sections are approximately 500m apart. It is clear from Figure 5.50 that RCMM is a useful tool for describing the three-dimensional form of river channel at regional scale. 5.8 MODELING OF PREDICTION ERRORS As discussed in section 5.1, the channel bathymetry can be represented as a combination of deterministic and stochastic components (equation 5.1). RCMM describes a mean surface ( z ) for the channel bed thus contributing towards the description of deterministic component of equation 5.1. Although the primary objective of this research is to describe the deterministic component, to make the description of the channel bathymetry complete, a brief discussion on the stochastic component (?(x,y)) is covered in this section. ),(? yx 194 If RCMM captures all the spatial trends in the bathymetry, the prediction errors should have no spatial correlation (white noise). However, it is usual to observe some small spatial correlation even after removing the trends (Kanevski et. al, 2002). The first task, therefore, is to study the spatial correlation of the prediction errors. The geostatistical analyst extension in ArcGIS is used to fit the semivariograms for studying the spatial correlation of the prediction errors. To get a better sense of the errors, the data are analyzed in (s,n,z) coordinate system so that semivariograms can be fitted in both n- (across the flow) and s- (along the flow) directions. First, the errors from calibration data (Brazos River Site 2) are analyzed. Figures 5.51 and 5.52 show spherical semivariograms (yellow lines) for the errors in n- and s-directions, respectively. Model Data (m 2 ) (m) Figure 5.51 Semivariogram of prediction errors (across flow, n-direction) for Site 2 on Brazos River 195 196 Model Data (m) (m 2 ) Figure 5.52 Semivariogram of prediction errors (along flow, s-direction) for Site 2 on Brazos River Spherical semivariograms are also fitted to the prediction errors along the Brazos River on Site 2. The semivariogram parameters for both Site 1 and Site 2 along the Brazos River are shown in Table 5.6. Brazos River at Site 2 Brazos River at Site 1 Parameters n-direction s-direction n-direction s-direction Range 80 88 60 60 Partial Sill 4.5 3.2 1.8 1.2 Nugget 0 0.12 0 0.09 Table 5.6: Semivariogram parameters for prediction errors at Site 1 and Site 2 along Brazos River Similarly, spherical semivariogram parameters for residuals along the Guadalupe River and Sulphur River are shown in Table 5.7. 197 Guadalupe River Sulphur River Parameters n-direction s-direction n-direction s-direction Range 24 21 30 24 Partial Sill 0.6 0.4 1.6 0.8 Nugget 0 0 0 0.1 Table 5.7: Semivariogram parameters for prediction errors Guadalupe and Sulphur Rivers It is clear from Tables 5.6 and 5.7 that the prediction errors do have a spatial correlation. For the residuals at all sites, the range in both n- and s- directions is nearly the same, which suggests RCMM is doing well in modeling the deterministic trend. However, the partial sill in the n-direction is still a little higher than the semivariance in the s-direction. The stochastic small-scale spatial variations exhibited by the prediction errors can be handled in several ways: spectral methods, moving average processes, geostatistical methods, etc. Christakos (1992) provides a good overview of these different methods. Geostatistical methods would be the best way to deal with errors in a GIS environment. The current version of geostatistical analyst in ArcGIS does not offer any functions with regard to simulation of random fields. However, there are external libraries such as GSLIB (Deutch and Journel, 1992) and GSTAT (Pebesma and Wesseling, 1998) which offer a technique called Sequential Gaussian Simulation (SGS), and SGS can be used to simulate random fields. A SGS model requires two inputs: semivariogram (Figure 5.50 or 5.51) and the extent of area to be simulated. In the absence of any pre- defined condition or data, SGS starts with a single cell (data point) at any random 198 location in the area of interest, and employs the input semivariogram to predict a value at that point using the kriging technique. This first point is then used to condition the prediction value at the second point, the first two points are then used to condition the prediction value at the third point, and this process continues until all the points are visited only once. A single prediction for each point constitutes a single realization of the procedure. Multiple realizations are performed to get a distribution at each point, and a value is picked for each cell. The output from the SGS process can be added to the mean surface obtained from RCMM to get the final result, which will then have both the deterministic and a realization of the stochastic component of the channel bathymetry. As mentioned earlier, the SGS process in not available with ArcGIS, and it will require considerable time even to employ libraries such as GSLIB or GSTAT to perform SGS in ArcGIS. 5.9 SUMMARY This chapter answers the second question (Given the standardized representation of channel bathymetry, how can this information be used to describe the three-dimensional structure of river channels in areas with no bathymetric data?) posed in section 1.4. An analytical model (RCMM) is developed to relate the channel planform with cross-sectional form so that the channel planform can be used to describe the cross-sections in areas with no bathymetric data. The applicability of RCMM is demonstrated by using it to describe the river channel bathymetry along Brazos, Guadalupe and Sulphur Rivers. RCMM is a novel approach that incorporates four different concepts for 199 extrapolating the river channel bathymetry. These concepts are: 1) analysis of river channel bathymetry in a normalized domain; 2) relationship between the radius of curvature of the channel planform and the thalweg location; 3) beta distributions for fitting river channel cross-sections; and 4) use of hydraulic geometry relationships for scaling normalized cross-sections. Although the results of application of RCMM to Brazos, Guadalupe, and Sulphur Rivers in Texas are promising, RCMM is limited by its ability to model all the variables that are involved in the complex river geometry. This limitation may be overcome by incorporating geology, land use, climate and other factors that are responsible in the development of river channel cross-sections. 200 Chapter 6 Conclusions and recommendations This dissertation has two main components. The first component deals with a procedure to create a three-dimensional of river channel bathymetry in the form of cross-sections and profile-lines. The second component deals with development of an analytical model for describing the three-dimensional form of river channels. The conclusions drawn from each component are discussed in sections 6.1 and 6.2, followed by recommendations for future work in sections 6.3. 6.1 THREE-DIMENSIONAL DESCRIPTION OF RIVER CHANNELS USING CROSS- SECTIONS AND PROFILE-LINES The PolylineMZ geometry available in the ArcGIS system is used to describe the three-dimensional features of river channels, such as thalweg, cross- sections, and profile-lines. The thalweg is used as a reference to assign (s,n) coordinates to the bathymetry points. For a given dataset, the location of the thalweg is non-variable. Therefore, using the thalweg as a reference ensures unique (s,n) coordinates for any given point in the river channel. In addition, the thalweg is represented as a smooth Bezier curves rather than sequences of straight line segments. This overcomes the problem of calculating the (s,n) coordinates for the bathymetry points lying adjacent to the sharp corners at the intersection of polyline segments. The regular FishNet tool that is currently available in ArcGIS does not produce a network of cross-sections and profile-lines in the Cartesian coordinate system. However, when a FishNet is created in the (s,n,z) coordinates, it provides Comment: 3D what? 3D model? Framework? Comment: Either: ?a smooth Bezier curve? or ?smooth Bezier curves (no ?a?)? 201 a network of cross-sections and profile-lines. Since the FishNet is created from a surface, creating an accurate surface representation using the points is an important step in the overall process. The river channel bathymetry is anisotropic, which means it is less variable in the direction of flow compared to the direction across the flow. In the (s,n,z) domain, the channel is always straight, and this helps to account for the anisotropy in the bathymetry data. The elliptical inverse distance weighting interpolation scheme developed in this research takes into account the anisotropy in the data. The elliptical inverse distance weighting scheme performed better than the regular inverse distance weighting scheme, which does not account for anisotropy. However, anisotropic kriging which also accounts for anisotropy in the data, worked the best among different interpolation schemes tested in this research. Anisotropic kriging is followed by elliptical inverse distance weighting with a difference of only one cm in the RMSE. Considering the simple procedure and computing time, elliptical inverse distance weighting can be used with sufficient confidence in the absence of the geostatistical analyst extension in ArcGIS, which has the kriging functionality. In conclusion, the anisotropic nature of the river bathymetry should be taken into consideration while interpolating the bathymetry points to get a good interpolated surface. The representation of channel bathymetry as a network of cross-sections and profile-lines (FishNet) is useful in several ways. First, it is easy to store the PolylineMZ features in a geodatabase compared to raster grids and triangulated irregular networks (TIN). Second, the bathymetry is easily rendered in a three- 202 dimensional environment, in contrast to a raster or TIN representation which may take a long time to be displayed. Finally, the cross-sections and profile-lines can be stored as standard channel features in the Arc Hydro data model, and are also suitable to be used as input for hydrodynamic modeling in one, two or three- dimensions. 6.2 ANALYTICAL MODEL FOR EXTRAPOLATION OF RIVER CHANNEL BATHYMETRY RCMM (River Channel Morphology Model) is based on inter-relating the channel planform, the thalweg location and the cross-sectional form with each other. Therefore, by knowing only the channel planform, the three-dimensional form of the channel can be described. The resulting three-dimensional description is network of cross-sections and profile lines. RCMM is developed in a domain that is independent of location and scale. The concept of using the normalized domain overcomes the issues that are associated with different scales. For example, any river channel, irrespective of its location, shape or dimensions, is always straight with unit width and unit depth in the normalized domain. In this normalized domain, the thalweg location is a function of the meandering channel planform, which in turn is related to the cross- sectional form. The channel planform is characterized by using the radius of curvature, which has a logarithmic relationship with the thalweg location. The relationship between the curvature of the channel planform and the thalweg is dependent on the sinuosity of the channel, and this is demonstrated by application of RCMM to the study reaches along the Brazos River, Guadalupe River and the Sulphur River. The results from application of the non-dimensional form of radius 203 of curvature versus thalweg relationship proved that a single equation is inadequate in describing the river channel morphology. Among different analytical functions considered in this research, a combination of two beta distributions worked the best for describing the cross- sectional form of the river channels. The parameters of the beta distributions for describing the cross-sectional form are derived by using the data at Site 2 along the Brazos River. However, unlike the thalweg-radius of curvature relationship, the parameters of beta distributions are not changed while applying RCMM to the study reaches along the Guadalupe River and the Sulphur River. The beta distributions only describe the general shape of the cross-sections, and it is assumed that this general shape remains unchanged for different channels. This may be true for river channels in similar geologic and climatic settings. For example, all the river channels studied in this research are alluvial channels located in Texas. River channels in mountainous regions with hard bedrock may have a different cross-sectional form. Therefore, the parameters of the beta distributions used for describing the cross-sectional form may have to be changed while using RCMM under different geologic and climatic settings. The cross-sections described in normalized domain are re-scaled to fit real cross-sections by using the hydraulic geometry relationships, which relate the width and the depth of the channel to the flow (or the upstream watershed area). The data available from the USGS gaging stations are used to develop the hydraulic geometry relationships, and these relationships worked well for the study cases. The cross-sections are then used to generate profile lines, which are 204 described using Bezier curves. However, unlike cross-sections, the Bezier curves used to represent profile-lines are not related to any channel properties. In effect the parameters of Bezier curves are related to the cross-sections, but they do not have any physical meaning in RCMM. RCMM is developed (or calibrated) by using the data at Site 2 along the Brazos River is verified by applying it to three datasets. These three datasets are: Site 1 along the Brazos River (upstream of Site 2), the Guadalupe River, and the Sulphur River. Two measures, normalized root mean squared error (RMSE*) and coefficient of determination (R 2 ), are used to test the suitability RCMM. It is found that the R 2 estimates were significantly influenced in a negative manner by a small portion of the verifying data that are not well captured by RCMM. Therefore, the verifying data had to be separated into two groups. The first group (Group 1) contained the data with low prediction errors, and the second group (Group 2) contained the data with high prediction errors. Only Group 1 that contained more than 75 percent of the data is used to estimate the R 2 value for the verifying datasets. The spatial distribution of the prediction errors showed that the data from Group 2 (high prediction errors) lie in those areas where the thalweg described by RCMM does not match closely with the measured thalweg. Therefore, a precise description of thalweg is an important factor in the channel bathymetry described by RCMM. This may also be related to other factors such as geology, climate, land use, etc. that are not included in RCMM. The data on the channel planform are available as blue lines from the National Hydrography Dataset, and the flow data are also available from USGS. 205 RCMM thus provides a procedure to describe the river channels in three dimensions over large spatial domains. Since the bathymetry data on river channels are not available over large spatial domains, RCMM, is a useful tool for describing channel bathymetry that may be useful in large-scale regional studies. This is demonstrated by application of RCMM to generate the channel bathymetry along a 300 miles segment of Brazos River. However, it is cautioned that RCMM should be calibrated at least with some data before applying it in regional scale studies. 6.3 RECOMMENDATIONS FOR FUTURE WORK The research presented in this dissertation has two basic limitations. First, it is applicable to flow within the channel, but not on the floodplain; second, it is applicable to a single channel, not to a flow confluence, or to a braided river channel. Therefore, the first recommendation would be to enhance the procedures presented in this dissertation to overcome these limitations. Channel bathymetry is used for floodplain modeling, and it is important to develop a procedure to integrate the cross-sections and profile-lines with the rest of the terrain so as to create a seamless terrain model that also has detailed channel bathymetry. Tributaries and flow confluences are common in all river channels. Therefore, incorporation of these features is important for application of the research presented in this dissertation to a wide variety of river channels. Creating an accurate description of river channels is a complex problem. Using an analytical form to describe cross-sections and relating their parameters to channel planform is a promising approach. However, using channel planform 206 and flow data, the river channels can be described only to a limited extent. Therefore, RCMM needs to be modified to accommodate additional information. The incorporation of sinuosity to modify the relationship between the thalweg location and radius of curvature is one example. Similarly, the role of geology (bed material type, floodplain deposits), land use, flow regimes, etc. in the development of river cross-sections can also be studied to enhance RCMM. In addition, meandering river channels also have large-scale variability in terms of pools and riffles sequences, which can be identified from the channel planform. A pool generally occurs at a bend followed by a riffle. The model should be extended to account for these large-scale variable features. Finally, RCMM should also be modified to account for small-scale stochastic variability by using the concept of sequential gaussian simulation. Appendix A Tutorial on Geospatial representation of river channels Prepared by Venkatesh Merwade Center for Research in Water Resources University of Texas at Austin 1.0 Introduction This tutorial is an exercise that performs all the functions described Chapter 4, Development of a Geospatial Structure for River Channels. It is expected that the user is familiar with ArcGIS software and knows how to perform basic operations such as adding data, opening attribute tables, selecting features, and editing the data. 2.0 Data Requirement All the data necessary for this exercise are stored in a geodatabase called Channel.mdb in the Data folder available at: ?http://www.crwr.utexas.edu/gis/gishydro03/Channel/Data.zip.? The Data folder also contains a sub-folder named Worked Data, which contains all the results for this exercise. The user may use some of these results in the exercise. The channel dataset in the Channel.mdb contains four feature classes as shown below: The BathymetryPoints contains the bathymetry data with elevation as one of its attributes, and the BoundaryBolygon is the boundary for the bathymetry data. The other two feature classes, FishNetXY and Thalweg, are empty and will be populated in this exercise. 207 3.0 Getting started Please copy the Data folder on the hard drive. Also make sure that a folder named ?c:\temp? exists on the drive. Open an empty ArcMap document and save it as channel.mxd. The first step is to add the channel tool to the map document by loading the ChannelTool.dll. To load the dll click on Tools-->Customize? as shown below: In the customize window, click Add from file? button to browse the ChannelTool.dll in the Data folder. Click open and then OK to add the Channel Tool to the toolbars list. Click on Channel Tool as shown below: 208 Close the window and you should see the Channel Tool as shown below: Save the ArcMap document. The river channel toolbar has three command buttons as below: (1) Locate Thalweg: to locate the thalweg using the bathymetry data, channel banklines, and an arbitrary centerline digitized over the bathymetry data. (2) Assing (s,n) coordinates: to assign (s,n) coordinates to the bathymetry data using the thalweg as the reference line and 209 (3) Create FishNet: to transfer the FishNet from (s,n) coordinates to (x,y) coordinates. These functions are illustrated one by one in the following sections. 4.0 Locating the Thalweg There are two different approaches to use the Locate Thalweg tool. First, the user can create a bathymetry grid using the bathymetry data and the channel banklines (boundary). Second, the user can let the program create the bathymetry grid. Both approaches are explained here. 4.1 First approach (using the initial line and the bathymetry grid) Add BathymetryPoints, and BoundaryPolygon from Channel.mdb to the ArcMap document. Using the Spatial Analyst toolbar in ArcMap, create a raster grid by interpolating the ?ELEV? attribute of the bathymetry points. The bmetrygrid stored in the Workded Data folder is the result. The user can either create this grid or add the bmetrygrid from the Worked Data folder. If a new grid is created, please name it as bmetrygrid. Add the bmetrygrid to the map. Turn off the bmetrygrid layer, and add the Thalweg feature class from Channel.mdb to the map. Start the edit session, and click on the sketch tool with Create New Feature task. Make sure the target is set to Thalweg layer as shown below: 210 Digitize an initial centerline on top of the bathymetry points as shown below. While digitizing, avoid sharp angles (try to keep it as smooth as possible) and the length of individual segments should be more than 100m. The initial centerline should be digitized in the direction of flow (north to south in this case). A sample initial centerline with its segments is shown below: Starting point Ending point Please make sure that the starting and ending points are outside the bathymetry data as shown in the above figure, and the starting and ending segments are perpendicular to the channel boundary. In the channel tool, click on River Channel ---> Locate Thalweg. The program will display the following message: Since bmetrygrid represents river bathymetry, click YES to get the following interface. 211 In this approach, the first two boxes (select point layer and select field) are disabled because they are not used if the bathymetry grid is provided as the input. Select the Thalweg feature class in the select thalweg layer box, BoundaryPolygon in the select boundary box, and select the bmetrygrid in the select raster layer box. Drag the interface to any corner of the map to see the program locating the thalweg. Click OK and the program will flash some points on the map. These points are the deepest points along the normal at each vertex of the input centerline, which are then joined to locate the thalweg. After the thalweg is located, the program will display a message as shown below: Click OK to see the result as shown below: Since the study area is located along a curve, the thalweg is pushed outwards towards the right bank (looking downstream). 4.2 Second approach (using bathymetry points, channel banklines, and the initial centerline) This approach is the same as the first one except that the program creates the bathymetry grid using the input data. Remove the bmetrygrid from the map. Digitize one more intial centerline and save the edits. Select this line using the select feature button. Selecting the new initial centerline is necessary because the 212 Thalweg layer already has one feature in it and if there are more than one features in the thalweg layer, the program works with the selected feature, if selected. If the new initial line is not selected then the program considers the first feature in the Thalweg feature class as the initial centerline. Since the first feature in the thalweg feature class is the thalweg located by the first approach, it is necessary to select the initial centerline digitized in the second approach. In the channel toolbar, click on River Channel ---> Locate Thalweg. The program will display following interface: Disabled In this approach, the Select raster layer box is disabled because the program will create the raster using the ELEV attribute of the BathymetryPoints. Select BathymetryPoints for the point layer, select ELEV field from the point layer, BoundaryPolygon for the boundary, and Thalweg for the thalweg layer. Click OK. In the first approach, the program immediately displayed the flashing points to locate the thalweg, but with this approach, the program has to first interpolate the raster and then locate the thalweg. It will take a while (a minute or two) before the program starts flashing the deepest points along the river channel. After the thalweg is located, the program will display the same message box that was displayed with the first approach. Click OK to see the results. Save the edits. Change the editor task to Modify Feature and click on the modify button as shown below: Select the thalweg feature to see the vertices of the polyine as shown below: 213 Vertices Bring the cursor over one of the vertices and right click to see the Edit Sketch window. Click on Properties? as shown below: ArcMap will then display the Edit Sketch Properties of the selected feature as shown below: 214 Notice that the thalweg is a three-dimensional polyline (PolylineZM) with two extra coordinates, m (measure) and z (elevation), besides (x,y), at all its vertices. Close the Edit Sketch Properties window, save the edits, and stop the edit session. 5.0 Assigning (s,n) coordinates to the bathymetry data This section deals with assigning (s,n) coordinates to the bathymetry data. These coordinates are stored as attributes in the feature class. Before using the tool, it is necessary to add new fields to the feature class to store the (s,n) coordinates. Open the attribute table of the BathymetryPoints feature class. On the bottom-right corner of the table, click on Options--> Add Field? as shown below: 215 In the Add Field window, add a new field named sCoordinate of type Double as shown below: 216 Repeat the same process to add one more field named nCoordinate of type Double. Make sure that the names are spelled correctly otherwise the program will display a message to add the field again. Start the edit session and click on River Channel--> Assign (s,n) coordinates as shown below: The program will display the Assign (s,n) coordinates interface as shown below: Select BathymetryPoints for the input point layer and Thalweg for the centerline layer. Click OK, and the program will display the following message after storing (s,n) coordinates in the attribute table of the BathymetryPoints feature class. Click OK. Save the edits and stop the edit session. Open the attribute table of the BathymetryPoints feature class to see the coordinates. Select all the points with negative n coordinates and all the points lying to the left hand side (looking 217 downstream) of the thalweg are selected. Similarly all the points to the right hand side of the thalweg have positive n coordinates. 6.0 Creating a FishNet This section, besides the Channel Tool, makes use of the FishNet tool from ESRI. The description of the FishNet tool, and how to use it can be found at the following Internet link: http://arcobjectsonline.esri.com/ArcObjectsOnline/Samples/3D%20Analyst/3D% 20Visualization/Surface%20Fishnet/Surface%20Fishnet.htm The tool itself is located in the Data folder as FisnNetCommand.dll. Add the FishNetCommand.dll to ArcMap by following the same instructions given in the Getting started section at the beginning. The FishNet command button is shown below: The procedure for creating a fishnet that will give a network of profile lines and cross-sections is illustrated below in four steps. Step 1 - plotting the bathymetry data in (s,n,z) coordinate system Open the attribute table of the BathymetryPoints feature class. Click on Options in the bottom-right corner and then click on Export? as shown below: 218 This will save the attribute data into a dbf table. Save the attribute data as Bathymetrydata.dbf in the Data folder. ArcMap will display the following message: Click YES to add the table to the current ArcMap document. Right click on Bathymetrydata.dbf and then click on Display XY Data? as shown below: 219 In the Display XY Data window, specify sCoordinate and nCoordinate for XField and Yfield, respectively as shown below: A layer named BathymetryDataEvents will be added to the ArcMap document. Right click on this layer and click on Zoom To Layer to see the new layer as shown below: The bathymetry data are now transformed into a new coordinate system (s,n,z) where the river channel is represented straight in the direction of flow. Right click on the layer and click on Data-->Export Data? to export the data to a shapefile as shown below: 220 Save the data as StBathymetry.shp and click YES to add this new file to the map when a message shown below is displayed: Remove the BathymetryDataEvent layer and the BathymetryData.dbf table from the map. StBathymetrydata.shp stores the bathymetry data in the new (s,n) coordinate system. Save the ArcMap document. Step 2 ? creating the bathymetry grid in the (s,n,z) coordinate system Add the StBoundary.shp from the Worked Data folder or create a new shapefile and digitize the new boundary for the bathymetry data in the (s,n) coordinate system as shown below: Using the spatial analyst toolbar in ArcMap, create a raster grid for the bathymetry data in (s,n,z) coordinates by interpolating the ELEV attribute of the StBathymetry shapefile. Use a cell size of 2.5m and use the new boundary as the mask for interpolation. Name the output grid as StGrid. The resulting raster in the new coordinates is shown below: 221 Save the ArcMap document Step 3 ? creating a FishNet in (s,n,z) This step involves using the ESRI FishNet tool to create a FishNet in (s,n,z) coordinates. Click on the FishNet? command button to get the Create Fishnet From Surface interface as shown below: Select the input surface as StGrid. Since this area is about 700m long, input the approximate number of mesh lines on the longest side (along the flow) to be 50 to get an approximate spacing of 15m between the lines. Do not check the option for grouping the surface and the FishNet. Save the output feature class as FishNetSN.shp. Click OK and the tool will create a FishNet as shown below: 222 (Note: If this tool gives any error message, please register the FishNetCommand.dll in ArcScene and create the FishNet in ArcScene by adding the stgrid as the input surface) Save the ArcMap document. Step 4 ? Transferring the FishNetSN to the original coordinate system Start the edit session. Add the FishNetXY.shp shapefile from the Channel.mdb to the map. Click on River Channel-->Create FishNet as shown below: The program will display the Transfer FishNet interface as shown below: Use FishNetSN.shp for the Input FishNet(s,n) Layer and Thalweg for the Centerline Layer. Select the empty feature class, FishNetXY, for the Output FishNet(x,y) Layer. The program will transfer the FishNet in (s,n) coordinates to 223 (x,y) coordinates to get a network of lines that are parallel and transverse to the flow direction as shown below: Save the edits. Change the editor task to Modify Feature and click on the modify button. See the edit sketch properties to notice that the FishNet lines are three- dimensional lines of type PolylineMZ. This result may be useful to use at finite element mesh in hydraulic modeling studies. OK, you are done! 224 225 References Ackerman, Cameron, T., Evans, Thomas, A., and Brunner, Gary, W., 1999. HEC- GeoRAS: Linking GIS to Hydraulic Analysis Using ARC/INFO and HEC-RAS. 1999 International ESRI User Conference. ESRI, Redlands, CA. Addiscott, T., M., and Mirza, N. A., 1998. Modelling contaminant transport at catchment or regional scale. Agriculture, Ecosystems and Environment. Volume 67, No. 2-3, pp. 211-221. Anderson, David, J., 2000. GIS based hydrologic and hydraulic modeling for floodplain delineation at highway river crossings. M.S. thesis. The University of Texas at Austin. Center for Research in Water Resources online report 00-11. Austin Barney, and Wentzel Mark, 2001. Two-dimensional fish habitat modeling for assessing instream flow requirements. Proceedings of the International Symposium on Integrated Water Resources Management. University of California, Davis (LAWR), April 2000. IAHS Publ. no. 272, 2001, pp. 393?399. Azagra, Esteban, Olivera, Francisco, and Maidment, David, 1999. Floodplain Visualization using TINs. Center for Research in Water Resources online report 99-5. B?zier, Pierre, E., 1993. The First Years of CAD/CAM and the UNISURF CAD System, in Fundamental Developments of Computer-Aided Geometric Modeling, edited by Les Piegl, Academic Press, pp. 13-26. Biftu, G. F., and Gan T. Y., 2001. Semi-distributed, physically based, hydrologic modeling of the Paddle River Basin, Alberta, using remotely sensed data. Journal of Hydrology. Volume 244, No. 3-4, pp. 137-156. Blench, Thomas, 1970. Regime theory design of canals with sand beds. Journal of the Irrigation and Drainage Division, ASCE Proceedings. Volume 96, pp. 205-213. 226 Bovee, Ken, D., 1996. Perspectives on two-dimensional river habitat models: the PHABSIM experience. Proceedings of Second International Symposium on Habitat Hydraulics Ecohydraulics 2000. Quebec, Canada, June 1996, Volume B, 149-162. Burroughes, Janet, and George, Ken, 2001. Automation of interpolation by zoned inverse distance weighting for linearly distribution of soundings. GeoCoast. Volume 2, No. 1, pp. 16-35. Callander, R., A., 1978. River Meandering. Annual Review of Fluid Mechanics. Volume 10, pp. 129-58. Carter, Glenn, S., and Shankar, Ude, 1997. Creating rectangular bathymetry grids for environmental numerical modeling of gravel-bed rivers. Applied Mathematical Modeling. Volume 20, No. 11, pp. 699-708. Chang, H., H. and Hill, J., C., 1976. Computer modeling of erodible flood channels and deltas. Journal of the Hydraulics Division, ASCE Proceedings. Volume 102, pp. 1461-77. Chang, Howard, H., 1980. Geometry of gravel streams. Journal of the hydraulics division, ASCE Proceedings. Volume 106, pp. 1443-1455. Chau, K. W., and Jiang, Y. W, 2001.3D numerical model for pearl river estuary. Journal of Hydraulic Engineering. Volume 127, No. 1, pp. 72-82. Christakos, George, 1992. Random field models in earth sciences. Academic Press, San Diego, CA, 474 p. Cody, Ronald, P., and Smith, Jeffrey, Smith, K., 1997. Applied Statistics and the SAS Programming Language. Prentice Hall, NJ, pp. 445. Colla, V., Reyneri, L., M., and Sgarbi, M., 1999. Orthogonal least squares algorithm applied to the initialization of multi-layer perceptrons. Proceedings of the European Symposium on Artificial Neural Networks, Bruges, Belgium, April 1999. Crowder, D., and Diplas, P., 2000. Using two-dimensional hydrodynamic models at the scales of ecological importance. Journal of Hydrology. Volume 230, Issues 3-4, pp. 172-191. Danish Hydraulic Institute (DHI), 2000. MIKE 11 Reference Manual and User Guide, 2000. 227 Danish Hydraulic Institute (DHI), 2001. MIKE 21 Reference Manual and User Guide, 2001. Deutch, C., V., and Journel, A., G., 1992. GSLIB: geostatistical software library and user?s guide. Oxford University Press New York, 340 p. Devore, Jay, L., 1995. Probability and Statistics for Engineering and the Sciences: Fourth Edition. Duxbury Press, Pacific Grove, CA. Dierckx, Paul, 1993. Curve and Surface Fitting with Splines. Clarendon Press, New York. Dietrich, W., E., 1987. Mechanics of flow and sediment transport in river bends. In River Channels: Environment and Process edited by K. S. Richards. Institute of British Geographers Special Publication No. 18, Basil Blackwell Inc., p. 179-227. Djokic, D., Beavers, M., A., and Deshakulakarni, K., 1994. Arc/HEC2- An ArcInfo-Hec2 Interface. Proceedings of the 21 st annual conference on water policy and management. American Water Resources Association, ASCE, pp. 41-44, May 1994. Dodov, B., and Efi-Foufoula-Georgiou, 2003. Generalized hydraulic geometry based on a multiscaling formalism, submitted, Water Resources Research, November 2003. Donnell, Barbara , P., Letter, Joseph, V., McAnally, W., H., and Thomas, W., A., 2001. Users Guide to RMA2 WES Version 4.5. Eaton, Brett, C., Church, Michael, and Millar, Robert, G., 2004. Rational regime model of alluvial channel morphology and response. Earth Surface Processes and Landforms. Volume 29, pp. 511-529. Einstein, H., A., and Shen, W., A., 1964. A study of meandering in straight alluvial channels. Journal of Geophysical Research. Volume 69, pp. 5239- 5247. Eriksson, M., and Siska, P. P., 2000. Understanding anisotropy computations. Mathematical Geology. Volume 32, No. 6, pp. 683-700. 228 Fissel, David, Birch, Rick, and Jiang, Jianhua, 2002. Three-dimensional computational flow modeling and high resolution flow surveys for fisheries environmental studies on the upper Columbia River. Proceedings of Hydro Vision 2002 conference. Portland, Oregon, July 2002. Fread, D., L., and Lewis, J., M., 1988. FLDWAV: A Generalized Flood Routing Model. Proceedings of National Conference on Hydraulic Engineering. Colorado Springs, Colorado. French, J., R., and Clifford, N., J., 2000. Hydrodynamic modeling as a basis for explaining estuarine environmental dynamics: some computational and methodological issues. Hydrological Processes. Volume 14, Issues 11-12, pp. 2089-2108. Froehlich, D., C., 1992. Finite element surface-water modeling system: two- dimensional flow in a horizontal plane. U. S. Department of Transportation Report No. FHWA-RD-92-057. Fukoka, S., and Sayre, W., W., 1973. Longitudinal dispersion in sinuous channels. Journal of Hydraulics Division, ASCE Proceedings. Volume 99, pp. 195-217. Ghanem, Ashraf, Steffler, Peter, Hicks, Faye, and Katapodis, Chris, 1996. Two- dimensional hydraulic simulation of physical habitat conditions in flowing streams. Regulated Rivers: Research & Management. Volume 12, No. 2-3, pp. 185-200. Gilvear, D., J., Bryant, R., and Hardy, T., 1999. Remote sensing of channel morphology and instream fluvial processes. Progress in Environmental Science. Volume 1, No. 3, pp. 257-284. Goff, John, A., and Nordfjord, Sylvia, 2004. Interpolation of Fluvial Morphology Using channel-oriented coordinate transformation: A case study from the New Jersey Shelf. Mathematical Geology (in press). Gregory, K. J., Davis, R. J., and Downs, P. W., 1992. Identification of river channel change due to urbanization. Applied Geography. Volume 12, No. 4, pp. 299-318. Harman, W. A., Jennings, G. D, Patterson, J. M., Clinton, D. R, Slate, L. O., Jessup, A. G., Everhart, J. R., and Smith, R. E., 1999. Bankfull Hydraulic Geometry Relationships for North Carolina Streams. AWRA Wildland 229 Hydrology Symposium Proceedings. Edited By: D.S. Olsen and J.P. Potyondy. AWRA Summer Symposium. Bozeman, MT. Heinzer, Tom, Sehbat, Michael, Feinberg, Bruce, and Kerper, Dale, 2000. The Use of GIS to Manage LIDAR Elevation Data and Facilitate Integration with the MIKE21 2-D Hydraulic Model in a Flood Inundation Decision Support System. Proceedings of the 2000 ESRI users international conference. San Diego, CA, June 2000. Hey, R., D., 1978. Determinant hydraulic geometry of river channels. Journal of Hydraulics Division, ASCE Proceedings. Volume 104, pp. 869-885. Hodges, B., R., and Imberger, J., 2001. Simple curvilinear method for numerical methods of open channels. Journal of Hydraulic Engineering. Volume 127, No. 11, pp. 949-958. Holley, E., R., and Jirka, G., H., 1986. Mixing in rivers. Technical Report E-86- 11, U.S. Army Engineer Waterways Experiment Station, Vicksburg, Mississippi. Isaaks, Edward, H., and Srivastava, Mohan, R., 1989. Applied Geostatistics. Oxford University Press, New York. James, L. A., 1996. Polynomial and Power Functions for Glacial Valley Cross- Section Morphology. Earth Surface Processes and Landforms. Volume 21, No., pp. 413-432. Jennings, Aaron, A., 2003. Modeling sedimentation and scour in small urban lakes. Environmental Modelling & Software. Volume 18, Issue 3, pp. 281- 291. Jenson, Jerry, L., Lake, Larry, W., Corbett, Patrick W. M., and Goggin, David, J., 1997. Statistics for Petroleum Engineers and Geoscientists. Prentice Hall PTR, New Jersey. Johannesson, H., and Parker G., 1989. Linear theory of river mechanics. In River Meandering edited by Ikeda Syunsuke and Parker Gary. Water Resources Monograph, Washington D.C., Volume 12, pp. 181-213. Johnston, Kevin, Ver Hoef, Jay M., Krivoruchko, Konstantin, and Lucas, Neil, 2001. Using ArcGIS Geostatistical Analyst. ESRI Press, Redlands, CA. 230 Jowett I. G., 1997. Instream flow methods: A comparison of approaches. Regulated Rivers: Research & Management. Voume 13, No. 2, pp. 115- 127. Kanevski, M., R. Parkin, A. Pozdnukhov, V. Timonin, M. Maignan, B. Yatsalo, and S. Canu (2002). Environmental Data Mining and Modelling Based on Machine Learning Algorithms and Geostatistics. Proceedings of the First Biennial Meeting of the International Modelling and Software Society, IEMSS?2002, June 2002, Lugano, Switzerland. A. E. Rizzoli & A. J. Jakeman (Eds), Vol. 3., pp. 414-419. Karim, Khalid, Gubbels, Maureen, E., and Goulter, Ian, C., 1995. Review of determination of instream flow requirements with special application to Australia. Water Resources Bulletin. Volume 31, No. 6, pp. 1063-1077. Knighton, A. D., 1981. Asymmetry of river channel cross-sections: Part I. Quantitative indices. Earth Surface Processes and Landforms. Volume 6, pp. 581-588. Knighton, A. D., 1982. Asymmetry of river channel cross-sections: Part II. Mode of development and local variation. Earth Surface Processes and Landforms. Volume 7, pp. 117-131. Knighton, A. D., 1984. Indices of flow asymmetry in natural streams: definition and performance. Journal of Hydrology. Volume 73, pp. 1-19. Knighton, A. D., 1998. Fluvial Forms and Processes: A new perspective. Oxford University Press Inc., New York. Knighton, A.D., 1985. Channel form adjustment in supraglacial streams, Austre Okstindbreen, Norway. Arctic and Alpine Research. Volume 17, pp. 451- 466. Lacey, G., 1930. Stable channels in alluvium. Minutes Proc., Institute of Civil Engineers, London. Volume 229, pp. 259-292. Leclerc, M., Bobee, A., Boudreault A., Shooner, G., and Corfa, G., 1991. Instream flow incremental methodology and 2-D hydrodynamic modeling: efficient tools to determine guaranteed minimum flow for biological purposes. In: Computer methods in water resources-II, proceedings of the second international conference. Marrakesh, Morocco, 20-22 February 1991, pp. 289-300. 231 Leclerc, Michel, and Lafeur, Julie, 1997. The fish habitat modeling with two- dimensional hydraulic tools: a worthwhile approach for setting minimum flow requirements? Presented at the Instream and Environmental Flow Symposium. Houston, TX, December 1997. Leopold, L., B., and Maddock, T., 1953. The hydraulic geometry of stream channels and some physiographic implications. U.S. Geological Survey Professional Paper, 252. Leopold, L.B., Wolman, M. G., and Miller, J. P., 1964. Fluvial Processes in Geomorphology. W. H. Freeman and Company, San Francisco. Lisle, Thomas, E., 1987. Using residual depths to monitor pool depths independently of discharge. Research Note PSW-394, Pacific Southwest Forest and Range Experiment Station, Forest Service, U.S. Department of Agriculture, Berkeley, CA. Maidment, D., 2002. Arc Hydro-GIS for water resources. ESRI Press, Redlands, California. McCoy, Jill, and Johnston, Kevin, 2001. Using ArcGIS Spatial Analyst. ESRI Press, Redlands, CA. Mertes Leal A. K., 2002. Remote sensing of riverine landscapes. Freshwater Biology. Volume 47, No. 6, pp. 799-816. Milhous, R. T., M. A. Updike, and D. M. Schneider. 1989. Physical Habitat Simulation System Reference Manual - Version II. Instream Flow Information Paper No. 26. U.S. Fish and Wildlife Service Biological Report 89(16). Miller, S. N., Guertin, D. P., and Goodrich, D. C., 1996. Investigating Stream Channel Morphology Using a Geographic Information System. Proceedings of the 1996 ESRI International User Conference. Palm Springs, CA. Mitas, L., and Mitasova, H., 1988. General variational approach to the interpolation problem. Computers & Mathematics with Applications. Volume 16, Issue 12, pp. 983-992. Mitas, L., and Mitasova, H., 1999. Spatial Interpolation. In: P.Longley, M.F. Goodchild, D.J. Maguire, D.W.Rhind (Eds.). Geographical Information 232 Systems: Principles, Techniques, Management and Applications. GeoInformation International, Wiley, 481-492. Moody, J. A., and Troutman, B. M., 2002. Characterization of the Spatial Variability of Channel Morphology. Earth Surface Processes and Landforms. Volume 27, No. 12, pp. 1251-1266. Nelson, J. M., and Dungan J. S., 1989. In River Meandering edited by Ikeda Syunsuke and Parker Gary. Water Resources Monograph, Washington D.C., Volume 12, pp. 69-102. Nelson, Johathan, M., and Smith, Dungan, J., 1989. Flow in meandering channels with natural topography. In River Meandering edited by Ikeda Syunsuke and Parker Gary. Water Resources Monograph, Washington D.C., Volume 12, pp. 69-102. Oetter, Doug, R., Ashkenas, Linda, R., Gregory, Stanley, V., and Minear, Paula, J., 2004. GIS Methodology for Characterizing Historical Conditions of the Willamette River Flood Plain, Oregon. Transactions in GIS. Volume 8, No. 3, pp. 367-383. Parasiewicz, P., and Dunbar M., J., 2001. Physical habitat modelling for fish: a developing approach. Large Rivers. Volume 12. No. 2-4, pp. 239-268. Park, C., C., 1977. World-wide variations in hydraulic geometry exponents of stream channel: An analysis and some observations. Journal of Hydrology. Volume 33, pp.133-146. Parker G., and Johannesson, H.,1989. Observations on several recent theories of resonance and overdeepening in meandering channels. In River Meandering edited by Ikeda Syunsuke and Parker Gary. Water Resources Monograph, Washington D.C., Volume 12, pp. 379-416. Pebesma, Edzer, J., and Wesseling, Cees, G., 1998. Gstat: a program for geostatistical modelling, prediction and simulation. Computers & Geosciences. Volume 24(1), pp. 17-31. Rathburn, Sara, and Wohl, Ellen, 2003. Predicting fine sediment dynamics along a pool-riffle mountain channel. Geomorphology. Volume 55, Issues 1-4, pp. 111-124. Renschler, C. S., and Harbor, J., 2002. Soil erosion assessment tools from point to regional scales ? the role of geomorphologists in land management 233 research and implementation. Geomorphology. Volume 47, No. 2-4, pp. 189-209. Richardson, Barbara, A., 1986. Evaluation of Instream Flow Methodologies for Fresh Water Fish in New South Wales. In: Stream Protection, The Management of Rivers for Instream Uses, Ian C. Campbell (Editor). Water Studies Center, Chisholm Institute of Technology, East Caulfield, pp. 143-167. Rosgen, David, L., 1994. A classification of natural rivers. CATENA. Volume 22, pp. 169-199. Runge, Morten, and Olesen, Kim Wium, 2003. Combined one- and two- dimensional flood modeling. Fourth Iranian Hydraulic Conference, 21-23 October, Shiraz, Iran. Savenije, Hubert, H.G., 2003. The width of a bankfull channel; Lacey's formula explained. Journal of Hydrology. Volume 276, Issues 1-4, pp. 176-183. Schultz, G. A., 1993. Hydrological Modeling Based on Remote Sensing Information. Advances in Space Research. Volume 13, No. 5, pp. 149- 166. Schumm, S.A., 1960. The shape of alluvial channels in relation to sediment type. U.S. Geological Survey Professional Paper, 352-B: 17-30 Schumm, S.A., Mosley, M.P., and Weaver, W.E., 1987. Experimental Fluvial Geomorphology, John Wiley, New York. Seminara, G., and Tubino, M., 1989. Alternate bars and meandering: free, forced and mixed interactions. In River Meandering edited by Ikeda Syunsuke and Parker Gary. Water Resources Monograph, Washington D.C., Volume 12, pp. 267-320. Shen, H., W., and Komura, S., 1968. Meandering tendencies in straight alluvial channels. Journal of Hydraulics Division, ASCE Proceedings. Volume 94, pp. 997 - 1016. Simons, D., B., and Albertson, M., L., 1960. Uniform water conveyance channels in alluvial material. Journal of the Hydraulics Division, ASCE Proceedings. Volume 86, pp. 33-71. 234 Sinha S. K., Sotiropoulos F., and Odgaard A. J., 1998. Three-dimensional numerical model for flow through natural rivers. Journal of Hydraulic Engineering. Volume 124, No. 1, pp. 13-24. Snead, Daniel, B., 2000. Development and application of unsteady flood models using geographic information systems. Departmental Report. The University of Texas at Austin. Center for Research in Water Resources online report 00-12. Stalnaker Clair, Lamb Berton L., Henriksen Jim, Bovee Ken, and Bartholow John, 1995. The instream flow incremental methodology-A primer for IFIM. National Biological Service, U. S. Department of the Interior. Biological Report 29. Stein, J. L., Stein, J. A., and Nix, H. A., 2002. Spatial analysis of anthropogenic river disturbance at regional and continental scales: identifying the wild rivers of Australia. Landscape and Urban Planning. Volume 60, No. 1, pp. 1-25. Tarbet, Karl, and Hardy, Thomas, B.,1996. Evaluation of one-dimensional and two-dimensional hydraulic modeling in a natural river and implications in instream flow assessment methods. Proceedings of Second International Symposium on Habitat Hydraulics Ecohydraulics 2000. Quebec, Canada, June 1996, Volume B, 395-406. Tate, E. C., Maidment, D. R., Olivera, F., and Anderson D. J., 2002. Creraing a terrain model for floodplain mapping. Journal of Hydrologic Engineering. Volume 7, No 2, pp. 100-108. Tennant, D.L. (1976), ?Instream Flow Regiments for Fish, Wildlife, Recreation, and Related Environmental Resources,? in Orsborn, J.F. and C.H. Allman (Editors). Proceedings of the Symposium and Specialty Conference on Instream Flow Needs II. American Fisheries Society, Bethesda, MA, pages 359-373. Texas Instream Flow Studies: Programmatic Work Plan, 2002. Document prepared by Texas Parks and Wildlife Department, Texas Commission on Environmental Quality, and Texas Water Development Board. 235 Tomczak, Maciej, 1998. Spatial Interpolation and its Uncertainty Using Automated Anisotropic Inverse Distance Weighting (IDW) - Cross- Validation/Jackknife Approach. Journal of Geographic Information and Decision Analysis. Volume 2, No. 2, pp. 18-30. U. S. Army Corps of Engineers (USACE), 2002a. HEC-GeoRAS: An extension for support of HEC-RAS using ArcView. Users Manual. Version 3.1, October 2002, Hydrologic Engineering Center, Davis, California. U. S. Army Corps of Engineers (USACE), 2002b. HEC-RAS: River Analysis System. Hydraulic Reference Manual, Version 3.1, November, Hydrologic Engineering Center, Davis, California. Van Winkle., W., Jager, H. I., Railsback, S. F., Holcomb, B. D., Studley, T. K., and Baldrige, J. E., 1998. Individual-based model of sympatric populations of brown and rainbow trout for instream flow assessment: model description and calibration. Ecological Modeling. Volume 110, No. 2, pp. 175-207. Waddle, T., Bovee, K., Bowen, Z., 1997. Two-dimensional habitat modeling in the Yellowstone / Upper Missouri River System. North American Lake Management Society (NALMS) meeting in Houston, TX (December 2-4 1997). Wadzuk, B., and Hodges, B. R., 2001. Model bathymetry for sinuous, dendritic reservoirs. Proceedings of the 6th International Workshop on Physical Processes in Natural Waters, University of Girona, Catalonia, Spain, June 2001. Waltuch, M., Cameron, E., Lafframboise, A., and Zeiler, M., 2001. In Chapter 1, Introducing ArcObjects. Exploring ArcObjects Vol.1-Applications and Cartography edited by Michael Zeiler. ESRI Press, Redlands, California, pp. 1-3. Wentzel, Mark, W., 2001. Two-dimensional Fish Habitat Modeling for Instream Flow Assesment in Texas. PhD Dissertation, Department of Civil Engineering, New Mexico State University, Las Cruces, New Mexico. Werner, M. G. F., 2001. Impact of grid size in GIS based flood extent mapping using a 1D flow model. Physics and Chemistry of the Earth, Part B: Hydrology, Oceans and Atmosphere. Volume 26, Issues 7-8, pp. 517-522 236 Werritty, Alan, and Leys, Katherine, F., 1999. River channel planform change: software for historical analysis. Geomorphology. Volume 29, Issues 1-2 , pp. 107-120. Williams, G. P., 1978. Bankfull Discharge of Rivers. Water Resources Research. Volume 14, No. 6, pp. 1141-1154. Winterbottom, S. J., 2000. Medium and short-term channel planform changes on the Rivers Tay and Tummel, Scotland. Geomorphology. Volume 34, No. 3-4, pp. 127-132. Yalin, M., S., 1992. River Mechanics. Pergamon Press, Oxford. Yang, Qing-wen, Chen, Philip, and Tang, Zesheng, 1994. Efficient algorithm for the reconstruction of 3D objects from orthographic projections. Computer- Aided Design. Volume 26, Issue 9, pp. 699-717. Yang, X., Damen, M.C. J., and Van Zuidam R. A., 1999. Satellite remote sensing and GIS for the analysis of channel migration changes in the active Yellow River Delta, China. International Journal of Applied Earth Observation and Geoinformation. Volume 1, No. 2, pp. 146-157. Ye, J., and McCorquodale, J. A., 1997. Depth-averaged hydrodynamic model in curvilinear collocated grid. Journal of Hydraulic Engineering. Volume 123, No. 5, pp. 380-388. Zeiler, Michael, 1999. Modeling Our World: The ESRI Guide to Geodatabase Design. ESRI Press, Redlands, California. 237 Vita Venkatesh Mohan Merwade was born in Kolhapur, India on August 26, 1975, the son of Mohan Ramchandra Merwade and Ratna Mohan Merwade. After completing his schooling at Premier English High School and pre-university courses at Vivekanand College, he entered Kolhapur Institute of Technology?s College of Engineering in August, 1993. He received the degree of Bachelors of Engineering in Environmental Engineering from Shivaji University in May, 1997. Upon graduation, he worked as a project engineer for Montgomery Watson India from July, 1997 to August, 1999. In September, 1999, after winning an Irish Government Fellowship, he left India to pursue a post-graduate degree in Hydrology from The National University of Ireland, Galway, Ireland. In December, 2000, he graduated with Master of Science in Engineering Hydrology from The National University of Ireland. He entered The Graduate School at The University of Texas at Austin in January, 2001. Permanent address: 1735 Rutland Drive, # 225, Austin, TX 78758 This dissertation was typed by the author.