Numerical modeling of solute transport in tidal canal networks

Material Information

Numerical modeling of solute transport in tidal canal networks
Walton, Raymond, 1951-
Copyright Date:
Physical Description:
xxii, 476 leaves : ill. ; 28 cm.


Subjects / Keywords:
Canals ( jstor )
Diffusion coefficient ( jstor )
Dyes ( jstor )
Modeling ( jstor )
Saltwater ( jstor )
Simulations ( jstor )
Three dimensional modeling ( jstor )
Velocity ( jstor )
Velocity distribution ( jstor )
Wedge bodies ( jstor )
Canal ecology ( lcsh )
Canals -- Florida ( lcsh )
Civil Engineering thesis Ph. D
Dissertations, Academic -- Civil Engineering -- UF
City of Gainesville ( local )
bibliography ( marcgt )
non-fiction ( marcgt )


Thesis--University of Florida.
Bibliography: leaves 468-475.
General Note:
General Note:
Statement of Responsibility:
by Raymond Walton.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
000065654 ( alephbibnum )
04361986 ( oclc )
AAH0868 ( notis )


This item has the following downloads:

Full Text







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. Pop-Stojanovic 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 Ho-Shong 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-


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 Newcastle-Upon-Tyne, 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/OE-4, 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.




. ii

LIST OF TABLES . . . . . . vii

LIST OF FIGURES . . . . . . .viii

NOTATION. . . . . . . . xiv

ABSTRACT . . . . . . . xx


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.1 Development of One-Dimensional 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 One-Dimensional System . 18
2.1.4 Transport Equation and Dispersion Coefficient 23
2.1.5 Boundary Conditions . . . 27
2.2 Finite-Difference and Finite-Element Approaches. 29
2.3 Method-of-Characteristics 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 One-Dimensional
Mass Transport .... . . . 49

3 FIELD MEASUREMENTS. ....... . . 98

Introduction and Site Descriptions
Geometry and Tides . .
Velocity and Wind. . .
Salinity and Temperature . .
Dye Dispersion Studies . .


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.1 Transport Mechanisms. . . . . 175
5.2 Longitudinal and Lateral Diffusion Coefficients 180
5.3 Vertical Diffusion Coefficient. . . 184


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.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


9.1 Summary of One-Dimensional Modeling. . .. 301
9.2 Summary of Three-Dimensional Modeling. . .. 306
9.3 Future Research. . . ... . 310
9.4 The Numerical Model as a Design Tool . .. 312


A USER'S MANUAL . . . . . 314


C PROGRAM LISTING . . . . . 366


REFERENCES. . . . . . . 468




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 Three-Dimensional 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 Three-Dimensional Model Test Canal. .. 299

8.5 Variability Studies Using Three-Dimensional Model. .. 300


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 Cross-Section. .. .... 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--
Finite-Difference 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 Method-of-Characteristics. . .. 65

Figure Page

2.13 Low Tide Concentration Profile for Various Ax--Fixed
Grid Method-of-Characteristics. . . . 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 Cross-Sections at K and Y in 57 Acres System. . ... 109

3.6 Cross-Sections and Plan View of Loxahatchee North Canal,
June 1977 . . . 110

3.7 Electromagnetic Velocity Meter Set-Up . . 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 Set-Up. . . . ... 114

3.11 Strip Chart Recording of Wind Data--57 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 3--57 Acres System,
October, 1977 . . . . . ... 119

3.16 Salinity Profiles--Loxahatchee River System, June 1977. 120


. 121

Photograph of Fluorometer and Recorder Set-Up. .

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,

1977 .

1977 .

1977. .

1977. .

1977 .

1977. .

1977 ....

Profiles--57 Acres
. . . . 122

Profiles--57 Acres
. . . . 123

Profiles--57 Acres
. . . . 124

Profiles--57 Acres
. . . . 125

Profiles--57 Acres
. . . . 126

Profiles--57 Acres
. . . . 127

Profiles--57 Acres

4.1 Example of Typical Non-Gaussian 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 Theoretical-Wind-Induced 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 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 Three-Dimensional 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 Cross-Sections, 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 ppm--Ten Tidal Cycles. ... 280

8.9 Case 1W: Initial Concentrations, ci = 100 ppm, Back-
ground Concentration, cpy = 5 ppm--Fifty Tidal Cycles. 281

Figure Page

8.10 Case 2W: Initial Concentrations, ci 5 ppm, Back-
ground Concentration, CRW = 100 ppm--Ten Tidal Cycles. .. 282

8.11 Case 2W: Initial Concentrations, ci = 5 ppm, Back-
ground Concentration, cRW = 100 ppm--Fifty Tidal Cycles. 283

8.12 Case 3W: Lateral Inflow Distribution Along Length of
Canal--Ten Tidal Cycles. .. .. . . 284

