CHANNEL MIGRATION MODEL FOR MEANDERING RIVERS
by
Timothy J. Randle
B.S.C.E., University of Utah, 1981
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Masters of Science
Civil Engineering
2004
@2004 by Timothy J. Randle
All rights reserved.
This thesis for the Master of Science
degree by
Timothy J. Randle
has been approved
by
Randle, Timothy J. (M.S. Civil Engineering)
Channel Migration Model for Meandering Rivers
Thesis directed by Professor James Guo
ABSTRACT
A numerical channel migration model has been developed for meandering rivers to
predict the future river channel alignment. The model is based the hypothesis that the
magnitude of river bank erosion and channel migration is primarily a function of the
channel radius of curvature, stream discharge, river channel dimensions, sediment
transport capacity, and the bank material properties acting to resist the erosion. A
variable planformphase lag is used in the model to determine the proportions of
channel migration that are lateral (across the floodplain) and downstream (along the
valley axis). When the phase lag is zero, the meander bends migrate laterally across
the floodplain so that the meander amplitude and sinuosity both increase. As the
channel sinuosity increases, the variable phase lag increases and more of the channel
migration occurs along the downstream valley axis.
The model input requirements include the initial channel alignment; bank material
properties of the channel, floodplain and terraces; channel roughness and bed
material grain size; discharge hydrograph; and the upstream sediment supply. The
model was applied to a 4.2mile reach of the Hoh River near Forks, Washington to
replicate the changes that occurred over the 52year historic period from 1950 to
2002. The model predicted the magnitude and planform phase lag of channel
migration reasonably well. Although the model has no explicit criteria for the cutoff
of meander bends, a cutoff can be inferred from the model results.
This abstract accurately represents the content of the candidates thesis. I recommend
its publication.
Signed
IV
CONTENTS
Figures..................................................................vii
Tables.....................................................................x
Chapter
1. Introduction............................................................1
1.1 Background.............................................................3
1.2 Thesis Statement.......................................................4
2. Literature Review.......................................................5
2.1 Reasons for River Meandering...........................................6
2.2 Planform Geometry of River Meander Bends..............................9
2.2.1 Bankfull Channel Width...............................................9
2.2.2 Meander Bend Wavelength.............................................10
2.2.3 Planform Layout.....................................................12
2.2.4 Equations of Channel Hydraulics and Geometry........................15
3. Methodology............................................................17
3.1 Riverbank Erosion Rate................................................18
3.2 Planform Phase Lag....................................................21
3.3 Channel Hydraulic Equations...........................................22
3.4 Channel Curvature Equations..........................................27
3.5 Model Input Requirements.............................................31
3.5.1 Initial Conditions..................................................31
3.5.2 Boundary Conditions.................................................33
3.6 General Model Algorithm..............................................33
4. Example Results........................................................35
5. Model Calibration and Validation.......................................38
v
6. Discussion..............................................................50
7. Conclusions.............................................................51
Appendix
A. Channel Migration Model: Fortran Program Listing........................52
Bibliography.............................................................. 162
vi
FIGURES
Figure 1. Erosion of a terrace bank along the Hoh River, near Forks,
Washington, threatens the county road near the edge of the floodplain.........1
Figure 2. Polynomial regression equation relating the channel sinuosity
(0) to the maximum angle (/3) of a meander loop, relative to the valley
axis..........................................................................14
Figure 3. Comparison of a sinegenerated meander curve, sine curve, and
a halfcircle curve of constant radius........................................14
Figure 4. Planform phase lag results in at least a portion of the meander
bend migrating downstream along the valley axis...............................21
Figure 5. Possible cross section shapes considered in the model. The
trapezoidal shape is used for river crossings while the triangular shape is
used for meander bends........................................................24
Figure 6. Example computation results for the determining the minimum
unit stream power (US').......................................................26
Figure 7. Example of smoothing the initial channel alignment points
along the Hoh River, one point at a time, using a regression analysis of
nearby points. Seven consecutive points are selected for regression
analysis to smooth the middle point of the selection..........................28
Figure 8. Definition sketch for computing the radius of curvature along a
meander bend..................................................................29
Figure 9. Example grid used to describe the variation in material
properties throughout the river channel, floodplains, and terraces............32
Figure 10. Example model results assuming an initial channel alignment
from a sine generated curve. The flow is from left to right...................37
Figure 11. Closeup view of example model results shown in figure 10..........37
Figure 12. Initial conditions for model simulation were digitized from the
1950 aerial photograph........................................................39
vii
Figure 13. Correlation in meandaily discharge between two Hoh River
streamgage locations for the period 1961 to 1964..............................40
Figure 14. Correlation in annualpeak discharge between two Hoh River
streamgage locations for the period 1961 to 1964..............................41
Figure 15. Meandaily and annualpeak discharge of the Hoh River..............41
Figure 16. The Manning roughness coefficient, n, was specified as a
power function of the hydraulic depth, D......................................42
Figure 17. Computed sediment load and concentration for the upstream
supply to the Morgan's Crossing Reach of the Hoh River........................43
Figure 18. Areas of different material properties for the Hoh River
channel, floodplains, and terraces along the Morgan's Crossing Reach.
Flow is in the westerly direction.............................................43
Figure 19. Predicted channel alignments for the period 1950 to 2002...........45
Figure 20. Initial (1950) and predicted (1960) channel alignments plotted
on the 1960 aerial photograph..................................................46
Figure 21. Initial (1950) and predicted (1971) channel alignments plotted
on the 1971 aerial photograph..................................................47
Figure 22. Initial (1950) and predicted (1977) channel alignments plotted
on the 1977 aerial photograph. A meander cutoff occurred (circled area)
sometime between 1971 and 1977. This meander cutoff was not predicted
by the model...................................................................47
Figure 23. Initial (1950) and predicted (1981) channel alignments plotted
on the 1981 aerial photograph..................................................48
Figure 24. Initial (1950) and predicted (1994) channel alignments plotted
on the 1994 aerial photograph. The sharp channel bend (near 80 degrees)
predicted by the model could be interpreted as a meander bend cutoff..........48
Figure 25. Initial (1950) and predicted (2001) channel alignments plotted
on the 2001 aerial photograph. The sudden change in predicted channel
alignment indicates a meander cutoff................................
49
Figure 26. Initial (1950) and predicted (2002) channel alignments plotted
on the 2002 aerial photograph..............................................49
IX
TABLES
Table 1. Regime equations by Soar and Thome (2001) which include the
95 percent confidence limits on the mean response.............................11
Table 2. Example model input data for the hypothetical problem................36
Table 3. Available streamflow gage records for the Hoh River.................40
Table 4. Model input data for equation 17, calibrated to the 4.2mile
reach of the Hoh River near Forks, Washington.................................44
x
1. Introduction
The migration of river channels across their floodplains and the occasional erosion of
terrace banks is a natural process (Leopold et al. 1964, Yang 1971, Dunne and
Leopold 1978, Leopold 1994, and Thome 1992 and 2002). This process becomes
especially import to people living in or near the floodplain or to organizations
planning or maintaining infrastructure within or along the edge of the floodplain (see
figure 1). Natural rates of channel migration can be accelerated by human or natural
disturbance. For example, the clearing of floodplain vegetation can accelerate the
rate of channel migration (Beeson and Doyle 1995).
Figure 1. Erosion of a terrace bank along the Hoh River, near Forks, Washington,
threatens the county road near the edge of the floodplain.
1
The development and application of a channel migration model for meandering river
is described in this thesis. The Introduction (Chapter 1) describes the need for a
channel migration model and the hypotheses to be tested by this research. The
Literature Review (Chapter 2) describes the present state of knowledge about the
reasons for river meandering and the geomorphic relationships between river
discharge, bankfullchannel width, meander wavelength, sinuosity, and the energy
slope for meandering rivers. The Methodology (Chapter 3) describes the
mathematical equations, data input requirements, and the model algorithm. Example
Results (Chapter 4) describes the application of the model to an idealized problem.
The Model Calibration and Validation (Chapter 5) is described through the simulation
of a 4.2mile reach of the Hoh River near Forks, Washington for a 52year period.
The Discussion of model performance is presented in Chapter 6. The Conclusions
about the hypotheses, the usefulness of the model, and the needs for future research
are presented in Chapter 7. The Bibliography (Chapter 8) lists the references for all
literature cited in this report. The Appendix includes the complete listing of the
Fortran computer code developed for this model.
The purpose of this research is to develop a numerical model for predicting future
alignments of meandering river channels given the following information:
existing alignment of the river channel thalweg
river valley alignment and slope
discharge hydrograph for the period of simulation
upstream sediment supply
bedmaterial sediment size of the river channel
channel roughness
bank material properties
riparian vegetation characteristics
The importance of riparian vegetation is described by Beeson and Doyle (1995). The
importance of large woody debris is described by Collins and Montgomery (2002).
Current and historic data of the river alignment and overall vegetation densities can
be obtained from orothorectified aerial photographs and digital elevation models.
River discharge data can be obtained from stream gage records. Sometimes, sediment
load or concentration data also can be obtained from stream gage records. When
sediment load data are not available, the upstream sediment supply can be computed
as a function of the sediment transport capacity at a stable cross section. The bed
material sediment size data can be obtained from field measurements. Estimates of
the channel roughness, bank material properties, and riparian vegetation
characteristics can be obtained during field inspections.
2
The channel migration model could be used as a planning tool to help evaluate the
risk of bank erosion and the potential impacts to structures along a river channel, such
as roads, bridge embankments, and buildings. The model could also be used to assess
the effects of bank protection structures along the river channel on the downstream
channel alignment. In addition, the model could be used to evaluate the effects of
alternative land management practices.
1.1 Background
Numerous numerical models have been developed over the past few decades to
predict river hydraulics and the transport, erosion, and deposition of sediment. For
example, onedimensional hydraulic models, such as HECRAS (Brunner 2002), have
been commonly applied to natural river channels to predict flood stage and mean
river velocities. Twodimensional hydraulic models, such as MIKE21 (DHI Water
and Environment 2002), are gaining more popularity and threedimensional hydraulic
models are also available (Lai et al. 2003). Typically, hydraulic models assume that
the channel geometry is static and they ignore the possible erosion or deposition of
sediment along the channel bed and banks.
Onedimensional sediment transport models, such as HEC6 (U.S. Army Corps of
Engineers 1993), STARS (Orvis and Randle 1987) and GSTARS (Yang and Simoes
2003 and Yang et al. 2004), have been applied to natural river channels to predict
aggradation and degradation of the river bed. The GSTARS model also has some
ability to predict whether sediment deposition or erosion should be applied to the
channel bed or banks, but the initial river alignment is fixed. These sediment
transport models have been successfully applied to a wide variety of river channels
where the erosion or deposition of sediment is important. However, the existing one
dimensional sediment transport models cannot predict the rate and extent of channel
migration for meandering river channels.
Twodimensional hydraulic and sediment transport models have been developed
(Darby et al. 2002), but the prediction of bank erosion is difficult. This is because
secondary currents, the important threedimensional effects of the velocity field
through a meander bend (see section 2.0), are not considered in a twodimensional
model. In addition, the time scale that two and threedimensional models can be
applied is typically limited to weeks or months.
Therefore, there is a need for a numerical model to aid in the prediction of river
channel migration over a time scale of years or decades. The model would have to be
calibrated to historic data and could then be applied to future predictions based on
assumptions of river flow and sediment load.
3
1.2 Thesis Statement
The linked set of hypotheses, represented by a numerical model, can be used to
predict the future rate and extent of migration of a single meandering river channel
where there is no meander cutoff. The rate of river channel migration is primarily a
function of the hydraulic forces acting to erode the bank and the properties of the
bank acting to resist the erosion.
Hypothesis 1: The hydraulic forces acting to erode the bank can be predicted
by the river discharge, channel slope, mean velocity, width, depth, and radius of
curvature. Except for the radius of curvature, all these hydraulic parameters are
accounted for in the hydraulic capacity to transport bedmaterial load, which can be
expressed as a dimensionless sediment concentration (e.g., parts per million, ppm, by
weight).
Hypothesis 2: The bank material properties acting to resist erosion can be
predicted by the riparian vegetation root density and root length, bank height, the
diameter and density of trees and large woody debris along the river bank, the
cohesive properties of the bank, and the percentage of coarse sediment that is too
large for incipient motion.
Hypothesis 3: The minimum rate of energy dissipation theory; combined with
Mannings equation, the continuity equation, and a sediment transport equation; can
be used to predict the channel width, depth, velocity, and slope. The prediction of the
channel slope, associated with the minimum unit stream power, can be used to
determine the proportion of river channel migration that is across the floodplain and
the proportion that is down valley.
4
2. Literature Review
A tremendous volume of literature exists on the subject of river meanders. The
literature review for this project includes papers or text on the reasons for river
meandering and the planform geometry of natural river channels, including the
geomorphic relationships between river discharge, bankfullchannel width, meander
wavelength, sinuosity, and the energy slope for meandering rivers.
Meandering river bends are observed in natural river channels, at all scales, and in all
parts of the world (Leopold et al. 1964, Thome 1992, and Leopold 1994). The
upstream inputs of water and sediment interact with the downstream landscape, the
material through which the river is flowing, and the riparian vegetation. For example,
two different rivers, which have the same sequence of river flows, but with different
sediment loads, will have different meander and planform characteristics.
As river flow passes through a meander bend, a secondary current is induced that
causes flow to travel in a downstream helical motion, which is downward near the
outer bank, across the deep portion of the channel bottom toward the inner bank, and
then back across the channel near the water surface. At the upper portion of the outer
bank, a small cell of reverse circulation often exists that forces flow upward along the
bank. This reverse circulation is highly significant to the flow processes acting to
erode the outer bank (Thome and Hey 1979 and Thome 1992).
The meander bend flow circulation increases the flow velocities and shear stresses
along the outer bank and decreases them along the inner bank. Consequently,
sediment erosion occurs along the outer bank and deposition occurs along the inner
bank. This results in the deepest portion of the channel being near the outer bank and
a point bar deposit along the inner bank (Leopold et al. 1964, Thome 1992, and
Leopold 1994). The tendency to erode the outer bank and deposit sediment along the
inner bank causes the channel to migrate laterally. Typically, the maximum rate of
bank erosion along the outer bank is located downstream of the meander bend apex,
which causes the river channel to migrate both across the floodplain and downstream
along the river valley.
5
2.1 Reasons for River Meandering
Yang (1971) evaluated several published hypotheses for river meandering, which are
listed below:
rotation of the earth (Eaking 1910, and Neu 1967)
local disturbances (Werner 1951)
bank erosion and sedimentation (Friedkin 1945, Mattes 1941, and Ackers and
Charlton 1970)
maximization of energy loss (Jefferson 1902 and Inglis 1947)
minimization of energy loss (Langbein and Leopold 1966)
Yang (1971) concluded that most of these hypotheses emphasize some special
phenomena observed in meandering channels and neglect the basic physical
reasoning which creates them. The rotation of the earth does not explain the
fundamental reason for river meandering because rivers are known to meander at or
near the equator where the influence of the earths rotation on a river alignment is
zero. In addition, the earths rotation should have no effect on the meandering of
small scale streams.
The hypothesis that river meandering is initiated by random local disturbances does
not seem plausible given that meandering is a very consistent planform for streams of
all sizes throughout the world. River meandering has also been observed in streams
with a large variety of geologic settings and sediment transport rates. Leopold (1994)
observed meandering in streams cut in glacial ice with virtually no sediment transport
and he observed meandering currents in the ocean where there are no river banks.
Jefferson (1902) and Inglis (1947) hypothesized that a meandering river alignment
was a means of dissipating excess energy. The steeper the valley slope, the more
meandering or sinuous the river alignment will be, within the constraints of the valley
walls. This is based on the hypothesis that water flowing through a curved channel
has a higher energy loss than water flowing through a straight channel. However, the
total energy loss, or drop in water surface elevation, over a long valley length is
determined by the valley slope and not whether the river alignment is straight or
meandering.
Langbein and Leopold (1966) hypothesized that the geometry of a meander is that of
a random walk whose most frequent form minimizes the sum of the squares of the
changes in direction in each unit length. They also state that meandering tends to
increase in the downstream direction, which results in a decrease in river slope, in
order to offset the increases in river discharge from tributaries and maintain a more
uniform total stream power along the channel length. Total stream power (7 Q S) is
6
defined by Chang (1979 and 1992) as the product of the unit weight of water (7),
discharge (Q), longitudinal slope (S).
Yang (1971) applied the concept of entropy, introduced by Langbein and Leopold
(1966), to a stream system to obtain two basic laws which govern the formation of a
natural stream network. The first law states that, under a dynamic equilibrium
condition, the ratio of average fall between any two different order streams, in the
same river basin, is unity. The second law states that during the evolution toward
its equilibrium condition, a natural stream chooses its course of flow in such a manner
that the time rate of potential energy expenditure per unit mass of water along this
course is a minimum, dependent upon the external constraints applied to the system.
Yang (1971) expressed the second law with the following equation.
ATT Jy
= = 0{Q,Sv, Cs,G..) = minimum (1)
At At
where Mil At = the time rate of potential energy expenditure per unit mass of water
Y = the drop in water surface elevation along a stream reach
k a. unit conversion factor between energy and the drop in water
surface
$ = minimum rate of energy expenditure subject to external constraints
Q = water discharge
Sv = valley slope
Cs = sediment concentration
G = geologic constraints
The strength of geologic constraints depends at least on the credibility of bank
material (including vegetation), channel roughness, and the stream valley width.
Yang (1971) showed that a meandering river alignment will always have lower rate
of energy dissipation than a straight stream and that the law is applicable to all
streams of all sizes and sediment concentrations. In addition, the law can explain
how river meandering changes with variations in water discharge, sediment
concentration, channel slope, channel geometry, valley slope, and geologic
constraints.
Chang (1992) points out that the river valley slope is usually formed by episodes of
sediment deposition over a period of centuries under various hydrologic conditions
that are different from those shaping the present channel. However, the valley slope
is not easily reduced during periods of erosion because channel degradation is slowed
or limited by coarsening of the bed material. The morphological response of a river
to a valley slope that is steeper than necessary to transport the upstream sediment
7
supply is for the channel to reduce its slope by meandering across the valley
floodplain.
Yang and Song (1979) introduced the minimum rate of energy dissipation theory and
unit stream power (velocityslope product), which is equal to the average rate of
potential energy dissipation per unit weight of water.
dY dx dY
dt , VS minimum subiect to dt dx
where Y= water surface elevation
t = time
x = distance along the river channel
Y= water surface elevation
V= average flow velocity
S = average channel slope
VS = unit stream power
The theory of minimum unit stream power states that, given enough time, an alluvial
river channel will adjust its velocity, slope, roughness, and geometry in such a
manner that a minimum amount of unit stream power is used to convey the upstream
discharge and transport the upstream sediment load (Yang and Song 1979). Because
the upstream inputs of water and sediment are continually changing, the river channel
is continually evolving in an attempt to reach the conditions corresponding to the
minimum unit stream power.
Yang (1973 and 1984) also developed sediment transport equations for sand and
gravel based on unit stream power (VS). The equation parameters were based on
dimensional analysis of the controlling variables and take the following form:
Log Cs = /+ JLog
VS
co
Krl
CD
(3)
where Cs =
Vcr~
W 
the sediment concentration of the bedmaterial load, in ppm
the critical mean velocity required for incipient motion of the sediment
particle
the fall velocity of the sediment particle
The terms, I and J, in equation 3 are defined in equations 4 and 5. The basic form of
equation 3 can also be derived from basic turbulence theory (Yang and Molinas
1982).
8
J = bv + b2Log
(5)
where d = the diameter of the sediment particle
v = the fluid viscosity
U* = the shear velocity U,
t shear stress
p = fluid density
7 = the unit weight of the fluid
D = hydraulic depth
The coefficients aj, ci2, CI3, bj, b2, and 63 were determined from multiple regression
analysis of laboratory data. The set of coefficients are different for sand and gravel
size sediments.
2.2 Planform Geometry of River Meander Bends
Many methods have been proposed over the years to predict stable channel geometry.
Most calculate width and depth and use discharge as the independent variable. Many
equations neglect sediment entirely, others use the bedmaterial sediment size, and
most assume a constant sedimentdischarge relationship.
2.2.1 Bankfull Channel Width
The hydraulic geometry approach (Leopold and Maddock 1953) uses empirically
based power equations to define relationships among channel variables based on the
bankfull discharge.
where W = the channel width
Q = the reference channel discharge (e.g., mean annual or bankfull discharge)
a = empirical coefficient
W = aQb
(6)
b = empirical exponent, often set equal to 0.5
Several regime equations describing channel geometry have been developed based on
measured data of many rivers. Lacey (1929), Simons and Albertson (1963), Blench
(1957), and Williams (1978) have presented such equations for sand bed streams.
Gravelbed equations have been presented by Kellerhals (1967), Bray (1982), and
Hey and Thome (1986). Julien and Wargadalam (1995) use a semitheoretical
approach based on the principals of open channel flow. The equation coefficients do
vary between locations and flow regimes and must be checked for applicability to a
specific site.
Soar and Thome (2001) present regime equations for natural sand and gravelbed
rivers considered to be in dynamic equilibrium. Different equations are presented for
riverbanks that are resistant to erosion because of riparian trees and for riverbanks
that are more erosive because of a lack of riparian trees. For 58 sandbed rivers in the
United States, two regime equations explained at least 85 of the variance in bankfull
width as a function of the bankfull discharge (see table 1). For the 85 gravelbed
rivers in the United States and United Kingdom, the four regime equations explained
at least 92 of the variance in bankfull width as a function of the bankfull discharge.
Other approaches use minimum stream power, sediment concentration, or energy
dissipation to predict channel shape. Yang (1971 and 1976) and Chang (1979 and
1992) have used minimum stream power to provide a third relationship beyond
sediment transport and alluvial resistance relationships to calculate channel geometry.
Care must be taken to validate any of these approaches for a particular location and
time period.
2.2.2 Meander Bend Wavelength
Leopold and Wolman (1957 and 1960) developed a predictive equation for the
wavelength of meander bends, \ as a function of the bankfull channel width. Soar
and Thome (2001) determined the equation at the 95 percent confidence limits
A = (11.26 to 12.47) (7)
where X= meander wavelength
Wb = bankfull channel width
10
Table 1. Regime equations by Soar and Thome (2001) which include the 95 percent confidence limits on the mean response.
River Type Regime Equations in Metric Units1 Regime Equations in English Units2 Number of rivers Coefficient of Determination (R2)
Sandbed rivers in the United States with less than 50 percent treelined banks Wb = 5.32 Qb0'5 e0'082 Wb = 2.94 Qb0'5 e0'082 32 0.87
Sandbed rivers in the United States with more than 50 percent treelined banks Wb = 3.38 Qb0 5 e0'085 Wb = 1.87 Qb0'5 e0'085 26 0.85
Gravelbed rivers in the United States with thin bank vegetation Wb = 4.17 Qb0'5 e0'087 Wb = 2.30 Qb05 e0'087 9 0.95
Gravelbed rivers in the United States with thick bank vegetation Wb = 3.67 Qb e Wb = 2.03 Qb0'5 e0'065 14 0.96
Gravelbed rivers jn the United Kingdom with less than 5 percent tree and shrub vegetation along the banks Wb = 3.75 Qb05 e0'064 Wb = 2.07 Qb0'5 e0,064 29 0.92
Gravelbed rivers in the United Kingdom with at least 5 percent tree and shrub vegetation along the banks Wb = 2.48 Qb0'5 e0'051 Wb =1.37 Qb0'5 e0 051 33 0.93
1 The bankfull discharge (Qb) is in m3/s anc 2 The bankfull discharge (Qb) is in ft3/s and Note: The term, e00xx in these equations, the bankfull width (Wb) the bankfull width (Wb) accounts for the 95 percc is in meters, is in feet. :nt confidence limits on the mean response.
2.2.3 Planform Layout
Langbein and Leopold (1966) found that the use of a sine generated function, to
describe the horizontal planform curve, minimizes the sum of squares of the changes
in direction in each unit length.
(j) = /?sin
2 n
(8)
where =
P =
s =
\ =
Q
angle of the channel direction, relative to the central downstream
valley axis
maximum angle of the meander bend relative to the central
downstream valley axis
distance along the river channel
meander wavelength
channel sinuosity
The maximum angle of the meander bend, /3, determines the channel sinuosity, Q.
Soar and Thome (2001) presented the following equation for Q as a function of /3:
Q =
4.84
4.84 p
(9)
Soar and Thome (2001) presented the following equation for the amplitude of a
meander bend:
An =
m
3/?
f
1coslfJ+^sinGu^
(10)
where Am= the meander bend amplitude
aa = an adjustment factor such that the product (au ft) equals the average
planform angle through the nearstraight reach between bends.
In order to illustrate these planform equations, a theoretical meander bend alignment
was generated by specifying the meander bend wavelength and amplitude or
sinuosity. For comparison, three types of theoretical meander curves were generated:
1. circle curves or curves of constant radius,
2. sine curves, and
3. sine generated curves.
12
The sine generated curve minimizes the sum of the squares of the changes in direction
in each unit length along the meander curve. The xy coordinates of a sine generated
meandering channel were determined from an initial coordinate point and by using
equation 8, along small increments of channel length, s, to determine the local
direction of the alignment.
Since equation 8 produces the tangent or derivative of the meander bend, the integral
of equation 8 produces the distance, y, from the central downstream valley axis.
For this study, an empirical polynomial equation was derived for channel sinuosity,
Q, as a function of maximum angle of the meander bend relative to the valley axis, /3.
This equation correlates well with the data (see figure 2) and is applicable for the
range of /3 between 0.17 and 1.9 radians.
The initial coordinate point for the theoretical alignment was defined at x = 0 where
the direction of the meander curve, (j>, is 0 radians, and value of the ycoordinate was
computed from equation 11. Using these assumptions, a sine generated curve was
derived for a meander wavelength of 2,000 feet, a maximum angle of 70 degrees
(1.2217 radians), and a sinuosity from equation 12 of 1.514. The plot of this sine
generated curve, along with a sine curve and half circle curves of constant radius, is
presented in figure 3.
(11)
Q = 0.9218/f4 2.8025fi3 + 3.3359/J2 1.2941/?+1.159
(12)
13
Sinuosity and Maximum Angle
MaximumAngleofa MeanderLoop.p (degrees)
50 60 70 80 90 100 110 120
Figure 2. Polynomial regression equation relating the channel sinuosity (O) to the
maximum angle (/3) of a meander loop, relative to the valley axis.
River Meander Curves
Figure 3. Comparison of a sinegenerated meander curve, sine curve, and a half
circle curve of constant radius.
14
2.2.4 Equations of Channel Hydraulics and Geometry
The book entitled: Fluvial Processes in River Engineering, by Chang (1992),
includes specific chapters on Flow in Curved Channels, Analytical River
Morphology, and Plan Geometry and Processes of River Meanders.
Chang (1992) presents equations for the additional friction slope through curved
channels where the total friction slope (S) is the sum of longitudinal friction slope (S)
and the transverse friction slope (S).
S= S'+ S" (13)
The additional transverse friction slope is attributed to the following causes:
1. Internal fluid friction due to secondary currents.
2. Boundary resistance associated with transverse flow.
3. Eddy loss resulting from flow separation in sharp bends.
4. Eddy loss due to sudden jumps occurring at super critical flow.
Chang (1992) presented an equation to predict the transverse friction slope as
function of the square of the Froude number (F), the square of the ratio depth (D)
over the radius of curvature (rc), and channel roughness:
S"=
( n c\n r A no yl/2 i oo ^3/2^ f 2
2.07/ + 4.68/ 1.83/
0.565 + /
1/2
\RCJ
where S =
f =
D =
Rc =
F =
g =
V=
energy gradient due to secondary currents
DarcyWeisback friction factor
depth of flow
center radius of curvature
Froude number F = v f=
/
gravitational acceleration
flow velocity
(14)
Chang (1992) also presented equations for the transverse bed slope (6) across a
meander curve.
15
il/2
1 + m'
(15)
. ~ 3a D
sm o =
2 R
V2
A _i
8da
m
i(2 + nf)
m = K
l/2
A_i
\p ,
gdcr^c
(16)
where 5
a =
Dy =
Ry =
Ps
P =
g =
dcr
T* =
Tc =
K =
the projected area to volume ratio of a sediment particle, normalized to
that of a sphere
projected areatovolume ratio for a sediment particle, normalized by
the areatovolume ratio for a sphere
depth of flow at some point across the channel
radius of curvature at some distance (y) across the channel
density of the sediment particle
density of the river flow
acceleration of gravity
diameter of a sediment particle whose motion is impending
critical Shields stress r, = zcj p)g
critical shear stress
calibrated constant
Chang (1992) states that the phenomenon of river meanders are related to the stream
wise variation of the helicalmotion strength of secondary currents. The changing
curvature of the channel flow path, and its phase lag with the channel curvature, have
import effects on lateral and downstream valley channel migration and should be
related to the variation in helicalmotion strength.
16
3. Methodology
The channel migration model is based (1) on a new equation to predict the rate of
bank erosion and (2) on a new application of minimum unit stream power (minimum
VS) to determine the planformphase lag between the changing curvature of flow and
the changing curvature of the river channel. The equations for the rate of bank
erosion and channel migration are described in section 3.1. A description of the
planformphase lag and the equations to compute it are presented in section 3.2. The
equations that relate the river discharge, sediment load, and channel dimensions are
presented in section 3.3. The equations for channel curvature are presented in section
3.4. The model input requirements are described in section 3.5. The general model
algorithm is presented in section 3.6.
The model simulates channel migration by computing bank erosion as a function of
the radius of channel curvature and the bank material properties resisting erosion. In
general, the bank erosion is applied along the outside river bank of each meander
bend with a corresponding volume of sediment deposition along the inside bend, thus
forming a point bar. A variable phase lag is incorporated so that bank erosion can be
applied at some distance downstream from where the short radius of channel
curvature induces the bank erosion. Through the use of the variable phase lag, the
model can simulate channel migration that is both across the floodplain and down
valley.
The new equation to predict the rate of bank erosion is also assumed to predict the
rate of channel migration. The bank erosion rate is computed as function of the unit
stream power and sediment transport capacity, the radius of channel curvature, and
the bank material properties acting to resist the erosion: vegetation, large woody
debris, cohesion, and armoring. The equation for bank erosion is based on
dimensional analysis of the controlling variables and the use of empirical coefficients.
The variable phase lag controls the amount of channel migration that is laterally
across the floodplain and the amount that is along the downstream valley axis. The
variable phase lag is determined by evaluating the river channel slope relative to the
valley slope and the limiting channel slope associated with the minimum unit stream
power. The phase lag is zero when the channel slope is equal to the valley slope,
which causes the meander bend to migrate entirely in the lateral direction. The extent
of lateral migration can be limited in the model by the specification of nonerodible
bank materials. The phase lag reaches its maximum value (onequarter of a meander
17
bend wave length) when the channel slope is equal to the slope associated with the
minimum unit stream power. When the phase lag is at the maximum value, the
channel migrates entirely in the downstreamvalley direction.
The model for this thesis is focusing on river channel migration and assumes that the
river is naturally capable of transporting the upstream sediment load through a reach
without deposition or erosion along the channel bed.
3.1 Riverbank Erosion Rate
The bank erosion rate is a function of the channel hydraulic properties acting to erode
the bank and the material properties of the river bank acting to resist the erosion. The
following hydraulic parameters are considered related to river bank erosion:
channel velocity (V),
energy slope (S),
hydraulic depth (D),
fluid viscosity (v),
bankfull channel width (Wb), and
meander bend radius of curvature (Rc).
The computation of the hydraulic capacity to transport sediment (using equations 3,4,
and 5) yields a dimensionless sediment concentration (Cs). The computation of the
sediment concentration already incorporates the mean channel velocity (V), energy
slope (S), hydraulic depth (D), median particle gain size (dso), fluid viscosity (v), and
particle fall velocity (w). In addition, the dimensionless ratio of bankfull channel
width over the meander bend radius of curvature (Wb /Rc) is used as an index to the
strength of secondary currents acting to cause bank erosion.
The material properties of the river bank acting to resist erosion include the presence
of riparian vegetation and large woody debris, cohesive properties of the bank
material, and the presence of coarse material that can armor the bank. The following
material properties of the river bank are considered to be indicators of erosion
resistance:
riparian vegetation
 fraction of bank area covered by vegetation roots (r7)
 vegetation root depth (r
 bank height (hb)
large woody debris
fraction of bank area covered by trees or large woody debris (LWD)
18
 average height of large woody debris j ams (dw)
bank cohesion
 plastic index (PI)
bank armoring
 portion of bank sediment too coarse for incipient motion (dc)
The bankfull channel width (Wb) and hydraulic depth (D) were used to normalize the
properties resisting bank erosion. Equation 17 is proposed to predict the rate of bank
erosion and channel migration.
B=\
a,Cs
UJ
at S
\r\j
+ al\LWD^\+aA(Pi) + al
>
dc \ >
\ c D)
at S+Lag^
a6V
(17)
where Be = rate of bank erosion [L/T],
Cs = bedmaterial sediment concentration [ppm]
Wb = bankfull channel width [L]
Rc channel radius of curvature [L]
S = distance along the channel
Lag = planform phase lag along the channel
rw = fraction of bank area covered by vegetation roots [%]
= vegetation root depth [L]
= bank height [L]
LWD = fraction of bank area covered by trees or large woody debris [%]
dw = average height of large woody debris jams [L]
= hydraulic depth of the channel [L]
= plastic index
= portion of bank sediment too coarse for incipient motion [%]
= mean channel velocity [L/T]
'7
n
hb
D
PI
dc
V
ci], ci2, CI3, CI4, as, and ab are empirical coefficients
All the parameters on the righthand side of equation 17 produce dimensionless terms
except for the average channel velocity, which is a primary variable controlling the
rate of bank erosion and is needed to provide dimensions to the bank erosion rate
[L/T]. The first set of terms in equation 17 [a 1 Cs (Wb / 7?c)] are computed at some
distance S along the channel. The following terms of equation 17 act to resist the
bank erosion and are computed at a phase lag distance (Lag) downstream from point
S. The velocity term of equation 17 (ab V) is computed at a distance S along the
19
channel. The computed rate of bank erosion from equation 17 is multiplied by a time
step (At) and then applied to the channel at a distance equal to S + Lag.
Only positive values of the bank erosion rate (Be) are considered, while negative
values are set to zero. The coefficients ai, a2, as, a4, and as must be determined
empirically through calibration. However, only some of these coefficients may be
important for a particular reach of river. For example, the coefficient as wouldnt be
important if there were no large woody debris or the coefficient a4 wouldnt be
important if there were no cohesive soils in the river bank. A river could also have a
reach with high cohesive banks with no protection from vegetation, large woody
debris, or coarse material. The calibration of coefficients is simplified for rivers
where the resistance to riverbank erosion, along a particular reach, is dominated by
only one type of material (vegetation, large woody debris, cohesion, or armoring).
The material properties of the river bank acting to resist erosion are discussed below.
The root structure of riparian vegetation can add cohesive properties to the river bank
and reduce the rates of bank erosion if the roots are as deep as the bank is high.
Otherwise, the bank can erode underneath the vegetation roots causing the top portion
of the bank to collapse. Trees that fall into the river from an eroding bank can add
large woody debris to the bank. Large woody debris present along the river bank
tends to increase the channel roughness, slow velocity, and also reduces the rates of
bank erosion. However, the woody debris may be floated away if the flow depths are
larger than the cumulative height of the large woody debris jam.
The presence of riparian vegetation can be expressed as the root density along the
local bank area (r7). The effectiveness of these roots can be expressed as the ratio of
the root length to the bank height (r^ Ihb). A ratio much less than unity would mean
that the roots are not deep enough to prevent erosion at the bank toe from
undercutting the vegetation. The presence of riparian forests and large woody debris
along the eroding bank can be expressed as a percentage of the local bank area
occupied by trees and large woody debris (L WD). The ability of large woody debris
to remain along the river bank can be expressed as the average ratio of the height of
large woody debris jams to the hydraulic depth of the river channel (dw ID). A ratio
much less than unity would indicate that the woody debris would be floated away.
The plasticity index (PI) (Das 1985) is used to characterize the cohesive properties of
the bank soil. The percentage of bank sediment that is too coarse for incipient motion
in the river channel (dc) is used to characterize the ability of coarse bank material to
armor the toe of an eroding river bank. The ratio of bank height to hydraulic depth
ihtJD) accounts for the greater mass of armoring material when a river bank is higher
than normal. For alluvial banks within the floodplain, this ratio is normally equal to
20
unity, but the ratio can be significantly greater than unity in the case of high terrace
banks.
3.2 Planform Phase Lag
The outward migration of a meander bend causes the river channel to migrate across
the floodplain. However, the planform phase lag between the curving flow path the
channel curvature also causes the meander bends to migrate downstream (see figure
4). If there were no phase lag, the meander bends would only migrate laterally across
the floodplain and increase in amplitude. If the phase lag were equal to the channel
distance along onefourth of a meander wavelength, then the meander bends would
only migrate downstream without any increase in amplitude.
Figure 4. Planform phase lag results in at least a portion of the meander bend
migrating downstream along the valley axis.
As stated in section 2.1, an alluvial river channel, given enough time, will adjust its
velocity, slope, roughness, and geometry in such a manner that a minimum amount of
unit stream power is used to convey the upstream discharge and transport the
21
upstream sediment load (Yang and Song 1979). This means that if the initial slope of
a river channel is much greater than the slope associated with the minimum unit
stream power (slope at minimum VS), then the meander bends will tend to migrate
laterally across the floodplain, the meander amplitude will increase, and the channel
slope will decrease. As the channel slope begins to approach the slope at minimum
VS, the migration of the meander bends will be primarily down valley. Since the
channel slope is limited by both the valley slope and the slope at minimum VS, the
planform phase lag should be proportional to the differences in these slopes.
Phase Lag =
te sj
(S, .S.ra)
Q /T
, 4 ,
(18)
where Sv = valley slope
Sc = channel slope
SminVS = channel slope at minimum VS
X = meander wavelength (X/4 is half the distance between river crossings)
0 = channel sinuosity
Cpl = a phaselag coefficient, normally equal to 1.0
Both the valley slope and channel slope can be computed from the floodplain
topography and the channel sinuosity. The slope at minimum VS is dependent on the
upstream river discharge (Q) and sediment concentration (Cs). The product of the
meander wavelength and sinuosity equal the channel length. Onefourth of this
product equals the channel length through onehalf of a meander bend. The phaselag
coefficient is constant with time, but can be set independently for each meander bend
to account for local conditions.
3.3 Channel Hydraulic Equations
The following hydraulic variables can be determined for the channel slope and cross
section by solving equations 19 through 25 and by minimizing the unit stream power
(minimum VS):
longitudinal channel slope (5)
average velocity (V)
wetted top width (W)
wetted perimeter (P)
maximum depth (Dmax)
cross sectional area (A)
22
hydraulic depth (D)
hydraulic radius (R)
The solution of these equations requires that the following variables be known or
assumed:
water discharge (Qw)
upstream sediment concentration (Cs)
median bedmaterial size (dso)
channel roughness (n)
acceleration of gravity (g)
sediment particle fall velocity (co)
water viscosity (v).
Mannings equation (Chow 1959) predicts the mean flow velocity as a function of the
hydraulic radius, longitudinal slope, and channel roughness:
V=RinSV2 (19)
n
where n = channel roughness coefficient
R = hydraulic radius
c = constant equal to 1.000 for SI units or 1.486 for English units
The hydraulic radius is computed as the ratio of the crosssectional flow area over the
wetted perimeter (equation 20). The hydraulic depth is similar to the hydraulic
radius, but is computed as the ratio of the crosssectional flow area over the wetted
top width (equation 21)
R = A P (20)
D = A W (21)
where A = cross sectional area
P = perimeter of the wetted channel
D = hydraulic depth
W= channel top width
Equations 22 and 23 are used to compute the crosssectional area and wetted
perimeter for a trapezoidal channel (see figure 5). These equations can also be
23
applied to a triangular cross section when the bottom width of the trapezoid is
assumed to be zero and the maximum depth is twice the hydraulic depth.
A = WDmax z(d)2
max \ max /
P= W 2zD + 2J(zD )2 + (d f
max \\ max / \ max /
(22)
(23)
Cross Section Geometries
Rectangular
Triangular
Trapezoid
W.S.
Figure 5. Possible cross section shapes considered in the model. The trapezoidal
shape is used for river crossings while the triangular shape is used for meander bends.
where Dmax maximum depth at a channel cross section
z = channel size slope (1 vertical to z horizontal)
The continuity equation relates the flow rate to the mean velocity and the cross
sectional flow area (equation 24).
Q = VA= VWD (24)
A sediment transport equation (such as equation 25) is used to match the hydraulic
capacity of the cross section geometry and channel slope to the upstream sediment
supply.
24
Cs = f(V,S,d50,U*,y,o))
(25)
where U* =
dso =
v =
GO
shear velocity = [g/?S]1/2
median bedmaterial particle size
kinematic viscosity of the fluid
sediment particle fall velocity
The steps to solving equations (19 to 25) are listed below:
1. Assume a range of wetted channel top widths (W,j and a width increment
{W^).
2. Assume a cross section shape (e.g., trapezoid or triangular) and bank side
slope (z).
3. For each channel width (Wi), assume an initial hydraulic depth (D) and
compute the cross sectional area (A) using equation 21.
4. Using an iterative procedure, solve equation 22 for the maximum water
depth (Dmax) so that the computed crosssectional area (A) matches the
area from step 3 within a specified tolerance.
5. Compute the wetted perimeter (P) of the channel bed and banks using
equation 23.
6. Compute the hydraulic radius (R) using equation 20.
7. Compute the average velocity (V) using the continuity equation (24).
8. Compute the longitudinal channel slope (S) assuming normal depth from
Mannings equation (19).
9. Compute the sediment transport capacity concentration (Cs) using
equations 3, 4, and 5 or some other transport function (equation 25).
10. Compare the computed sediment concentration (Cs) with the given
upstream supply concentration. If the computed concentration does not
match the upstream supply within a specified tolerance, then adjust the
assumed hydraulic depth (D) and go to step 3. If the computed
concentration is too low, then decrease the assumed hydraulic depth. If
the computed concentration is too high, then increase the assumed
hydraulic depth.
11. Compute the unit stream power (VS).
12. Increment the assumed channel width (Wi = Winc) and repeat steps 2
through 11.
13. For range of assumed channel widths, determine which width has the
minimum unit stream power (VS) and record the corresponding hydraulic
properties (V, S, Dmax, P, A, D, and R).
25
An example set of results from this computational procedure are shown in figure 6.
For each point along the curves shown in figure 6, the river discharge (Q), sediment
concentrationCs), channel roughness (n), viscosity (v), and the median sediment
particle diameter (dso) and fall velocity (u) are all constant.
Several patterns can be identified from this figure: As the channel width (W)
increases, the flow depth and velocity both decrease. The channel slope (S) initially
decreases, reaches a minimum and then steadily increases. The decreasing depth and
velocity and the increasing slope produce a minimum unit stream power (VS). The
channel width corresponding to the minimum VS is close to, but wider than, the width
corresponding to the minimum slope. This is because the unit stream power (VS) can
only minimize after the slope begins to increase with increasing channel width. The
hydraulic properties corresponding to the minimum VS are used in this model rather
than those corresponding to the minimum slope.
Minimum Unit Stream Power (Minimum VS)
Figure 6. Example computation results for the determining the minimum unit stream
power (VS).
26
3.4 Channel Curvature Equations
The channel thalweg alignment of the meandering river channel is initially defined by
the user. The points defining the channel alignment are smoothed before the first
time step and before subsequent time steps at a certain userspecified interval. The
smoothing procedure is applied to prevent small local irregularities in the channel
alignment from evolving, through the simulation of channel migration, into an erratic
channel alignment. Local irregularities along the channel bank frequently occur in
natural channels. These local irregularities produce local eddy currents along the
channel banks that affect the rate of energy loss. However, these local eddy currents
are not the cause of river meandering or channel migration (see Chapter 2).
Therefore, the local irregularities of the channel banks must be smoothed in the model
to prevent local and erratic channel migration.
The river channel alignment points are smoothed, one point at a time, using a second
order polynomial regression analysis (Devore 1991) of the nearby points. The
polynomial regression is used to help preserve the shape of meandering alignment.
Seven consecutive channelalignment points are used to smooth the one point in the
middle of that group (see figure 7). Three of the seven points are immediately
upstream, and three are immediately downstream, of the point to be smoothed. A
secondorder polynomial regression equation is produced from the coordinates of the
seven selected points. The regression equation is then used to compute the smoothed
horizontal coordinates of the middle point. The procedure is then repeated to smooth
the next downstream point. A new set of seven points are selected that surround the
next downstream point to be smoothed. Only one of the seven points is different
from the previous selection. The polynomial regression analysis repeated for the new
set of seven points and the regression equation is used to smooth the middle point.
The model keeps track of the changing alignment throughout the simulation. The
radius of curvature continually changes along the meander bend alignment. At the
beginning of each time step, the model computes the radius of curvature for each
point along the channel alignment. In addition, the model computes the azimuth of
the local direction of flow and the azimuth and origin for the local radius of curvature.
The changing radius of curvature along the meander bend is a key variable for
determining the rate of channel migration. The azimuth of the radius is important for
determining the direction of channel migration.
The local radius of curvature along a meander curve can be defined for any three
points (Xi,Yi), (X2, Y2), (A3, Y3) along the curve (see figure 8). The know values are the
coordinates of the three points along the meander curve: (Xj, Yi), (X2, Y2), and (X3,
7j). The unknown values are the radius of curvature (R) and coordinates for the
27
center of curvature (Xc, Yc). These three unknown values can be solved using the
following three equations:
314,500
5 314,000
.0
o>
c
Â£ 313,500
313,000
312,500
Hoh River near Forks, WA
Valley Axis
Initial Points
7 Regression Points
X Point to be Smoothed
Smoothed Points
.. I .... I .... I i .
825,500 826,000 826,500 827,000 827,500 828,000 828,500 829,000 829,500
x Easting (feet)
Figure 7. Example of smoothing the initial channel alignment points along the Hoh
River, one point at a time, using a regression analysis of nearby points. Seven
consecutive points are selected for regression analysis to smooth the middle point of
the selection.
I rj + l (26)
^l(X,X2f+(Y,Y1f (27)
M.Xj+fcY,)2 (28)
Equations 27 and 28 can be combined to produce the following equation.
(x, X, y +  ys Y = (x, x,)! + (r, r,)! (29)
Equation 29 can be rearranged to put the X and Y terms on different sides of the
equation.
28
(X, Xj(xc x2)2 = (r, y2)2(rc r,)! (30)
Radius of Curvature
Figure 8. Definition sketch for computing the radius of curvature along a meander
bend.
Equation 30 can be rearranged to combine the terms for Xc and Yc.
(2X,2X2)X, + X,2X2 =(2Y,2Y,)Yr +Y,1 Y,2 (31)
Equation 31 can be rearranged to solve for Yc.
(2X2 2X3 )X'+X12X22+Ys2Y22
(2Y,2Y,)
Equations 26 and 27 can be combined to produce the following equation.
(Xt X,f (.X, X, f = (n y,)! fc r2 f (33)
Equation 33 can be rearranged in a manner similar to equation 31.
29
(2X,2X2)X2 +X 2 X,2 =(2y22yX+y,2rf
(34)
Equation 34 can be rearranged to solve for Xc.
X.
hv (2X22X,)XC +x,2 x22 +r2 r2 ,y2 y2 + X 2 X 2
\zi2 zil) (2Y,2Y2) 1 jÂ£ 1 jÂ£ 2 1  yi 2
(2X,2X2)
(35)
Equation 35 can be used to solve for Xc by trial and error iterations. An initial value
is assumed for Xc and then used to compute Xc. The difference between the initial
value, Xc, and the computed value, Xc, is error term.
Error A A < tolerance
C C
(36)
Through successive iterations, the computed values of Ac can be substituted back into
the right side of equation 35 until the error term is within a certain tolerance.
The azimuth of the radius, 6, can be computed with the following equation.
9 = tan"1
ffcr.n
yxxc))
+ c
(37)
If (Y2 Yc) >0 and (X2 Xc) >0, (northeast quadrant), then c = 0
If (Y2 Yc) >0 and (X2 Xc) < 0, (northwest quadrant), then c = 7r
If (Y2 Yc) <0 and (X2Xc) < 0, (southwest quadrant), then c = it
If (Y2 Yc)< 0 and (X2 Xc) >0, (southeast quadrant), then c = 2ir
The tangent of the meander curve, $, at point 2 is perpendicular to the azimuth of the
radius, 6, and can be computed with one of the following two equations:
If 0 <6 <7r, then 0 = 6* 2 (38)
If 7T < 6 < 2ir, then 0 = 6' + 2 (39)
30
A solution to equation 35 will not exist if the coordinates of points 1,2, and 3
{{Xj, Yi), (X2, Y2), and (X3, Y3)} are in a straight line. In this case, the radius R is
infinite and set equal to 100 times the channel width, the tangent has the same
alignment as points 1, 2, and 3, and the azimuth of the radius is perpendicular to the
tangent at point 2.
3.5 Model Input Requirements
The meander model input requirements include both initial conditions and boundary
conditions. The initial conditions describe the river channel alignment and the
material properties (vegetation, large woody debris, cohesion, and armor material) of
the floodplain and terraces, including any lateral limits of channel migration. The
boundary conditions describe the upstream inputs of water and sediment.
3.5.1 Initial Conditions
The initial alignment of the river channel thalweg is digitally measured from rectified
aerial photographs or from channel surveys. The model requires that the digital
coordinates (northing and easting) of the channel thalweg alignment be specified in
the upstream to downstream order.
The river valley alignment and longitudinal slope are also measured from oroth
rectified aerial photography and topographic surveys.
A horizontal grid is used to describe the variation of material properties throughout
the river channel, floodplains, and terraces. The grid is rectangular and of uniform
size. Horizontal coordinates for the southwest most and northeast most comers of the
grid area are specified and the grid spacing is determined by dividing the total grid
length by a fixed number of grid nodes. An example grid is presented in figure 9.
The following material properties are defined for each node of the horizontal grid:
fraction of bank area covered by vegetation roots [%]
vegetation root depth [L]
surface elevation [L]
fraction of bank area covered by trees or large woody debris [%]
average height of large woody debris jams [L]
plastic index (PI) of bank material
portion of bank sediment too coarse for incipient motion [%].
31
Channel Migration After 200 Time Steps
Variable Phase Lag
Down Valley axis (feet)
Figure 9. Example grid used to describe the variation in material properties
throughout the river channel, floodplains, and terraces.
The material properties at each grid node are applied to the local area centered about
the grid node, which extends in four directions (each halfway to the next nearest grid
node).
The floodplain and terrace properties can be specified uniformly throughout the
horizontal grid to represent the normal condition. Different floodplain and terrace
properties can be specified for certain areas of the grid to locally override the normal
condition.
The rates of channel migration may be slowed or stopped by logjams, consolidated
terrace banks, bedrock outcrops, canyon walls, and bank protection structures. The
locations of these features can be digitally measured from orthorectified aerial
photographs, maps, or field surveys. The reduced erosion from these features can be
modeled by specifying resistant material properties in localized areas. Erosion can be
completely stopped in localized areas by specifying in the model that 100 percent of
the bank material is too coarse for incipient motion.
32
As the model determines that a portion of river bank is eroded, the model will assume
that riparian vegetation and roots along the eroded bank dies. However, trees that fall
into the river along the eroding bank may persist as large woody debris.
At the same time the channel is eroding the bank along the outside of the meander
curve, a similar volume of deposition is assumed to occur along the inside of the
meander curve. Therefore, the net volume or mass of sediment erosion is assumed to
be zero.
3.5.2 Boundary Conditions
The upstream boundary conditions include the discharge hydrograph (specified as a
time series of dates, hours, and river discharge) and sediment supply. The upstream
sediment supply can be specified as a time series of sediment concentrations or loads,
specified as an equation from a sedimentdischarge rating curve or computed by the
model from a specified cross section shape and channel slope.
3.6 General Model Algorithm
The general algorithm for the model is described below:
1. Read the coordinate data defining the alignment of the river channel thalweg.
2. Read the coordinates of the valley axis, the grid area, and read the material
properties of the river channel, floodplains, and terraces.
3. Read the discharge hydrograph data (date, time, and discharge) describing the
upstream boundary.
4. Read the sediment grain size and water temperature data.
5. Compute the fluid viscosity and sediment particle fall velocity from the grain
size and temperature data.
6. Begin the loop for each time step. If the river discharge is less than the user
defined threshold for channel migration, then skip to the next discharge.
7. Compute the upstream sediment supply (equations 3, 4, and 5) from the river
discharge and an assumed trapezoidal cross section with specified bottom
width, bank angle, and longitudinal slope.
8. Smooth the channel alignment coordinates to eliminate irregularities using a
secondorder polynomial regression applied to a moving group of seven
points. Smooth the channel coordinates at the beginning of the first time step
and before every time step at a specified interval.
33
9. For each point describing the river channel alignment, compute the radius of
curvature, azimuth of the radius, and the tangent angle of the meander curve
(equations 32 and 35 through 38).
10. Compute the direction of channel curvature (+1 = left turn; 1 = right turn).
11. Determine which channel alignment points cross the valley axis (river
crossings), compute the river channel slope between these crossings, and
compute the planform phase lag using equation (18).
12. Compute the rate of bank erosion and channel migration for each point along
the channel alignment (equation 17).
13. Apply the time step to the bank erosion rate and adjust the coordinates
describing the channel alignment.
14. Repeat steps 5 through 13 for each discharge of the hydro graph.
The detailed procedure for all computations is documented by the Fortran program
named: Channel Migration Model.for (see Appendix A).
34
4. Example Results
In order to demonstrate the model capabilities, the migration of an idealized channel
alignment is simulated with a constant discharge and uniform bank material
properties. The idealized channel alignment is from a sine generated curve computed
using equations 11 and 12. A plot of the initial channel alignment, along with future
channel migration alignments, is presented in figure 10. A constant discharge of
8,000 ft3/s was assumed for 200 time steps and uniform bank material properties were
specified throughout the hypothetical floodplain (see table 2). The upstream sediment
supply was computed using the sediment transport capacity equation by Yang (1984)
and an assumed trapezoidal cross section with a bottom width of 140 feet, a bank
angle of 30 degrees, and a longitudinal slope of 0.0005. The minimum unit stream
power and the corresponding hydraulic properties were computed for the upstream
discharge and sediment supply (see figure 6).
The valley slope was set equal to 0.000967. With an initial channel sinuosity of 1.10,
the initial channel slope was 0.000883 and the slope associated with the minimum
unit stream power was 0.00486. From equation 18, the difference between the valley
slope and the initial channel slope is 17 percent of the difference between the valley
slope and the slope associated with the minimum unit stream power. The relatively
low sinuosity results in an initially small planformphase lag, which causes most of
the initial channel migration to be in the lateral direction (across the floodplain). As
the channel migration continues, the sinuosity increases and the channel slope
decreases. The decrease in channel slope causes the planformphase lag to increase,
which results in more of the channel migrating is in the downvalley direction. The
results from this example model simulation are presented in figure 10, with a closeup
view presented in figure 11.
35
Table 2. Example model input data for the hypothetical problem.
Model Input Parameters for Equation 17 Variable Value
Empirical Weighting Coefficients:
hydraulic coefficient ai 30.
vegetation coefficient a2 100.
large woody debris coefficient 100.
cohesion coefficient a4 2.7778
armoring coefficient as 100.
velocity scaling coefficient a<5 0.002
Vegetation Properties Resisting Bank Erosion
vegetation root density rs 0.3
vegetation root depth [L] rd 6.0
bank height [L] hb 10
Large Woody Debris Resisting Bank Erosion
fraction of trees or large woody debris along the LWD 0.2
average height of large woody debris jams [L] dw 2.0
Cohesion Resisting Bank Erosion
plastic index PI 6.0
Armoring Resisting Bank Erosion fraction of bank material coarser than incipient
motion dc 0.10
36
Channel Migration After 200 Time Steps
Variable Phase Lag
Figure 10. Example model results assuming an initial channel alignment from a sine
generated curve. The flow is from left to right.
Channel Migration After 200 Time Steps
Variable Phase Lag
Initial Conditions Time Step 1 Time Step 20 Time Step 40
Time Step 60 Time Step 80 Time Step 100 Time Step 120
Time Step 140 Time Step 160 Time Step 180 Time Step 200
O
w
o
1,200
2,000 2,400 2,800 3,200 3,600 4,000 4,400 4,800 5,200 5,600
Down Valley axis (feet)
Figure 11. Closeup view of example model results shown in figure 10.
37
5. Model Calibration and Validation
Measurements of historical channel migration should be used to calibrate and then
validate the model results before the model is used for future predictions. Typically,
channel surveys and rectified aerial photographs provide the best data from which to
measure the historic channel alignments. Channel alignment data from historic maps
can also be used if the maps are accurate and can be rectified.
The historic measurements of channel alignment should be divided into two
independent time periods: one for model calibration and one for model validation.
Model input parameters should be calibrated so that the predicted channel alignment
best matches the measured channel alignment. The calibrated input parameters
should then be used to simulate a different historical time period to validate the model
results. The difference between the predicted and measured channel alignments can
then be used to describe the accuracy of the model predictions.
The data necessary to demonstrate the calibration and validation of the model are
available for the Hoh River near Forks, Washington. The U.S. Department of the
Interior, Bureau of Reclamation (Piety et al. 2004) is completing a geomorphic
assessment of a 23mile reach of this gravelbed river from river miles 17 to 40. The
channel migration model was applied to a 4.2mile reach of the Hoh River between
river miles 21.4 and 25.6. This reach is known as the Morgans Crossing Reach,
which is approximately 5.5 miles downstream from the boundary of Olympic
National Park.
Historic aerial photographs of the Hoh River are available about every decade for the
following years:
1939 1971 1994
1950 1977 2001
1960 1981 2002
The mosaic of aerial photographs taken in 2001 was orthorectified by the
Washington Department of Natural Resources (Piety et al. 2004). The aerial
photographs that were taken in the other years were registered to the same scale and
coordinate system as the 2001 orthorectified aerial photographs (1983 North
American Datum, Washington State Plane coordinates, North Zone). A geographic
information system (GIS) data base was then created that included all of these aerial
photographs (Piety et al. 2004).
38
The 1939 aerial photographs were taken at a higher altitude and lack the detail of the
later aerial photographs. The river alignment in 1939 was fairly straight and there
were only small adjustments in the river alignment by 1950. Therefore, 1950 was
selected and the starting point for model simulation. The initial channel alignment
was digitized from the 1950 aerial photograph (see figure 12).
Figure 12. Initial conditions for model simulation were digitized from the 1950 aerial
photograph.
The 52year period from 1950 to 2002 was used for model calibration and validation.
The 27year time period, from 1950 to 1977, was selected for model calibration. The
25year time period, from 1977 to 2002, was selected for model validation.
The U.S. Geological Survey has produced discharge records for the Hoh River at two
different locations (see table 3). The stream gage for the first time period (gage
number 12041000) is just upstream from the Morgans Crossing Reach while the
stream gage for the second time period (gage number 12041200) is about 6 miles
downstream. Since there is a fouryear overlap in the discharge records, a correlation
was developed between the two streamgage locations for meandaily and annual
peak discharge (see figures 13 and 14). Based on these correlations, the meandaily
discharge at the second location is 24 percent larger than at the first location and the
39
annualpeak discharge is 27 percent larger. Therefore, the regression equations
shown in figures 13 and 14 were used to estimate the meandaily and annualpeak
discharges at the first steamgage location (gage number 12041000) for water years
1965 to 2002 (see figure 15).
Table 3. Available streamflow gage records for the Hoh River .
Discharge Type Stream Gage Location Gage Number Available Water Years
Meandaily near Forks, WA 12041000 1927 to 1964
Meandaily at US Highway 101 near Forks, WA 12041200 1961 to 2002
Annual peak near Forks, WA 12041000 1927 to 1964
Annual peak at US Highway 101 near Forks, WA 12041200 1961 to 2002
Correlator! in Meandaily Discharge
O)
.c
o
5
_>*
*5
TJ
C
m
o
Qn=1.24 Q0128
R2 = 0.978
Meandaily Discharge, Q0 (cfs) at original upstream gage location
Figure 13. Correlation in meandaily discharge between two Hoh River streamgage
locations for the period 1961 to 1964.
40
Correlation in Peak Discharge
(A
o
5
a>
o>
w
(C
Â£
O
(A
b
C3
0)
o.
Qn=1.27 Qo+1009
R2 = 0.846
Peak Discharge, Q0 (cfs) at original upstream gage location
Figure 14. Correlation in annualpeak discharge between two Hoh River streamgage
locations for the period 1961 to 1964.
50.000
40.000
Hoh River Discharge near Forks, WA
Annual Peak Discharge Meandaily Discharge
:i   _ i 1 J
1950 1955 1960 1964 1970 1975 1980 1984 1990 1995 2000
Figure 15. Meandaily and annualpeak discharge of the Hoh River.
The annualpeak discharge is higher than the meandaily discharge. An hourly
discharge hydrograph would be more accurate for model simulations than a daily
41
discharge hydrograph. However, the data for an hourly discharge hydrograph is not
readily available for a 52year period. Therefore, the meandaily discharge
hydrograph was used in the model to represent the flows through the Morgans
Crossing Reach.
Since there are no measurements of sediment load, the upstream sediment supply to
the Morgans Crossing Reach of the Hoh River was computed based on the hydraulic
capacity of the river. A trapezoidal channel was specified with a bank slope of 30
degrees, a bottom width of 400 feet, and a longitudinal slope of 0.0026. The
Mannings roughness coefficient (see equation 19) was specified as a function of the
hydraulic depth (see figure 16). The computed sediment load and concentration for
this trapezoidal cross section is presented in figure 17.
Hoh River Channel Roughness Versus Depth
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
Hydraulic Depth, D (feet)
Figure 16. The Manning roughness coefficient, n, was specified as a power function
of the hydraulic depth, D.
The channel migration model was calibrated by adjusting empirical weighting
coefficients (see equation 17) and the phaselag coefficients for each meander bend
(see equation 18). The bank material properties were globally defined for the
modeled area and locally defined for seven smaller areas of the channel, floodplain,
and terraces (see figure 18). These seven areas were identified based on field
observations and inspection of the aerial photographs. The material properties for
vegetation, large woody debris, cohesion, and armoring (see equation 17 and table 4)
were determined for each area based on field inspection and were not subsequently
adjusted during the calibration process. These parameters control the rate and extent
42
of bank erosion and channel migration. The local values override the global values
and the local parameters that are specified last, override previously specified values.
Hoh River Sediment Supply
T3
ra
o
c
Q>
E
a>
V)
c
0)
E
T3
0)
C
,o
5 E
? o
0) Q_
O ^
c
o
o
Discharge (cfs)
Figure 17. Computed sediment load and concentration for the upstream supply to the
Morgan's Crossing Reach of the Hoh River.
Hoh River Model Grid
318,000
310,000 !'''l'''''''111*
818,000 822,000 826,000 830,000 834,000
* 1950 River Alignment
Debris Fan
Floodplain Area 1
O Floodplain Area 2
B Terrace Bank
Dense Woods
County Road Armoring
Glacial Bluff
OGrid Boundary
838,000
Easting (feet)
Figure 18. Areas of different material properties for the Hoh River channel,
floodplains, and terraces along the Morgan's Crossing Reach. Flow is in the westerly
direction.
43
The weighting coefficients (a2 through as) of equation 17 were set to equally weight
the parameters resisting bank erosion. Coefficients a2, as, and as, were set equal to
100. The coefficient a4 was set equal to the ratio 100/60 (a4 1.67) because a
maximum plasticity index (PI) of 60 was specified.
The calibration strategy was to first assume an initial value (0.002) for the coefficient
as The coefficient, a2, was then adjusted until the predicted increase in meander
amplitude was close to the actual increase measured from the aerial photographs. The
coefficient a&, was then adjusted until the predicted increase in meander amplitude
best matched the actual increase. The planform phaselag coefficients of equation 18
were initially set equal to 1.0. However, the phaselag coefficient for the upstream
most meander bend was reduced to 0.7 to prevent migration through a very cohesive
bank and debris fan. The final calibration parameters are listed in table 4.
Table 4. Model input data for equation 17, calibrated to the 4.2mile reach of the Hoh River near Forks, Washington.
Input data Global para meters Local Parameters
Debris fan & canyon Flood plain area 1 Flood plain area 2 Terrace bank Dense woods County Road rip rap Glacial bluff
Empirical Weigh ting Coefl icients (only applied globa ly)
ai 6.0
a2 100.
100.
a4 1.67
as 100.
ae 0.00194
Vegetation Properties Resisting Bank Erosion
re 0.5 0.3 0.4 0.4 0.5 0.7 0.2 0.2
rd 8.0 6. 6. 6. 14. 14. 6. 6.
h 12.0 20. 8. 8. 12. 12. 20. 100.
Large Woody De :>ris Resisting Ban c Erosion
LWD 0.4 0.3 0.2 0.2 0.4 0.8 0.2 0.2
dw 8.0 5. 5. 5. 8. 8. 5. 5,
Cohesion Resisting Bank Erosion
PI 15. 60. 5. 5. 20. 10. 36. 36.
Armoring Resisting Bank Erosion
dc 0.3 1.0 0.2 0.2 0.4 0.4 1.0 0.3
44
Figure 19 shows a series of predicted channel alignments over the period 1950 to
2002. The predicted alignment of the Hoh River in 1960, 1971, and 1977 is plotted in
figures 20, 21, and 22. For each figure, the predicted river alignment is plotted on top
of the aerial photographs for visual comparison. The model was not able to predict
the cutoff of the third meander bend that occurred sometime between 1971 and 1977
(see figure 22). Therefore, the predicted channel alignment was calibrated for the
first two meander bends. The combination of coefficients a\ and and the phaselag
coefficient that provided the best match between the predicted and measured channel
alignments was determined by visual inspection of the results (plotted on top of the
aerial photographs).
Hoh River near Forks, WA
320,000
ai
c
!c
t
o
z
a, = 6.0 i
a6 = 0.001941
Â£ 316,000
312,000
308,000
Data Input
8/15/1950
10/7/1960
10/22/1971
10/25/1977
10/19/1985
8/22/2001
Valley Axis
Initial Conditions
8/1/1955
10/5/1965
11/12/1974
9/28/1981
10/26/1994
6/29/2002
818,000 822,000 826,000 830,000 834,000 838,000
Easting (feet)
Figure 19. Predicted channel alignments for the period 1950 to 2002.
Once the model calibration process was complete, the model simulation of the
historic calibration period (1950 to 1977) was continued to 2002 for the validation
period. Figures 23 through 26 show the predicted Hoh River alignment for a series of
years during the validation period from 1981 through 2002. The model could have
been restarted for the validation period using the measured channel alignment in
1977. However, the model simulation was continued for the entire period (1950 to
2002) to see if the model would eventually predict a cutoff of the third meander bend.
45
In fact, the model results do indicate that the channel may cutoff the third meander
bend by 1994 and at least by 2001 (see figures 24 and 25). In addition, the predicted
channel alignments in 2001 and 2002 closely matched the measured alignments
through the first two meander bends.
Figure 20. Initial (1950) and predicted (1960) channel alignments plotted on the 1960
aerial photograph.
46
Figure 21. Initial (1950) and predicted (1971) channel alignments plotted on the 1971
aerial photograph.
Figure 22. Initial (1950) and predicted (1977) channel alignments plotted on the 1977
aerial photograph. A meander cutoff occurred (circled area) sometime between 1971
and 1977. This meander cutoff was not predicted by the model.
47
1981
Figure 23. Initial (1950) and predicted (1981) channel alignments plotted on the 1981
aerial photograph.
Figure 24. Initial (1950) and predicted (1994) channel alignments plotted on the 1994
aerial photograph. The sharp channel bend (near 80 degrees) predicted by the model
could be interpreted as a meander bend cutoff.
48
Figure 25. Initial (1950) and predicted (2001) channel alignments plotted on the 2001
aerial photograph. The sudden change in predicted channel alignment indicates a
meander cutoff.
Figure 26. Initial (1950) and predicted (2002) channel alignments plotted on the 2002
aerial photograph.
49
6. Discussion
In spite of all the parameters of equations 17 and 18, only three coefficients (fly, fl
and Cpi) required calibration. Except of the upstream most meander bend, the phase
lag coefficients, Cpl, could be left equal to unity. For the calibration period, the
model was able to correctly predict the direction and approximate magnitude of
channel migration for a least 10 years, until the cutoff of meander bends that first
occurred between 1960 and 1971 (see figures 20 and 21) and again between 1971 and
1977 (see figures 21 and 22). Although, the model was not able to predict the
meander bend cutoff during the calibration period, a continuation of the simulation
did eventually predict that a meander bend cutoff would occur in the downstream half
of the reach.
The model correctly predicted that initial direction of channel migration was in the
lateral direction, across the floodplain. As the meanderbend amplitude increased,
more of the meanderbend migration was in the downvalley direction (see figure 19).
For a river channel with a longitudinal slope that is much steeper than the slope
associated with the minimum unit stream power, the channel migration model does
predict that the initial channel migration will initially tend to be in the lateral direction
across the floodplain. As the longitudinal channel slope approaches the minimum
slope, the model predicts that the more of the channel migration is in the direction of
the downstream valley axis.
The successes of the model simulation would tend to confirm the three hypotheses
represented by equations 17 and 18. The radius of channel curvature, along with the
sediment transport capacity, does appear to correlate with the magnitude of channel
migration. However, the planform phase lag is needed to determine the direction of
channel migration.
Although the model does not yet have explicit criteria to predict channel avulsions,
the occurrence of channel avulsions can be implied by the prediction of a sharp curve
of sudden change of in direction of the channel alignment.
50
7. Conclusions
The channel migration model described in this thesis does show promise as a tool for
predicting future channel alignments. The model could be used to compare different
hydrologic scenarios and alternative land use plans. The model could also be used to
evaluate proposed highway and bridge alignments in the vicinity of the river channel
and floodplain.
The model could continue to be improved with additional research. The model
should be tested and validated on other rivers to determine if it can work for other
cases. By testing the model on a variety of rivers, it may be possible to determine
coefficient values that are applicable to a wide range of rivers. Additional research is
needed to develop and implement criteria for the occurrence of channel avulsions and
meanderbend cutoffs.
Future versions of the model could track the age of floodplain and terrace vegetation
over time. As the model simulates the creation of point bars along the inside curves
of meander bends, the model could also simulate the growth of vegetation on these
point bars. Young vegetation could either be scoured by subsequent floods or grow
to maturity and potentially protect the bank from future erosion. A future version of
the model could also into account the increase in sediment supply caused by the
erosion of terrace banks because they are higher in elevation than the floodplain
surface.
51
Appendix
A. Channel Migration Model: Fortran Program Listing
PROGRAM Channel_Migration_Model
C
C CHANNEL MIGRATION MODEL FOR MEANDERING RIVERS
C
C This program is designed to predict the future alignment of a
C meandering river channel.
C
C Timothy J. Randle
C University of Colorado, Denver Campus
C Masters Thesis, Summer 2004
C
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc c
c n List of subroutines in the order they appera in the program
c SUBROUTINE READGRID
c SUBROUTINE READXY(NUMPTS)
c SUBROUTINE SMOOTH(NUMPTS)
c SUBROUTINE SMOOTH_Bank_Ero s ion(NUMPTS)
c SUBROUTINE READFLOW(NUMQ)
c ri SUBROUTINE JULIAN(MM,DD,YYYY,TIME,JDATE)
L. c Planform geometry, bank erosion, and channel migration
c ri routines:
c SUBROUTINE RADIUS(NUMPTS,IPRINT)
c SUBROUTINE BANK_EROSION(NUMPTS,dt,Q,Lag,IPRINT)
c SUBROUTINE CURVE_TURN(NUMPTS)
c SUBROUTINE THALWEG_LOCATON(NUMPTS,I,Tloc,WIDTH,IPRINT)
c SUBROUTINE SLOPE_CAL(NUMPTS,IPRINT)
c ri SUBROUTINE MIGRATION(NUMPTS)
c n Sediment transport routines:
c SUBROUTINE UPSEDSUP(Q)
c SUBROUTINE MINVS
c SUBROUTINE SolveW(Tloc,WIDTH2,Hdepth2,Dmax,VELAV2)
c SUBROUTINE READSED
c SUBROUTINE SDTR(PX,PXN,NOBS)
c SUBROUTINE DIAM(PXN,SIZE,DSIZE,PCTN,NRGS)
c SUBROUTINE HEADING
c SUBROUTINE YANG(WIDTH)
c SUBROUTINE SETVEL(TDF)
52
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
c
c
Cslope(20),
, VSmin, Smin, Vmin,
DOUBLE PRECISION X, Y, CHANLNG, Xc, Yc, Rc, THETA, PHI,
+ B, SM, PIE, dSO, dSl, dS2, X2, Y2, Xvl, Yvl
COMMON / PLAN_DATA / X(1000), Y(1000), CHANLNG(1000),
+ Xc(1000), Yc(1000), Rc(1000), THETA(IOOO), PHI(1000),
+ SEIN(IOOO), SEIN2(1000), XLT(IOOO), YLT(IOOO),
XRT(1000),
+ YRT(IOOO), XMIN, YMIN
INTEGER SEIN, SEIN2
COMMON / XY2 / X2(1000), Y2(1000)
COMMON / SLOPES / Lmax, NC, CROSSING(20),
+ PLagFactor(20), PhaseLag(20)
COMMON / MINIMUMS / AMPmax(20), MaxAMP(20)
+ DmaxM, SHAPE
COMMON / HYDROLOGY / QTHRESHOLD, JDAY(20000), FLOW(20000)
COMMON / EROSION / CsWRc(lOOO), ECOEF(1000,50) BE(1000)
COMMON / WIDTHS / CHwidth(1000), CHwidthMax(1000),
CHWtime(1000),
+ WIDTHM, WIDTHL, WIDTHU, WIDTHI
COMMON / HYDRAULICS / VEL(1000), Dh(1000), Dmax(lOOO),
+ TlocL(1000), TlocR(1000)
COMMON / BESmooth / BE2(1000)
COMMON / HYDR / DISCH, AREA, Hdepth, HdepthM, SLOPE, TDF,
NVALUE
COMMON / MANNING / aMan, bMan, Rlow, Rhigh
COMMON / HYDRC / D35, D50, D65, D90, H_Radius, VELAV
COMMON / SDFLOW / QBI(15) QSI(15), QTI(15) ,
+ SETV(15), BEDLD, SUSLD, TOTLD, FBP(15), FSP(15)
COMMON / SEDDAT / DRL(15), DRU(15), GMD(15), SEDSUP, Cs, NRGS
COMMON / SIZER / DR(15), PPB(15), PSS(15)
COMMON / SNDREG / AX, BX
INTEGER CROSSING, DD, YYYY, TIME
REAL Lag, PLagFactor, LagRatio, JDAY, JDATE, NVALUE
C
COMMON / PTIMES / IPT, IPMM(20), IPDD(20), IPYYYY(20)
COMMON / OPTION / SWITCH(12), NTITLE, KPRINT(200), NPRINT(7),
+ LPRINT(4), LTSWTC, LCSWTC
INTEGER RIGHT, SWITCH
COMMON / CVAR / METHOD
COMMON / NAMES / TITLEXY(5) TITLEQ(5) TITLEG(5)
CHARACTER TITLEXY*80, TITLEQ*8 0 TITLEG*80
COMMON / NAMES2 / TITLE(5), HYDNME, STGNME, TMEINT, UNITS,
+ MILES, CSORDR
CHARACTER LABEL*2, TITLE*70, HYDNME*60, STGNME*60, STATION*70,
+TMEINT*10, UNITS*6, MILES*3, CSORDR*7, METHOD*2, XSHAPE*6,
ANS*1
53
o n
COMMON /
COMMON /
dyGrid,
+
+
REAL LWD
COMMON /
COMMON /
+
COMMON /
CONSTANTS / PIE
GRID / XLL, YLL, XUR, YUR, IMAXg, JMAXg, dxGrid,
rg(101,51), rd(101,51), hb(101,51),
LWD(101,51), dw(101,51), PI(101,51), dc(101,51)
BEcoef / al, a2, a3, a4, a5, a6
V_Axis / NumVpts, Vx(20), Vy(20), ALPHAv(20),
vdist(20), Vslope(20), CLAG(20)
UPCHAN / UPSLOPE, BANK_ANG, BottomW, QU1, QU2, IQU
c
c Open the output files
c 6 Meander.out
c 7 Bank Erosion.out
c 8 SEDVS .OUT
c 11 Debugger.out
c
Predicted channel coordinate data.
Predicted band erosion rates.
Minimum unit stream power computations.
Time steps and error messages.
OPEN(6,FILE='Meander.out',STATUS='UNKNOWN')
OPEN(7,FILE='Bank Erosion.out',STATUS='UNKNOWN')
OPEN(8,FILE=1SEDVS.OUT',STATUS='UNKNOWN')
OPEN(9,FILE=1SEDMINVS.OUT',STATUS='UNKNOWN')
OPEN(10,FILE='Smooth.out',STATUS='UNKNOWN')
OPEN(11,FILE='Debugger.out',STATUS='UNKNOWN')
WRITE(*,1)
1 FORMAT(/,'Welcome to the channel migration model for ',
+'meandering river',/,'channels.')
C
C X = East coordinate or downstream valley axis.
C Y = North coordinate or cross valley axis.
C CHANLNG = Cumulative channel length, increasing in the
downstream
C
C Xc =
C Yc =
C Rc =
C THETA
C PHI =
C
C Call <
C CALL ]
C
C Call i
c CALL ]
c
direction.
Xcoordinate of the point of curvature.
Ycoordinate of the point of curvature.
Radius of curvature.
Azimuth angle, in radians, of the radius.
Tangent angle, in radians, of the meander curve
Call a routnie to read the discharge data.
CALL READFLOW(NUMQ)
54
c
C Call a routine to read the sediment data.
C
CALL READSED
C
C Call a routine to compute the sediment particle settling
C velocities.
C
CALL SETVEL(TDF)
C
WRITE(6,10)
10 FORMAT('123459(1 234567890 '), 2345678901211 ( '
234567890'),
+/,5X,' Initial Conditions')
C
WRITE(8,12)
12 FORMAT(/,' Discharge Sed Load Sed Sup Ratio
+'Converge Hdepth Hdepth2 Width Ave Vel ',
+'Slope UnitVS')
C
WRITE(9,14)
14 FORMAT(' Discharge Sed Load Sed Cone Min VS Slope',
+' Velocity Width Hyd Depth Max Depth Thai Pos')
C
C Determine the printing interval
C NPNT = Number of time steps to print
C INTP = Time step printing interval
c
NPNT =10
INTP = NUMQ/NPNT
IPRINT = 0
ITP = 1
C
C Begin a loop for each time step.
C
IPRINT = 0
JSET = 8
JPRINT = 0
DO 100 N = 1,NUMQ
C
C Determine the calander date.
C
JDATE = JDAY(N)
CALL CDATE(MM,DD,YYYY,TIME,JDATE)
C
C Determine if the print switch should be turned on.
C
IF (MM .EQ. IPMM(ITP) .AND. DD .EQ. IPDD(ITP) .AND.
+ YYYY .EQ. IPYYYY(ITP)) THEN
IPRINT = 1
55
n o
ITP = ITP + 1
IF (ITP .GT. IPT) ITP = IPT
ENDIF
Assign the river discharge (Q).
Q = FLOW (N)
IF (N .EQ. 200) THEN
XXX = 1.
rt ENDIF
L. c If the flow is between the user defined limits,
c c routine to compute the upstream sediment supply
IF (Q .GE. QU1 .AND. Q .LE. QU2) THEN
C
C
C
C
C
c
c
20
C
C
C
CALL UPSEDSUP(Q)
WRITE(8,*) ' '
ELSE
GO TO 100
ENDIF
IF (IQU .EQ. 1) GO TO 100
Compute the time step as the difference between consecutive
dates and times.
JDAY is the julian date, including the time of day
expressed as a decimal fraction of the day.
IF (N .EQ. 1) THEN
dt = 1.
ELSE
dt = JDAY(N)  JDAY(Nl)
ENDIF
IF (dt .LT. 0.) THEN
WRITE(*,*) 'Negative time step:',dt
PAUSE ''
ENDIF WRITE(*,20) N, MM, DD, YYYY, TIME, dt, Q
WRITE(7,20) N, MM, DD, YYYY, TIME, dt, Q
WRITE(10,20) N, MM, DD, YYYY, TIME, dt, Q
WRITE(11,20) N, MM, DD, YYYY, TIME, dt. Q
FORMAT(15,2X ,12 , '/' ,12, '/',14 ,2X,14 , dt :
1 days, Q = 1,F10.1 ,' cfs')
Determine if the print switch should be turned on.
JPRINT = JPRINT + 1
C IF (JPRINT .EQ. 1 .OR. JPRINT/INTP*INTP .EQ. JPRINT) THEN
C IPRINT = 1
C ELSE
56
o n n n n
C IPRINT = 0
C ENDIF
C
C Call a routine to smooth the xy data, which define the
C channel thalweg, at the first time step and at a certain
C time step interval after that.
C
IF (JPRINT .LE. 2 .OR. JPRINT/JSET*JSET .EQ. JPRINT) THEN
CALL SMOOTH(NUMPTS)
ENDIF
C
C Call a routine to compute the radius of curvature, azimuth
C of the radius, and the tangent angle of the meander curve.
C
IF (JPRINT .EQ. 1) THEN
IPNT = 1
ELSE
IPNT = 0
ENDIF
CALL RADIUS(NUMPTS,IPNT)
Call a routine to determine the direction of channel
L1.L V CL L. LLI. w f
left
(SEIN = +1) or right (SEIN
1)
CALL CURVE_TURN(NUMPTS)
C
C Call a routine to compute channel slope.
C
CALL SLOPE_CAL(NUMPTS,IPRINT)
C
C Call a routine to compute the bank erosion rate.
C
CALL BANK_EROSION(NUMPTS,dt,Q,Lag,IPRINT)
C
C Call a routine to smooth the amount of bank erosion.
C
CALL SMOOTH_BANK_EROSION(NUMPTS)
C
C Print the time step and flow rate.
C
IF (IPRINT .EQ. 1) THEN
WRITE(6,30) N, MM, DD, YYYY, TIME, dt, Q
30 FORMAT(15,12,'/',12,'/',14,2X,14,4X,F10.3,
+ F10.1,' Cfs')
ENDIF
C
C Call a routine to compute the new coordinates of the
C migrating river channel.
C
57
CALL MIGRATION(NUMPTS,IPRINT)
IF (IPRINT EQ. 1) IPRINT = 0
100 CONTINUE
C
C Call a routine to smooth the ending xy data.
C
CALL SMOOTH(NUMPTS)
CLOSE (6)
CLOSE (8)
CLOSE(9)
CLOSE(10)
CLOSE(11)
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE READGRID
C
C Routine to read the grid data.
C
DOUBLE PRECISION X, Y, CHANLNG, Xc, Yc, Rc, THETA, PHI, PIE
COMMON / PLAN_DATA / X(1000), Y(1000), CHANLNG(1000),
+ Xc(1000), Yc(1000), Rc(1000), THETA(IOOO), PHI(1000),
+ SEIN(1000), SEIN2(1000), XLT(1000), YLT(1000),
XRT(1000),
+ YRT(1000), XMIN, YMIN
INTEGER SEIN, SEIN2
COMMON / GRID / XLL, YLL, XUR, YUR, IMAXg, JMAXg, dxGrid,
dyGrid,
+ rg(101,51), rd(101,51), hb(101,51),
+ LWD(101,51), dw(101,51), PI(101,51), dc(101,51)
REAL LWD
COMMON / BEcoef / al, a2, a3, a4, a5, a6
COMMON / V_Axis / NumVptS, Vx(20), Vy(20), ALPHAv(20),
+ Vdist(20), Vslope(20), CLAG(20)
COMMON / PTIMES / IPT, IPMM(20), IPDD(20), IPYYYY(20)
COMMON / NAMES / TITLEXY(5), TITLEQ(5), TITLEG(5)
CHARACTER TITLEXY*80, TITLEQ*80, TITLEG*80
CHARACTER LABEL*3, UNITS*3, FILEIN*40
COMMON / UPCHAN / UPSLOPE, BANK_ANG, BottomW, QU1, QU2, IQU
COMMON / MANNING / aMan, bMan, Rlow, Rhigh
COMMON / CONSTANTS / PIE
DATA IMAXg / 101 /
DATA JMAXg / 51 /
DATA PIE / 3.14159265359 /
DATA IQU / 0 /
C
C Initialize variables.
C
NTITLE = 0
NumVpts = 0
58
nonnnn no
nGlobal = 0
nLocal = 0
IPT = 0
NC = 0
C
C Open the input file.
C
C FILEIN = 'Grid.txt
FILEIN = 'Hoh Grid.txt
OPEN(5,FILE=FILEIN,STATUS='OLD')
OPEN(30,FILE=1 GRID.OUT 1 ,STATUS='UNKNOWN')
C
Initialize the grid values.
DO 20 N = 1,20
CLAG(N) =1.0
20 CONTINUE
DO 40 J = 1,JMAXg
DO 30 I = l,IMAXg
rg(I,J) = 0.
rd (I, J) = 0
hb (I, J) = 0
LWD (I, J) = 0
dw (I, J) = 0
PI(I,J) = 0
dc(I, J) = 0
30 CONTINUE
40 CONTINUE
50 READ(5,60,END=450) LABEL
60 FORMAT(A)
Skip comment lines.
IF (LABEL(1:2) .EQ. 'C .OR. LABEL(1:2) .EQ. 'C ') GO TO 50
Read the title lines.
IF (LABEL(1:2) .EQ. 'TI' .OR. LABEL(1:2) .EQ. 1ti' OR.
+ LABEL(1:2) .EQ. 'Ti') THEN
NTITLE = NTITLE + 1
IF (NTITLE .GT. 5) GO TO 50
BACKSPACE 5
READ (5,70) TITLEG(NTITLE)
70 FORMAT(10X,A)
WRITE(*,80) TITLEG(NTITLE)
80 FORMAT(5X,A)
GO TO 50
ENDIF
59
no no
Read the system of units, either English or S.I.
IP (LABEL(1:2) .EQ. 'UN' .OR. LABEL(1:2) .EQ. 'un1 .OR.
+ LABEL(1:2) .EQ. 'Un') THEN
BACKSPACE 5
READ(5,90) UNITS
90 FORMAT(1OX,A)
IF (UNITS .EQ. 'SI ' .OR. UNITS .EQ. ' SI' .OR.
+ UNITS .EQ. 'si ' .OR. UNITS .EQ. ' si') THEN
UNITS = 'SI '
nu = 1.3E06
g = 9.81
b = 1.
ELSE
UNITS = 'ENG'
nu = 1.4E05
g = 32.2
b = 1.486
ENDIF
WRITE(*,*) 'Units = ',UNITS
WRITE(*,*) 'nu = 1 1 nu, g = ''/ A II Â£1
WRITE(7,*) 'Units = ',UNITS
WRITE(7,*) 'nu = 1 ,nu,', g = b = b
GO TO 50
ENDIF
C
Read the upstream channel dimensions.
IF (LABEL(1:2) .EQ. 'UP' .OR. LABEL(1:2) .EQ. 'Up' .OR.
+ LABEL(1:2) .EQ. 'up') THEN
BACKSPACE 5
READ(5,100) UPSLOPE, BANK_ANG, BottomW, QU1, QU2
100 FORMAT(10X,5F10.0)
GO TO 50
ENDIF
C
C Read the equation for the Manning's n roughness coefficient
C expressed as a function of hydraulic depth (n = a*D**b)
C
IF (LABEL(1:2) .EQ. 'MA' .OR. LABEL(1:2)
+ LABEL (1:2) EQ. 'ma') THEN
BACKSPACE 5 READ(5,100) aMan, bMan, Rlow, Rhigh
GO TO 50
ENDIF
C
C Read the switch that instructs the program to just compute
C minimum unit stream power for the upstream river channel.
C
60
.OR.
IF (LABEL(1:2) .EQ. 'IN' .OR. LABEL(1:2) .EQ. 'In'
+ LABEL(1:2) .EQ. 'in') THEN
IQU = 1
GO TO 50
END IF
C
C Read the valley axis coordinates and compute the angles and
C distances.
C
IF (LABEL(1:2) .EQ. 'VA' .OR. LABEL(1:2) .EQ. 'Va' .OR.
+ LABEL(1:2) .EQ. 'va') THEN
NumVpts = NumVpts + 1
BACKSPACE 5
READ(5,100) Vx(NumVpts), Vy(NumVpts), Vslope(NumVpts)
Vx(NumVpts) = Vx(NumVpts) XMIN
Vy(NumVpts) = Vy(NumVpts) YMIN
C
C Compute the angle of the valley axis (Alphal) and the
C cumulative distance along the valley axis (Vdist).
C
IF (NumVpts .GT. 1) THEN
dX = Vx(NumVpts) Vx(NumVptS1)
dY = Vy(NumVpts) Vy(NumVpts1)
IF (dY GE. 0 . .AND. dX GE. o.) C = 0.
IF (dY GE. 0 . .AND. dX .LT. o.) C = PIE
IF (dY LT. 0 . .AND. dX .LT. o.) C = PIE
IF (dY LT. 0. .AND. dX .GE. 0.) C = 2.*PIE
IF (ABS(dX) .GT. 1.0E 05) THEN
ALPHAv(NumVpts) = ATAN(dY/dX) + C
ELSE
IF (dY .GE. 0.) THEN
ALPHAv(NumVpts) = PIE/2.
ELSE
ALPHAv(NumVpts) = l.*PIE/2.
ENDIF
ENDIF
IF (ALPHAv(NumVpts) .GT. 2.*PIE) THEN
ALPHAv(NumVpts) = ALPHAv(NumVpts) 2.*PIE
ENDIF
AlphaDeg = ALPHAv(NumVpts)*180/PIE
Vdist(NumVpts) = (DX* *2 + dy* *2)**0.5 + Vdist(NumVpts1)
ENDIF
GO TO 50
ENDIF
C
C Read the phase lag coefficients (CLag) that can amplify or
C reduce the phase lag.
C
IF (LABEL(1:2) .EQ. 'PH' .OR. LABEL(1:2) .EQ. 'Ph' .OR.
+ LABEL(1:2) .EQ. 'ph') THEN
61
BACKSPACE 5
NC = NC + 1
READ(5,100) CLAG(NC)
GO TO 50
ENDIF
c Read the global parameters for the channel and floodplain.
c Check to see if the parameters will be globally applied to
c the grid.
IF (LABEL (1:2) , .EQ. 'GL' .OR. LABEL(1:2) EQ. . 'Gl' .OR.
+ LABEL (1:2) , EQ. 'gi') THEN
nGlobal = 1 GO TO 50
r ENDIF
L. c Check to see if the parameters will be locally applied to
c n the grid.
IF (LABEL(1:2) . .EQ. 'LO' .OR. LABEL(1:2) EQ. . 'Lo1 .OR.
+ LABEL(1:2) . EQ. 'lo') THEN
nLocal = 1
GO TO 50
ENDIF
C
C Read the grid coordinates for the lower left and upper right
C (southwest and northeast) corners of the grid.
C
IF (LABEL(1:2) .EQ. 'GR' .OR. LABEL(1:2) .EQ. 'Gr1 .OR.
+ LABEL(1:2) .EQ. 'gr1) THEN
BACKSPACE 5
READ(5,110) XLL0, YLL0, XUR0, YUR0
110 FORMAT(10X,2F10.0,/,10X,2F10.0)
XLL0 = XLL0 XMIN
YLL0 = YLL0 YMIN
XUR0 = XUR0 XMIN
YUR0 = YUR0 YMIN
C
C Determine if the grid coordinates are to be applied
C globally or locally. The beginning and ending indicies
C of the grid are defined by Ibeg and lend for the xdirection
C and Jbeg and Jend for the ydirection.
C
IF (nLocal .EQ. 0 .OR. nGlobal .EQ. 0) THEN
C Apply globally
XLL = XLL0
YLL = YLL0
XUR = XUR0
YUR = YUR0
dxGrid = (XUR XLL)/FLOAT(IMAXg 1)
62
dyGrid = (YUR YLL)/FLOAT(JMAXg 1)
Ibeg = 1
lend = IMAXg
Jbeg = 1
Jend = JMAXg
WRITE(30,115) XLL+XMIN, YLL+YMIN, XUR+XMIN, YUR+YMIN,
+ dxGrid, dyGrid, Ibeg, lend, Jbeg, Jend
115 FORMAT(4F10.1,2F6.1,416)
ELSE
C Apply locally
XLL2 = XLLO
YLL2 = YLLO
XUR2 = XURO
YUR2 = YURO
Ibeg = INT((XLL2  XLL)/dxGrid + 0.4999) + 1
lend = INT((XUR2  XLL)/dxGrid + 0.4999) + 1
Jbeg = INT((YLL2  YLL)/dyGrid + 0.4999) + 1
Jend = INT((YUR2  YLL)/dyGrid + 0.4999) + 1
WRITE(30,115) XLL2+XMIN, YLL2+YMIN, XUR2+XMIN, YUR2+YMIN,
+ dxGrid, dyGrid, Ibeg, lend, Jbeg, Jend
ENDIF
GO TO 50
ENDIF
C
C Read the hydraulic coefficient for bank erosion (al).
C
IF (LABEL(1:2) .EQ. 'HY' .OR. LABEL(1:2) .EQ. 'Hy' .OR.
+ LABEL(1:2) .EQ. 'hy') THEN
BACKSPACE 5
READ(5,100) al
GO TO 50
ENDIF
C
C Read the vegetation coefficient (a2) and parameters for bank
C erosion: root density (rg), root depth (rd), and bank height
C (hb). The root density (rg) is the fraction of bank covered
C by roots. The root depth (rd) and bank height (hb) must be
C given in the same units (e.g., feet, inches, meters, or
C centimeters)
C
IF (LABEL .EQ. 'VEG' .OR. LABEL .EQ. 'Veg' .OR.
+ LABEL .EQ. 'veg') THEN
BACKSPACE 5
IF (nLocal .EQ. 0) THEN
READ(5,120) a2, rg(Ibeg,Jbeg), rd(Ibeg,Jbeg),
hb (Ibeg, Jbeg)
120 FORMAT(10X,4F10.0)
ELSE
READ(5,125) rg(Ibeg,Jbeg), rd(Ibeg,Jbeg), hb(Ibeg,Jbeg)
125 FORMAT(20X,4F10.0)
63
n n
ENDIF
Assign the vegetation parameters to the grid, either
globally or locally.
DO 140 I = Ibeg,lend
DO 130 J = Jbeg,Jend
rg(I,J) = rg(Ibeg,Jbeg)
rd(I,J) = rd(Ibeg,Jbeg)
hb(I,J) = hb(Ibeg,Jbeg)
CONTINUE
CONTINUE
GO TO 50
ENDIF
C
C Read the large woody debris coefficient (a3) and parameters
C for bank erosion: large woody debris density (LWD) and the
C average height of woody debris jams (dw). The large woody
C debris density (LWD) is the fraction of bank covered by large
C woody debris. The average height of woody debris jams (dw)
C must be given in the same units as the hydraulic depth.
C
IF (LABEL(1:2) .EQ. 'LW' .OR. LABEL(1:2) .EQ. 'Lw' .OR.
+ LABEL(1:2) .EQ. 'lw') THEN
BACKSPACE 5
IF (nLocal .EQ. 0) THEN
READ(5,120) a3, LWD(Ibeg,Jbeg), dw(Ibeg,Jbeg)
ELSE
READ(5,125) LWD(Ibeg,Jbeg), dw(Ibeg,Jbeg)
ENDIF
C
C Assign the large woody debris jam heigth to the grid, either
C globally or locally.
C
DO 160 I = Ibeg,lend
DO 150 J = Jbeg,Jend
LWD(I,J) = LWD(Ibeg,Jbeg)
dw(I,J) = dw(Ibeg,Jbeg)
150 CONTINUE
160 CONTINUE
GO TO 50
ENDIF
C
C Read the cohesion coefficient (a4) and the parameter for bank
C erosion: Plasticity index (PI), which is the liquid limit
minus
the plastic limit of the soil (PI = LL PL). A Pivalue of
zero would indicate no cohesion, where as a PIvalue between
C
C
C
C
130
140
64
C about 30 and 70 would be for inorganic clays of high
plasticity.
C
IF (LABEL(1:2) .EQ. 'CO' .OR. LABEL(1:2) .EQ. 'Co' .OR.
+ LABEL(1:2) .EQ. 'CO') THEN
BACKSPACE 5
IF (nLocal EQ. 0) THEN
READ(5,120) a4, PI(Ibeg,Jbeg)
ELSE
READ(5,125) PI(Ibeg,Jbeg)
ENDIF
C
C Assign the plasticity index to the grid, either
C globally or locally.
C
DO 180 I = Ibeg,lend
DO 170 J = Jbeg,Jend
PI(I,J) = PI(Ibeg,Jbeg)
170 CONTINUE
180 CONTINUE
GO TO 50
ENDIF
C
C Read the armoring coefficient (a5) and the parameter for bank
C erosion: Fraction of bank material particles that are too
C large for sediment transport (dc). A dcvalue of zero would
C mean that all particles sizes within the bank can be
transported
C by the streamflow. A dcvalue of 1 would mean than that all
C particles sizes are too large for transport and the bank is
C considered armored.
C
IF (LABEL(1:2) .EQ. 'AR' .OR. LABEL(1:2) .EQ. 'Ar' .OR.
+ LABEL(1:2) .EQ. 'ar') THEN
BACKSPACE 5
IF (nLocal .EQ. 0) THEN
READ(5,120) a5, dc(Ibeg,Jbeg)
ELSE
READ(5,125) dc(Ibeg,Jbeg)
ENDIF
C
C Assign the armoring parameters to the grid, either
C globally or locally.
C
DO 200 I = Ibeg,lend
DO 190 J = Jbeg,Jend
dc(I,J) = dc(Ibeg,Jbeg)
190 CONTINUE
200 CONTINUE
GO TO 50
65
CJ L)
ENDIF
C
C Read the velocity coefficient (a6).
C
IF (LABEL .EQ. 'VEL' .OR. LABEL .EQ. 'Vel' .OR.
+ LABEL .EQ. 'vel') THEN
BACKSPACE 5
READ(5,120) a6
GO TO 50
ENDIF
C
C Read the times that model results are to be printed.
C
IF (LABEL(1:2) .EQ. 'PT' .OR. LABEL(1:2) .EQ. 'Pt'
+ LABEL(1:2) .EQ. 'pt') THEN
BACKSPACE 5
IPT = IPT + 1
IF (IPT .GT. 20) IPT =20
READ(5,210) IPMM(IPT), IPDD(IPT), IPYYYY(IPT)
210 FORMAT(10X,12,IX,12,IX,14)
GO TO 50
ENDIF
C
C Read the Manning's n roughness coefficients.
C
IF (LABEL .EQ. 'NV' .OR. LABEL .EQ. 'Nv' .OR.
+ LABEL EQ. 'nV' .OR. LABEL .EQ. 'nv') THEN
BACKSPACE 5
READ(5,120) NVALUE
GO TO 50
ENDIF
GO TO 50
450 CLOSE(4)
505 CONTINUE
DO 510 J = JMAXg,1,1
WRITE(30,500) (rg(I,J), I = l,IMAXg)
500 FORMAT(4F10.2,101F6.2)
510 CONTINUE
DO 520 J = JMAXg,1,1
WRITE(30,500) (rd(I,J), I = l,IMAXg)
520 CONTINUE
DO 530 J = JMAXg,1,1
WRITE(30,500) (hb(I,J), I = l,IMAXg)
530 CONTINUE
DO 540 J = JMAXg,1,1
WRITE(30,500) (LWD(I,J), I = l,IMAXg)
540 CONTINUE
DO 550 J = JMAXg,1,1
WRITE(30,500) (dw(I,J), I = l,IMAXg)
.OR.
66
550 CONTINUE
DO 560 J = JMAXg,l,l
WRITE(30,500) (PI(I,J), I = l,IMAXg)
560 CONTINUE
DO 570 J = JMAXg,1,1
WRITE(30,500) (dc(I,J), I = l,IMAXg)
570 CONTINUE
CLOSE (3 0)
RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE READXY (NTJMPTS)
c
C Routine to read the XY coordinate data defining the alingment
C of the meandering river channel.
C
DOUBLE PRECISION X, Y, CHANLNG, Xc, Yc, Rc, THETA, PHI
COMMON / PLAN_DATA / X(1000), Y(1000), CHANLNG(1000),
+Xc(1000), Yc(1000), Rc(1000), THETA(IOOO), PHI(1000),
+SEIN(1000), SEIN2(1000), XLT(1000), YLT(1000), XRT(1000)
+YRT(1000), XMIN, YMIN
INTEGER SEIN, SEIN2
COMMON / NAMES / TITLEXY(5), TITLEQ(5), TITLEG(5)
CHARACTER FNAME5*100, TITLEXY*80, TITLEQ*80, TITLEG*80
C FNAME5(1:40) =
C FNAME5(41:80) =
C FNAME5(81:90) =
C FNAME5(1:40) =
C FNAME5(41:80) =
C FNAME5(81:90) =
C FNAME5(1:40)
C FNAME5(41:80)
C FNAME5(81:100) =
C FNAME5(1:40)
C FNAME5(1:40)
FNAME5(1:40)
'C:\Tim Randle\Graduate School\Thesis Cha'
'nnel Migration Model\Sine Generated Curv'
'e.prn '
'F:\Graduate School\Thesis\Sine Generated'
1 Curve.prn '
I I
'C:\Tim Randle\Graduate School\Thesis Cha'
'nnel Migration Model\Sine Generated curv'
'e 1.1 sinuosity.prn '
'Sine Generated curve 1.1 sinuosity.prn '
'1950 dense Alignment.prn '
'1950 dense Alignment.prn
C
C
c
FNAME5(41:80) = '
FNAME5(81:100) = '
OPEN(UNIT=5,FILE=FNAME5,STATUS=
Read the first 5 title lines.
'OLD'
DO 40 N = 1,5
READ(5,30) TITLEXY(N)
IF (N .LE. 3) WRITE(6,30) TITLEXY(N)
3 0 FORMAT (A)
40 CONTINUE
67
o n
Read the XY coordinate data of the meandering river channel.
XMIN = 9999999.
YMIN = 9999999.
1 = 1
100 READ(5,*,END=200) X(I), Y(I)
C
C Determine the minimum coordinates.
C
IF (X(I) .LT. XMIN) XMIN = X(I)
IF (Y(I) .LT. YMIN) YMIN = Y(I)
1 = 1 + 1
GO TO 100
200 NUMPTS =11
C
C Subtract the minimum values from the X and Y coordinates.
C
INT(XMIN/1000. :
INT(YMIN/1000. :
*1000.
*1000.
IF (XMIN .GT. 1000.) XMIN
IF (YMIN .GT. 1000.) YMIN
DO 210 I = 1, NUMPTS
X(I) = X(I) XMIN
Y(I) = Y(I)  YMIN
210 CONTINUE
CLOSE (5)
RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE SMOOTH(NUMPTS)
C
C Routine to smooth the XY coordinate data using a quadradic
C regression equation (y = bo + bl x + b2 x*2) that is locally
C applied to a moving group of seven points.
C
DOUBLE PRECISION X, Y, CHANLNG, Xc, Yc, Rc, THETA, PHI,
+X2, Y2
COMMON / PLAN_DATA / X(1000), Y(1000), CHANLNG(1000),
+Xc(1000), Yc(1000), Rc(1000), THETA(IOOO), PHI(1000),
+SEIN(1000), SEIN2(1000), XLT(1000), YLT(1000), XRT(1000),
+YRT(1000), XMIN, YMIN
INTEGER SEIN, SEIN2
COMMON / XY2 / X2(1000), Y2(1000)
COMMON / REGDAT / XR(10), YR(10), NUMREG
DOUBLE PRECISION XR, YR
REAL BETA0, BETA1, BETA2, R2, AR2
DATA NUMREG / 7 /
C
C NUMPTS is the number points describing the river alingment.
C PCTpts is the percentage of the total number of points.
C NUMREG is the number of points in the moving average.
68
c
Ibeg = NUMREG/2 + 1
lend = NUMPTS NUMREG/2
Iend2 = lend
DO 40 I = Ibeg,lend
Kbeg = I NUMREG/2
Kend = I + NUMREG/2
K = Kbeg
C
C Compute incremental distance in X and Y (dX and dY) for the
C local group of points.
C
dX = X(Kend) X(Kbeg)
dY = Y(Kend) Y(Kbeg)
c
C If there is a greater increment in the Xdirection, then
C determine the regression equation in terms of X (dX > dY).
C Otherwise, determine the regression equation in terms of Y.
c
DO 10 J = 1,NUMREG
IF (ABS(dX) .GE. 0.001*ABS(dY)) THEN
XR (J) = X (K)
YR(J) = Y (K)
ELSE
YR (J) = X (K)
XR (J) = Y (K)
ENDIF
. K = K + 1
10 CONTINUE
C
C Call a routine to compute a quadratic regression equation
C based on the local group of NUMREG points.
C
CALL POLY2REG(BETAO,BETA1,BETA2,R2,AR2)
WRITE(10,20) I, X(I), Y(I), R2, AR2, BETAO,BETA1,BETA2
20 FORMAT(I5,2F10.1,2F6.3,3E14.6)
IF (I .EQ. 330 .OR. I .EQ. 31S) THEN
XXX = 1.
ENDIF
C
C Compute a new x and y coordinate.
C
C IF (R2 .GT. 0.85) THEN
IF (R2 .GT. 0.50) THEN
C
C If the correlation is high enough (RA2 > 0.95), then
C compute a new x and y coordinate based on the regression
C equation.
C
IF (ABS(dX) .GE. 0.001*ABS(dY)) THEN
69
c
C For each x(i) position, compute a new y(i) position
C based on the regression equation.
C
X2 (I) = X (I)
Y2(I) = BETAO + BETA1*X(I) + BETA2*X(I)**2
IF (I .EQ. lend) THEN
DO 22 IK = Iend+1,NUMPTS
X2 (IK) = X (IK)
Y2(IK) = BETAO + BETA1*X(IK) + BETA2*X(IK)**2
22 CONTINUE
Iend2 = NUMPTS
ENDIF
ELSE
C
C For each y(i) position, compute a new x(i) position
C based on the regression equation.
C
Y2 (I) = Y (I)
X2(I) = BETAO + BETA1*Y(I) + BETA2*Y(I)**2
IF (I .EQ. lend) THEN
DO 24 IK = Iend+1,NUMPTS
X2(IK) = BETAO + BETA1*Y(IK) + BETA2*Y(IK)**2
Y2 (IK) = Y (IK)
24 CONTINUE
Iend2 = NUMPTS
ENDIF
ENDIF
ELSE
C
C Compute new coordinates using a 5point moving average.
C
SUMX = 0.
SUMY = 0.
Mbeg =12
Mend =1+2
DO 30 M = Mbeg,Mend
SUMX = SUMX + X(M)
SUMY = SUMY + Y(M)
3 0 CONTINUE
X2(I) = SUMX/5.
Y2(I) = SUMY/5.
XXX = 1.
ENDIF
40 CONTINUE
C
C Replace the points describing the river alignment with the
C smoothed points from the moving average.
C
DO 60 I = Ibeg,Iend2
70
X (I) = X2 (I)
Y(I) = Y2(I)
60 CONTINUE
RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE POLY2REG(BETAO,BETA1,BETA2,R2,AR2)
COMMON / REGDAT / XR(10), YR(10), NUMREG
DOUBLE PRECISION XR, YR, SUMX, SUMY, SUMXY, SUMX2, SUMY2,
SUMX3,
+SUMX4, SUMX2Y, AveX, AveY, AveX2, Sly, S2y, Sll, S12, S22,
+SSE, SST
REAL BETAO, BETA1, BETA2, R2, AR2
C
C Initialize the sums to zero.
C
SUMX = 0.
SUMY = 0.
SUMXY = 0.
SUMX2 = 0.
SUMY2 = 0.
SUMX3 = 0.
SUMX4 = 0.
SUMX2Y = 0.
C
C Compute the sums of x, y, xy, x2, x3, and x2y
C
DO 10 N = 1, NUMREG
SUMX = SUMX + XR (N)
SUMY = SUMY + YR (N)
SUMXY = SUMXY + XR (N) *YR (N)
SUMX2 = SUMX2 + XR (N) **2
SUMY2 = SUMY2 + YR(N)**2
SUMX3 = SUMX3 + XR(N)**3
SUMX4 = SUMX4 + XR(N)**4
SUMX2Y = SUMX2Y + (XR(N)* *2)*YR(N)
10 CONTINUE
C
C Compute the statistical properties: average x, y, x2
C
dNUM = FLOAT(NUMREG)
AveX = SUMX/dNUM
AveY = SUMY/dNUM
AveX2 = SUMX2/dNUM
C
C Compute the equation coefficients (BetaO, Betal, Beta2)
C
Sly = SUMXY dNUM*AveX*AveY
S2y = SUMX2Y dNUM*AveX2*AveY
71
511 = SUMX2 dNUM*AveX**2
512 = SUMX3 dNUM*AveX*AveX2
S22 = SUMX4 dNUM*AveX2**2
BETA1 = (Sly*S22 S2y*S12)/(S11*S22 S12**2)
BETA2 = (S2Y*S11 Sly*S12)/(S11*S22 S12**2)
BETAO = AveY BETAl*AveX BETA2 *AveX2
C
C Compute the sum of squared residuals (SSE)
C
SSE = SUMY2 BETAO*SUMY BETA1*SUMXY BETA2*SUMX2Y
C
C Compute the coefficient of multiple determination
C
SST =0.
dK = 2.
DO 20 N = 1,NUMREG
SST = SST + (YR(N)  AveY)**2
20 CONTINUE
R2 = 1. SSE/SST
AR2 = ((dNUM 1.)*R2 dK)/(dNUM 1. dK)
RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE SMOOTH_Bank_Erosion(NUMPTS)
c
C Routine to smooth the XY coordinate data using a 5point
C moving average.
C
DOUBLE PRECISION X, Y, CHANLNG, Xc, Yc, Rc, THETA, PHI,
+X2, Y2, XR, YR
COMMON / EROSION / CsWRc(lOOO), ECOEF(1000,50), BE(1000)
COMMON / BESmooth / BE2 (1000)
COMMON / PLAN_DATA / X(1000), Y(1000), CHANLNG(1000),
+Xc(1000), Yc(1000), Rc(1000), THETA(IOOO), PHI(1000),
+SEIN(1000), SEIN2(1000), XLT(1000), YLT(1000), XRT(1000),
+YRT(1000), XMIN, YMIN
INTEGER SEIN, SEIN2
COMMON / REGDAT / XR(10), YR(10), NUMREG
REAL BETAO, BETA1, BETA2, R2, AR2
DATA NUMREG / 7 /
C
C Compute the 7point moving average for the bank erosion.
C
Npts = 7
Ibeg = Npts/2 + 1
lend = NUMPTS Npts/2
DO 20 I = Ibeg,lend
SUMBE = 0.
Kbeg = I Npts/2
72
n n p, n n
Kend = I + Npts/2
DO 10 K = Kbeg,Kend
SUMBE = SUMBE + BE(K)
10 CONTINUE
BE2(I) = SUMBE/FLOAT(Npts)
20 CONTINUE
C
C Replace the computed bank erosion amounts with the smoothed
C values.
C
DO 30 I = Ibeg,lend
BE(I) = BE2(I)
30 CONTINUE
Smooth the bank erosion amounts near the upstream and
ends.
(1./3.)*BE2(Ibeg)
(2 ./3 .) *BE2 (Ibeg)
(2./3.)*BE2(lend)
(1./3.)*BE2(lend)
BE(Ibeg2) =
BE(Ibeg1) =
BE (Iend+1) =
BE(Iend+2) =
RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE READFLOW(NUMQ)
C
C Routine to read the discharge or river flow data defining
C the alignment of the meandering river channel.
C
C
C
C
C
C
C
C
DOUBLE PRECISION X, Y, CHANLNG, Xc, Yc, Rc, THETA, PHI
COMMON / HYDROLOGY / QTHRESHOLD, JDAY(20000), FLOW(20000)
INTEGER DD, YYYY, TIME
REAL JDAY, JDATE
COMMON / NAMES / TITLEXY(5) TITLEQ(5) TITLEG(5)
CHARACTER FNAME5*90, TITLEXY*80, TITLEQ*80, TITLEG*80
FNAME5(1:40) = 'C:\Tim Randle\Graduate School\Thesis Cha'
FNAME5(41:80) = 'nnel Migration Model\Flow Hydrograph.prn'
FNAME5(81:90) = ' '
FNAME5(1:40) = 'Flow Hydrograph.prn '
FNAME5(1:40) = 'Hoh River Discharge.prn '
FNAME5(41:80) = ' 1
FNAME5(81:90) = 1 '
OPEN (UNIT=5,FILE=FNAME5,STATUS='OLD')
Read the first 5 title lines.
DO 40 N = 1,5
READ(5,30) TITLEQ (N)
73
IF (N .LE. 3) WRITE(6,30) TITLEQ(N)
3 0 FORMAT(A)
40 CONTINUE
C
C Read the flow threshold where only flows above this threshold
C need to be considered by the channel migration model,
c
READ (5,50) QTHRESHOLD
50 FORMAT(1SX,F10.0)
C
C Read the discharge or river flow data of the meandering river
C channel.
C
C OPEN(20,FILE='DATES.OUT',STATUS='UNKNOWN')
N = 1
100 READ(5,110,END=200) MM, DD, YYYY, TIME, FLOW(N)
110 FORMAT(12,IX,12,IX,14,2X,14,F10.0)
CALL JULIAN(MM,DD,YYYY,TIME,JDATE)
JDAY(N) = JDATE
CALL CDATE(iMM,iDD,iYYYY,TIME,JDATE)
JM = iMM MM
JD = iDD DD
JY = iYYYY YYYY
C WRITE(20,120) MM, DD, YYYY, TIME, JDATE, iMM,
C +JM, JD, JY
C 120 FORMAT(12,'/',12,'/',14,2X,14,F10.2,4X,12,1/
C +4X,315)
N = N + 1
GO TO 100
200 NUMQ = N 1
CLOSE(5)
C CLOSE(20)
RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE JULIAN(MM,DD,YYYY,TIME,JDATE)
DIMENSION MDAYS(12)
INTEGER DD, YYYY, TIME, HOUR
REAL JDATE
DATA MDAYS / 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
C
C Compute the julian date from January 1, 1900.
C
C Test for a leap year. If so, change the number of days in
C February to 29 rather than 28.
C
IF (YYYY/4*4 .EQ. YYYY) THEN
MDAYS(2) =29
ELSE
iDD, iYYYY,
',12,'/',14,
74
n n
MDAYS(2) = 28
ENDIF
C
C Start by computing the number of days from 1/1/1900 to the
C end of the year preceeding MM/DD/YYYY.
C
LEAPYEARS = (YYYY 1900 1)/4 + 1
JDATE = (YYYY 1900)*365 + LEAPYEARS
IF (MM .EQ. 1) THEN
JDATE = JDATE + DD
ELSE
DO 20 M = 1,MM1
JDATE = JDATE + MDAYS(M)
2 0 CONTINUE
JDATE = JDATE + DD
ENDIF
C
Convert the time to a decimal fraction of a day.
HOUR = TIME/100
MIN = TIME HOUR*100
JDATE = JDATE + FLOAT(HOUR)/24. + FLOAT(MIN)/60.
RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE CDATE(MM,DD,YYYY,TIME,JDATE)
DIMENSION MDAYS(12)
INTEGER DD, YYYY, TIME, HOUR
REAL JDATE
DATA MDAYS / 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
C
C Compute the number of years (IYY) and the number of leap years
C (LY) since 1900.
C
IF (JDATE .GT. 18628.) THEN
XXX = 1.
ENDIF
IYY = INT(JDATE/365)
LY = INT(IYY/4) + 1
C IF(INT(JDATEIYY*365) .LT. LY) THEN
C IYY = IYY 1
C ENDIF
JTEST = INT(JDATEIYY*365)
IF (IYY/4*4 .NE. IYY) THEN
IF (JTEST .LE. LY) IYY = IYY 1
ELSE
IF (JTEST .LT. LY) IYY = IYY 1
ENDIF
LY = INT(IYY/4) + 1
75
c
C Test for a leap year. If so, change the number of days in
C February to 29 rather than 28.
C
IF (IYY/4*4 .EQ. IYY) THEN
MDAYS(2) = 29
JDAYS2 = 366
LY = LY 1
ELSE
MDAYS(2) = 28
JDAYS2 = 365
ENDIF
YYYY = 1900 + IYY
C
C Compute the julian date from January 1, 1900.
C Compute the number of julian days within the year (JDAYS)
C
JDAYS = INT(JDATE) IYY*365 LY
C
C Determine the month based on the number of julian days within
C the year.
C
DO 20 M = 12,1,1
JDAYS1 = JDAYS2 MDAYS(M)
IF (JDAYS .LE. JDAYS2 .AND. JDAYS .GT. JDAYS1) THEN
MM = M
DD = JDAYS JDAYS1
DAY = (JDATE FLOAT(IYY*365 + LY + JDAYS))
HOUR = INT(DAY*24.)
MIN = INT(DAY FLOAT(HOUR)/24.)
TIME = HOUR*100 + MIN
GO TO 30
ENDIF
JDAYS2 = JDAYS1
20 CONTINUE
30 RETURN
C
C Convert the time to a decimal fraction of a day.
C
C HOUR = TIME/100
C MIN = TIME HOUR* 100
C JDATE = JDATE + FLOAT(HOUR)/24. + FLOAT(MIN)/60.
C* RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE RADIUS(NUMPTS,IPRINT)
c
C Routine to compute the radius of curvature, the point of
C curvature, the radius azimuth, and the tangent angle of the
76
n o
meandering river channel.
DOUBLE PRECISION X, Y, CHANLNG, Xc, Yc, Rc, THETA, PHI,
+B, SM, PIE, dSO, dSl, dS2, XW, YW, XcW, YcW
COMMON / PLAN_DATA / X(1000}, Y(1000), CHANLNG(1000),
+Xc(1000), Yc(1000) Rc(100 0) THETA(IOOO), PHI (1000),
+SEIN(1000), SEIN2(1000), XLT(1000), YLT(IOOO), XRT(IOOO),
+YRT(1000), XMIN, YMIN
INTEGER SEIN, SEIN2
COMMON / HYDR / DISCH, AREA, Hdepth, HdepthM, SLOPE, TDF,
NVALUE
COMMON / WIDTHS / CHwidth(1000), CHwidthMax(1000),
CHWtime(1000),
+ WIDTHM, WIDTHL, WIDTHU, WIDTHI
REAL NVALUE
COMMON / NAMES / TITLEXY(5), TITLEQ(5), TITLEG(5)
CHARACTER TITLEXY*80, TITLEQ*80, TITLEG*80
COMMON / CONSTANTS / PIE
C
IF (IPRINT .EQ. 1) WRITE(6,10)
10 FORMAT(' I chanlng X Y XLT YLT
I
/
+ ' XRT YRT Xc Yc Rc ,
+ ' Theta Phi Turn RcW TlocL
I
/
+ 'TlocR Bank Eros Phase Lag CH Width HydDepth Max
Depth' )
C
C Compute the geometric parameters for each point along the
C meandering river channel.
C
C
C
C
C
C
C
C
+
CHANLNG(1) = 0.
DO 100 I = 2,NUMPTS1
Test to see if the three X points (11, I, 1+1) or if any of
the three Y points are identical.
60 dXl = (X(I) X(Il))
dX2 = (X(I + 1)  X(I) )
dYl = (Y(I)  Y(Il))
dY2 = (Y(1+1) Y(I))
Check to see if X(i) and X(il) are the same.
IF (ABS(dXl) .LT. 0.001) THEN
Yc(I) = (Y(I)**2 Y(Il)**2)/(2.*Y(I)  2.*Y(I1))
Xc(I) = (2.*Y(1 + 1)  2.*Y(I))*Yc(I) + Y(I)**2 Y(I + 1)**2
+ (X(I)**2 X(I+1)**2)/(2.*X(I)  2.*X(I+1))
77
o n
Rc(I) = ((X(I) Xc(I)) **2 + (y(I) Ye(I))**2)**0.5
GO TO 75
ENDIF
C
C Check to see if X(i) and X(i+1) are the same.
C
IF (ABS(dX2) .LT. 0.001) THEN
Yc (I) = (Y(I + 1)**2 Y(l)**2)/(2. *Y (1 + 1) 2.*Y(I))
Xc(I) = (2.*Y(I) 2.*Y(11))*Yc(I) + Y(11)**2 Y(I)**2
+
+ (X(11)**2 X(I)**2)/(2.*X(I1)  2.*X(I))
Rc(I) = ((X(I)  Xc(I))**2 + (Y(I)  Yc (I))**2)**0.5
GO TO 75
ENDIF
C
C Check to see if Y(i) and Y(il) are the same.
C
IF (ABS(dYl) .LT. 0.001) THEN
Xc(I) = (X(I)**2 X(Il)**2)/(2.*X(11) 2.*X(I))
Yc(I) = ((2.*X(I)  2.*X(1+1))*Xc(I) + X(1+1)**2 X(I)**2
+
+ Y(1+1)**2 Y(I)**2)/(2.*Y(I+1)  2.*Y(I))
Rc(I) = ((X(I) Xc(I))**2 + (Y(I) Yc(I))**2)**0.5
GO TO 75
ENDIF
C
C Check to see if Y(i) and Y(i+1) are the same.
C
IF (ABS(dY2) .LT. 0.001) THEN
Xc(I) = (X(I)**2 X(I+1)**2)/(2.*X(I)  2.*X(I+1))
Yc(I) = ((2.*X(11)  2.*X(I))*Xc(I) + X(I)**2 X (11)**2
+
+ Y(I)**2 Y(Il)**2)/(2.*Y(I)  2.*Y(11))
Rc(I) = ((X(I)  Xc(I))**2 + (Y(I)  Yc(I))**2)**0.5
GO TO 75
ENDIF
Test to see if any of the three points (11, I, and 1+1)
form
C a straight line by comparing the xy slopes from 11 to I
and
C I to 1+1.
C
Slopel = dYl/dXl
Slope2 = dY2/dX2
C
C If the two slopes are the same, then a straight line exists
C and the radius of curvature is infinite (set to a default of
C Rc = 100*Width).
C
78
on noon oono non onnnnn no
c
IF (ABS(slopel Slope2) .LT. 1.0E5) THEN
Compute the azimuth of the channel alingment (PHI).
PHI(I) = ATAN(Slope1)
PHIDEG = PHI(I)*180./PIE
Test for the NW and SW quadrants.
IF ((X(I+1) X(Il)) .LT. 0.) PHI(I)
PIE + PHI(I)
Test for the SE quadrant.
IF ((Y(1+1) Y(Il)) .LT. 0. .AND.
+ (X(I+1) X(Il)) .GE. 0.) PHI(I) = 2.*PIE + PHI(I)
PHIDEG = PHI(I)*180./PIE
Compute the azimuth of the radius (THETA).
THETA(I) = PHI(I)  PIE/2.
IF (THETA(I) .LT. 0.) THETA(I) = THETA(I) + 2.*PIE
THETADEG = THETA(I)*180./PIE
IF (I .GT. 1) THETA(I) = THETA(I1)
Assign the radius to 100 times the channel width and
compute the coordinates of the radius point (Xc and Yc)
Rc(I) = 100.*WIDTHM
Xc(I) = X(I) Rc(I)*C0S(THETA(I))
Yc(I) = Y(I) Rc(I)*SIN(THETA(I))
GO TO 80
ENDIF
Compute the xcoordinate (Xc) for the center point of
curvature.
+
+
+
+
SM = 1.  (2.*Y(I)  2.*Y(I1))*(2.*X(I)  2.*X(I+1))/
((2.*Y(1+1)  2.*Y(I))*(2.*X(11)  2.*X(I)))
IF (ABS(SM) .GT. 1.0E5) THEN
Xcl = ((2.*Y(I) 2.*Y(I1))*
(X(1 + 1)**2 X(I)**2 + Y(1 + 1)**2 Y(I)**2)/
(2.*Y(1+1)  2.*Y (I)) + Y(11)**2 Y(I)**2 +
X (11)**2 X(I)**2)/(2.*X(I1)  2.*X(I))/SM
ELSE
Xc(I) = Xc(11)
Yc(I) = Yc(11)
Rc(I) = Rc(Il)
C
Compute the azimuth (THETA) for the radius of curvature
79
oooooo noon n noon
dX = X(I)  Xc (I)
dY = Y (I) Yc (I)
IF (dY .GE. 0 . .AND. dX .GE. 0.) C = 0.
IF (dY .GE. 0 . .AND. dX .LT. 0.) C = PIE
IF (dY .LT. 0 . .AND. dX .LT. 0.) C = PIE
IF (dY .LT. 0 . .AND. dX .GE. 0. ) C = 2.*PIE
IF (ABS (dX) .GT. 1.0E 05) THEN
THETA(I) = ATAN(dY/dX) + C
ELSE
IF (dY .GE. 0.) THEN
THETA(I) = PIE/2.
ELSE
THETA(I) = 1.5 *PIE
ENDIF
ENDIF
IF (THETA(I) .GT. 2.*PIE) THETA(I) = THETA(I)
2.*PIE
Compute the azimuth (PHI) for the alignment of the
PHI(I) = PHI(11)
GO TO 80
ENDIF
Xc2 = ((2.*Y(I)  2.*Y(11))*((2.*X(I)  2.*X(1+1))*Xcl +
+ X(1 + 1)**2 X(I)**2 + Y(1 + 1)**2 Y(I)**2)/
+ (2.*Y(1 + 1) 2.*Y(I) ) + Y(I1)**2 Y(I)**2 +
+ X(I1)**2 X(I)**2)/(2.*X(11)  2.*X(I))
XERR = Xcl Xc2
Xc (I) = Xc2
Compute the ycoordinate (Yc) for the center point of
curvature.
Yc (I) = ( (2.*X(I) 2.*X(I + 1))*Xc(I) +
+ X(1 + 1)**2 X(I)**2 + Y(1 + 1)**2 Y(I)**2)/
+ (2.*Y(1 + 1) 2.*Y(I))
Compute the radius of curvature (Rc).
Rc(I) = ((Xc(I) X(I))**2 + (Yc(I) Y(I))**2)**0.5
Test for proper convergance
IF (ABS(XERR)/Rc(I) .GT. 0.01) THEN
WRITE(11,70) I, Xcl, Xc2, XERR, Rc(I)
70 FORMAT('I =',15,' Xcl =',F10.2,', Xc2 =',F10.2,
+ ', XERR =',F10.2,/,'Radius of Curvature is',F10.1)
ENDIF
80
o o
Compute the azimuth (THETA) for the radius of curvature.
dX = XI [I)  Xc (I)
dY = Y 1 [I)  Yc (I)
IF (dY .GE. 0. .AND. dX .GE. 0.) C = 0.
IF (dY .GE. 0. .AND. dX LT. 0.) C = PIE
IF (dY .LT. 0. .AND. dX .LT. 0.) C = PIE
IF (dY .LT. 0. .AND. dX .GE. 0.) C = 2 *PIE
IF (ABS(dX) .GT. 1.0E 05) THEN
THETA(I) = ATAN(dY/dX) + C
ELSE
IF (dY .GE. 0.) THEN
THETA(I) = PIE/2.
ELSE
THETA(I) = l.*PIE/2.
ENDIF
ENDIF
IF (THETA(I) .GT. 2.*PIE) THETA(I) = THETA(I)  2.*PIE
C
C Compute the tamgent angle (PHI) for the menader curve.
C PHIch is the approximate angle along the direction of the
C channel alignment.
C
dYch = Y(1+1)  Y(Il)
dXch = X(I+1)  X(Il)
IF (dYch .GE. 0 . .AND. dXch .GE. 0.) C = 0.
IF (dYch .GE. 0 . .AND. dXch .LT. 0.) C = PIE
IF (dYch .LT. 0 . .AND. dXch .LT. 0.) C = PIE
IF (dYch .LT. 0 . .AND. dXch .GE. 0.) C = 2.*PIE
IF (ABS(dXch) .GT. 1.0E 05) THEN
PHIch = ATAN(dYch/dXch) + C
ELSE
IF (dYch .GE. 0.) THEN
PHIch = PIE/2.
ELSE
PHIch = 1.5*PIE
ENDIF
ENDIF
C
C Compute the angle difference between the azimuths of the
C channel alingment (PHIch) and the radius of curvaure (THETA)
C
dANGLE = (PHIch THETA(I))*180./PIE
IF (dANGLE .LT. 0.) dANGLE = 360. + dANGLE
IF (ABS(dANGLE 90.) .LT. 45.) THEN
C Left turn, right bank erosion, positive sign.
PHI(I) = THETA(I) + PIE/2.
ELSE
C Right turn, left bank erosion, negative sign.
PHI (I) = THETA(I)  PIE/2.
81
c
c
c
c
c
c
c
c
c
c
ENDIF
IF (PHI(I) .GT. 2.*PIE) PHI(I) = PHI(I)  2.*PIE
IF (PHI(I) .LT. 0.) PHI(I) = PHI(I) + 2.*PIE
Compute the cumulative channel length (CHANLNG).
Compute the cord lenghts (CD1 and CD2) between points
11 to I and I to 1+1
80 CD1 = (((X(I) X(11))**2 + (Y(I) Y(Il))**2)**0.5)/2.
CD2 = (((X(I+1) X(I))**2 + (Y(1+1) Y(I))**2)**0.5)/2.
Compute the angles (ANG1 and ANG2) between points
11 to I and I to 1+1
ANG1 = ATAN(CDl/Rc(I))
ANG2 = ATAN(CD2/RC(I))
dSl = ANGl*Rc(I)
dS2 = ANG2*Rc(I)
IF (I .EQ. 2) dSO = ANG1/2.*Rc(I)
CHANLNG(I) = CHANLNG(I1) + dSO + dSl
dSO = dS2
C
C
C
IF (IPRINT .EQ. 1) THEN
XW = X(I) + XMIN
YW = Y(I) + YMIN
XcW = Xc(I) + XMIN
YcW = Yc(I) + YMIN
WRITE(6,90) I, CHANLNG(I), XW, YW, 0., 0., 0., 0.,
+ XcW, YcW, Rc(I), THETA(I)*180./PIE, PHI(I)*180./PIE,
+ 9., 0., 0., 0., 0., 0., 0., 0., 0.
90 FORMAT(I5,9F10.1,F12.0,2F10.1,F10.0,3F10.3,2F10.2,F10.0,
+ 2F10.2)
ENDIF
100 CONTINUE
Assign the radius for the last point (I = NUMPTS).
Xc(NUMPTS)
Yc(NUMPTS)
Rc(NUMPTS)
Xc(NUMPTS1)
Yc(NUMPTS1)
Rc(NUMPTS1)
C
C
C
Compute the azimuth (THETA) for the radius of curvature.
dX
dY
X(NUMPTS)
Y(NUMPTS)
 Xc(NUMPTS)
 Yc(NUMPTS)
IF (dY .GE. 0 . .AND. dX .GE. 0. ) C = 0.
IF (dY .GE. 0 . .AND. dX LT. 0. ) C = PIE
IF (dY .LT. 0 . .AND. dX .LT. 0. ) C = PIE
IF (dY .LT. 0 . .AND. dX .GE. 0. ) C = 2.*PIE
82
c
c
c
IF (ABS(dX) .GT. 1.0E05) THEN
THETA(NUMPTS) = ATAN(dY/dX) + C
ELSE
IF (dY .GE. 0.) THEN
THETA(NUMPTS) = PIE/2.
ELSE
THETA(NUMPTS) = l.*PIE/2.
ENDIF
ENDIF
IF (THETA(NUMPTS) .GT. 2.*PIE)
+THETA(NUMPTS) = THETA(NUMPTS) 2.*PIE
Compute the azimuth (PHI) for the alignment of the channel.
PHI(NUMPTS) = PHI(NUMPTS1)
RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE BANK_EROSION(NUMPTS,dt,Q,Lag,IPRINT)
C
C Routine to compute the new coordinates of the migrating
channel.
C
DOUBLE PRECISION X, Y, CHANLNG, Xc, Yc, Rc, THETA, PHI,
+B, SM, PIE, dSl, dS2
COMMON / PLAN_DATA / X(1000), Y(1000), CHANLNG(1000),
+ Xc(1000), Yc(1000), Rc(1000), THETA(IOOO), PHI(1000)
+ SEIN(1000), SEIN2(1000), XLT(1000), YLT(1000),
XRT(1000),
+ YRT(1000), XMIN, YMIN
INTEGER SEIN, SEIN2
COMMON / HYDR / DISCH, AREA, Hdepth, HdepthM, SLOPE, TDF,
NVALUE
COMMON / WIDTHS / CHwidth(1000), CHwidthMax(1000),
CHWtime(1000),
+ WIDTHM, WIDTHL, WIDTHU, WIDTHI
COMMON / HYDRAULICS / VEL(1000), Dh(1000), Dmax(lOOO),
+ TlocL(1000), TlocR(1000)
COMMON / MINIMUMS / AMPmax(20), MaxAMP(20), VSmin, Smin, Vmin
+ DmaxM, SHAPE
COMMON / SLOPES / Lmax, NC, CROSSING(20), Cslope(20),
+ PLagFactor(20), PhaseLag(20)
COMMON / EROSION / CsWRc(1000) ECOEF(1000,50) BE (1000)
COMMON / SEDDAT / DRL(15), DRU(15), GMD(15), SEDSUP, Cs, NRGS
COMMON / NAMES / TITLEXY(5), TITLEQ(5), TITLEG(5)
CHARACTER TITLEXY*80, TITLEQ*80, TITLEG*80
INTEGER CROSSING
REAL NVALUE, Lag, PLagFactor, LWD
83
COMMON / GRID / XLL, YLL, XUR, YUR, IMAXg, JMAXg, dxGrid,
dyGrid,
+ rg(101,51), rd(101,51), hb(101,51),
+ LWD(101,51) dw(101,51) PI (101,51) dc(101,51)
COMMON / BEcoef / al, a2, a3, a4, a5, a6
COMMON / CONSTANTS /PIE
C
C Hydraulic Parameters Causing Channel bank erosion:
C al = hydraulic coefficient
C Cs = Sediment concentration of the bed material load [ppm]
C Rc = radius of curvature [L]
C D = hydraulic depth [L]
C Vegetation Resisting Bank Erosion
C a2 = vegetation coefficient
C rg = vegetation root density
C rd = vegetation root depth [L]
C hb = bank height [L]
C Large Woody Debris Resisting Bank Erosion
C a3 = large woody debris coefficient
C LWD= fraction of trees or large woody debris along the
C eroding bank
C dw = average height of large woody debris jams [L]
C D hydraulic depth
C Cohesion Resisting Bank Erosion
C a4 = cohesion coefficient
C PI = plastic index
C Armoring Resisting Bank Erosion
C a5 = armoring coefficient
C dc = fraction of bank material coaser than incipient motion
C Velocity scalling factor
C a6 = velocity scalling coefficient
C V mean channel velocity
C
DATA LMAX / 50 /
C
IF (IPRINT .EQ. 1) WRITE (6,10)
10 FORMAT(' I Chanlng X Y XLT YLT
+ ' XRT YRT Xc Yc Rc ,
+ ' Theta Phi Turn RcW TlocL
+ 'TlocR Bank Eros Phase Lag CH Width HydDepth Max
Depth ')
IF (IPRINT .EQ. 1) WRITE(7,11)
11 FORMAT(' I Chanlng CsWRc al Cs
Width',
+' Radius')
C
C Initialize the arrays to zero.
C
84
no non
DO 30 I = 1,NUMPTS
BE(I) = 0.0
CsWRc(I) =0.0
VEL(I) = 0.0
Dh(I) = 0.0
Dmax(I) = 0.0
TlocL(I) =0.0
TlocR(I) =0.0
DO 20 L = 1,LMAX
ECOEF(I,L) = 0.0
20 CONTINUE
30 CONTINUE
C
C Initialize the channel slope (SLOPE) and lag distance (LAG).
C
NCi = 2
SLOPE = Cslope(NCi)
Lag = PhaseLag(NCi)
C
C Initialize the beginning and ending channel array indicies
C (Ibeg and lend).
C
C lend = NUMPTS + 1 Ibeg
D = HdepthM
Ibeg = CROSSING(1)
lend = CROSSING(NC)
JLOW = 2
C
C Compute the rate and amount of bank erosion and river channel
C migration across the floodplain.
C
DO 100 I = 2,NUMPTS1
IF (I .EQ. 108 .OR. I .EQ. 155) THEN
XXX = 1.
ENDIF
Reset the channel slope (SLOPE) and lag distance (Lag).
IF (I .EQ. CROSSING(NCi)) THEN
NCi = NCi + 1
IF (NCi .GT. NC) NCi = NC
SLOPE = Cslope(NCi)
Lag = PhaseLag(NCi)
xxx = 1.
ENDIF
Set the channel width equal to the width of the previous
time
C step or, if this is the first time step, equal to the width
C associated with the minimum unit stream power (VS).
85
nnnnnonnnnn on
C
c c c c c IF (CHwidth(I) .LE. 0.) CHwidth(I) = WIDTHM Call a routine to compute the distance from the channel thalweg to the outside river bank by use of an empirical formula. CALL THALWEG_LOCATON(NUMPTS,I,JLOW,LAG,Tloc,IPRINT)
non Determine if the channel slope is less than the minimum. IF (SLOPE .GT. SMIN) THEN
C C C supply, C p Call a routine to compute the channel width, depth, and velocity given the river discharge, upstream sediment channel slope, roughness, and bedmaterial grain size.
L. 40 CALL SolveW(Tloc,WIDTH,D,Dmx,V) ELSE IF (I .EQ. CROSSING(NCi)) THEN WRITE(*,40) NCi, SLOPE, SMIN WRITE(11,40) NCi, SLOPE, SMIN FORMAT(/,13,' Channel slope is ',E10.4,/,
+ c c c c c c c c c c 5X,'Minimum slope is ',E10.4) ENDIF ENDIF CHwidth(I) = WIDTH Dmax(I) = Dmx Dh(I) = D VEL(I) = V If the maximum channel width (CHwidthMAX) is zero, then initialize it by setting it equal to the channel width (CHwidth). IF (CHwidthMAX(I) .LT. 1.) CHwidthMAX(I) = CHwidth(I) Check to see if the channel width could be limited. IF (CHwidth(I) .GT. CHwidthMAX(I)) THEN Check to see if the channel width needs to be limited.
c c c IF (CHwidth(I) .GT. 1.10*CHwidthMAX(I)) THEN CHwidth(I) = 1.10*CHwidthMAX(I) ENDIF Update the maximum channel width.
C
86
CHwidthMAX(I) = CHwidth(I)
ENDIF
C
C If point I, along the channel path, is not within the
C computatonal reach, then skip to the end of the loop,
c
IF (I .LT. Ibeg .OR. I .GT. lend) THEN
C* IF (I .GT. lend) THEN
C* IF (I .LT. Ibeg) THEN
CsWRc(I) = 0.
GO TO 100
ENDIF
C
C Compute the bank erosion rate (BER)
C
IF (CHwidth(I)/Rc(I) .LT. 3.) THEN
CsWRc(I) = al*Cs*CHwidth(I)/Rc(I)
ELSE
C WRITE(*,*) 'WIDTH =1,CHwidth(I),', Rc(',I,') =',Rc(I)
CsWRc(I) = al*Cs*3.
ENDIF
C
C From the direction of channel curvature (left or right
turn),
C determine if the left or right bank is to be eroded. A
right
C turn means that the left bank is subject to erosion and a
C negative sign is used. A left turn means that the right
bank
C is subject to erosion and a positive sign is used.
C
C Left turn, right bank erosion, positive sign {SEIN2(I) = +1}
C Right turn, left bank erosion, negative sign {SEIN2(I) = l}
C
CsWRc(l) = SEIN2(I)*CsWRc(I)
C
IF (IPRINT .EQ. 1) THEN
WRITE(7,50) I, CHANLNG(I), CsWRc(I), al, Cs, CHwidth(I),
Rc(I)
50 FORMAT(I5,F10.2,3F10.3,2F10.2)
ENDIF
C
C Determine the coefficients of influence (ECOEF) for a
C particular radius of curvature at point, I, on the next
C several downstream points.
C
DO 70 L = 1,LMAX
IF ( (Il+L) .GT. NUMPTS) THEN
ECOEF(I,L) = 0.0
GO TO 70
87
ENDIF
C
C
C
C
C
C
C
C
C
C
C
+
c
c
c
c
(and
C
c
back
C
distance.
C
C
C
distance
C
C
2 .
C
C
C
the
C
C
C
Compute the cumulative distance (DS) from the point, I,
to the next several downstream points (I 1 + L). When
L equals 1, the point is the same as locaton I and dS = 0.
dS2 is equal to the distance between 1+1 and I.
dS = CHANLNG(Il+L) CHANLNG(I)
dS2 = CHANLNG(1+1) CHANLNG(I)
Test for an Iarray index (Il+L) that might be too large
and for points past the end of the modeled reach were the
channel length is not defined (dS < 0.)
IF ((Il+L) .GT. 1000 .OR. dS .LT. 0. .OR.
dS2 .LT. 0.) THEN
ECOEF(I,L) =0.0
GO TO 70
ENDIF
Compute the coefficient of influence (ECOEF) as the ratio
of the channel length to the lag distance (dS/Lag). The
coefficient will linearly increase from zero at point I
L = 1) to a value of 1 at channel length equal to the lag
distance. The coefficient will then linearly decrease
to zero at a channel length equal to twice the Lag
IF (Lag GE. dS2) THEN
ECOEF(I,L) = dS/Lag
When the lag distance is at least as much as the
between the points at I and 1+1, the coefficient of
influence at L = 1 will be set to half the value at L =
IF (L EQ. 1) THEN
ECOEF (1,1) = dS2/Lag/2.
ENDIF
ELSE
If the lag distance is less than the distance between
points at I and 1+1, then the coefficient of influence
(ECOEF) will equal 1 at L = 1 and zero everywhere else.
IF (L .EQ. 1) THEN
88
