NUMERICAL MODELING OF SOLUTE
TRANSPORT IN TIDAL CANAL NETWORKS
BY
RAYMOND WALTON
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
ACKNOWLEDGEMENTS
The completion of a study of this magnitude necessarily involves
the combined efforts of a large number of people. The author wishes to
acknowledge collectively the assistance of the many people who have
contributed substantially to this dissertation.
The author would like to express his appreciation to his committee
chairman, Dr. B. A. Christensen, for his guidance and support during the
course of this study. He would also like to thank other committee
members, Dr. Wayne C. Huber of the Environmental Engineering Sciences,
and Dr. Z. R. PopStojanovic of the Mathematics Department for their
assistance and input to the project.
Robert Snyder, of Snyder Oceanography Services in Jupiter, Florida,
deserves special recognition for his many fundamental contributions to
the canal research project. Bob and his wife, Bea, shared the author's
enthusiasm for canal work and provided support and encouragement many
times during the course of the investigation.
Jon Lee, presently with USGS in Bay St. Louis, Mississippi,contri
buted substantially to both the field surveys and the mathematical part
of the modeling.
Rachel Christensen prepared the computer graphics for concentration
displays and contributed substantially to the success of several field
surveys. Among the many other students who volunteered to participate in
field work, were HoShong Hou, Clark Clement, Ron Chernik, and Jim Dlubac.
A special mention for assistance far beyond that which any rea
sonable person would ask is due for the drafting accomplished by Carol
Dillard and Dave Bloomquist in the last few weeks of preparation of
this dissertation. Without their good nature and highly professional
skill this dissertation would have taken much longer to prepare.
The Northeast Regional Data Center (NERDC) computer was used for
the development of the model and all of the canal design simulations.
The author particularly appreciates and wishes to acknowledge the
friendly assistance provided by the NERDC operators, when computer work
extended late into the night and early morning hours.
The author wishes to thank Alice Moreau, Ilona Weber, and Adele
Koehler for all the typing, retyping, and many hours of tedious proof
reading.
Finally, and most importantly, the author would like to acknowledge
the love and encouragement of family and friends, particular Pat who
got him here, Rysie who got him through his qualifying examinations, and
Sue who carried him through the late night programming and Writing. In
this category, he would also like to thank Dr. Roland K. Price of the
Hydraulics Research Station, England (formerly a visiting Assistant
Professor at the University of Florida) and Dr. P. Novak of the Univer
sity of NewcastleUponTyne, England, for their support in applying for
a graduate position in the United States.
The work described was funded by the Office of Sea Grant, National
Oceanic and Atmospheric Administration under Project No. R/OE4, through
the State University System of Florida Sea Grant College Program, the
Board of Regents of the State University System, and the Board of
Commissioners of Palm Beach County.
w
TABLE OF CONTENTS
ACKNOWLEDGEMENTS. . .
. ii
LIST OF TABLES . . . . . . vii
LIST OF FIGURES . . . . . . .viii
NOTATION. . . . . . . . xiv
ABSTRACT . . . . . . . xx
CHAPTER
1 INTRODUCTION. . . . . . 1
1.1 Background Information . . 1
1.2 Research Aim . . . . 3
1.3 Philosophy of Numerical Modeling . . 4
1.4 Scope of Report. ..... . . 6
2 ANALYSIS OF NUMERICAL MODELING TECHNIQUES . . 15
2.1 Development of OneDimensional Model . . 15
2.1.1 Classification and Definitions of Canal
Network Geometry. ... . . 15
2.1.2 Model Canal Network Geometry. . .. 17
2.7.3 Hydrodynamics of OneDimensional System . 18
2.1.4 Transport Equation and Dispersion Coefficient 23
2.1.5 Boundary Conditions . . . 27
2.2 FiniteDifference and FiniteElement Approaches. 29
2.3 MethodofCharacteristics Approach . . 32
2.4 Hybrid Computer Approach . . . ... 34
2.5 Second Upwind Differencing Method. . . .. 39
2.6 Method of Second Moments .. . . 42
2.7 Comparison of Techniques . . . 46
2.8 Effect of Varying Model Parameters on OneDimensional
Mass Transport .... . . . 49
3 FIELD MEASUREMENTS. ....... . . 98
Introduction and Site Descriptions
Geometry and Tides . .
Velocity and Wind. . .
Salinity and Temperature . .
Dye Dispersion Studies . .
Page
4 THREEDIMENSIONAL HYDRODYNAMICS IN CANAL NETWORKS. ... 133
4.1 Introduction. ........ . . 133
4.2 Tidal Velocities. .. . . . 137
4.3 Wind Induced Circulation. . . . 138
4.4 Secondary Currents. .... . . 146
4.5 Density Induced Currents. . . ... 154
5 THREEDIMENSIONAL MASSTRANSPORT EQUATION. . ... 175
5.1 Transport Mechanisms. . . . . 175
5.2 Longitudinal and Lateral Diffusion Coefficients 180
5.3 Vertical Diffusion Coefficient. . . 184
6 DEVELOPMENT OF THREE DIMENSIONAL NUMERICAL MODEL . 191
6.1 Layout of Geometry. .... . . 193
6.2 Treatment of the Velocity Field . . .. 199
6.2.1 Tidal Velocities . . ... 199
6.2.2 Wind Induced Circulation . . 204
6.2.3 Secondary Currents . . ... 206
6.2.4 Density Currents . . . 208
6.3 Treatment of Dispersive Terms of the Transport
Equation. . . . . 210
6.3.1 Longitudinal Dispersion Term . . 212
6.3.2 Lateral and Vertical Dispersion Terms. .. 213
6.4 Lateral Inflows .. . . 215
6.5 Decay Coefficients . . . 217
6.6 Boundary Conditions .. . . 220
6.7 Results Presentation. .. .. . . 222
7 MODEL STABILITY AND CONVERGENCE CRITERIA . . 231
7.1 Conservative Property, Order of Accuracy, and
Transportiveness ..... . . 231
7.2 Stability Criteria .... . . 233
7.2.1 Velocity Criteria. . . . 235
7.2.2 Dispersion Criteria. . . . 239
7.3 Convergence Criteria. ... . . 242
8 MODEL TESTS AND RESULTS. ... . . . 251
8.1 Introduction. . . . . . 251
8.2 Case History #1: Results of Big Pine Key Canal
III Runs . . . . . 252
8.3 Case History #2: Results of 57 Acres Model Runs. 256
8.4 Case History #3: Results of Loxahatchee River Runs 262
8.5 Effect of Varying Model Parameters on Three
Dimensional Mass Transport. . . .. 264
8.5.1 Effect of Wind Induced Circulation . 265
8.5.2 Effect of Density Current. . . .. 269
Page
9 SUMMARY OF NUMERICAL MODELING . . . 301
9.1 Summary of OneDimensional Modeling. . .. 301
9.2 Summary of ThreeDimensional Modeling. . .. 306
9.3 Future Research. . . ... . 310
9.4 The Numerical Model as a Design Tool . .. 312
APPENDIX
A USER'S MANUAL . . . . . 314
B FLOW CHART FOR THE COMPUTER PROGRAM . . 323
C PROGRAM LISTING . . . . . 366
D COMPUTER RESULTS FOR 57 ACRES OCTOBER CASE HISTORY. .. 416
REFERENCES. . . . . . . 468
BIOGRAPHICAL SKETCH ... . . . 476
I
LIST OF TABLES
Table Page
2.1 Typical Measured Canal Parameters. . . .. 93
2.2 Comparison Between Horizontal Water Surface Assumption
and Harleman and Lee's Hydrodynamic Model. . ... 94
2.3 Standard Data Set for First Test Canal . . 95
2.4 Standard Data Set for Second Test Canal. . .. 96
2.5 Comparison of Numerical Techniques . . . 97
3.1 Field Surveys Conducted by Hydraulic Laboratory During
the Canal Design Research Project. . . .. 129
3.2 Summary of Data for South Loop, 57 Acres . . 131
3.3 Summary of Data for Loxahatchee River North Canal. .. 132
7.1 Constant Parameters for ThreeDimensional Model Test
Canal . . . . . . 250
8.1 Parameters for Big Pine Key Canal III Case History . 295
8.2 Parameters for 57 Acres Case History . . 296
8.3 Parameters for Loxahatchee North Canal Case History. 298
8.4 Parameters for ThreeDimensional Model Test Canal. .. 299
8.5 Variability Studies Using ThreeDimensional Model. .. 300
LIST OF FIGURES
Figure Page
1.1 Undeveloped Canal Network in Florida. . . 9
1.2 Developed Canal Network Takeover of Natural System
in Florida . . . . . .. 10
1.3 Example of Polluted Canal Network in Florida. . .. 11
1.4 Example of Polluted Canal Network in Florida. . 12
1.5 Criteria for Canal Design Without Consideration of the
Environment (after Snyder, 1976). . . ... 13
1.6 Suggested Criteria for Canal Design With Consideration
of the Environment (after Snyder, 1976) . . 14
2.1 Definition Drawing of Canal Network . . 54
2.2 Schematic Drawing of Trapezoidal CrossSection. .. .... 55
2.3 First Test Canal Network. ... . . 56
2.4 Second Test Canal Network ... . . 57
2.5 Typical Atlantic Coast Tide Curve . . . 58
2.6 Typical Gulf of Mexico Tide Curve . . . 59
2.7 Typical Logarithmic Velocity Profile. . . 60
2.8 Schematic Drawing of Horizontal Water Surface Assumption. 61
2.9 Schematic Representation of Effect of Numerical Dis
persion . . . . . ..... 62
2.10 Low Tide Concentration Profiles for Various Ax and At
FiniteDifference Method. .... . . 63
2.11 Schematic Representation of Convective Steps of Method
of Characteristic Methods . . . ... 64
2.12 Comparison of Concentration Profiles for Fixed Grid and
Movable Grid MethodofCharacteristics. . .. 65
Figure Page
2.13 Low Tide Concentration Profile for Various AxFixed
Grid MethodofCharacteristics. . . . 66
2.14 Schematic Drawing of Hybrid Method Test Canal . 67
2.15 Generalized Analog Diagram of Single Node of Hybrid
Model . . . . . . 68
2.16 Typical Analog Output .. . . 69
2.17 Low Tide Concentration Profiles for Various Analog
Methods . . . . . . 70
2.18 Schematic Representation of Second Upwind Difference
Method . . . . . . 71
2.19 Low Tide Concentration Profiles for Various Ax and At
Second Upwind Difference Method . . ... 72
2.20 Rectangular Distribution Approximation to Actual
Distribution. . . . . . .. 73
2.21 Schematic Representation of Conservation of First Moment. 74
2.22 Low Tide Concentration Profiles for Various Ax and At
Method of Second Moments. . . . .. 75
2.23 Comparison of Techniques' Accuracy in Modeling Pure
Convection ............... . . . 76
2.24 Comparison of Techniques' Accuracy in Conserving Mass 77
2.25 Concentration Profiles for Second Test Canal Network
Second Upwind Difference Method . . . 78
2.26 Concentration Profiles for Second Test Canal Network
Method of Second Moments. . . . . 79
2.27 Variability of Tidal Entrance Time Decay Coefficient, T 80
2.28 Variability of Tidal Amplitude, a . . . 81
2.29 Variability of Canal Length, L. . . . 82
2.30 Variability of Dimensionless Canal Length . . 83
2.31 Variability of Bottom Width, b. .. . . ... 84
2.32 Variability of Inverse Side Slope, s. . . .. 85
2.33 Variability of Lateral Inflow Rate, q, . .. 86
Figure Page
2.34 Variability of Lateral Inflow Concentration, c . 87
2.35 Variability of Mean Tidal Depth, d ........... 88
2.36 Variability of Nikuradse's Equivalent Sand Roughness, k 89
2.37 Variability of Dimensionless Dispersion Coefficient, K. 90
2.38 Variability of Low Tide Concentration Profiles for
Various Branch Canal Locations. . . . 91
2.39 Variability of High Tide Concentration Profiles for
Various Branch Canal Locations. ... . . 92
3.1 Location Map for Canal Studies. . . ... 105
3.2 57 Acres Site Plan Showing Measurement Stations . 106
3.3 Topographic (7 1/2" quadrangle) Map for Loxahatchee
River Canals ... ... . . . 107
3.4 Center Line Depths in 57 Acres System . . 108
3.5 CrossSections at K and Y in 57 Acres System. . ... 109
3.6 CrossSections and Plan View of Loxahatchee North Canal,
June 1977 . . . 110
3.7 Electromagnetic Velocity Meter SetUp . . Ill
3.8 Typical Vertical Velocity Profile in 57 Acres System. 112
3.9 Typical Vertical Velocity Profile in Loxahatchee River
System, 13 June, 1977 . . 113
3.10 Wind Recorder SetUp. . . . ... 114
3.11 Strip Chart Recording of Wind Data57 Acres System,
October, 1977 . . . . ... 115
3.12 Measured Wind Velocity, 57 Acres Canal System,
October, 1977 . . 116
3.13 Station Salinities Versus Time in 57 Acres System . 117
3.14 Station Salinities Versus Time in 57 Acres System . 118
3.15 Salinity Versus Time at Station 357 Acres System,
October, 1977 . . . . . ... 119
3.16 Salinity ProfilesLoxahatchee River System, June 1977. 120
Page
. 121
Photograph of Fluorometer and Recorder SetUp. .
3.18 Longitudinal Dye
System, October,
3.19 Longitudinal Dye
System, October,
3.20 Longitudinal Dye
System, October,
3.21 Longitudinal Dye
System, October,
3.22 Longitudinal Dye
System, October,
3.23 Longitudinal Dye
System, October,
3.24 Longitudinal Dye
System, October,
Concentration
1977 .
Concentration
1977 .
Concentration
1977. .
Concentration
1977. .
Concentration
1977 .
Concentration
1977. .
Concentration
1977 ....
Profiles57 Acres
. . . . 122
Profiles57 Acres
. . . . 123
Profiles57 Acres
. . . . 124
Profiles57 Acres
. . . . 125
Profiles57 Acres
. . . . 126
Profiles57 Acres
. . . . 127
Profiles57 Acres
4.1 Example of Typical NonGaussian Concentration Profile. 165
4.2 Definition of Coordinate System. . . ... 166
4.3 Schematic Drawing of Mean Depths in Lateral Cells of a
Vertical Layer . . . . . 167
4.4 TheoreticalWindInduced Vertical Velocity Profile . 168
4.5 Comparison Between Observed and Theoretical Vertical
Velocity, Profiles, With and Without Wdith Correction
(Nz = 0.002 ft /sec)57 Acres . . . 169
4.6 Schematic Drawing of Helical Current Induced by a Bend. 170
4.7 Comparison Between Observed and Computed Lateral Veloci
ties Induced by Bend in South Loop of 57 Acres System. 171
4.8 Schematic Drawing of a Salt Wedge Entering a Canal
Showing Definitions. ..... . . 172
4.9 Typical Salinity Profile in Loxahatchee River North Canal,
Showing Presence of a Saltwater Wedge. . . .. 173
4.10 Comparison of Observed and Computed Velocity Profiles
for Loxahatchee River Site . . . ... 174
5.1 Lateral Variation of Velocity. .. .. . . 189
Figure
3.17
Figure Page
5.2 Comparison of Saltwater Interface Stability Functions. 190
6.1a Astronomical Tides Meeting at Null Point . . 225
6.1b Astronomical Tides Recombining at Upstream Junction. 225
6.2 Schematic Layout of Canal Network Showing Features . 226
6.3 Cell Structure in Reach. .. .. . . 227
6.4 Cell Structure in Junction .. . . 228
6.5 Buffer Cell Structure. .... . . 229
6.6 Schematic Layout of Bend ..... . . 230
7.1 Schematic Canal for ThreeDimensional Model Tests. .. 246
7.2 Variation of Ax and At after 50 Tidal Cycles . 247
7.3 Variation of NLAYZ (NLAYY = 1) . . 248
7.4 Variation of NLAYY (NLAYZ = 3) . . . 249
8.1 Model Layout of Reaches and Junctions in 57 Acres System 273
8.2 Typical Longitudinal Sections and CrossSections, Big
Pine Key III, Florida [EPA, May, 1975] . . 274
8.3 Observed and Predicted Concentration Profiles for Big
Pine Key Canal III Case History. . . .. 275
8.4 Observed and Predicted Concentration Profiles for 57
Acres July, 1977, Case History . . . 276
8.5 Observed and Predicted Concentration Profiles for 57
Acres October, 1977, Case History. . . .. 277
8.6 Variation of Vertical Dispersion Coefficient With Time
and Wind Speed at Mid Point of Reach 1 (Figure 8.1),
57 Acres October, 1977, Case History . . 278
8.7 Comparison Between Observed and Computed Concentrations
for Loxahatchee River North Canal, June, 1977. . 279
8.8 Case 1W: Initial Concentrations, ci = 100 ppm, Back
ground Concentration, CRW = 5 ppmTen Tidal Cycles. ... 280
8.9 Case 1W: Initial Concentrations, ci = 100 ppm, Back
ground Concentration, cpy = 5 ppmFifty Tidal Cycles. 281
Figure Page
8.10 Case 2W: Initial Concentrations, ci 5 ppm, Back
ground Concentration, CRW = 100 ppmTen Tidal Cycles. .. 282
8.11 Case 2W: Initial Concentrations, ci = 5 ppm, Back
ground Concentration, cRW = 100 ppmFifty Tidal Cycles. 283
8.12 Case 3W: Lateral Inflow Distribution Along Length of
CanalTen Tidal Cycles. .. .. . . 284
8.13 Case 3W: Lateral Inflow Distribution Along Length of
CanalFifty Tidal Cycles . . . ... 285
8.14 Case 4W: Lateral Inflow Distribution Along Upper 1/2
of CanalTen Tidal Cycles . . . ... 286
8.15 Case 4W: Lateral Inflow Distribution Along Upper 1/2
of CanalFifty Tidal Cycles . . . 287
8.16 Case 5W: Lateral Inflow at DeadEndTen Tidal Cycles 288
8.17 Case 5W: Lateral Inflow at DeadEndFifty Tidal
Cycles . . . . . . 289
8.18 Case IS: Effect of Salt Wedge With Initial Concentrations,
c 0 = 100 ppn, Background Concentration, cRW = 5 ppm
Fifty Tidal Cycles . . . . 290
8.19 Case 2S: Effects of Salt Wedge With Initial Concentrations,
ci = 5 ppm, Background Concentration, cRW = 100 ppm
Fifty Tidal Cycles . . . . . 291
8.20 Case 3S: Effect of Salt Wedge With Lateral Inflow
Distribution Along Length of CanalFifty Tidal Cycles 292
8.21 Case 4S: Effect of Salt Wedge With Lateral Inflow
Distribution Along Upper 1/2 of CanalFifty Tidal
Cycles . . . . . . 293
8.22 Case 5S: Effect of Salt Wedge With Lateral Inflow
Distribution at DeadEndFifty Tidal Cycles . 294
NOTATION
a tidal amplitude (L)
a a4 constants of integration
A crosssectional area (L2)
A area of water surface upstream of section (L2
A(k) function defined by Equation (4.31)
b bottom width (L)
B top width (L)
B(z) function defined by Equation (4.29)
cI concentration of lateral inflow dimensionlesss)
CRW background concentration in receiving waterbody dimensionlesss)
cm concentration of rectangular distribution dimensionlesss)
c turbulent time mean concentration dimensionlesss)
c' fluctuation from turbulent time mean concentration
dimensionlesss)
c initial concentration dimensionlesss)
C Chezy coefficient (L /T)
d depth (L)
d mean tidal depth (L)
e exponential constant = 2.718
e.ijk error in cell ijk
E dispersion coefficient (L2 /T)
Ek layer averaged vertical dispersion coefficient (L 2/T)
Ek layer averaged, Richardson number dependent, vertical dis
persion coefficient (L /T)
E longitudinal dispersion coefficient (L2/T)
E numerical dispersion coefficient (L2/T)
E longitudinal diffusion/dispersion coefficient (L2/T)
E lateral diffusion/dispersion coefficient (L 2/T)
E lnvertical diffusion/dispersion coefficient (L2/T)
Z
E background dispersion coefficient (L2/T)
f Coriolis parameter (1/T)
F( ) exact solution of partial differential equation
F( ) numerical approximation to partial differential equation
F dimensionless location of center of mass in cell
m
F (n) function defined by Equation (4.38)
F2(n) smooth bed form of F4(n)
F4(n) function defined by Equation (4.39)
FI(k) layer averaged form of FI(n)
F4(k) layer averaged form of F4(k)
g acceleration due to gravity (L/T2
i number of reach (Chapter 2)
number of segment
j number of reach (Chapter 2)
number of lateral layer
k Nikuradse's equivalent sand roughness (L)
constant of proportionality (Equation (5.7))
number of vertical layer
K dimensionless dispersion coefficient (Chapter 2)
decay coefficient (I/T)
KR reach uniform decay coefficient (1/T)
K wind drag coefficient dimensionlesss)
W
length scale of turbulent eddy (L)
a characteristic length of the crosssection of a canal (L)
L length of reach (L)
length scale of convective period (Section 4.1) (L)
Ld length of decay of secondary current (L)
L length of saltwater wedge (L)
N. number of concentration distributions in cell (Section 2.6)
Nu number of upstream reaches (Equation (2.9))
N vertical momentum transfer coefficient (L 2/T)
Nz constant defined by Equation (4.33)
p variable used in Section 4.4
permissible deviation from background velocity (Equation (4.51))
Pa atmospheric pressure (M/LT2 )
q, lateral inflow per unit length (L2/T)
Q discharge (L3/T)
Q discharge at upstream section of reach (L3/T)
Q discharge defined by Equation (6.3) (L3/T)
r radius of bend (L)
r rate of production or loss of substance (1/T)
R hydraulic radius (L)
Ri Richardson number dimensionlesss)
R dimensionless width of distribution in cell
m
s inverse side slope dimensionlesss)
sL inverse side slope of left bank dimensionlesss)
sR inverse side slope of right bank dimensionlesss)
S slope of energy grade line dimensionlesss)
t time (T)
t' time since low tide (T)
T tidal period (T)
u crosssectional mean velocity (L/T)
uD dispersion velocity in xdirection (L/T)
x
upD dispersion velocity in ydirection (L/T)
y
upD dispersion velocity in zdirection (L/T)
Z
uF velocity of front of saltwater wedge (L/T)
u crosssectional mean velocity at upstream section of reach (L/T)
u* bed shear velocity (L/T)
u' turbulent fluctuation from time mean velocity in xdirection
(L/T)
u densimetric velocity dimensionlesss)
u constant defined by Equation (4.54)
u2 constant defined by Equation (4.56)
u3 constant defined by Equation (4.59)
u4 constant defined by Equation (4.64)
v lateral velocity component (L/T)
V volume of tidal prism upstream of section (L3
V transfer volume due to wind (L3 )
w vertical velocity component (L/T)
w wind speed (L/T)
w' turbulent fluctuation from time mean velocity in zdirection
(L/T)
x longitudinal distance from upstream section of reach (L)
x' distance from tidal entrance (L)
y lateral coordinate direction
z vertical coordinate direction
Greek Letters
S partial derivative operator
A dimensionless distance traveled by characteristic velocity to
reach node (Equation (2.34))
AA increase in crosssectional area (L2
AF( ) truncation terms of numerical approximation
AH head loss (L)
3
AQ change in discharge along length of reach (L /T)
At time increment (T)
AV increase in volume (L3)
Ax longitudinal spatial increment (L)
Ay lateral spatial increment (L)
AZ vertical spatial increment (L)
0AAP incremental density (M/LT2 )
E small number (Equation (7.11))
n elevation of water surface from the mean depth (L)
e angle between wind and positive xdirection or reach (degrees)
K von Karman's constant = 0.4
u kinematic viscosity (L 2/T)
S dimensionless longitudinal variable
R universal constant = 3.141593
p,ro density (M/L3)
T time decay coefficient at tidal entrance (1/T)
Txx shear stress in xdirection with respect to xdirection (M/LT2
IT shear stress in ydirection with respect to xdirection (M/LT2)
Txz shear stress in zdirection with respect to xdirection (M/LT2
S angle between interface and positive xdirection (degrees)
xviii
'i(Ri) function of Richardson number dimensionlesss)
h,, tidal frequency (1/T)
Subscripts
av average over two time layers
b bottom layer
c center of cell
f variables in freshwater layer above saltwater wedge
P vertical layer
LT low tide
m spatial mean value
max maximum value
p previous time level
s variables in saltwater wedge
t tidal variables
TE tidal entrance
w wind variables
0 node at upstream section of reach
1 note adjacent to upstream node of reach
2 node, two away from upstream node of reach
Superscripts
n time level
turbulent time mean value
value of variable at intermediate step
Abstract of Dissertation Presented to the Graduate Council
of the University of Florida in Partial Fulfillment
of the Requirements for the Degree
of Doctor of Philosophy
NUMERICAL MODELING OF SOLUTE
TRANSPORT IN TIDAL CANAL NETWORKS
By
Raymond Walton
June 1977
Chairman: B. A. Christensen
Major Department: Civil Engineering
Florida's population is growing at a rate of twice that of the
world and three times that of the United States. This rapid growth
has been accompanied by a large demand for residential dwellings in
the state's coastal zones. As a result, many finger canal networks
have been dredged.
In the early days of canal development, there were few regula
tions or design criteria, construction was often indiscriminate and
frequently led to ecosystem degradation. Recently, public opinion
has created a demand for research into environmentally acceptable
canal design.
A three year study has recently been completed by the Hydraulic
Laboratory, University of Florida, to develop a manual for the eval
uation and design of canal networks. The study included field
measurements and physical and numerical modeling of the transport
processes.
Low energy canals are commonly found along the majority of the
Gulf of Mexico and the Eastern Atlantic seaboard. They are typified
by small tidal amplitides and small water surface slopes. Measured
maximum water surface slopes are such that in most practical cases a
horizontal water surface slope can be assumed. This allows the tidal
velocity field and the dispersion coefficients to be calculated as
closed form functions of the longitudinal displacement from the up
stream section of a canal.
Traditional numerical analysis of proposed or existing canal
networks used onedimensional digital models to study the transport
mechanism. Such models can only simulate net flux through any cross
section of the flow, such as that produced by the tide or by lateral
inflow within the network. Field measurements indicate that these
mechanisms are often of minor importance when compared with wind
induced circulation, or density currents.
A further drawback to these models has often been the severe
damping of longitudinal concentration profiles due to the numerical
dispersion inherent in the approximation schemes used. A number of
solution techniques are compared to determine their effectiveness in
reducing this error, while still producing an economic model. The
effect of varying a canal network's geometric parameters is investi
gated, and a solution technique selected for an extension to a three
dimensional model.
A threedimensional masstransport model is presented based on
a modified method of second moments technique. The model incorporates
the horizontal water surface assumption to simplify the tidal component
to produce a composite velocity field. Again, dispersion coefficients
are explicitly determined.
The model has been used to simulate the results measured from dye
dispersion studies on the East Coast of Florida, and can explain the
welldocumented inability of onedimensional models to simulate observed
conditions in this type of canal network. Furthermore, in calibrating
the model on a variety of canal networks operating under widely dif
fering conditions, the parameters were remarkably uniform, indicating
that the transport mechanism is well formed.
CHAPTER 1
INTRODUCTION
1.1 Background Information
Florida's population is growing at a rate of about twice that of the
world and three times that of the United States. This rapid growth,
caused mainly by people moving into the state, has been accompanied by
a large demand for residential dwellings in the coastal zone. In par
ticular, the demand for waterfront property has led to the dredging of
many finger canal networks (Figures 1.1 and 1.2).
The state has 9,000 miles of coastline [Carter, 1974, p. 59], with
over 3,000 miles along the Atlantic Ocean, and over 5,000 miles on the
Gulf of Mexico (these figures include estuaries, bays, islands, and
sounds). This far exceeds all the other states with the single exception
of Alaska. Today, most of the residential finger canal development is on
the southwest and southeast coasts, but as the demand continues, develop
ment sites are chosen ever further northward.
In the early days of canal development, there were few regulations
or design criteria. Construction was often indiscriminate and frequently
led to ecosystem degradation, as evidenced by poor water quality and
loss of marine life and vegetation due to loss of a varied habitat
(Figures 1.31.4).
Within the last few years the public has increasingly focused its
attention on environmental issues, a point reflected in the President of
the United States' opening message in the First Annual Report of the
Council on Environmental Quality presented to Congress in August, 1970
[Council on Environmental Quality, 1970, p. v]
The recent upsurge of public concern over environmental
questions reflects a belated recognition that man has
been too cavalier in his relations with nature. Unless
we arrest the depredations that have been inflicted so
carelessly on our nature systemswhich exist in an
intricate set of balanceswe face the prospect of
ecological disaster.
The hopeful side is that such a prospect can be avoided.
Although recognition of the danger has come late, it has
come forcefully. There are large gaps in our environ
mental knowledge, but a great deal of what needs to be
done can be identified. Much of this has already been
begun, and much more can be started quickly if we act
now.
In Florida this "awareness" has focused on the coastal finger canal
networks. The demand is now for environmentally compatible designs which
either maintain the original diversity of the ecosystem, or if possible,
improve upon it. To achieve this, the need has arisen for a comprehen
sive study of the hydrodynamic and transport properties of these canal
networks to develop realistic design criteria. Such criteria should be
based on providing a varied habitat within the canal network for the
natural wildlife, and on ensuring optimum mixing and dispersion in these
low energy systems. Snyder [1976] examined present design criteria
(Figure 1.5) and suggested improvements to the procedure (Figure 1.6)
that might lead to environmentally acceptable designs. The work of this
report is intended to aid the canal designer and research scientist to
better analyze existing and proposed network to determine if such
criteria are, or will be, met.
1.2 Research Aim
The current research effort is a three year program sponsored by the
Office of Sea Grant, National Oceanic and Atmospheric Administration
(NOAA), and the Board of Regents (BOR) of the State of Florida University
System to study the hydrodynamics and transport properties of residential
finger canal networks in the coastal zone. The research has been divided
into three major areas:
a) Physical modeling
b) Numerical and analog modeling
c) Field measurements.
The aim of the research effort is to produce a canal design manual
which will assist engineers, planners, and research scientists analyze
existing or proposed canal networks with a view to meeting a reformed set
of design criteria. To this end a number of related reports and hydro
dynamic surveys have already been published [Walton, Morris, and
Christensen, 1975; Walton, Morris, Evans, and Christensen, 1975; Morris
and Christensen, 1976; Morris, Walton, Dlubac, and Christensen, 1977;
Morris, Walton, and Christensen, 1978]. The canal design manual will be
a summary of all the work undertaken on this project.
The material presented in this dissertation is really the first half
of a twopart publication. The first half of the work is an examination
of numerical methods and modeling techniques to predict the hydrodynamics
and transport phenomenon of canal networks. From this study a predictive
threedimensional model has been developed and will be outlined in the
proceeding chapters. The second half of the work [Morris, 1978] is to
take this model and, together with knowledge obtained from the physical
modeling research and field measurements, develop a design procedure.
1.3 Philosophy of Numerical Modeling
Today, there is an abundance of numerical models to predict hydro
dynamics, transport phenomenon, and water quality for almost every con
ceivable type of waterbody imaginable. A number of these models have
been used so often that the results obtained from them are considered to
be reliable; that is to say no obvious discontinuities appear in the
solution, and computer operation does not uncover faculty programming.
In fact, there are so many models that some research effect is frequently
spent in either cataloging them or making comparisons between them
[Lombardo, 1973; Grimsrud et al., 1976; Graham, 1977].
Unfortunately, however, there is some tendency to accept a model as
proven and almostly blindly, as it were, apply it to one's own, and
frequently very different, situation. Added to this, there has also been
some tendency to develop oversimplified models, or at the other extreme,
elaborate models in which complex techniques are used to approximate
terms which have only a small effect on the results, and which can usually
be neglected when one considers the quality of the input data or of the
other approximations and assumptions made.
In a recent final draft report from the Vertex Corporation, Virginia
[Horowitz and Bazel, 1977], a study was made of the investigative and
legislative procedures used in the development of what the U.S. Environ
mental Protection Agency (U.S. EPA) considered to be the six best examples
of Advanced Wastewater Treatment (AWT) plants in the United States. Their
report not only criticized the legal aspects as embedded in the Federal
Water Pollution Control Act Ammendments of 1972 and the local laws of
many states, but also the data sampling procedures and the numerical
analysis procedures,
WATERQUALITY SURVEYS are generally suspect on technical
grounds, beset with irregularities in sampling and
analysis, and naive in matters pertaining to hydraulics,
sediments, and water chemistry.
MATHEMATICAL MODELS are oversimplified and filled with
elaborate guesswork. They are intricate, abstruse
fictions. They rarely account for the principle feature
of the waterways they claim to represent, and they are
usually built from inadequate data on hydrodynamics and
water quality. [Horowitz and Bazel, 1977, p. 12]
Now, this view may be somewhat overstated, but it seems that many
people are not aware of the old adage G. I .G.O.garbage in, garbage out,
or as it could be modified, garbage in initially, always garbage out.
This statement can be looked at from another point of view. Snyder [1977]
examined the potential accuracy of a certain water quality model. He
assumed that if each constant "to be read in" was known to be within
95 percent confidence limits, and that the reaeration coefficient can
vary over two orders of magnitude, then the results from the model,
reflecting the combined additive and multiplicative effects of errors in
the constants would be accurate to within less than 0.04 percent con
fidence limitsnot a very reassuring prospect!
The problem may be, for example, as it certainly is in some models,
that the water quality module is well defined but that the hydrodynamic
module is oversimplified, or else too elaborate to be physically meaning
ful. A frequent oversimplification is to assume that the hydrodynamics
can be simulated with a onedimensional model. This is fatal when used
for the analysis of coastal finger canal networks in Florida because they
are low energy systems with small velocities and dispersion coefficients.
In these low energy systems, external influences such as the wind and
salinity gradients give rise to multilayered flows with flow reversals.
Similarly the elaborate hydrodynamic models tend to concentrate on
accurate solutions of a complex set of equations, the complexity of which
becomes redundant when one looks at how such models handle averaging
processes over confined reaches with highly variable roughness and input
data.
The aim of this research is to develop a predictive threedimensional
model which incorporates the physically important factors in canal hydro
dynamics and transport. Simple closed form solutions for the velocity
field are used instead of the full governing equations by simplifying the
latter equations based on available field data and by considering the
relative importance of terms when compared with the probable effect of
unknowns, accuracy of input data, and discretization techniques.
1.4 Scope of Report
Florida's canal networks are low energy systems, that is low veloci
ties are commonly measured and the energy available for mixing processes
is small. Because the diffusion and dispersion coefficients associated
with the mixing processes are small, the numerical dispersion produced
by the discretization techniques, used in many models for the convection
and time derivative terms of the mass transport equation, frequently
dominates the natural dispersion being modeled. The result is that con
centration profiles can be underestimated unless extremely small time and
spatial increments are used, which is computationally and economically
unrealistic.
The accuracy and stability of various types of models and numerical
techniques are examined in Chapter 2 to determine their applicability in
modeling the hydrodynamics of Florida's coastal finger canal systems. A
onedimensional mass transport model is developed, using the assumption
that the water surface in the system is horizontal to derive closed form
expressions for the velocity field and dispersion coefficients. The most
favorable scheme was chosen as a basis for extending the model to three
dimensions.
The velocity field in coastal finger canal networks is rarely one
dimensional. In fact observations, discussed in Chapter 3, often show
multilayered systems in which flow reversals occur due to wind circula
tion, density gradients, and secondary currents. However, approximate
solutions of the full governing equations of momentum and continuity are
unrealistic because of the averaging procedures used in determining
roughness effects, wind circulation, and the effects of density gradients
In Chapter 4, the governing equations are examined and simplifying
assumptions made consistent with the expected accuracy of input data and
model results. Closed form solutions of the time and spatial variables
for the various component parts of the hydrodynamics of a canal network
are derived, which may be superimposed to determine the velocity field.
The threedimensional masstransport equation, developed in Chapter
5, describes the transport of a substance in a canal network under the
influence of the velocity field and dispersion processes. Knowing the
velocity field as closed form functions of the time increment and spatial
variables, the dispersion coefficients are developed from classical
relationships giving similar functions. These coefficients are modified,
however, to take into account the fact that mixing still occurs during
periods of slack tide when the tidal velocities are theoretically zero.
Using the form of the hydrodynamics and transport processes dis
cussed in the proceeding two chapters, a predictive threedimensional mass
transport model is developed in Chapter 6 using a massincell technique.
The model forms of the boundary conditions at solid boundaries, tidal
entrances, and of the lateral inflows are presented, and are shown to be
realistic using a comparison with the onedimensional mass transport
model developed in Chapter 2, and a simple twodimensional finite element
model of flow in the receiving waterbody.
In Chapter 7, the stability criteria for the model are developed and
analyzed using a series of test runs which vary the time and spatial in
crements. The model is transportive, conservative, and is shown to be
stable and accurate. In calibrating the model, field data, obtained on
a canal network in southeast Florida, were simulated numerically by varying
the model parameters. The model was then run to determine its accuracy
by comparing the results with a second set of data in the same canal
network. These comparisons are presented in Chapter 8. A summary and
conclusions are given in Chapter 9.
Appendix A is a users' manual for the model. A flow chart of the
program is given in Appendix B, and Appendix C is a program lising.
Appendix D presents a typical model output for one of the simulations
discussed in Chapter 8.
ton:.,
mow
Figure 1.1 Undeveloped Canal Network in Florida.
V.'
Developed Canal Network Takeover of Natural System in
Florida.
t
4. V
.4
Figure 1.2 
Figure 1.3 Example of Polluted Canal Network in Florida.
;AA
iIf IN1? II
Figure 1.4 Example of Polluted Canal Network in Florida.
I=3 B
". r,,4
m z
2 = 
52 35> ^52>
< 20~ 202
u 0u
aI 2^si i
> 3itj ?S 1
F0 1 o, = z
0 z
5,j 1s0 ^ l l
2 o 2 10 00
2'z
"1^ 32i a ~ 2 5 2c 93
4 ~ ~ c 0c j m.3? 0 " ;< 
z I z
_o 0 
21ai< 22
222<
V
U
o 5002
24
02 CCI Co
<
a
0
z
u
0
0:
0
0 5
** 0 2
221,I
22 :t s23
2CC .022 20
41
C:
C
a
4
0
V
C
a
:3
5
cl
u,,
C
41
sa
U r
S oz
0< IE
i A
I1
00 :
0<
O 0 
<0 r
<0 C 0
P0< Or
5 0< 0
Qr~ ~ 5 0
0<4<0 0
S
A 
O >0
A 
O 0
01Or
0<
O 0
0)40m
0000 a
<0
I t Q
O z
0 I
0r 0 'r'
u o
0 1 r o 0 0 s 
u< o> ; < 0i 0
s i o
w C
A 0) 0
= 0 0
 i 4 C