8.13 Case 3W: Lateral Inflow Distribution Along Length of
Canal--Fifty Tidal Cycles . . . ... 285

8.14 Case 4W: Lateral Inflow Distribution Along Upper 1/2
of Canal--Ten Tidal Cycles . . . ... 286

8.15 Case 4W: Lateral Inflow Distribution Along Upper 1/2
of Canal--Fifty Tidal Cycles . . . 287

8.16 Case 5W: Lateral Inflow at Dead-End--Ten Tidal Cycles 288

8.17 Case 5W: Lateral Inflow at Dead-End--Fifty 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 Canal--Fifty Tidal Cycles 292

8.21 Case 4S: Effect of Salt Wedge With Lateral Inflow
Distribution Along Upper 1/2 of Canal--Fifty Tidal
Cycles . . . . . . 293

8.22 Case 5S: Effect of Salt Wedge With Lateral Inflow
Distribution at Dead-End--Fifty Tidal Cycles . 294


a tidal amplitude (L)

a -a4 constants of integration

A cross-sectional 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

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)

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
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)

length scale of turbulent eddy (L)

a characteristic length of the cross-section 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
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 cross-sectional mean velocity (L/T)

uD dispersion velocity in x-direction (L/T)
upD dispersion velocity in y-direction (L/T)
upD dispersion velocity in z-direction (L/T)
uF velocity of front of saltwater wedge (L/T)

u cross-sectional mean velocity at upstream section of reach (L/T)

u* bed shear velocity (L/T)

u' turbulent fluctuation from time mean velocity in x-direction

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 z-direction

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 cross-sectional area (L2

AF( ) truncation terms of numerical approximation

AH head loss (L)
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 x-direction 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 x-direction with respect to x-direction (M/LT2
IT shear stress in y-direction with respect to x-direction (M/LT2)

Txz shear stress in z-direction with respect to x-direction (M/LT2

S angle between interface and positive x-direction (degrees)


'i(Ri) function of Richardson number dimensionlesss)

h,, tidal frequency (1/T)


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


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



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


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 one-dimensional 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 three-dimensional mass-transport 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

well-documented inability of one-dimensional 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.



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.3-1.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 systems--which exist in an
intricate set of balances--we 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

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 two-part 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

three-dimensional 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,

WATER-QUALITY 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. 1-2]

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 limits--not 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 one-dimensional 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 multi-layered 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


The aim of this research is to develop a predictive three-dimensional

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


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

one-dimensional 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-


The velocity field in coastal finger canal networks is rarely one-

dimensional. In fact observations, discussed in Chapter 3, often show

multi-layered 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 three-dimensional mass-transport 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 three-dimensional mass

transport model is developed in Chapter 6 using a mass-in-cell 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 one-dimensional mass transport

model developed in Chapter 2, and a simple two-dimensional 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.



Figure 1.1 Undeveloped Canal Network in Florida.


Developed Canal Network Takeover of Natural System in


4. V


Figure 1.2 -

Figure 1.3 Example of Polluted Canal Network in Florida.

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-
F-0 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 -

21-a-i< -22


-o 5002
02 CCI C-o



0 5

** 0 2
22 :t s23
2CC .022 20











U- r

S oz

0< IE

i- A
00 :


O 0 -
<0 r
<0 C 0
P0< Or
5 0< 0
Qr~ ~ 5 0
0<4<0 0

A -
O >0

A -
O 0



O 0
0000- a

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 ?
















re f~

U- (

't 0

0- )- 0

OJC fl 00

C z
4 >

Z z


O <
4I 0

0< ^
> u
> Z I



2.1 Development of One-Dimensional Model

Before the development of a three-dimensional model was begun,

it was decided to develop a one-dimensional mass-transport 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 three-dimensional model.

In the remainder of this section, the geometry, hydrodynamics,

transport equation and dispersion coefficients, and boundary conditions

for the one-dimensional 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 cross-sections 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 dead-ends. 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 x-axis 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

x-axis. 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 x-axis 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 cross-section with

bottom width, b, and with the same inverse side slope, s, on each

side of the canal (Figure 2.2). The cross-sectional 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 dead-end 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 dead-end 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 One-Dimensional 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 cross-section is zero at any

time, and hence which cannot be modeled using a one-dimensional

numerical scheme. To compensate, many modelers try to reproduce the

three-dimensional effects by adjusting the one-dimensional 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 one-dimensional 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 Saint-Venant 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


t = time (T)

Q = discharge (L3/T)

x = longitudinal displacement (L)

q- = lateral inflow per unit length of reach (L2/T)
u = cross-sectional 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.,


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 one-dimensional 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 10-5

and 10-6.

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 one-dimensional 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)


d = mean tidal depth in reach i, (L)


2 2i (2.5)


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

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

Nu = number of upstream reaches.


d = B dd (2.10)
dt dt


A = BL
A = area of water surface in canal network upstream of

a section (L 2),


B = top width of cross-section (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 one-dimen-

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 one-dimensional canal network.

2.1.4 Transport Equation and Dispersion Coefficient

The three-dimensional mass-transport 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. 576-578; 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


x, y, z = coordinate directions (L)

u, v, w = velocities in x, y, z directions respectively

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 one-dimensional mass transport equation can then be obtained

by cross-sectionally 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)


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 one-dimensional

mass-transport 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. 28-33]. A

