i
CRWR Online Report 0408
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 787124497
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 crosssections and
profilelines ...............................................................................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 Onedimensional Hydrodynamic Modeling ..............................26
2.3.2 Twodimensional Hydrodynamic Modeling..............................30
2.3.3 Threedimensional 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 EchoSounder Data .................................................60
3.6 Data Description ................................................................................61
Chapter 4 Development of a Geospatial Structure for River Channels ................65
4.1 Introduction........................................................................................65
4.2 Threedimensional 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 Crosssectional
Form 143
5.3.4 Rescaling the normalized crosssections .................................161
5.3.5 Creating profilelines using crosssections ..............................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 crosssections...................................................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 Threedimensional Description of River Channels Using Cross
sections and Profilelines .................................................................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
crosssection.....................................................................................156
Table 5.2: Summary of regression analysis for fitting a combination of
beta distributions to the crosssection 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 2year 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, crosssections and a finiteelement
mesh.....................................................................................................8
Figure 1.4 A wireframe model representation of sphere ...............................9
Figure 1.5 Threedimensional description of a river channel using cross
sections and profilelines....................................................................11
Figure 1.6 Threedimensional description of a river channel using cross
sections, profilelines and hydrovolume ...........................................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) Crosssections and
profilelines........................................................................................21
Figure 2.4 Basic components of instream flow studies ...............................23
Figure 2.5 River system schematic in HECRAS........................................28
Figure 2.6 A typical HECRAS crosssection .............................................29
Figure 2.7 Finite element mesh for twodimensional 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 crosssection with a steep slope; b) Crosssection
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) crosssectional forms for different thalweg locations...129
Figure 5.2 Flowchart for the development of RCMM...............................132
Figure 5.3 A typical crosssection with (s,n,z) coordinates. P(n
i
,z
i
) is any
point in a crosssection; 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 crosssection. (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 crosssection ................................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 crosssection............146
Figure 5.13 Thirddegree polynomial fitted to the channel crosssection..147
Figure 5.14 Interpolating and smoothing splines ......................................148
Figure 5.15 Smoothing spline fitted to the crosssection...........................149
Figure 5.16 Beta distribution.....................................................................150
Figure 5.17 Beta distribution fitted to the crosssection ............................151
Figure 5.18 Combination of two beta distributions ...................................152
Figure 5.19 Combination of Beta distributions fitted to the channel
crosssection.....................................................................................153
Figure 5.20 Gamma density functions for different (a,b) pairs .................154
Figure 5.21 Gamma distribution fitted to the crosssection.......................155
Figure 5.22 Combination of gamma distributions fitted to the channel
crosssection.....................................................................................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 2year return period flow (Q
2
) and
drainage area (A) for the Brazos River.............................................166
Figure 5.27 Generation of profilelines using crosssections.....................168
xiii
Figure 5.28 Thalweg prediction using the radius of curvaturethalweg
location relationship (equations 5.5 and 5.6). Circled areas show
locations where the thalweg prediction is imprecise. Sections AA
and BB are plotted in Figures 5.30 and 31 ......................................170
Figure 5.29 Thalweg prediction using the nondimensional radius of
curvaturethalweg location relationship (equations 5.7 and 5.8) ......171
Figure 5.30 Crosssection described by RCMM at section AA shown in
Figure 5.28.......................................................................................172
Figure 5.31 Crosssection described by RCMM at section BB shown in
Figure 5.28.......................................................................................173
Figure 5.32 Threedimensional description of a river channel in the form
of crosssections 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 crosssections 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, sdirection)
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 threedimensional structure of river channels. The research
presented in this dissertation provides a way of representing the threedimensional
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 threedimensional structure of river channels in the
form of a mesh of crosssections and profilelines. The crosssections transverse
the flow, and the profilelines run along the flow. A network of crosssections and
profilelines thus provides a threedimensional 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 smallscale (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 planview), the thalweg location, and the
crosssections. In addition, River Channel Morphology Model characterizes the
channel planform and crosssections by analytical functions. Therefore, by
knowing only the planform of the channel which is available for large spatial
domains, the crosssections and profilelines 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 onedimensional
model such as HECRAS requires input bathymetry data in the form of cross
sections. RMA2, on the other hand, is a twodimensional 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 threedimensional representation of river
bathymetry is a primary motivation for this research.
7
Bathymetry Points Crosssections Finite element mesh
Figure 1.3 Bathymetry data as points, crosssections and a finiteelement 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 threedimensional 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
threedimensional lines to represent surfaces. Therefore, an efficient way of
storing and representing channel bathymetry is by using threedimensional lines,
crosssections and profilelines.
9
10
Channel planform refers to the twodimensional (x,y) planview of a river
channel in space. Each river crosssection 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 crosssections 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 crosssections and the channel planform are interrelated using analytical
functions, it may be possible to describe the channel crosssections 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, crosssections 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 crosssection and profilelines. A network of cross
sections and profilelines forms a threedimensional mesh to represent the
channel bathymetry (Figure 1.5).
Figure 1.5 Threedimensional description of a river channel using crosssections
and profilelines
2. Develop a model in which the channel crosssections and channel planform
are represented using analytical functions. Such a model will allow inter
relating crosssections with channel planform. Therefore, by knowing only the
channel planform, crosssections 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 crosssections and profilelines?
11
12
2. Given the standardized representation of channel bathymetry, how can this
information be used to predict the threedimensional 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 crosssectional 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 crosssections. Therefore, the change in
channel planform and the crosssectional shape are interrelated. In addition to
change in the shape, the size (crosssectional area) of crosssections 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 crosssections is related to the channel planform.
ii. The size of river crosssections 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 crosssections.
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 crosssections 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 crosssections
and profilelines. The channel bathymetry in the form of crosssections and
profilelines 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, crosssections and profilelines can serve as input to one, two,
or threedimensional 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 crosssections and profilelines.
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 threedimensions.
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 crosssections and
profilelines (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 (HydroVolume), and this concept of
HydroVolume opens doors for developing hydrodynamic models and water
quality models (mass balance models) for river channels using GIS.
Figure 1.6 Threedimensional description of a river channel using crosssections,
profilelines and hydrovolume
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 threedimensional description of river bathymetry using cross
sections and profilelines. The fifth chapter describes the River Channel
Morphology Model for creating a threedimensional description of river channel
bathymetry by interrelating the crosssectional 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, timeseries records of discharge
measurements stored in a table is an object class.
? Geodatabase: A geodatabase is the toplevel 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 twodimensional (x,y)
planview of a river channel in space (Figure 2.3a).
? Channel crosssection: A river channel crosssection is a ground profile
transverse to the flow direction (Figure 2.3b).
? Channel profileline: A river channel profileline 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 straightline valley length is
the sinuosity of the meandering bend.
Valley length
ML
x
y
A
Flow direction
a)
b)
Crosssection
Profile lines
Thalweg
Threedimensional view at A
Figure 2.3 (a) Channel planform definition; (b) Crosssections and profilelines
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 nonconsumptive 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 flowrates. 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 crosssection surveys and stagedischarge
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 onedimensional
hydrodynamic model, several habitat rating methods employing twodimensional
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 twodimensional 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 onedimensional
habitat model to a twodimensional 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 twodimensional models. While others,
26
Leclerc and Lafleur (1997), for example, used the bathymetry data collected for a
twodimensional model to compare the two approaches. The conclusion, however,
was the same in all the studies: the twodimensional habitat model performed
better than the onedimensional model. Advantages provided by twodimensional
models include fewer velocity measurements, easier calibration, and flexibility
with respect to modeling under different scenarios.
Besides advantages, twodimensional models were also reported to be
very sensitive to channel bathymetry (Bovee, 1996; Tarbet and Hardy, 1996). The
accuracy of any twodimensional 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) twodimensional, and 3) threedimensional. Each of these groups
is briefly discussed in the following sections.
2.3.1 Onedimensional Hydrodynamic Modeling
Onedimensional 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 onedimensional models solve onedimensional
energy equations to compute water surface elevations at each crosssection 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 onedimensional models simulate unsteady
state flow conditions in river channels and solve crosssectional averaged Saint
Venant equations to route the flow hydrograph and compute water surface
elevation at each crosssection. For example, HECRAS (USACE, 2002b), a
commonly used onedimensional hydrodynamic model, has the capability to
perform both steady and unsteady state simulations.
A typical onedimensional model, HECRAS in this case, requires
geometric and flow data. Basic geometric data consist of the river system
schematic and the crosssection 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
Crosssection
Crosssection
Crosssection
Mcs
M
s
Reach 3
Figure 2.5 River system schematic in HECRAS
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 crosssection for HECRAS. Each reach on the river system
schematic has a unique identifier. Each crosssection on the river system
schematic is associated with a reach identifier, and a station identifier (M
s
), which
specifies distance along the reach. The crossreach distance along the cross
section is specified by station numbers (M
cs
). Looking downstream, the left end of
a crosssection has the lowest station number and the right end has the highest
station number.
28
M
s
M
cs
Figure 2.6 A typical HECRAS crosssection
The flow data for HECRAS 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 dischargerating curve.
Besides geometry and flow data, additional information such as substrate type to
calculate Manning?s n is also required. Besides HECRAS, MIKE 11 (DHI, 2000)
from the Danish Hydraulic Institute and FLDWAV (Fread and Lewis, 1988) from
the National Weather Service are also used for onedimensional modeling.
29
30
Onedimensional models have traditionally been used for flood plain
mapping, and are still commonly used for this purpose. The model setup is easy
and the computations are fast. However, they are limited by their ability to model
twodimensional characteristics such as channel meander, velocity currents,
flooding in flat urban environments, lateral distribution of flows, etc. Some
limitations of onedimensional models are overcome by using twodimensional
models.
2.3.2 Twodimensional Hydrodynamic Modeling
Onedimensional 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. Twodimensional
hydrodynamic models solve depthaveraged 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 twodimensional 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 twodimensional 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 twodimensional models use either the finitedifference method or the
finiteelement method to solve the energy and momentum equations. The choice
of solution method depends on several factors such as the computational speed,
model setup, 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 twodimensional 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 finitedifference model. Typical data
requirements for twodimensional 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 onedimensional models. The collection
of geometry data for twodimensional models, however, demands more resources
compared to onedimensional models.
Finite element nodes
Finite element mesh
Bathymetry points
Figure 2.7 Finite element mesh for twodimensional hydrodynamic modeling
The most commonly used boundary conditions for twodimensional
models are downstream water surface elevation and upstream flowrate.
Calibration data for twodimensional 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 twodimensional models consist of eddy viscosity values and roughness
coefficients for various substrate types.
Twodimensional models have both advantages and disadvantages
compared to onedimensional models. Advantages include twodimensional
distribution of flow and velocity, and the ability to simulate complicated flow
patterns. The main disadvantage of twodimensional models is that they are not as
powerful as onedimensional models for simulating flow across control structures
such as weirs, pumps, tidal gates, etc. This disadvantage of twodimensional
models can be overcome by combining them with onedimensional models to use
the best of both models (Runge and Oslen, 2003). The crucial factor while using
twodimensional 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 twodimensional
models have been made both by scientists and engineers. Some of these attempts
include refining the mesh resolution to account for smallscale 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 twodimensional hydrodynamic
modeling.
2.3.3 Threedimensional Hydrodynamic Modeling
Threedimensional models are similar to twodimensional models, except
that the governing equations are not depthaveraged. 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 threedimensional models
include their ability to model the vertical distribution of flow and velocity, which
include stratification, diffusion and dispersion processes. Threedimensional
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 threedimensional models in the future.
The data requirements for threedimensional 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 threedimensional 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 HECRAS (then
HEC2), and then import the output results from HECRAS to GIS for
visualization. This effort was carried forward by Ackerman et. al., (1999), which
led to the development of HECGeoRAS (USACE, 2002a). HECGeoRAS is a
GIS extension that allows GIS users to create and export HECRAS geometry
files such as channel centerline and crosssections using a digital terrain model.
The output from HECRAS is also available in GIS format for visualization.
Development of HECGeoRAS 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
HECRAS 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 threedimensional 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 preprocess raw data
for two and threedimensional models. A typical preprocessing step in GIS to
prepare the input data for a two or threedimensional 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 preprocessing
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 threedimensional 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 threedimensional description of river channels in the form
of crosssections 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 geoprocessing or editing tools available in GIS can be used
to preprocess 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 crosssections and
profilelines 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 flowlength, L
pq
, which is greater than
d
pq
. Using d
pq
for calculating traveltime, 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 floworiented 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 floworiented 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 singlebeam 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 multibeam 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 singlebeam depth
sounder. The minimum cost associated with singlebeam depth sounder and
multibeam 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 mainstem 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 crosssections 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 satellitebased 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 crosssection 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 autocorrelated, and similar flow perturbations are approximately
separated by equal distance ( ?2 times width), leading to alternate regions of high
speed and lowspeed 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 maximumvelocity 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 crossstream variations in
the boundary shear stress and velocity fields (Dietrich, 1987), which in turn cause
and propagate the asymmetries in crosssectional 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 crosssections 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 crosssectional 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 crosssectional 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 crosssectional 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 (atastation hydraulic geometry) or can be developed across several
crosssections 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 crosssectional 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 threedimensional models instead of onedimensional
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 preprocessed 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 threedimensional description of river
channel bathymetry in GIS. Therefore, a procedure that will take into account the
orientation of channel morphology to create a standardized threedimensional
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, crosssections that are input to onedimensional
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 preprocessing. Arc Hydro already
has a channel component to store the data in vector form using crosssections 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
crosssectional 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 rivermiles 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 boatmounted 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 setup 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 setup.
3.5.1 Multibeam EchoSounder 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 echosounder. Multibeam ecosounder 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
echosounder is performed at several different locations on the channel bottom at
once (Geyer, 1992). The multibeam echosounder 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 highresolution
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 twodimensional 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 saltwater 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 crosssection with a steep slope; b) Crosssection measured
using a single beam depth sounder
Figure 3.9(a) shows a crosssection 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 crosssections.
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 crosssections and profilelines. 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 THREEDIMENSIONAL 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 channelfitted
(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
threedimensional 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 ncoordinates, and all the points
lying to the right hand side of the centerline have positive ncoordinates. For
example, in Figure 4.3, point P has a negative ncoordinate, and point Q has a
positive ncoordinate. 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 flowindependent 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 userdefined 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 ncoordinates (Nelson and Dungan, 1989; Johannesson and Parker,
1989). A thalweg or a centerline or a bank can introduce singularities (intersecting
of crosssections) 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 threedimensional lines (PolylineZM) that are
parallel to the x and yaxes. 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 threedimensional surface.
4.5 METHODOLOGY
The desired result is a threedimensional vector dataset for storing river
channel bathymetry in the form of crosssections and profilelines. The surface
FishNet tool available in ArcGIS provides a network of threedimensional 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 crosssections and the lines parallel to the yaxis are not profilelines.
However, the concept of (s,n) coordinate system can be used to convert a regular
FishNet into a network of crosssections 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 saxis and transverse to the naxis.
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 naxes,
which are actually parallel and transverse to the flow. Therefore, a FishNet in
(s,n) coordinates gives a network of crosssections (lines parallel to naxis) and
profilelines (lines parallel to saxis). 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 profilelines.
The basic methodology for creating a threedimensional 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 profilelines.
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 0135 degrees or between 225360 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 crosssections covering the entire width of the channel.
Each normal line is a threedimensional 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 threedimensional 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 scoordinate at each vertex is equal to the
measure value at that point. The zcoordinate 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 (ncoordinate) 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 ncoordinate of P, while the scoordinate 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 scoordinate of P is equal to the scoordinate of point
B, and the ncoordinate 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
bedslope along the flow (?x in Figure 4.16) is smaller than the bedslope across
the flow (?y in Figure 4.16).
A A
B
B
x
z
y
Plan view of a river channel
Section AA Section BB
?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 naxis to the variance in the bathymetry along saxis.
)(
)(
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
m1
x
m1
+
a
m2
x
m2
+?+ a
1
x + a
0
, where x is the independent variable, m is a nonnegative
integer (degree of the polynomial function), and a
m1
, a
m2
,?, 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 twodimensional Euclidian space (area of interest), and is the
square of the secondorder 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
semivariogram model, which provides a measure of variance as a function of
distance between the observation points. The semivariance (?
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 leastsquares 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
Semivariogram
model
Empirical
Semivariogram
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 yaxis is called as the sill. The value along the y
axis at a point where the semivariogram intercepts with the yaxis is called as
the nugget. At zero separation distance (d
ij
= 0), the semivariance 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 microscales 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 twodimensional
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.
Semivariogram envelope
Semivariogram 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 subdatasets: 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 worstcase
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 worstcase scenario (0.52m) is better than even the bestcase
scenario of the regular inverse distance weighting approach (0.65m). The RMSE
for the bestcase scenario (20 points) is 0.49m, which is only (0.520.49) 0.03m
lower than the worstcase 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 bestcase scenario and
the worstcase scenario is only (0.530.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.
Semivariogram envelope
Semivariogram 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 thirddegree 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 yaxes (Figure 4.4). Similarly, the FishNet created in the
(s,n) coordinate system has the lines parallel to the s and naxes. 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 crosssections and profilelines.
4.5.6 Transforming the FishNet
To get a threedimensional description of river bathymetry using cross
sections and profilelines in the Cartesian coordinate system, the FishNet created
in (s,n) coordinates is transformed back to the (x,y) coordinate system. Each
crosssection and profileline 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 crosssections while the lines parallel to
the centerline are profilelines (Figure 4.42).
2D Hydraulic FishNet 3D Hydraulic FishNet
Figure 4.42 Hydraulic FishNet
In addition, when the crosssections and profilelines 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 threedimensional
models, the result can be used as it is to define a finite element mesh. For one
126
127
dimensional models the crosssections can be separated, and then used as an input
to define the onedimensional 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 crosssections and profilelines?) 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 threedimensional channeloriented mesh
comprising of crosssections and profilelines. The procedure uses the concept of
(s,n) coordinate system to transform a regular ArcGIS FishNet into a channel
oriented network of crosssections and profilelines. 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 threedimensional description
of a river channel. The procedure uses a set of bathymetry points as an input, and
the end result is a network of crosssections and profilelines. 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 threedimensional
features of the FishNet (crosssections and profilelines) 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 threedimensional 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 crosssectional 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)
crosssectional forms for different thalweg locations
129
130
As mentioned in section 2.7, in meandering river channels, the crosssections 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 crosssectional shapes by moving the thalweg towards the outer
banks. Therefore, the development of crosssectional asymmetry in river channels
can be related to the change in the channel planform (Leopold et. al., 1964;
Knighton, 1998). RCMM assumes the crosssectional 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 crosssections have asymmetric or
symmetric forms (Figure 5.1b). For example, when the thalweg is close to the
bank, the crosssections have asymmetric form. If the thalweg is
approximately at the center of the channel, the crosssections have a
symmetric form.
In summary, the direction of channel meander relates to the thalweg
location, which in turn helps to define the crosssectional asymmetry. Therefore,
if the knowledge about the channel planform is available, it may be possible to
predict the crosssectional 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 crosssectional 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 stepbystep 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 crosssections
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 crosssections
(9)
in the crosssections with profilelines to create
a threedimensional 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 nondimensional form,
in which the width and the depth of the channel are unity. For example, Figure 5.3
shows a typical crosssection 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 crosssection. 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 crosssection with (s,n,z) coordinates. P(n
i
,z
i
) is any point in
a crosssection; 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 nondimensional 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 crosssection 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 ncoordinate (n
L
) and maximum ncoordinate
(n
R
) are used to calculate the width (W = abs[n
L
] + n
R
). Similarly the points with
minimum zcoordinate (Z
d
) and maximum zcoordinate (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 nondimensional form. The conversion only changes the
scale while preserving the original shape of the crosssection.
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 crosssection. (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 crosssection 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 crosssectional 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 crosssection 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 ncoordinate). 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 crosssection. If the radius of curvature is large, the
thalweg is approximately at the center of the channel with a nearly symmetric
crosssection. 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 1014 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 subsegments 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
nondimensional 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 nondimensional 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 Crosssectional 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 crosssectional shape has to be quantified with an analytical form.
The concept of fitting an analytical form to crosssections is not new. Attempts
have been made to fit analytical forms to glacial valley crosssections. For
example, James (1996) compared power functions and secondorder polynomials
for glacial valley crosssections 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 crosssections. 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 crosssections is found, the significance of each
parameter is tested by doing a nonlinear regression analysis using SAS (Cody
and Smith, 1997).
A brief description of all the functions that are considered for fitting
channel crosssections is given below by using one sample crosssection 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 crosssections such as glacial valley crosssections 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
crosssection from the thalweg. Therefore, the crosssection has to be split at the
thalweg, and each side of the crosssection has to be fit with an individual power
function. Figure 5.9 shows bathymetry data for a crosssection.
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 crosssection
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 crosssection 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 crosssection 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
m1
x
m1
+ a
m2
x
m2
+?+ a
1
x +
a
0
, where x is the independent variable, m is a nonnegative integer, and a
m1
,
a
m2
,?, a
0
are the coefficients of the polynomial function. The nonnegative
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
crosssection. A quadratic function does produce one bend, but it does not provide
enough flexibility to fit the crosssectional 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 crosssection
A thirddegree 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*
Thirddegree Polynomial
Bathymetry Points
Extra Bend
Figure 5.13 Thirddegree polynomial fitted to the channel crosssection
A thirddegree 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 Bspline function that is fitted to the observed data. A Bspline 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 nonnegative 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 Bspline with ten
knots fitted to the crosssection 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 crosssection
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 crosssectional area towards the end of the river channel cross
sections. Second, in the normalized domain, z* < 1 for channel crosssections.
This condition (z* < 1) is violated with a single beta distribution. In other words,
it is difficult to fit a pdf to a channel crosssection 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 crosssection 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 crosssection
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 crosssection 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 crosssection
Similar to beta distribution, the parameters of gamma distribution are also
determined by using the NewtonRhapson 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 crosssection. 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(xa
?
?b
?
) + f(xa
2
?b
2
)}*k, where a
??
?
?
b
??
, a
2?
=
?
b
2
, and k < 1. (5.19a)
Which for fitting a crosssection 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 crosssection.
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 crosssection
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 nonlinear 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 crosssection 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 crosssections so that the cross
sectional form can be related to the thalweg location.
Relating Channel Crosssections to Thalweg locations
This section describes step 6 of the flowchart shown in Figure 5.2. Once
the channel crosssections 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 crosssections) 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 nonlinear 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 crosssection 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 crosssection 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 crosssection 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 pvalue 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 crosssections
This section describes steps 7 and 8 of the flowchart shown in Figure 5.2.
The crosssections 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 crosssectional 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
channelwidth 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 channelwidth,
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 timeseries of stream
levels, steamflow, and crosssection measurements for more than 850,000 gaging
stations. To develop the hydraulic geometry relationships for the Brazos River,
flow data and crosssection measurements are downloaded from the USGS
website (http://nwis.waterdata.usgs.gov/usa/nwis/measurements). The USGS
performs crosssection measurements, which include channelwidth, 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
crosssections surveyed at the gaging stations. For example, the crosssections 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
(twoyear 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 2year 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
(channelwidth 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 crosssections. 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 2year return period flow is used as a standard
to develop hydraulic geometry relationships and rescale the normalized cross
sections. In addition, the 2year 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 2year 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 crosssections.
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 profilelines using crosssections
This section describes step 9 of the flowchart shown in Figure 5.2. To
create a complete threedimensional description of river channels, the cross
sections must be joined by profilelines. Profilelines can be generated using two
different approaches: 1) depthbased approach and 2) areabased approach. The
depthbased approach divides the crosssection into regions with equal depths. If a
crosssection that is 10m deep (from banks to the thalweg) has to be divided into
three regions, the depthbased approach will introduce the profile lines at d/3 on
each side (Figure 5.27). On the other hand, the areabased approach will divide
the crosssection into regions with equal areas. After the profilelines are
generated, they are converted into Bezier curves for creating a smooth
representation of the threedimensional channel form.
a) Depthbased profilelines b) Areabased profilelines
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 profilelines using crosssections
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 threedimensional 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 2year 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 2year flow. The 2year 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 curvaturethalweg 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 crosssections. The parameters in
Figure 5.23 can be used to generate beta distributions for any thalweg
location.
5. Rescale the crosssections 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 crosssections 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 curvaturethalweg location
relationship (equations 5.5 and 5.6). Circled areas show locations where the
thalweg prediction is imprecise. Sections AA and BB 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 nondimensional 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 nondimensional 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 crosssections using a combination of two beta distributions
(equation 5.16). The resulting crosssections are then rescaled to fit the observed
data. Figure 5.30 shows a crosssection described by RCMM at location AA
shown in Figure 5.28.
16
18
20
22
20 0 20 40 60 80
ncoordinate (m)
24
z