0 z 0 w
p 0 5
oz <
j C "^"
4 0 J !S
i.. S 2 ?
II I i II I OIH_
C)
0
rJ
E
c
0
4
C
LU
J
0
o
0
CD
c
40
C
r0
c,3
&
0
C_)
(c
C'
0
40
r0
4)
00*
Q)40^
0)4
re f~
U (
Z.
rz
't 0
0 ) 0
OJC fl 00
C z
4 >
Z z
I
0I
O <
4I 0
0< ^
> u
> Z I
CHAPTER 2
ANALYSIS OF NUMERICAL MODELING TECHNIQUES
2.1 Development of OneDimensional Model
Before the development of a threedimensional model was begun,
it was decided to develop a onedimensional masstransport model to
predict the transient concentration profiles of a passive, conserva
tive substance in a simple canal network. A number of numerical
techniques were tested and the results analyzed for accuracy, stability,
and transportiveness. The various techniques used are discussed in
the next five sections of this chapter. In Section 2.7 the tech
niques are compared, their results analyzed, and a method selected
for extension to a threedimensional model.
In the remainder of this section, the geometry, hydrodynamics,
transport equation and dispersion coefficients, and boundary conditions
for the onedimensional model are developed.
2.1.1 Classification and Definitions of Canal Network Geometry
The majority of coastal finger canal networks in Florida consist
of straight, prismatic reaches with trapezoidal crosssections and
relatively uniform depths, separated by junctions (see Figures 1.1 and
1.2). At the t'.,,l/ nt, c,'., a canal network has one or more hydraulic
connections with a receiving waterbody, usually either the Intracoastal
Waterway (ICW) or the ocean, or in some cases (particularly on the
Gulf Coast) to an estuary. From the tidal entrances, a canal network
consists of an upward branching system of canals, possibly with loops,
culminating in what will be called deadends. The itad(.', of a
canal network are the upward limits of the system at which there is
either an impervious boundary, a freshwater flow over a salinity
control structure, or else at a distance sufficiently far up a river
that the hydrodynamic effects of tides are negligible. A loop is
defined as that part of a canal network in which a closed continuous
line can be drawn along the longitudinal center lines of component
reaches. The area of the water surface associated with each loop
includes not only the surface areas of the reaches and junctions
joined by the continuous line, but also the surface areas of canals
and junctions which are hydraulically further away from that point
of the continuous line which has the shortest hydraulic connection
to a tidal entrance, and which have hydraulic connections to reaches
and junctions along the continuous line (see Figure 2.1).
In this report, the positive xaxis for each component reach will
be defined from the point of the reach which is hydraulically farthest
away from the tidal entrances, along the longitudinal center line in
the direction of the ebb tide (or towards the point closest to the
tidal entrances). This will also be termed the downstream direction.
Similarly the uqtream, direction will be defined along the negative
xaxis. Facing downstream, the left bank will be the bank to the left
of the longitudinal center line of the reach, and the riaht bank
will be the bank to the right.
At a junction, which is the meeting point of two to four dis
tinct reaches, the downstrra,' reach for that junction is defined as
the reach which lies hydraulically closer to the tidal entrance.
One of the remaining reaches coming together at the junction is then
defined to be the ,q/rtion cc for that junction. When four reaches
meet at a junction, this is the middle of the three upstream reaches.
When only three reaches meet, the upstream reach is designated as that
which is a more continuous extension of the downstream reach in the
modeler's judgement. Once the upstream and downstream reaches are
defined, the positive xaxis through the junction is from the former
to the latter. Facing in this direction, a canal joining the junction
from the left is termed the left branch for that junction, and a
canal joining from the right is called the right branch.
All the above definitions are embodied in Figure 2.1.
2.1.2 Model Canal Network Geometry
To analyze some of the various numerical solution techniques
available, two simple canal networks were defined. Each reach of
the two networks was assumed to be trapezoidal in crosssection with
bottom width, b, and with the same inverse side slope, s, on each
side of the canal (Figure 2.2). The crosssectional area, A, for
any depth of water, d, is then given by,
A = d(b + sd) (2.1)
The first network (Figure 2.3) consisted of a single straight
prismatic canal of length, L, with a deadend at the upstream end and
a tidal entrance downstream. The second network (Figure 2.4) consisted
of two canals meeting at a junction with a single tidal entrance at
the downstream end, dividing the system into three reaches of lengths
L where i refers to the number of the reach. In Figures 2.3 and
2.4 the reach number is the uncircled number and the junction number
the circled number. For computational purposes, it was easier to
refer to each deadend as junction number 1, to number interior junc
tions consecutively from 2 to (NJUNC + 1), where NJUNC is the number
of interior junctions, and to number tidal entrances consecutively
from (NJUNC + 2) upward.
2.1.3 Hydrodynamics of OneDimensional System
There are several mechanisms which contribute to the hydrodynamics
of canal networks, and these can be classified into two groups which
have an internal or an external effect on the system. The internal
mechanisms are those which cause a variation within the canal network
by acting directly on the interior of the system. The external
mechanisms are those which operate on the receiving waterbody, and
which transfer their effects through the tidal entrances.
Examples of internal mechanisms are lateral inflow, secondary
currents, wind induced circulation and density driven currents, in
duced by a saline wedge entering the canal network during the flood
tide and retreating during the ebb tide, or by thermal gradients.
However, all these effects except lateral inflow, produce circulation
patterns whose net mass flux through a crosssection is zero at any
time, and hence which cannot be modeled using a onedimensional
numerical scheme. To compensate, many modelers try to reproduce the
threedimensional effects by adjusting the onedimensional model's
coefficients, but the results are rarely satisfactory as the circulation
pattern is not followed accurately enough. Only the lateral inflow, of
the above list, can be thus modeled, because it has positive net flux
in the system.
The main external influences are the astronomical tides, baro
clinically induced tides and wind tides. All these phenomena have
the effect of raising or lowering the water elevation at the tidal
entrances to the canal network, and thus directly inducing a response
within the system. In the onedimensional model, these effects can
be modeled in one of two ways. If the record of tidal elevations at
the tidal entrances to the canal network follows a simple harmonic
function, then the amplitude and frequency (or period) of this distri
bution can be read into the model, and the tidal elevations generated.
Otherwise, the record can be digitized and read directly into the
model as an input tidal entrance boundary condition.
Most hydrodynamic models are based on the de SaintVenant Equa
tions of continuity and momentum (or the dynamic equation as it is
frequently termed). In one such model [Harleman and Lee, 1969],
which will be used as a comparison with the hydrodynamic model presented
in this section, the equations are written as follows, continuity,
A q = 0 (2.2)
and momentum,
Q+ u Q + Q u + g d A + g Q = 0 (2.3)
Xt ')x +X X AC2R
where
t = time (T)
Q = discharge (L3/T)
x = longitudinal displacement (L)
q = lateral inflow per unit length of reach (L2/T)
u = crosssectional mean velocity (L/T)
g = acceleration due to gravity (L/T2
C = Chezy coefficient (L1/2/T)
R = hydraulic radius (L)
This model was used to study the hydrodynamics of a single canal with
a variety of end boundary conditions. However, a number of investi
gators have used similar models to study canal networks with multiple
branches, tidal entrances, and loops [Vreugdenhil, 1973; Abbott et al.,
1975].
A study of the hydrodynamic variables measured in some of Florida's
coastal finger canal networks indicates that it is usually not neces
sary to develop a onedimensional hydrodynamic model based on the full
set of equations. Canal networks in Florida are mostly fairly short,
less than 10,000 ft long for example. They usually consist of straight
prismatic reaches with negligible bed slopes separated by junctions.
Typically, tidal ranges on the Atlantic coast are about 2 to 4 ft
(Figure 2.5), and along the Gulf coast are 2 to 3 ft (Figure 2.6)
[NOAA, 1974]. These small tidal ranges induce very small velocities
in the networks (Figure 2.7), and is not uncommon to measure maximum
velocities of less than 0.5 ft/sec at the tidal entrances (Figure
2.7). These and other typically measured hydrodynamic parameters are
listed in Table 2.1. As can be seen from this table, maximum
measured water surface slopes have been found to be between 105
and 106.
This last observation in particular led Christensen [1975] to
suggest that the acceleration terms in the momentum Equation (2.3)
might be negligible in this case, and that a good approximation to the
flow would be to consider a horizontal water surface rising and falling
with the tidal frequency and amplitude. Thus, the continuity equation
alone would be sufficient to uniquely determine the velocity field
[Walton, 1976a]. To test this assumption, a simple hydrodynamic
model based on the conservation of mass principle was developed for
the first test canal network (Figure 2.3), and the results compared
with those obtained from Harleman and Lee's model for the same system.
For all the onedimensional modeling, it was assumed that the
depth varied sinusoidally with the tidal amplitude, a, and frequency,
w (Figure 2.8). Thus the depth, d., in any reach, i, of the canal
network, is given by,
d. = do + a cos mt (2.4)
where
d = mean tidal depth in reach i, (L)
and
2 2i (2.5)
T
where
T = tidal period = 12.42 hrs.
Such a representation is consistent with tides on the Atlantic coast
Figure (2.5).
As the discharge for reach i, Qi, can be written,
Q. = A.u. (2.6)
where the area A. is not a function of x, Equation (2.2) becomes,
dA. du.
1 + A. I q 0 (2.7)
Integrating with respect to x. from the upstream section of the reach,
xi] dA.
u (x,t) = (f ql. dx xi dTi)/A + u (2.8)
0 1 i
where
u = velocity at upstream section of the reach (L/T).
The velocity, u can be calculated in a similar manner in the upward
branches of the reach considered, giving,
Nu Lj dA. x dA.
u (x,t) [ ( q dx L dt) + q dx x dt]/Ai
j 0 J 0 1
(2.9)
where
Nu = number of upstream reaches.
Now
dA.
d = B dd (2.10)
dt dt
and
A = BL
Ws
where
A = area of water surface in canal network upstream of
a section (L 2),
and
B = top width of crosssection (L),
thus Equation (2.9) can be simplified to,
L. x. A
Nu 1 1 WSi dd
u.(x,t) = [ f J q dx + J qi dx]/Ai A dt (2.12)
0j 0 1 o i
Equations (2.4) and (2.12) are then sufficient to determine the depth
of flow, and the velocity field at any point in the model onedimen
sional canal network as a closed function of x and t.
To test the accuracy of this model, the depths and velocities
were computed at various points inside the first test canal, the
single prismatic canal (Figure 2.3), and the results compared with
the results obtained from Harleman and Lee's model for the same system.
Many of the parameters were varied, in particular, the canal length
up to 11,000 ft. In all cases the depths and velocities calculated
from the two models were within 2 percent of each other. A typical
set of comparative results for a 9,500 ft canal is listed in Table 2.2.
This was considered to be sufficiently accurate considering the
errors inherent in numerical modeling and in the practical measure
ment of field data. Equations (2.4) and (2.12) were therefore used
to describe the hydrodynamics in a onedimensional canal network.
2.1.4 Transport Equation and Dispersion Coefficient
The threedimensional masstransport equation is developed by
considering the conservation of a substance in an elemental volume
of the flow. Using Fick's first law for molecular diffusion and an
analogy for turbulent diffusion, the equation can be written
[Harleman, 1966, pp. 576578; Pritchard, 1971, p. 16],
ac + (cu) + (cv) + (cw) (E )c
Tt x y y +z xx x x y
+ (E 2c) + (E ac + r (2.13)
+y y y 2z z )z p
where
x, y, z = coordinate directions (L)
u, v, w = velocities in x, y, z directions respectively
(L/T)
Ex, Ey Ez = turbulent diffusion coefficients in x, y, z
directions respectively (L2/T)
c = concentration dimensionlesss)
r = rate of production or loss of substance (1/T).
The onedimensional mass transport equation can then be obtained
by crosssectionally averaging Equation (2.13) [Harleman, 1971, pp. 37
38] and by using a Fick's first law analogy for longitudinal dispersion
[Taylor, 1954; Aris, 1956],
S(Ac) + (Auc) = (A(Ex + E) D) + Ar (2.14)
where
E, = longitudinal dispersion coefficient (L 2/T).
Elder [1959] examined the orders of magnitude of the longitudinal tur
bulent diffusion coefficient, E x, and the longitudinal dispersion
coefficient, E and found that the latter was about ten times the
former. Thus, Equation (2.14) is often written,
S(Ac) + (Auc) (AE c) + Ar (2.15)
7t px ax x p
although measurement techniques for the longitudinal dispersion coef
ficient are integral methods and thus automatically include the effect
of the longitudinal diffusion coefficient.
As the substance being modeled was considered to be passive and
conservative, and as there were no sinks for the flow except at the
tidal entrances, the rate of production term, r includes only the
lateral inflow, q If the concentration of the lateral inflow is cI,
then the rate of production is given by,
rp = qIc1/A (2.16)
and Equation (2.15) becomes
S(Ac) + a (Auc) (AE a'c + q (2.17)
Tt ax ax x IC
The above equation is called the conservative form of the onedimensional
masstransport equation. The term conservative is used to imply that
a finite difference approximation of the equation preserves the integral
relationship of the continuity equation [Roache, 1972, pp. 2833]. A
noncoznsZCrvativrc form of the equation does not have this property.
Roache states that a conservative formulation does not always give
more accurate results, but that the experience of many researchers has
indicated that such formulations generally are more accurate. A non
conservative form of Equation (2.17), using the continuity Equation (2.3)
and the fact that the area, A, is not a function of the longitudinal
displacement, is,
cc + u a I (E +c) c)/A (2.18)
Dt ix x E ax qIc
It will be seen later that it is sometimes necessary to use this form
of the equation.
Classically, the longitudinal dispersion coefficient, E is
written [Taylor, 1954; Elder, 1959],
E = KRu* (2.19)
where
K = dimensionless dispersion coefficient
R = hydraulic radius (L)
u* = bed shear velocity (L/T).
The dimensionless dispersion coefficient, K, is a function of the
geometry of the system, depending on radii of bends, bottom roughness,
wall Reynolds' number, etc. Although classically considered to be
about twenty (based on the hydraulic radius), many researchers have
measured values of K in various types of waterbodies, from less than
ten to the order of several hundreds.
In many Floridian canals, logarithmic velocity profiles have been
observed when the effects of wind induced circulation and density
driven currents have been small (Figure 2.7). Such a profile gives
a unique relationship between the spatial mean velocity, u, and the
bed shear velocity, u* [Nikuradse, 1933],
u = 2.5u* ln(10.9d/k) (2.20)
where
k = Nikuradse's equivalent sand roughness (L).
Thus, from Equation (2.19)
F 0.4KRu (2.21)
n (10n.9gd/k)Y
As the spatial mean velocity, u, is a closed form function of
the longitudinal displacement and time, the longitudinal dispersion
coefficient is also a closed form function of x and t. Therefore,
the only remaining unknown in the onedimensional masstransport
Equation (2.17) or (2.18) is the concentration, c.
2.1.5 Boundary Conditions
The onedimensional masstransport Equation (2.17) or (2.18) is
a secondorder parabolic equation, being secondorder in the spatial
variable, x, and firstorder in the time variable, t. Thus, to close
the system, two spatial boundary conditions and a set of initial con
ditions are required.
For both the test canal networks, the initial conditions were
that the concentrations at time t = 0 in the canal either equaled the
background concentration, c RI, or else that they were read in.
At the deadends of the canal system, a zeroflux condition was
initially used in the form,
ac 0 (2.22)
However, it soon became apparent that simple numerical approximations
to Equation (2.22), such as a forward difference operator,
o 0 (2.23)
where
c = concentration at deadend dimensionlesss)
cI = concentration at adjacent node dimensionlesss)
implied that the gradient of the concentration profile in the segment
of the reach adjacent to the deadend is zero, which led to erroneous
results. Higherorder schemes were discounted because it was felt that
they used information which was too far away from the deadend to be
meaningful.
To overcome this problem, a simple mass balance expression was
developed based on the mean velocity between the two nodes at the
deadend. However, this condition was eventually simplified when it
was observed that the concentration profiles resulting from this model
were linear at low tide, and that the mass balance condition maintained
this linearity all the way to the deadend. Thus, the final form of
the deadend boundary condition for this test model was that the
concentration at the deadend was a linear extrapolation of the con
centrations at the twoadjacent nodes,
cQ = 2c1 c2 (2.24)
At the tidal entrances to the system, a dual condition was estab
lished. During the ebb tide a "floating" type boundary condition was
used in which the concentration at the entrance was calculated using
a backward difference operator. Dispersion was neglected at this
point because the concentration would either be convected out of the
system at the next time step, or else would be lost in the description
of the flood tide boundary condition. Thus during the ebb tide, the
n+l
concentration at the tidal entrance, cTE was given from the
equation
n+l n n n n n
TE t TE uTE CTE TE1 TEl (2.25)
At Ax
where
At = time increment (T)
Ax = spatial increment (L)
n = time level
Assuming that the concentration at the tidal entrance reaches a
values of cLT at low tide, a first order decay was used to describe the
time history of the concentration at the tidal entrance during the
flood tide,
cTE = cRW + (cLT cRW) exp (3t'/T) (2.26)
where
t' = time since low tide (T)
T= time decay coefficient (1/T).
The form of this boundary condition results in the concentration at
the tidal entrance being within 2 percent of the background concen
tration, cRW, after T units of time after low tide.
To test the accuracy of this assumption, a simple twodimensional
finiteelement model was developed to describe the transport in a
receiving waterbody such as the ICW under a point source/sink
loading. Results indicated that a first order decay type condition
was a fair approximation to the resulting flood tide concentration
profile at the tidal entrance. Furthermore, as will be seen later in
Section 2.8, variation of the time decay coefficient, T, does not
dramatically alter the resulting concentration profiles in the canal
network.
2.2 FiniteDifference and FiniteElement Approaches
Measured longitudinal dispersion coefficients in Floridian canals
are usually very small, less than 5 ft 2/sec being common. This means
that the dispersion term in the onedimensional mass transport
Equation (2.17) or (2.18) is small. If this term is dropped altogether,
then the remaining firstorder hyperbolic equation is the onedimen
sional convection equation,
c + u 0 (2.27)
A number of people have studied this equation [von Neumann and
Richtmyer, 1950; Stone and Brian, 1963], and have concluded that model
studies of this equation using the more common low order finitedif
ference approximations usually lead to inaccurate results, and possibly
instabilities, depending on the scheme used. A survey of the more
common finitedifference schemes such as the forward, backward, central
explicit or implicit schemes, will not be given here as they are ex
tensively covered in the literature [Ames, 1969; Milne, 1970; Roache,
1972; Smith, 1975].
The problem, termed numerical dispersion, arises out of a Taylor
series approximation to the terms of Equation (2.27). If, for
example, a constant velocity, u, is assumed in the positive x
direction and a backward difference scheme is used to approximate the
equation to the second order [Molenkamp, 1968],
t ,+ u c : lu (1 1ul ) xA (2.28)
ux
the term on the right hand side of Equation (2.28) is the numerical
dispersion term with an associated numerical dispersion coefficient
En r u (1 u ) x (2.29)
n lul ( lul Ax
The effect of this pseudodispersion is to smooth the profile for positive
values of En, and to produce instabilities in the results for negative
values. This same problem can be shown diagramatically (Figure 2.9)
for the case [Pedersen and Prahm, 1973],
u 2 ax (2.30)
which gives a positive numerical dispersion coefficient.
From Equation (2.28), it can be seen that one way of reducing
the error is to decrease the spatial increment. However, this tends
to increase computer costs to astronomical levels as the time inter
val, at, usually has to be reduced also to meet stability criteria.
A second method would be to choose Ax and At such that
u = x (2.31)
At
and this has been done successfully for some steady, uniform flow
problems. In the tidal canal problem, however, the flow is unsteady
and nonuniform, and the condition of Equation (2.31) cannot be met
everywhere. In fact unless this condition is met at the tidal entrance,
E will be negative somewhere in the canal network in this example.
The above analysis was based on the backward difference scheme,
but serves to illustrate the problem. Similar numerical dispersion
expressions can be developed for other schemes for idealized flow
conditions. The aim of this study then, is to develop a scheme in
which the inherent numerical dispersion present in the numerical
solution is much smaller than the natural dispersion being modeled,
and also one which does not require excessive computer time because
of the smallness of the spatial and time increments, Ax and At,
respectively.
A number of schemes were developed by this author [Walton, 1976a],
forward, backward and central explicit, and backward and central
implicit, as well as by Langley [1976] who used an explicit central
difference scheme to model flow in a canal with a boat basin. For the
first test canal network (Figure 2.3), using the standard data set
listed in Table 2.3, it was found in every case that the only way to
improve the accuracy of the schemes was to reduce Ax, and therefore
at (Figure 2.10). For the 1,000 ft canal (Table 2.3), it was necessary
to make Ax = 25 ft and at = 100 secs.
In conclusion it was decided that this type of restriction on
the size of the spatial and time increments would serve to produce
an uneconomic model for much larger canal networks and for a three
dimensional model. A similar finiteelement model developed by Leimkuh
ler et al [1975] was also tested, and was found to suffer from the
same drawbacks as the finitedifference models. This latter method
was therefore also discounted.
2.3 MethodofCharacteristics Approach
As was seen in Section 2.2, the numerical dispersion inherent in
the more common low order finitedifference and finiteelements
schemes is sufficiently large to swamp the natural dispersion being
modeled, for economic choices of the spatial and time increments.
The problem is then to model the convective part, Equation (2.27), of
the onedimension masstransport equation economically and accurately,
and minimize the numerical dispersion produced.
Equation (2.27) is a firstorder hyperbolic equation. A tradi
tional way to solve this equation is the MethodofCharacteristics.
If the onedimension masstransport equation in its nonconservative
form, Equation (2.18), is considered, the left hand side of the
equation can be expressed as a material derivative,
Dc Ac uc I Ac c
= t + u x x (E ( c)/A (2.32)
A two step solution scheme can then be developed by firstly con
vecting the concentration field along the characteristic lines of the
velocity forming new nodal concentrations cn where i represents a
spatial node and n+1l is the new time level. Equation (2.32) is now a
secondorder parabolic equation with no dominating firstorder hyper
bolic term, and can readily be approximated by a centraldifference
scheme, in which the time derivatives are approximated along the
characteristic line of the velocity, and the spatial derivatives
are estimated by averaging between the two time levels,
n+1 cl 1 + + n+l n+
c 1 c. + 2 [(Ei+l + Ei)(c cn )
c 2Ax2 i avi+l avi
(Ei + E i )(c av c av l
1 i av. av1
+ Atq (C cavn+l)/An+l (2.33)
i
where
c n+l = 1 n+l n (2.34)
av = ) (2.34)
in which
A = dimensionless distance traveled by characteristic
velocity to reach node i.
The convective step can be performed in two ways as illustrated
in Figure 2.11. The first method is to convect the concentrations at
the existing nodes to form a new nodal grid at the next time level.
This is called the movable grid formulation. The second method is to
interpolate between nodal values using the velocity characteristic
lines to estimate where new nodal concentrations originated. This
is called the fixed grid formulation.
The two variations were run with the standard data set of Table
2.3 for a variety of spatial and time increments. The results shown
in Figures 2.12 and 2.13, and others, are analysed in Section 2.7.
2.4 Hybrid Computer Approach
If the convective Equation (2.27) is again considered, a function
F(x,t) is defined as follows,
F(x,t) = c + u 0 = 0 (2.35)
for a uniform velocity field, u, and F(x,t) is defined to be the
finitedifference approximation to F(x,t), then
F(x,t) = F(x,t) + AF(x,t) (2.36)
where
AF(x,t) = truncation terms of finitedifference series
expansion.
For a central difference scheme of the spatial derivative, and a
forward difference scheme for the time derivate, using Taylor series
expansions up to the second order,
2
AF(x,t) = At (2.37)
Sat
as the central difference scheme is second order accurate. Differen
tiating Equation (2.35) with respect to t and noting that the dif
ferentiation is commutative,
=)x () (2.38)
at
Using Equation (2.35) again,
P2c 2 22c
u 2 (2.39)
at 'x
and substituting into Equation (2.37)
2
AF(x,t) = u2 (2.40)
at
Thus, for a central difference scheme, an expression has been derived
for the numerical dispersion coefficient,
it 2
En = u (2.41)
This form of the numerical dispersion coefficient illustrates
an intriguing possibility. If the time increment could be decreased
to zero then the numerical dispersion would vanish irrespective of the
spatial increment, Ax. Although Equation (2.41) does not accurately
represent the numerical error present in modeling a nonuniform velocity
field, it does at least suggest its basic form.
One way to decrease the time interval to zero, or make time con
tinuous, is to use an analog or hybrid computer. This consists of a
digital computer controlling the potentiometer settings of an analog
computer. A project was initiated with the Hybrid Computations
Laboratory at MartinMarietta Aerospace in Orlando, Florida [Morris,
Walton, and Christensen, 1977; Sorondo and Baldwin, 1977], to study
the feasibility of applying analog modeling techniques to the transport
in a onedimensional canal network. A finitedifference model was
also designed for the same system to verify the model.
An analog model has several advantages over the digital models
which at that time were being considered for canal networks. Firstly,
analog solutions are obtained very quickly, and many hundreds of tidal
cycles could be simulated in less than a minute. On a digital computer
this would take several minutes of central processing unit (CPU) time,
plus a much longer time to print out the results, Secondly, because
of the hybrid capability, many test runs can be run in a very short
time, allowing a much more complete look at the effects of varying
each parameter alone.
The objectives of the hybrid model study were
a) To determine whether the hybrid approach and the Martin
Marietta facilities were capable of implementing a two
dimensional simulation of pollution transport in a tidal
canal network.
b) To evaluate the hybrid model by comparison with digital
models being developed by the Hydraulic Laboratory.
c) To determine whether the nonconservative or the conserva
tive formulation of the masstransport equation was
preferred.
d) To determine the most consistent finitedifferencing
techniques.
e) To determine which of the boundary conditions provides
the most realistic results consistent with other models
using difference techniques.
f) To assess the relative influence or significance of
variations in the independent variables of the analog
model.
The analog model was set up as the first test canal (Figure 2.3)
with five segments in the canal and a sixth segment in the receiving
water to control the tidal entrance boundary condition (Figure 2.14).
The model, which included the capability of switching between four
different explicit finitedifference techniques, the forward, backward,
central, and first upwind differencing schemes (the latter being simply
a backwarddifference scheme for the ebb tide, and a forwarddifference
for the flood tide), required a large percentage of an analog
patchpanel, sometimes called a hybrid terminal. An example of a
circuit diagram for one node of the canal is given in Figure 2.15.
It was estimated that a total of twenty nodes could be programmed onto
one panel. The output was in hard copy form on an eight channel strip
chart recorder, the eighth channel giving the velocity at the tidal
entrance (Figure 2.16). The results for the first test canal case for
each of the four different finitedifference techniques are given in
Figure 2.17.
After the work was completed a more detailed report was written
[Morris, Walton, and Christensen, 1977] and submitted to the Office of
Sea Grant. The value of this modeling technique and its possible
extension to threedimensional modeling will be covered in a comparative
discussion of all the techniques investigated in Section 2.7. The
conclusions are given below.
1. The central difference form of the nonconservative
equation appears to give the most consistent results,
both in the analog and the digital models. The
capability of obtaining a set of concentration plots
over a thousand tidal cycles in a matter of seconds
(after the model has been initialized) is very
convenient.
2. The conservative form of the equation has potential,
but its results were obscurred by an unrealistic
boundary condition and some instabilities. Further
testing of this form would be required before a
choice between them can be made.
3. The boundary condition at the openend of the canal
is satisfactory.
4. The zeroflux (3c/;x =0) boundary condition at the
closedend of the canal appears to distort the results.
However, the linear boundary condition gives results
which are in very good agreement with other models
developed.
5. The mass balance boundary condition at the deadend
was not functioning properly, although it has been
shown to operate properly in digital models developed
by the Hydraulic Laboratory. Further testing of
this formulation would be desirable.
6. The fivesegment model is somewhat coarse spatially,
but better than was initially expected. A ten or
twentysegment model would be useful for verifying
the results of the fivesegment model, but could
not be expected to provide much additional
information.
7. Flushing time, the time for the model to reduce
concentrations at all sections to "steadystate"
values from some initial, nonequilibrium condition,
cannot be inferred from the present model results
due to inconsistencies in the variability portion
of the study. However, these inconsistencies may
be due in part to difficulty in accurately reading
plotted results.
8. Expansion of the model to about twenty sections is
possible using one analog programming board on the
hybrid computer. A twodimensional model incorporating
three nodes in the vertical direction would therefore
require three analog consoles. Multipleboard
simulations are difficult to schedule and expensive to
run. Therefore, it is concluded that for canal
networks and twodimensional models the analog/
hybrid approach is extremely limited by hardware
availability, not generally practical, and quite
expensive.
9. It typically would take several days to a week to have
a specific set of analog or digital computer runs
completed and mailed to the Hydraulic Laboratory,
although it is possible to greatly reduce this time
if expense is secondary.
10. Numerous refinements, as outlined in the initial proposal
were not tested as they required more hardware and
would have obscured the basic evaluation. These refine
ments included variable canal geometry and additional
water quality parameters. Incorporation of these
additional variables, which would ultimately be re
quired for a complete canal design capability, would
be prohibitive from an analog hardware standpoint.
11. The basic canal neocork' formulation was not tested due
to lack of a good closedend boundary condition and
lack of time. A single canal junction in onedimension
could be implemented without much difficulty on the
analog computer, but extension to a full canal network
is definitely limited by analog hardware availability.
2.5 Second Upwind Differencing Method
The "Second Upwind Differencing" method, as termed by Roache
[1972] and used by Bella and Dobbins [1968], or "the Donor Cell"
method as it was called by Gentry, Martin, and Daly [1966], is a
finitedifference technique used by Lee [1977] to solve the one
dimensional masstransport Equation (2.17). Written in its conser
vative form in terms of the discharge, Q, rather than the cross
sectional mean velocity, u,
(Ac) + x (Qc) = x (AE c) + qc (2.42)
Using the schematic canal of Figure 2.18 as a guide, the method uses
a backward finitedifference expression for the convective term (Qc)
during the flood tide, and a forward finitedifference term during the
ebb tide.
The Second Upwind Differencing method has several advantages over
other finitedifferencing techniques. Firstly, it is conservative,
that is the net accumulation of mass in a region equals the net flux
across the boundaries of the region during the convective process.
Secondly, it is transportive during the convective phase, that is "the
effect of a perturbation in a transportive property is advected only
in the direction of the velocity" [Roache, 1972, pp. 6772]. No
stationary finitedifference formulation of the oscillatory convective
process possesses this property.
The method suffered originally from the major disadvantage that
numerical dispersion was large, as mentioned by Roberts and Weiss
[1966]. This problem was examined for the First Upwind Differencing
method by Noh and Protter [1963] who derived an expression for
the numerical dispersion present in the solution. The theory was
extended here to the Second Upwind Differencing method by Lee [1977].
Bella and Dobbins [1968] had suggested that the numerical dispersion
might be subtracted from the natural dispersion to give a composite
dispersion term that accurately modeled prototype conditions. However,
in their studies on rivers and estuaries, the natural dispersion
was much larger than the numerical dispersion, giving positive values
for the composite dispersion. The problem is that this is not the
case in Floridian canals, in which the resulting composite dispersion
is negative for reasonable choices of Ax, giving rise to unstable solu
tions. As mentioned previously, a negative dispersion term has the
effect of augmenting instabilities present in the numerical scheme
rather than smoothing them as positive dispersion does.
Recently, Boris and Book [1973, 1976] and Book, Boris, and Hain
[1975] developed a correction technique termed "fluxcorrected trans
port." This technique extended the idea of Bella and Dobbins in sub
tracting the numerical dispersion to include the extra condition that
the antidiffusion stage should generate no new maxima or minima in the
solution, nor should it accentuate already existing extrema. To
accomplish this, the antidispersion flux terms are corrected or limited
so that no antidispersion transfer of mass can push the concentration
at any grid point beyond the concentration value at neighboring
points.
Using the schematic canal of Figure 2.18 as a guide, the Second
Upwind Differencing model with fluxcorrected transport may be de
scribed as operating as follows during one time step:
1. addition of pollutant through lateral inflow,
2. net tidal convection of pollutant into the segment
from adjacent segment with values of concentration,
c, defined at the center of the segment, and values
of the discharge, Q, defined at boundaries of the
segment,
3. limited antidispersion to correct for the numerical dis
persion errors introduced in the convective step but
controlled by the fluxcorrected transport criterion,
4. net transport of pollutant into the segment by natural
longitudinal dispersion.
This method is essentially secondorder accurate, as it has been
corrected for the secondorder error, the numerical dispersion. Even
with the correction term, the method remains both conservative in all
steps and transportive in the convection step, properties which
cannot be claimed in all cases by the other finitedifference and
finiteelement models so far considered. These properties alone do
not ensure that the desired solution will be obtained. The stability
and convergence of the solution must also be considered.
Although there is no universally accepted definition of stability,
this term is generally used to refer to the damping with time of
small perturbations (or errors) introduced into the finitedifference
solution. It can readily be seen from Figure 2.19 that the solution
shows no visible instabilities to be present. One is also concerned
about convergence, or whether the solution approaches the solution
of the partial differential Equation (2.17) as Ax and At approach
zero.
For equations such as the onedimensional masstransport Equation
(2.17), with variable coefficients, necessary and sufficient stability
and convergence conditions are not available. However, conditions
such as
uma < x (2.43)
max ,t
and
Ax2
E a xAt (2.44)
make sense physically and provide guidelines for choosing Ax and At.
In the final analysis, comparisons of the numerical solutions of
the difference models must be relied on in attempting to assess the
accuracy and, considering computing costs, the economics of any one
model. This comparison is made in Section 2.7.
2.6 Method of Second Moments
The Method of Second Moments is essentially an upwind finite
difference scheme. However, instead of considering a concentration
represented at a single point, or node, the concentration is considered
to be a block with the height being the magnitude of the concentration
and the width being the width of the cell, Ax, the spatial increment.
Thus, like the nodal finitedifference representation, this method
suffers equally from numerical dispersion.
The advantage of this way of representing the concentration is
that it conserves the mass, or the zeroth moment of the mass. One
method that has been used in meteorology to reduce the numerical
dispersion present in finitedifference schemes, is to also conserve
higher order moments of the distribution [Egan and Mahoney, 1972;
Pedersen and Prahm, 1973]. Any number of moments could be chosen to
represent the distribution, but commonly the first and second moments
are used. These physically represent the center of the mass, and
the width of mass with respect to its center. Higher order moments,
such as third order moments which represent the skew of the distribu
tion, do not have the same intuitive feel, but could be used if the
lower order moments fail to represent the distribution accurately.
Consider the concentration distribution, c(E), shown in Figure
2.20 in a cell of unit width, represented by a rectangle having the
same mass, cm, the same center of mass, Fm and the same second moment
squared, R2 [Pedersen and Prahm, 1973], then,
m
0.5
cm = f c(E)de (2.45)
0.5
where
x (iAx (2.46)
and
0.5
Fm c(t)td!/cm (2.47)
m 0.5
2 0.5 2
R2 = 12 c()( F )2 dt/cm (2.48)
0.5
Figure 2.9 illustrates what happens if only the convection of a square
wave as defined by Equations (2.27) and (2.30) is considered. This
form conserves the zeroth moment of the distribution, or its mass. If
the first moment of the distribution is also conserved (Figure 2.21),
it can be seen that the analytic distribution is much more closely
approximated. Inclusion of the second moment for this particular
problem reproduces the square wave exactly, and is another reason for
considering only moments up to the second order.
The problem defined in Equations (2.27) and (2.30) is a uniform
flow case. For an oscillatory flow in canals the method has to be
amended. A closed form solution for the velocity at any point in
the canal was developed in Section 2.1.3 (Equation (2.9)). Assuming
that the lateral inflow entering each cell has a uniform longitudinal
variation in that cell, then the velocity variation over the length
of the cell is linear in x and represents the distortion the cell
undergoes as the tidal elevation changes from the old crosssectional
area, A to the new area, A. As each reach of the canal network is
prismatic, the ratio A /A represents the amount of expansion or con
traction of each cell in the reach. Thus, the first step in the
solution procedure is to distort the cell,
c = c A /A (2.49)
mp p
F = F (2.50)
R = RAp/A (2.51)
where
p = previous time level.
The second step is to move this new distribution in the direction of
the flow a distance equal to the distance moved by the center of mass
in the time interval, At
x = Ax(1/2 + Fm) V/A (2.52)
where
V = incremental volume of tidal prism upstream of cell (V is
defined as positive for a flood tide). (L 3).
The new distribution and the lateral inflow are then allocated to
the various cells into which they fall, where for each cell the
2
square of the second moment, R of each composite distribution is
c.
calculated about the center, c, of cell i,
N. 0.5
R = 12 1 c .()2 dF/c
i j=l 0.5 i (2.53)
where
N. = number of distributions that lie inside cell i.
When all the cells of each reach have been readjusted, the square of
the second moment of the resulting distribution about the new center
of mass, F m., is found using the parallel axis theorem for each cell,
R 2 = 12 F2 (2.54)
mi ci mi
Finally, dispersion is modeled at the new time level using the
central finitedifference scheme as before.
n+1 nl At E I E n+l n+l
ci = c 2A 2 [(Ei+ i+l c )
(E + E 1)(c + c n 1)] (2.55)
where
ci = is the concentration value at node i after the
convective step.
As this method is an upwind differencing method, the method also
possesses the properties of being conservative and transportive as
defined in Section 2.5. Similarly, the same stability and convergence
conditions must be expected to apply to this method. The results of
varying the spatial and time variables, Ax and At, for the first test
canal (Figure 2.3) using the standard data set of Table 2.3 are shown
in Figure 2.22. A comparative analysis of this and the other methods
described in the proceeding sections follows.
2.7 Comparison of Techniques
In Sections 2.2 through 2.6, several numerical techniques to
evaluate the onedimensional masstransport Equations (2.17) or (2.18)
were discussed, and the results shown for the simulation of the first
test canal, Figure 2.3, using the standard data set in Table 2.3.
Several methods are immediately discounted. The finiteelement
and finitedifference methods up to second order were discounted
because of the restrictive choice of the spatial and time increments,
Ax and At, respectively, to achieve comparable results with the other
models investigated. An extension to a complex canal network would
require an excessive amount of computer time. Higher order methods
were not considered in this analysis because they require small Ax and
At increments so that the information being used to model conditions
at one node is not spread out over too great a distance. The movable
grid methodofcharacteristics approach was also discounted because it
was found that the net convection out of system due to the lateral
inflow, also convected the grid out of the system over time.
The other four methods, the hybrid method, the fixed grid method
ofcharacteristics, the second upwind differencing method, and the
method of second moments, produced very similar results for the first
test canal with the standard data set. In obtaining this result, the
hybrid model used seven nodes. At the MartinMarietta facility, it
was estimated that twenty nodes could be used on one patchpanel and
that three patchpanels could be linked together in series. However,
this would tie up the laboratory and provide only sixty nodes which are
clearly not enough for a model of any complexity. Several ideas were
considered in which the solution domain could be divided up, using the
digital computer to control the operation, and the solution matched at
boundary nodes. These ideas were dropped when it became apparent that
the complexity of such a setup would be unrealistic. Also, once a
patchpanel is set up, the geometry becomes fairly rigid, and it is
not easy to program for varying spatial increments and numbers of
nodes per reach. For all these reasons, it was decided that the dis
advantages outweighed the primary advantage of speed, and further
research was discontinued.
The three remaining models are cheap to run and gave very simi
lar results for the first test canal. Therefore, a second run
was designed in which there was initially a square wave of a substance
of concentration twenty units located for x in [400, 600] at low tide.
The background concentration was five units. After ten tidal cycles,
the resulting profiles are shown in Figure 2.23. From this figure,
two things can be seen. Firstly, the effect of the limited antidis
persion is shown for the two runs of the second upwind differencing
scheme. Secondly, the profile for the methodofcharacteristics is
severely attenuated and dispersed. The reason for this is that the
interpolation procedure for redefining nodal values cannot handle the
steep gradients of the leading edge found in the square wave.
This can be shown again if another test is performed with the
first test canal in which convection only is modeled. The initial
concentrations in the canal are background, and a constant lateral
inflow load is applied uniformly along the length of the canal. The
resulting profiles after 290 tidal cycles are shown in Figure 2.24.
Once again, it is clear that the interpolation procedure cannot handle
the sharp gradients of the distribution at the tidal entrance, and
thus the methodofcharacteristics solution technique was also dis
counted here. It can also be seen from Figure 2.24 that the method
of second moments is the most successful method in conserving mass.
From Figure 2.23, and theoretically also, as explained in Section
2.6, the method of second moments convects the square wave very
accurately. Added to this, the second upwind differencing method's
limited antidispersion step is somewhat artificial and has a fairly
complex form in junctions. To examine the effect of junctions, a
second data set, Table 2.4 was developed for the second test canal,
Figure 2.4. The results are shown in Figures 2.25 and 2.26. The
results show that the junction in the second upwind differencing model
is causing a change in the gradients of the profile at that point.
Taking all the above factors into account (Table 2.5), it was
decided that the method of second moments, because of its extreme
accuracy in modeling the square wave, its superior ability to conserve
mass, its intuitive relations to physical properties of the distribu
tion, its ease of computation, and its relative economy, was ideally
suited for an extension to a threedimensional masstransport model of
a coastal finger canal network.
2.8 Effect of Varying Model Parameters on OneDimensional Mass Transport
In the previous sections, a onedimensional masstransport model
has been presented. After an extensive investigation of some of the
numerical solution techniques available, it was decided that the
method of second moments gave the most accurate results, was stable,
convergent, transportive, and was also fairly economical to run.
It has been pointed out in the introductory chapter that the
physical phenomena associated with mass transport in a coastal finger
canal network are fully threedimensional processes under the influence
of wind, salinity gradients, tides, secondary currents, and lateral
inflows. Their features have been incorporated into a threedimensional
model and the results generated using the numerical scheme selected in
this chapter.
However, a study of the variability of the onedimensional model's
parameters is useful here because the effects will be similar in
a threedimensional model, the only difference being the actions
of all the other effects modeled. These may be studied independently
using the latter model. The variation of these parameters, then,
will give a qualitative guide to design engineer as to the relative
effects of design elements. Once a canal network is designed, its
quantitative functioning can be evaluated using the threedimensional
model. All the variability runs unless otherwise stated used the
first test canal network (Figure 2.3) and varied the standard data
set (Table 2.3).
The parameters may be divided into two groups called non
controllable and controllable parameters. The first category includes
the tidal amplitude, and the decay coefficient, T, associated with
the flood tide, tidal entrance boundary condition. From Figure 2.27,
it can be seen that the decay coefficient, T, has a surprisingly small
effect on the equilibrium concentration profile in the canal over the
range of half a tidal period. This means that the inaccuracies
inherent in the selection of such a flood tide boundary condition
do not have a significant effect of the mass transport in the canal
network. On the other hand, the results from varying the tidal ampli
tude, a (Figure 2.28), show a much more dramatic effect as would be
expected from a simple volumetric tidal prims analysis. Areas with
low tidal ranges, and thus automatically low energies, could expect
severe pollution problems at the deadends of canal systems, unless
either suitable flowthrough systems can be designed, or other sources
of energy, such as the wind, can be utilized to improve flushing.
The second category, called the controllable parameters, in
cludes the geometric parameters of length, width, mean tidal depth,
and inverse side slope, as well as the lateral inflow rate and its
concentration, Nikuradse's equivalent sand roughness, and the dimen
sionless dispersion coefficient. The geometric parameters, with the
exception of the mean tidal depth, give results that would appear to
be intuitively obvious. The effect of increasing the length of the
canal, L, was simply to linearly increase the concentration at the
deadend (Figures 2.29 and 2.30), whereas the effect of increase of
the bottom width, b (Figure 2.31), or the inverse side slope, s
(Figure 2.32), for constant lateral inflow and concentration, was to
proportionally decrease the equilibrium profile as a greater volume of
water is presented to dilute the inflow. Conversely, the effects of
increasing the lateral inflow, q,, while holding its concentration
constant (Figure 2.33) or holding the inflow rate constant and in
creasing its concentration, c, (Figure 2.34), for fixed canal geometry,
proportionally increased the resulting equilibrium profile.
The effect of increasing the mean tidal depth, do, has a dif
ferent effect (Figure 2.35). The other geometric parameters merely
increased the size of the tidal prism and thus the velocity, for a
fixed mean depth. The variation in the mean tidal depth alters
the velocity of the flow; however, it also alters the dispersion
coefficient, and the two effects partially cancel. From Figure 2.35
it can be seen that for increasing depth, the resulting equilibrium
concentration profile decreases slowly. This result, however, is some
what misleading because of some of the other physical phenomena not
modeled. It is usually concluded for many reasons, such as low flows
near the bed if the actually vertical velocity profile is considered, the
depth of penetration of light to provide photosynthesis, the effect of
rearation over a deep water column, and so on, that deep canals are
not desirable as anoxic conditions are known to commonly occur. Thus,
this result should not be included in an analysis of proposed canal
design except to provide an estimate of the equilibrium profile, but
rather the effect of other forcing functions should be studied using
the threedimensional model.
The effects of varying Nikuradse's equivalent sand roughness, k
(Figure 2.36), and the dimensionless dispersion coefficient, K
(Figure 2.37), have similar effects, as would be expected through their
association in the longitudinal dispersion coefficient. Equilibrium
concentration profiles increase with decreasing values of these para
meters, indicating the fact that rougher canals and bends provide more
shear and turbulent eddies that aid the mixing process.
It should be noted that as the magnitudes of the resulting equi
librium profile increase, the time to achieve equilibrium also in
creases. This could have a potentially serious effect on canal networks
which are subject to high loads from time to time. From a statistical
point of view, looking at return period relationships, an efficient
canal would be one which flushes out a design percentage of the pol
lutant before the next load was expected. In an inefficient system,
a certain amount of the pollutant above the design amount considered
acceptable would be retained when the next load arrived. The concen
trations would then build up in the system until an equilibrium
condition is reached, possibly much above that considered environ
mentally acceptable by design standards.
The final test used the second model canal network (Figure 2.4)
in which a main 1,000 ft canal had a branch at varying locations along
its length. The length of the branch was 500 ft. The resulting
53
concentration profiles at low tide and high tide for location of the
branch canal at x 200, 400, 500. 600 and 800 ft are shown in
Figures 2.38 and 2.39 respectively. The results indicate that as the
branch canal is placed nearer to the deadend, the resulting concen
tration profile in the main canal is reduced as the effective excursion
of the tidal prism is increased. However, the converse is true for
the branch canal and the problem becomes a trade off dependent on
acceptable design criteria.
Jpstream
reach
In'er'or
Unctionrs
encs
branch\
0
Figure 2.1 Definition Drawing of Canal Network.
Figure 2.2 Schematic Drawing of Trapezoidal CrossSection.
1
C
IT
"2
]i
57
2
C
c)2^
0;C I 
J
z 
 I