non-coznsZCrvativrc 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)


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)


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 one-dimensional mass-transport

Equation (2.17) or (2.18) is the concentration, c.

2.1.5 Boundary Conditions

The one-dimensional mass-transport Equation (2.17) or (2.18) is

a second-order parabolic equation, being second-order in the spatial

variable, x, and first-order 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 dead-ends of the canal system, a zero-flux 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)


c = concentration at dead-end dimensionlesss)

cI = concentration at adjacent node dimensionlesss)

implied that the gradient of the concentration profile in the segment

of the reach adjacent to the dead-end is zero, which led to erroneous

results. Higher-order schemes were discounted because it was felt that

they used information which was too far away from the dead-end to be


To overcome this problem, a simple mass balance expression was

developed based on the mean velocity between the two nodes at the

dead-end. 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 dead-end. Thus, the final form of

the dead-end boundary condition for this test model was that the

concentration at the dead-end was a linear extrapolation of the con-

centrations at the two-adjacent 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
concentration at the tidal entrance, cTE was given from the


n+l n n n n n
TE t- TE uTE CTE TE-1 TE-l (2.25)
At Ax


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)


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 two-dimensional

finite-element 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


2.2 Finite-Difference and Finite-Element 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 one-dimensional mass transport

Equation (2.17) or (2.18) is small. If this term is dropped altogether,

then the remaining first-order hyperbolic equation is the one-dimen-

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 finite-dif-

ference approximations usually lead to inaccurate results, and possibly

instabilities, depending on the scheme used. A survey of the more

common finite-difference 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 1u|l ) xA (2.28)
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)

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,


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 finite-element model developed by Leimkuh-

ler et al [1975] was also tested, and was found to suffer from the

same drawbacks as the finite-difference models. This latter method

was therefore also discounted.

2.3 Method-of-Characteristics Approach

As was seen in Section 2.2, the numerical dispersion inherent in

the more common low order finite-difference and finite-elements

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 one-dimension mass-transport equation economically and accurately,

and minimize the numerical dispersion produced.

Equation (2.27) is a first-order hyperbolic equation. A tradi-

tional way to solve this equation is the Method-of-Characteristics.

If the one-dimension mass-transport 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

second-order parabolic equation with no dominating first-order hyper-

bolic term, and can readily be approximated by a central-difference

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)

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

finite-difference approximation to F(x,t), then

F(x,t) = F(x,t) + AF(x,t) (2.36)


AF(x,t) = truncation terms of finite-difference series


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,

AF(x,t) = -At (2.37)

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)
Using Equation (2.35) again,

P2c 2 22c
u 2 (2.39)
at 'x
and substituting into Equation (2.37)

AF(x,t) = u2 (2.40)
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 Martin-Marietta 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 one-dimensional canal network. A finite-difference 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 mass-transport equation was


d) To determine the most consistent finite-differencing


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


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 finite-difference techniques, the forward, backward,

central, and first upwind differencing schemes (the latter being simply

a backward-difference scheme for the ebb tide, and a forward-difference

for the flood tide), required a large percentage of an analog

patch-panel, 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 finite-difference 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 three-dimensional 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

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 open-end of the canal
is satisfactory.

4. The zero-flux (3c/;x =0) boundary condition at the
closed-end of the canal appears to distort the results.
However, the linear boundary condition gives results
which are in very good agreement with other models

5. The mass balance boundary condition at the dead-end
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 five-segment model is somewhat coarse spatially,
but better than was initially expected. A ten or
twenty-segment model would be useful for verifying
the results of the five-segment model, but could
not be expected to provide much additional

7. Flushing time, the time for the model to reduce
concentrations at all sections to "steady-state"
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 two-dimensional model incorporating
three nodes in the vertical direction would therefore
require three analog consoles. Multiple-board
simulations are difficult to schedule and expensive to
run. Therefore, it is concluded that for canal
networks and two-dimensional models the analog/
hybrid approach is extremely limited by hardware
availability, not generally practical, and quite

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 closed-end boundary condition and
lack of time. A single canal junction in one-dimension
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

finite-difference technique used by Lee [1977] to solve the one-

dimensional mass-transport 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 finite-difference expression for the convective term (Qc)

during the flood tide, and a forward finite-difference term during the

ebb tide.

The Second Upwind Differencing method has several advantages over