c
o
or
di
n
a
t
e
(
m
)
Observed
Model
Figure 5.30 Crosssection described by RCMM at section AA shown in Figure
5.28
The crosssection 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
crosssection BB shown in Figure 5.28, the crosssection 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 Crosssection described by RCMM at section BB shown in Figure
5.28
The disagreement between the model and the data at crosssection BB is a
result of an abrupt change in the bathymetry just upstream. Crosssection BB 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 crosssections are generated, profile lines are
generated to get a threedimensional mesh as shown in Figure 5.32. The profile
lines are generated using the depthbased approach to divide the crosssections
into ten regions.
173
Figure 5.32 Threedimensional 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 crosssections and profilelines 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 crosssections
are related to the thalweg location. If the thalweg is not located precisely, RCMM
is unable to create a reasonable description of crosssections 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 nondimensional 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 nondimensional 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 = 1014 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 straightline 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 crosssections
The beta distributions describe the general shape of the crosssections
based on the thalweg location, and the general crosssectional 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 threedimensional 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 crosssections are approximately 500m apart.
It is clear from Figure 5.50 that RCMM is a useful tool for describing the
threedimensional 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 sdirections, respectively.
Model
Data
(m
2
)
(m)
Figure 5.51 Semivariogram of prediction errors (across flow, ndirection) for Site
2 on Brazos River
195
196
Model
Data
(m)
(m
2
)
Figure 5.52 Semivariogram of prediction errors (along flow, sdirection) 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
ndirection sdirection ndirection sdirection
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
ndirection sdirection ndirection sdirection
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 ndirection is still a little
higher than the semivariance in the sdirection.
The stochastic smallscale 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 threedimensional 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 crosssectional form so that the
channel planform can be used to describe the crosssections 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 crosssections; and 4) use of hydraulic
geometry relationships for scaling normalized crosssections.
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 crosssections.
200
Chapter 6 Conclusions and recommendations
This dissertation has two main components. The first component deals
with a procedure to create a threedimensional of river channel bathymetry in the
form of crosssections and profilelines. The second component deals with
development of an analytical model for describing the threedimensional 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 THREEDIMENSIONAL DESCRIPTION OF RIVER CHANNELS USING CROSS
SECTIONS AND PROFILELINES
The PolylineMZ geometry available in the ArcGIS system is used to
describe the threedimensional features of river channels, such as thalweg, cross
sections, and profilelines. 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 nonvariable. 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 crosssections and profilelines 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 crosssections and profilelines. 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 crosssections
and profilelines (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 crosssections and profilelines 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 interrelating the
channel planform, the thalweg location and the crosssectional form with each
other. Therefore, by knowing only the channel planform, the threedimensional
form of the channel can be described. The resulting threedimensional description
is network of crosssections 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 nondimensional 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 crosssectional form are derived by using the data at Site 2 along
the Brazos River. However, unlike the thalwegradius 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 crosssections, 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 crosssectional form. Therefore, the parameters of the beta
distributions used for describing the crosssectional form may have to be changed
while using RCMM under different geologic and climatic settings.
The crosssections described in normalized domain are rescaled to fit real
crosssections 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 crosssections are then used to generate profile lines, which are
204
described using Bezier curves. However, unlike crosssections, the Bezier curves
used to represent profilelines are not related to any channel properties. In effect
the parameters of Bezier curves are related to the crosssections, 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 largescale 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 crosssections and profilelines 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 crosssections 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 crosssections can also be studied to enhance RCMM.
In addition, meandering river channels also have largescale 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 largescale variable features. Finally,
RCMM should also be modified to account for smallscale 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 subfolder 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 threedimensional 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 bottomright
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
crosssections 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 bottomright 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
HECRAS. 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. 23, pp. 211221.
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 0011.
Austin Barney, and Wentzel Mark, 2001. Twodimensional 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 995.
B?zier, Pierre, E., 1993. The First Years of CAD/CAM and the UNISURF CAD
System, in Fundamental Developments of ComputerAided Geometric
Modeling, edited by Les Piegl, Academic Press, pp. 1326.
Biftu, G. F., and Gan T. Y., 2001. Semidistributed, physically based, hydrologic
modeling of the Paddle River Basin, Alberta, using remotely sensed data.
Journal of Hydrology. Volume 244, No. 34, pp. 137156.
Blench, Thomas, 1970. Regime theory design of canals with sand beds. Journal of
the Irrigation and Drainage Division, ASCE Proceedings. Volume 96, pp.
205213.
226
Bovee, Ken, D., 1996. Perspectives on twodimensional river habitat models: the
PHABSIM experience. Proceedings of Second International Symposium
on Habitat Hydraulics Ecohydraulics 2000. Quebec, Canada, June 1996,
Volume B, 149162.
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. 1635.
Callander, R., A., 1978. River Meandering. Annual Review of Fluid Mechanics.
Volume 10, pp. 12958.
Carter, Glenn, S., and Shankar, Ude, 1997. Creating rectangular bathymetry grids
for environmental numerical modeling of gravelbed rivers. Applied
Mathematical Modeling. Volume 20, No. 11, pp. 699708.
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. 146177.
Chang, Howard, H., 1980. Geometry of gravel streams. Journal of the hydraulics
division, ASCE Proceedings. Volume 106, pp. 14431455.
Chau, K. W., and Jiang, Y. W, 2001.3D numerical model for pearl river estuary.
Journal of Hydraulic Engineering. Volume 127, No. 1, pp. 7282.
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 multilayer perceptrons.
Proceedings of the European Symposium on Artificial Neural Networks,
Bruges, Belgium, April 1999.
Crowder, D., and Diplas, P., 2000. Using twodimensional hydrodynamic models
at the scales of ecological importance. Journal of Hydrology. Volume 230,
Issues 34, pp. 172191.
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. 179227.
Djokic, D., Beavers, M., A., and Deshakulakarni, K., 1994. Arc/HEC2 An
ArcInfoHec2 Interface. Proceedings of the 21
st
annual conference on
water policy and management. American Water Resources Association,
ASCE, pp. 4144, May 1994.
Dodov, B., and EfiFoufoulaGeorgiou, 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. 511529.
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. 683700.
228
Fissel, David, Birch, Rick, and Jiang, Jianhua, 2002. Threedimensional
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 1112,
pp. 20892108.
Froehlich, D., C., 1992. Finite element surfacewater modeling system: two
dimensional flow in a horizontal plane. U. S. Department of
Transportation Report No. FHWARD92057.
Fukoka, S., and Sayre, W., W., 1973. Longitudinal dispersion in sinuous
channels. Journal of Hydraulics Division, ASCE Proceedings. Volume 99,
pp. 195217.
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. 23,
pp. 185200.
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. 257284.
Goff, John, A., and Nordfjord, Sylvia, 2004. Interpolation of Fluvial Morphology
Using channeloriented 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. 299318.
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 2D 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. 869885.
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. 949958.
Holley, E., R., and Jirka, G., H., 1986. Mixing in rivers. Technical Report E86
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. 413432.
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. 181213.
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. 414419.
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. 10631077.
Knighton, A. D., 1981. Asymmetry of river channel crosssections: Part I.
Quantitative indices. Earth Surface Processes and Landforms. Volume 6,
pp. 581588.
Knighton, A. D., 1982. Asymmetry of river channel crosssections: Part II. Mode
of development and local variation. Earth Surface Processes and
Landforms. Volume 7, pp. 117131.
Knighton, A. D., 1984. Indices of flow asymmetry in natural streams: definition
and performance. Journal of Hydrology. Volume 73, pp. 119.
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. 259292.
Leclerc, M., Bobee, A., Boudreault A., Shooner, G., and Corfa, G., 1991.
Instream flow incremental methodology and 2D hydrodynamic modeling:
efficient tools to determine guaranteed minimum flow for biological
purposes. In: Computer methods in water resourcesII, proceedings of the
second international conference. Marrakesh, Morocco, 2022 February
1991, pp. 289300.
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 PSW394, Pacific Southwest
Forest and Range Experiment Station, Forest Service, U.S. Department of
Agriculture, Berkeley, CA.
Maidment, D., 2002. Arc HydroGIS 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. 799816.
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. 983992.
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, 481492.
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. 12511266.
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. 69102.
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. 69102.
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. 367383.
Parasiewicz, P., and Dunbar M., J., 2001. Physical habitat modelling for fish: a
developing approach. Large Rivers. Volume 12. No. 24, pp. 239268.
Park, C., C., 1977. Worldwide variations in hydraulic geometry exponents of
stream channel: An analysis and some observations. Journal of
Hydrology. Volume 33, pp.133146.
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. 379416.
Pebesma, Edzer, J., and Wesseling, Cees, G., 1998. Gstat: a program for
geostatistical modelling, prediction and simulation. Computers &
Geosciences. Volume 24(1), pp. 1731.
Rathburn, Sara, and Wohl, Ellen, 2003. Predicting fine sediment dynamics along
a poolriffle mountain channel. Geomorphology. Volume 55, Issues 14,
pp. 111124.
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. 24, pp.
189209.
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. 143167.
Rosgen, David, L., 1994. A classification of natural rivers. CATENA. Volume 22,
pp. 169199.
Runge, Morten, and Olesen, Kim Wium, 2003. Combined one and two
dimensional flood modeling. Fourth Iranian Hydraulic Conference, 2123
October, Shiraz, Iran.
Savenije, Hubert, H.G., 2003. The width of a bankfull channel; Lacey's formula
explained. Journal of Hydrology. Volume 276, Issues 14, pp. 176183.
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, 352B: 1730
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. 267320.
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. 3371.
234
Sinha S. K., Sotiropoulos F., and Odgaard A. J., 1998. Threedimensional
numerical model for flow through natural rivers. Journal of Hydraulic
Engineering. Volume 124, No. 1, pp. 1324.
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 0012.
Stalnaker Clair, Lamb Berton L., Henriksen Jim, Bovee Ken, and Bartholow
John, 1995. The instream flow incremental methodologyA 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. 125.
Tarbet, Karl, and Hardy, Thomas, B.,1996. Evaluation of onedimensional and
twodimensional 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, 395406.
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. 100108.
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 359373.
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. 1830.
U. S. Army Corps of Engineers (USACE), 2002a. HECGeoRAS: An extension
for support of HECRAS using ArcView. Users Manual. Version 3.1,
October 2002, Hydrologic Engineering Center, Davis, California.
U. S. Army Corps of Engineers (USACE), 2002b. HECRAS: 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. Individualbased model of sympatric
populations of brown and rainbow trout for instream flow assessment:
model description and calibration. Ecological Modeling. Volume 110, No.
2, pp. 175207.
Waddle, T., Bovee, K., Bowen, Z., 1997. Twodimensional habitat modeling in
the Yellowstone / Upper Missouri River System. North American Lake
Management Society (NALMS) meeting in Houston, TX (December 24
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.1Applications and
Cartography edited by Michael Zeiler. ESRI Press, Redlands, California,
pp. 13.
Wentzel, Mark, W., 2001. Twodimensional 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 78, pp. 517522
236
Werritty, Alan, and Leys, Katherine, F., 1999. River channel planform change:
software for historical analysis. Geomorphology. Volume 29, Issues 12 ,
pp. 107120.
Williams, G. P., 1978. Bankfull Discharge of Rivers. Water Resources Research.
Volume 14, No. 6, pp. 11411154.
Winterbottom, S. J., 2000. Medium and shortterm channel planform changes on
the Rivers Tay and Tummel, Scotland. Geomorphology. Volume 34, No.
34, pp. 127132.
Yalin, M., S., 1992. River Mechanics. Pergamon Press, Oxford.
Yang, Qingwen, 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. 699717.
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. 146157.
Ye, J., and McCorquodale, J. A., 1997. Depthaveraged hydrodynamic model in
curvilinear collocated grid. Journal of Hydraulic Engineering. Volume
123, No. 5, pp. 380388.
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 preuniversity
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 postgraduate 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.