0.5
6 12 18 6 12 18 6
'NOVEMBER 16, 1975 NOVEMBER 16. 1975
20:30 0855
Figure 2,5 Typical Atlantic Coast Tide Curve.
i]13A NI NOIiVA311 1101i
Ll_
4 0
0_ 0
LU
6 0
0
80
.6 .4 .2 0 .2 .4 .6
LONGITUDINAL VELOCITY, FPS
.8
o ,
.4 O
L
S.2
/
.6 .4 2 0 .2 .4 .6
LONGITUDINAL VELOCITy, FPS
Figure 2.7 Typical Lzgarithmic Velocity Profile.
di= d, +a Cos t
Reach i L
7/7/'.'
Figure 2.3 Schematic Drawing of Horizontal Water Surface Assumption.
Upstreomr
junction
Downstream
junction
m m+1 m +2
m Mi 1M~2
T= I
m MAImt
T=2
I i7
Figure 2.9 Schematic Representation of Effect of Numerical Dispersion.
C
U
C
C
44.
1/ C
V C
5. I
U
"5 I
4~4
~
C
/ 
x
'A
~
'' C
C
"5 C
CC
C
1/ 
U C
C
U
C
C
~l// 5C
4
I & 'C 4
C
I! CU
5,
C
C
'C
'C u5 'C
C
J
(Lndd) uoL;~;u~uoD
'C
C"
C
U
CHARACTERISTIC METHODS
VELOCITY
nx
1.MOVEABLE GRID
n
i+1 Distance, x.
VELOCITY
lnterpola:edl Values
i i i 1 Di sance ,ix
2. FIXED GRID
Figure 2.11 Schematic Representation of Convective Steps of Method
ofCharacteristics Methods.
3 8 Low Tide