other finite-differencing 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. 67-72]. No

stationary finite-difference 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 "flux-corrected 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


Using the schematic canal of Figure 2.18 as a guide, the Second

Upwind Differencing model with flux-corrected 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


3. limited antidispersion to correct for the numerical dis-

persion errors introduced in the convective step but

controlled by the flux-corrected transport criterion,

4. net transport of pollutant into the segment by natural

longitudinal dispersion.

This method is essentially second-order accurate, as it has been

corrected for the second-order 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 finite-difference and

finite-element 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 finite-difference

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


For equations such as the one-dimensional mass-transport 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

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 finite-difference 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 finite-difference 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,
cm = f c(E)de (2.45)


x (iAx (2.46)


Fm- c(t)td!/cm (2.47)
m -0.5

2 0.5 2
R2 = 12 c()( F )2 dt/cm (2.48)

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 cross-sectional

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)


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)


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
square of the second moment, R of each composite distribution is
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)


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 finite-difference 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)


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 one-dimensional mass-transport 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 finite-element

and finite-difference 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 method-of-characteristics 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-

of-characteristics, 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 Martin-Marietta facility, it

was estimated that twenty nodes could be used on one patch-panel and

that three patch-panels 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 set-up would be unrealistic. Also, once a

patch-panel 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 method-of-characteristics 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 method-of-characteristics 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 three-dimensional mass-transport model of

a coastal finger canal network.

2.8 Effect of Varying Model Parameters on One-Dimensional Mass Transport

In the previous sections, a one-dimensional mass-transport 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 three-dimensional processes under the influence

of wind, salinity gradients, tides, secondary currents, and lateral

inflows. Their features have been incorporated into a three-dimensional

model and the results generated using the numerical scheme selected in

this chapter.

However, a study of the variability of the one-dimensional model's

parameters is useful here because the effects will be similar in

a three-dimensional 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 three-dimensional

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 dead-ends of canal systems, unless

either suitable flow-through 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

dead-end (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 three-dimensional 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


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 dead-end, 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.






Figure 2.1 Definition Drawing of Canal Network.

Figure 2.2 Schematic Drawing of Trapezoidal Cross-Section.







0;C I -


z -

- I


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 NOIiVA3-11 -1101i

-4 0

0_ 0
6 0


-.6 -.4 -.2 0 .2 .4 .6


o ,

.4 O




-.6 -.4 -2 0 .2 .4 .6

Figure 2.7 Typical Lzgar-ithmic Velocity Profile.

di= d, +-a Cos t

Reach i L


Figure 2.3 Schematic Drawing of Horizontal Water Surface Assumption.



m m+1 m +2

m Mi 1M~-2

T= I

m MA-Imt


I i7

Figure 2.9 Schematic Representation of Effect of Numerical Dispersion.



1/ C

5-. I
"-5 I

/ -

'-' C
"-5 C


1/ -
U- C

~l// 5-C


I & 'C -4


'C u-5 'C
(Lndd) uoL;~;u~uoD








i+1 Distance, x.


lnterpola:edl Values

-i i i- -1 Di s-ance ,ix

Figure 2.11 Schematic Representation of Convective Steps of Method-
of-Characteristics Methods.

3 8 Low Tide

High Tide

0 200 400 600 800 1000

Distance From Deco End (feet)

Figure 2.12 Comparison of Concentration Profiles for Fixed Grids and
MIovable Grid Metnod-of-Characteristics.








Ii-d N ,1.04 N : '







- a)



r~1 U U















________________________ 0

p ____________________________________________

p ------------

"' ".. : .., ,:'' .'" : .':' ,' ....... "" ': 1,. ,' -. ,.,



* L6~,s.t2LL.A.. ~.&.*. ~ P -


Figure 2.16 Typical Analog Output.


o -Q














200 400 600 800
DistancR from the ocean boundary in it


Figure 2.19 Low Tide Concentration Profiles for Various 4x and At -
Second Upwind Difference Method.



-0.5 Fm 0.0 0.5
1--.-- Rm -

Figure 2.20 Rectangular Distribution Approximation to Actual

T=O m*0+


m rn+ I M+2-

T 2-

Figure 2.21 Scnematic Representation of Conservation of First Moment.

___ ___ ___ o
4 0
__ __ :1

4 C

LI~ ,~
_____ 0
z =

< .~
0 U) C
0 -
I C'\J ~


(uidd) NOIiV~LLN~NO0

Distance Froml Dead end ft)

Figure 2.23 Comsarinc of Techniques' Accuracy in Modeling
Pure Corvectior.





C~j 0 o 0 N CD (

(Ludd) N I i JiN-- NO

0 200 400 600 800 ICC

Distance From Dead-End (feet)

Figure 2.25 Concentration Profiles for Second Test Canal Network -
Secord Upw.ind Difference Method.