High Tide
4
0 200 400 600 800 1000
Distance From Deco End (feet)
Figure 2.12 Comparison of Concentration Profiles for Fixed Grids and
MIovable Grid MetnodofCharacteristics.
66
U
0
0
co
Ll,
LL
41
Iid N ,1.04 N : '
~JY,
67
EE
no
C
0
0
a)
I
a
 a)
a
C
r~1 U U
C
U
U)
U
0
68
a
a
0
0)
a
U)
L
C)
a
0)
0
a
~0
QJ
N
a
w
a
a
(C
U)
a)
C
0)
U
________________________ 0
p ____________________________________________
J
p
p 
"' ".. : .., ,:'' .'" : .':' ,' ....... "" ': 1,. ,' . ,.,
***A..Ll^l;';";;l;J.;^^^^^^^.^
~
* L6~,s.t2LL.A.. ~.&.*. ~ P 
P
Figure 2.16 Typical Analog Output.
(Vdd) NOIiVHIN3NO0D
o Q
o0
0)
0
0)
.0
o
0
Li.
C)
z
w
0'
a)
LLi
Sz
aj
v)
cL
wa
uo
200 400 600 800
DistancR from the ocean boundary in it
1000
Figure 2.19 Low Tide Concentration Profiles for Various 4x and At 
Second Upwind Difference Method.
c(,)
/
0.5 Fm 0.0 0.5
1. Rm 
Figure 2.20 Rectangular Distribution Approximation to Actual
Distribution.
T=O m*0+
TIm
m rn+ I M+2
T 2
Figure 2.21 Scnematic Representation of Conservation of First Moment.
o
___ ___ ___ o
4 0
a
C
U
C
__ __ :1
4 C
I
A I
C
0
A I
0
LI~ ,~
LU
C
_____ 0
z =
< .~
I
0 U) C
0 
I C'\J ~
U
C
C
C,
ci
3
C
J
(D
c~J
(uidd) NOIiV~LLN~NO0
C)
C
C)
U
Distance Froml Dead end ft)
Figure 2.23 Comsarinc of Techniques' Accuracy in Modeling
Pure Corvectior.
0
0
0
CD
CC
C~j 0 o 0 N CD (
(Ludd) N I i JiN NO
0 200 400 600 800 ICC
Distance From DeadEnd (feet)
Figure 2.25 Concentration Profiles for Second Test Canal Network 
Secord Upw.ind Difference Method.