Front Cover
 Front Matter
 Title Page
 Table of Contents
 Selection of models
 Specific features of models...
 Model verification
 Determination of pollutant concentration...
 List of symbols
 A. Hydrodynamic and water quality...
 B. User's manual and complete verified...

Group Title: Technical paper - Florida Sea Grant College Program ; no. 23
Title: Modeling the Apalachicola system
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00074940/00001
 Material Information
Title: Modeling the Apalachicola system a hydrodynamic and water quality model with a hydrodynamic and water quality atlas of Apalachicola Bay
Series Title: Report Hydraulic Laboratory, Dept. of Civil Engineering, University of Florida
Physical Description: 1 v. (various pagings) : ill., maps, charts ; 29 cm.
Language: English
Creator: Conner, Carol
University of Florida -- Dept. of Civil Engineering. -- Hydraulic Laboratory
Publisher: Marine Advisory Program, Florida Cooperative Extension Service, University of Florida
Place of Publication: Gainesville
Publication Date: 1982
Subject: Estuaries -- Mathematical models -- Florida -- Apalachicola Bay   ( lcsh )
Hydrodynamics -- Mathematical models   ( lcsh )
Water quality -- Mathematical models -- Florida -- Apalachicola Bay   ( lcsh )
Estuaries -- Charts, diagrams, etc -- Florida -- Apalachicola Bay   ( lcsh )
Water quality -- Charts, diagrams, etc -- Florida -- Apalachicola Bay   ( lcsh )
Genre: government publication (state, provincial, terriorial, dependent)   ( marcgt )
non-fiction   ( marcgt )
Bibliography: Includes bibliographical references.
Statement of Responsibility: by Carol Conner ... et al..
General Note: "May 1982."
Funding: Technical paper (Florida Sea Grant College) ;
 Record Information
Bibliographic ID: UF00074940
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved, Board of Trustees of the University of Florida
Resource Identifier: aleph - 000990249
oclc - 08701826
notis - AEW7161

Table of Contents
    Front Cover
        Front Cover
    Front Matter
        Front Matter
    Title Page
        Title Page 1
        Title Page 2
    Table of Contents
        Table of Contents 1
        Table of Contents 2
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
    Selection of models
        Page 22
        Page 23
        Page 24
        Page 25
    Specific features of models used
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
    Model verification
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
    Determination of pollutant concentration and salinity...
        Page 80
        Page 81
        Page 82
        Page 83
    List of symbols
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
    A. Hydrodynamic and water quality atlas of the Apalachicola bay System, Franklin County, Florida
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
        Page 173
        Page 174
        Page 175
        Page 176
        Page 177
        Page 178
        Page 179
        Page 180
        Page 181
        Page 182
        Page 183
        Page 184
        Page 185
        Page 186
        Page 187
        Page 188
        Page 189
        Page 190
        Page 191
        Page 192
        Page 193
        Page 194
        Page 195
        Page 196
        Page 197
        Page 198
        Page 199
        Page 200
        Page 201
        Page 202
        Page 203
        Page 204
        Page 205
        Page 206
        Page 207
        Page 208
        Page 209
        Page 210
        Page 211
        Page 212
        Page 213
        Page 214
        Page 215
        Page 216
        Page 217
        Page 218
        Page 219
        Page 220
        Page 221
        Page 222
        Page 223
        Page 224
        Page 225
        Page 226
        Page 227
        Page 228
        Page 229
        Page 230
        Page 231
        Page 232
        Page 233
        Page 234
        Page 235
        Page 236
        Page 237
        Page 238
        Page 239
        Page 240
        Page 241
        Page 242
        Page 243
        Page 244
        Page 245
        Page 246
        Page 247
        Page 248
        Page 249
        Page 250
        Page 251
        Page 252
        Page 253
        Page 254
        Page 255
        Page 256
        Page 257
        Page 258
        Page 259
        Page 260
        Page 261
        Page 262
        Page 263
        Page 264
        Page 265
        Page 266
        Page 267
        Page 268
        Page 269
        Page 270
        Page 271
        Page 272
        Page 273
        Page 274
        Page 275
        Page 276
        Page 277
        Page 278
        Page 279
        Page 280
        Page 281
        Page 282
        Page 283
        Page 284
        Page 285
        Page 286
        Page 287
        Page 288
        Page 289
        Page 290
        Page 291
        Page 292
        Page 293
        Page 294
        Page 295
        Page 296
        Page 297
        Page 298
        Page 299
        Page 300
        Page 301
        Page 302
        Page 303
        Page 304
        Page 305
        Page 306
        Page 307
        Page 308
        Page 309
        Page 310
    B. User's manual and complete verified model on computer tape
        Page 311
        Page 312
        Page 313
        Page 314
        Page 315
        Page 316
        Page 317
        Page 318
        Page 319
        Page 320
        Page 321
        Page 322
        Page 323
        Page 324
        Page 325
        Page 326
        Page 327
        Page 328
        Page 329
        Page 330
        Page 331
        Page 332
        Page 333
        Page 334
        Page 335
        Page 336
        Page 337
        Page 338
        Page 339
        Page 340
        Page 341
        Page 342
Full Text



.. S 0 SOUTH

S\ ^Atlanta


S. *'. Savannah
Chattahoochee Flint Rive .
.River :::::. ATLANTIC
I- *.:-- OCEAN
-.. :. -. ._. . '.
Mobile :. .. Jacksonville
Mobile Ny^^
: *.*.: .
Apalachicola Bay Apalachicola River -.' '-.*
0 F
MEXICO :. .. .
Tampa L*.* ..
St. Pete. -

S":.*:. ;:: Miami
'.. :.. ...:
By Carol Conner
Ann Conway
Barry A. Benedict V*
B. A. Christensen




A Hydrodynamic and Water Quality Model
With a Hydrodynamic and Water
Quality Atlas of Apalachicola Bay
Project No. R/EM-13


Carol Conner, Ann Conway
Barry A. Benedict and B. A. Christensen

Technical Paper No. 23
May 1982

Hydraulic Laboratory
Department of Civil Engineering
University of Florida

Technical Papers are duplicated in limited quantities for specialized
audiences requiring rapid access to information and may recieve only limited
editing. This paper was compiled by the Florida Sea Grant College with support
from NOAA Office of Sea Grant, U.S. Department of Commerce, grant number
NA80AA-D-00038. It was published by the Marine Advisory Program which functions
as a component of the Florida Cooperative Extension Service, John T. Woeste,
Dean, in conducting Cooperative Extension work in Agriculture Home Economics,
and Marine Sciences, State of Florida, U.S. Department of Agriculture, U.S.
Department of Commerce, and Boards of County Commissioners, cooperating.
Printed and distributed in furtherance of the Acts of Congress of May d and
June 14, 1914. The Florida Sea Grant Collge is an Equal Employment Opportunity-
Affirmative Action employer authorized to provide research, educational
information and other services only to individuals and institutions that
function without regard to race, color, sex, or national origin.

*Also published as Report No. 8107, Hydraulic Laboratory, Department of Civil
Engineering, University of Florida.

1. INTRODUCTION . . . ... . . .. 1

1.1 Logistics of Model Approach . . . . 1
1.2 The Apalachicola Basin/Estuary System . . . 7
1.3 Classification of Apalachicola Bay . .. .18
1.4 The Apalachicola Bay Numerical Model and Atlas of
Hydrodynamics and Water Quality . . .... 19


2.1 Criteria for Selection . .... .. . . 22
2.2 Models Considered. . . .... . 23
2.3 Models Selected . . .... . . . 24


3.1 The CAFE-1 Model . . . . . 26
3.1.1 General Description . . . . 26
3.1.2 Equations for CAFE-1 . . . .. 31
3.1.3 Turbulent Eddy Viscosity Coefficient Eij . 38
3.1.4 Application and Debugging of CAFE-1 Model for
Apalachicola Bay . . . . ... 42
3.2 The DISPER-1 Model . . .. ..... 44
3.2.1 Equations and Boundary Conditions for DISPER-1 .. 45
3.2.2 Stability and Convergence of DISPER-1 . .. 46

4. MODEL VERIFICATION . . .... . . 50

4.1 Verification of CAFE-1 . . .... . 50
4.1.1 Field Data . . .... . . 50 Tide Recordings . . . 50 Velocity Measurements . . .... 54
4.1.2 Model Input for Verification Runs . . .. 54 Grid Configuration . . .. 54 Turbulent Eddy Viscosity Coefficients .58 Manning's n . . . . 59 River Flow . . . . 59 Wind Velocity . . . . 61 Tide Information at Boundaries ...... .61
4.1.3 Comparison of Model Results with Field Data .. 62
4.1.4 Comparison with Wang's Verification . .. 66
4.1.5 Further Calibration . . . . 66
4.2 Verification of DISPER-1 . . . . .. 68
4.2.1 Field Data for Quality Verification . .. 68
4.2.2 Model Input for Verification Runs . . .. 71
4.2.3 Results of DISPER-1 Verification Runs . .. 72
4.3 Satellite Verification of DISPER-1 . .... . 74


5.1 River Flows . . . ... .... .. 76
5.2 Tides . . . . .. .. . 76
5.3 Winds . . . . ....... 78
5.4 Generated Results . . . . ... 79

CONCENTRATIONS . . . .. . . 80

7. ACKNOWLEDGEMENTS . . . . ... ...... 83

8. LIST OF SYMBOLS . . . . ... ...... 84

9. REFERENCES . . . . ... .. .. .. 87


TAPE. (1 reel)


The present report discusses the establishment of a numerical model

of the Apalachicola system consisting of the River, its tributaries and

the biologically highly productive Apalachicola Bay. Much of the material

has been published elsewhere in the form of reports (e.g., Graham, DeCosta,

and Christensen, 1978), scientific papers (e.g., Christensen, 1979, Graham

and Christensen, 1978) and bulletins ofa more general nature (e.g., Hill

and Graham, 1980).

The present report should therefore be considered as a review of the

progress made with special emphasis on the final results that are presented

in the form of an atlas depicting the hydrodynamics and pollutant migration

in the Apalachicola Bay during an "average" year taking tidal flushing as

well as wind and river flow into consideration. This atlas will allow the

reader to determine water velocity, orientation of water velocity, water

quality corresponding to any known water quality ih the river, and salinity

of the bay waters. The atlas is attached to the present report.as Appendix

A. The complete verified numerical model is attached as Appendix B on

computer tape.

1.1 Logistics of Model Approach

Man's development of drainage basins draining to estuarine waters

will always have an impact on the quality of these waters and therefore

also on the biological system of the estuary in question. In the case of

Apalachicola Bay this has clearly been demonstrated by Livingston

(Livingston etal., 1981). A further impact on the economy depending on

the estuarine biosystem and the ocean biosystemits supports must therefore

be expected and should be considered by the responsible developer.

Such a consideration is certainly justified when it is recalled that

our estuarine systems are among the most productive areas as far as biomass

is concerned. Figure 1.1 shows a comparison of biomass production in ton

per acre per year for typical areas. It should be noted that the estuarine

production is two to four times that of our best agricultural land areas



Figure 1.1.

Comparison of Biomass Production Rates. (Source: Teal and
Teal, 1969).

that,by the way,only can sustain their high production rates by extensive

use of fossil fuel based fertilizers.

It is therefore recommended that the impact of development is evalu-
ated and that the results of the evaluation is fed back to the developers
or developing agencies to insure optimum utilization of our coastal
Figure 1.2 is a flow chart showing the logistics of optimum land
utilization by use of computer models. It is noted how the development's
impact on the riverine system is evaluated by a predictive model and the
output from this model serves as input to the estuarine model. The output
from the estuarine model will again provide information enough in the form
of water velocities, pollutant concentrations and salinities to enable the
marine biologist to evaluate the impact on the estuary's biosystem. From
these the economists may draw their conclusions as to the impact on the
region's economical system which will serve as a foundation for further
action in the political system providing the feedback to the planners and
developers. A slightly more detailed flow chart is given in Figure 1.3.






Figure 1.2. Logistics of Model Approach.

(calibrated by existing data)




Detailed Flow Chart of the Model Approach to Rational Utilization of Coastal Drainage Basins.



Figure 1.3.

The utilization of computer models in this scheme is of course quite

obvious since physical models or full scale experimentation would be

economically unrealistic if not completely impossible.
Since this work deals with large areas of land and sea,it is logical to

utilize satellite imagery as input to the models and for verification of

the results obtained by the models. This is indicated in Figure .2. Input

from the LANDSAT system is suggested in the river model. For instance the
LANDSAT images will give information about the present utilization of the
land by colors as shown in the example of Figure 1.4 indicating how the

drainage basin is used for silviculture in the neighborhood of Apalachicola


10 5 0 10 km
Figure 1.4. Computer Enhanced LANDSAT Image of Area North of Apalachicola
Bay Showing Silvicultural Activities. August 19, 1976.

In Figure 1.4 white indicates sand or urban development while the

color green shows marsh, swamp or natural forest. Purple areas are clearcut

forest that has been harvested within one year of the date of the image.

Orange is showing revegetated forest areas.

LANDSAT images enhanced for color (quality) of surface waters are

used for verification by pattern recognition of water quality as discussed

later in this report and by Graham, Hill and Christensen (1978).

Figure 1.5. LANDSAT Image of Apalachicola Bay Computer Enhanced for Water
Quality. Outgoing Tide.

Verification data are of course also obtained by direct observation

of velocities and salinities in the bay.

1.2 The Apalachicola Basin/Estuary System

The Apalachicola Basin and Estuary system was selected for the

modeling effort described in this report because of its relative simplicity

as far as pollution and pollution sources go. Pollution of the river and

bay waters is at the present time generally limited to a lowering of the pH

probably caused by clearcutting of areas around East Bay. This phenomenon

has been extensively studied and recorded by Livingston (Livingston, 1981)

during the last decade. Information relating the increased acidity of the

waters to the silvicultural activities is scarce and has not been available

to the writers.

Although the model developed and verified in this report is specifi-

cally dealing with the Apalachicola system it should be emphasized that

the general methodology applied may be used in other river-bay systems.

As shown in Figure 1.6 the water entering the Apalachicola Bay may

originate as far away as from the Atlanta, Georgia, region from where it

enters the Apalachicola River primarily through the Chattahoochee and

Flint rivers. Briefly, the basin is about 50,500 km2 in area and the

Apalachicola River properis considered to begin at the Florida line and

flows about 170 km to the Gulf of Mexico. The river has been improved by

the U.S. Army Corps of Engineers to Bainbridge, Georgia, and a large

reservoir has been established at the Florida line.

The Apalachicola River has very good hydrologic data. The mean

yearly hydrograph for the decade 1967 through 1976 is shown in Figure 1.7.

The yearly average discharge is about 700 m3s This is a substantial




S". S ":'. Savannah
Flint Rive .

.: Jacksonville

la River :.. .':


St. Pete


Figure 1.6. The Apalachicola River/Bay System with Tributaries.

discharge. However, compared to the rate of tidal exchange of waters in the

bay it is very small. The river is therefore having only a limited

influence on the hydrodynamics of the Apalachicola Bay.

Apalachicola Bay, shown in Figure 1.8, is a barrier island-contained

estuary on the Florida Panhandle. It is geomorphically typical of such

systems, but is significant because of the importance of its waters to

marine life (Livingston, 1981). At present, it suffers relatively low

levels of point-source pollution.

The bay is about 550 Km2 in area, depending upon where the boundaries

are drawn. Mean depth is about 2 meters and is everywhere less than 3

meters except in West Pass. As expected the system is dominated by the

I ,\ -


T 77/ ^ ------- --

"E /- QAVE = 705 m's '
S700 ------

U 600 //

X 500

400 -

Q: !

100 -O -- -- -- -



Figure 1.7. Mean Yearly Hydrograph for Apalachicola River for
1961 1976.

Apalachicola River whose mean annual discharge is about 700 m3/s. The

mean tidal range at Apalachicola is 0.40 m. The tidal prism is thus

220 million m3 on the average and the mean residence time per river water

about 17 days.

In general the bay is shallow, well-mixed and prone to being wind-

driven. Little else is known about its hydrography. It is quite remarkable

that the entire literature on the hydrography of one of the largest

estuaries on the Gulf Coast can be easily read in about fifteen minutes.

Most of the hydrographic data has been collected by biologists and is not

readily useful for engineering or oceanographic purposes, nor has it been

interpreted in a physically rigorous manner. Figure 1.8 shows the major

points of hydraulic interest in the modeling effort. These are,beginning

at the east side and proceeding in a clockwise manner, East Pass between

St. George Island and Dog Island, Sikes Cut penetrating St. George Island,

West Pass between Sand Island and St. Vincent Island, Indian Pass between

St. Vincent Island and the Florida Panhandle, the Apalachicola River

entering the bay north of the city of Apalachicola, and East Bay receiving

much of the runoff from the land areas used for silviculture in Franklin


The following Figures 1.9 through 1.21 are aerial photos of these

points taken from the flight path in Figure 1.8 at the locations marked

A through M.





Figure 1.8.

The Apalachicola Bay. Flight Path Used When Figures 1.9 Through 1.21 Were Taken

Figure 1.'9.

Apalachicola Bay.
Left: St. George

VIEW A: East Pass Seen from South.
Island. Right: Dog Island.

Figure 1.10.

Apalachicola Bay. VIEW B: St. George Island at Rattlesnake
Cove and Goose Island Seen from South. -

i I----L-I-
;~T~-~-TI~C ~-


~: i:

Figure 1.11.

Figure 1.12.

Apalachicola Bay.
VIEW C. St. George
Island at Causeway
Seen from South.

Apalachicola Bay.
VIEW D. St. George
Island between
Causeway and Sikes
Cut Seen from South.

Figure 1.13. Apalachicola Bay.

VIEW E. Sikes Cut Seen from the Mexican

3r. gIPtc'w4-Ci~ *c 53\.~j-
--- c~~r~--- .- -

h ~~e'4~z. -

Figure 1.14. Apalachicola Bay. VIEW F.
Left: St. Vincent Island.
(Sand Island).

West Pass Seen from Southwest.
Right: St. George Island

~r c~-~rar` -'-r n`r-p

rj -;

Figure 1.15. Apalachicola Bay. VIEW G. Indian Pass Seen from Southwest.
Right: St. Vincent Island.

Figure 1.16. Apalachicola Bay. VIEW H.
Left: St. Vincent Island.

Indian Pass Seen from Southwest.

Figure 1.17. Apalachicola Bay. V
Seen from Northwest.

IEW I. North Shore of Bay at Green Point

Figure 1.18.

Apalachicola Bay. VIEW J.
Bay at City of Apalachicola
John Gorrie Memorial Bridge.

The Apalachicola River Entering
(Right) Seen from Northwest.


li- ~ L -:'--.
r~rs~***~~r: ';l'-==~,iriiL
-~L ~-iT~



Figure 1.19. Apalachicola Bay.

VIEW K. East Bay Seen from Southeast.


Figure 1.20. Apalachicola Bay.
from North.

VIEW L. St. George Island Causeway Seen



Figure 1.21. Apalachicola Bay. VIEW M. East Pass Seen from West. Left:
Dog Island. Right: St. George Island.

1.3 Classification of Apalachicola Bay
Estuaries are often classified according to the degree and manner in
which they are stratified. Stratification is a function of river flow,
tidal prism, volume, depth-wide ratio, density difference, hydraulic rough-
ness and so forth. Classification schemes have been proposed (Pritchard,
1970, 1973) which are based on several dimensionless parameters:

P P2 1 Pl D R (

T = tidal period (12.42 hr)
P = tidal prism
Q- river discharge

P1,P2 = characteristic extreme densities
L = estuary length
D = average estuary depth
R = tidal range

If p, corresponds to s (salinity) = 0 and p2 to s = 32% (ocean value)

and R is (almost) constant, then any classification scheme must be a function

of at least 2 parameters one expressing inertial mixing ability and the

other reflecting geometry. For a given estuary, geometry and hydraulic

characteristics are fixed so a classification scheme is primarily based on

P (1.2)

and, since T is also fixed, it can be seen to be a simple function of

retention time providing the estuary depth is not greatly dependent upon

river discharge.

The foregoing analysis has not considered wind shear, which is a very

important mixing mechanism for small D/L.

Considering the above and Pritchard's estuarine classification scheme

shown in Figure 1.22 the Apalachicola Bay may be characterized as a width

dominant estuary controlled by tidal currents originating from astronomical

as well as wind tides. It is a type D estuary.

1.4 The Apalachicola Bay Numerical Model and Atlas of Hydrodynamics
and Water Quality

Since the Apalachicola Bay is only slightly influenced by the river

discharge as far as the bay's hydrodynamics goes it was decided to use the

mean yearly hydrograph for Apalachicola River shown in Figure 1.7 as river

model input to the hydrodynamic estuarine model.

Water velocities, their orientation and water surface'elevations are

modeled by use of the Wang-Connor CAFE-1 model (Wang and Connor, 1975).

This is a vertically integrated finite element model.



Figure 1.22.

Pritchard's Estuarine
(1970, 1973).

Classification Scheme.



Water quality is modeled by use of the Wang-Connor DISPER-1 model

(Wang and Connor, 1975) and the output from the CAFE-1 model. DISPER-1

is also a finite element model.

Since very little is known about the present and future water quality

of the river an arbitrary pollution concentration of 100 units have been

used in the quality model. This, together with a fairly simple formula for

conversion of the computer printout will allow the user of the reported re-

sults to compute the water quality at any location and time in the bay

corresponding to any pollutant concentration in the river. At the present

time the quality model is limited to conservative pollutants. However,

later introduction of decay and generation terms in the governing equation

will allow extention to nonconservative substances.

The final results are given for the 12 months of the typical year

during which average tides, winds and river flow are presumed to prevail.

The conditions are given for the typical tidalVcycTe for each month.

Velocity, net velocity as well as water quality corresponding to river

concentration 100 are given. Conditions are given for time increments

equal to one-eighth of a tidal cycle beginning at low tide. Thus each month

will be represented by 8 velocity maps, 1 net velocity map, and 8 water

quality maps. The total of these 204 maps makes up the atlas given at the

end of this report.


A short discussion of the selection of the model is given in this chapter.

For a more complete discussion the reader is referred to an earlier report

to Sea Grant (Graham, DeCosta and Christensen, 1978).

2.1 Criteria for Selection

The primary criterion for selection of any model must of course be

their ability to accurately reflect existing and projected behavior in the

water body of interest. The basic characteristics of behavior of

Apalachicola Bay are described briefly in Section 1 of this report and in

more detail by Graham, DeCosta and Christensen (1978). It was concluded

that a two-dimensional vertically, averaged numerical model was appropriate

for Apalachicola Bay due to its generally well-mixed body of water. It is

likely that Apalachicola Bay would represent the case in Florida where

this assumption-is most marginal, so that application to most other Florida

estuaries would be justified.

It was felt that a real-time hydrodynamic capability was required to

properly simulate transient velocities and quality in the bay, since tidal

flow dominates and stormwater inflow is transient by definition. Since

many of the concentration terms (of the form u c) in the equations are

nonlinear, it is not justified to use tidally-average net-flow models.

A primary goal was to reduce the sub-grid scale eddy diffusivity so

that the quality model will be as predictive as possible.

While the real-time 2-D hydrodynamic and dispersion models are

definitely state-of-the-art, 3-D and multi-layer models are not (in the

author's opinion). While promising results are being made with 2-layer

'and multi-layer models, they must still be considered to be in a development

phase at least in terms of the capabilities of most potential users in

Florida. Data requirements for loading a multilayer model are formidable,

so the likelihood of reliable applications is even further removed.

To summarize, the selection criteria used were:

1) model must be a least two-dimensional

2) model must simulate transient conditions

3) model must be sophisticated, yet readily available,

proven, and tested

4) model must be nonproprietary, or documented to the

extent that it is effectively nonproprietary

5) model must be "state-of-the-art", but past the

research phase

6) model must be able to be run by users with reasonable


7) model must be able to be run at a reasonable cost

8) model must be flexible, not site-specific, in its


9) model must be compatible with a water quality package

10) model must be able to simulate conditions typical of

the Apalachicola, and most other Florida estuaries, viz.:

a) shallow

b) influenced by freshwater inflow, tides and wind.

2.2 Models Considered

Earlier reports and papers describe the several models considered

for application to Apalachicola Bay according to the criteria in Section

2.1. The strategy was to select a model for the hydrodynamics first.

After this model was chosen, a compatible water quality model could be


Among models considered were the following, with more detail given

by Graham (1977), Graham, DeCosta and Christensen (1978), and Graham,

Daniels and Christensen (1979).

1. Hydroscience Estuary Model (Hydroscience, 1977).

2. Dynamic Estuary Model (also forms the receiving water model, RECEIV,

for the EPA Storm Water Management Model (SWMM) (Feigner and

Harris, 1970)). See also Huber et al. (1975).

3. Leendertse Finite-Difference Model (Leendertse, 1967).

4. Masch Gulf Coast Model (Masch et al., 1969).

5. M.I.T. Finite Element Models

CAFE-1 for hydrodynamics (Wang and Connor (1975), Pagenkopf
et al. (1976))

DISPER-1 for water quality (Leimkuhler, et al.(1975), and
Pagenkopf et al. (1976))

A number of other models were considered and discarded at an early

stage, especially those with no real time tidal characteristics (which is

also true of model 1 above). Extensive reviews are provided elsewhere and

will not be repeated here.

2.3 Models Selected

Based on a review of available models, coupled with the criteria

presented in Section 2.1, it was decided that one package, the CAFE-1,

DISPER-1 system developed at M.I.T. (also for Sea Grant) seemed overwhelmingly

superior. In addition to technical and verification advantages, it was

deemed important to select this finite element model for several reasons:

1) it provides more flexibility in discretizing the spatial


2) nodes can be placed at known measurement stations,

3) finite element models were recommended by Dr. Frank Masch,

formerly of WRE, and Dr. John Wang of the University of

Miami, both recognized experts in the field,

4) Wang and Connor's (1975) finite element model has

recently been applied by Swakon and Wang (1977) to

Biscayne Bay with good results,

and 5) tapes and advice are available from Dr. Wang at the

University of Miami.


In the following, only sufficientmodel detail will be presented to enable

understanding model features and limitations. More detail can be found in

the basic references and in earlier Sea Grant reports on this work. (See,

for example, Graham, Daniels and Christensen, 1979).

3.1 The CAFE-1 Model

3.1.1 General Description

The CAFE-1 (standing for 1-layer Circulation Analysis by Finite

Elements) model was developed and tested by M.I.T. (Wang and Conner, 1975).

It has the following general properties:

1) Real time, i.e., describing velocities and depths through

the tidal cycle.

2) Finite-element formulation.

3) Implicit time stepping.

4) One-layer, i.e., it is vertically integrated.

The model differs from other popular ones based on the Leendertse

finite-difference format in that a finite element computation scheme is

used. This computational scheme is more complicated, but has two very

significant advantages:

1) The elements can be of different size.

2) The elements can assume almost any configuration, providing

a triangular network is used.

Thus, the finite element grid can very closely approximate the geo-

metrical boundaries of the body of water being modeled, while a finite

difference scheme must approximate the boundary as the edges of squares.


Since the finite element grid sizes can be of arbitrary size (usually

finite difference meshes are of one size), the grid can be made finer in

narrow enclosed areas as well as in areas of interest. Conversely, the

grid size can be made more coarse in large open areas, and in areas outside

of the region of interest. The advantage of having finer resolution only

in the region of interest is that computing costs, which are proportional

to number of elements, are reduced. The advantage of closely approximating

the boundary geometry in a hydrodynamic model is also quite important for

obvious reasons. Since pressure terms dominate, and these are easy to

measure and/or obtain from published data, for instance from tide tables,

and since the results are not greatly sensitive to reasonable value

choices of bottom roughness (Manning's n) and internal stress coefficients

it follows that a very good approximation to the circulation can be com-

puted with very little (and inexpensively acquired) input data. Since

verification and calibration of numerical models can be quite expensive,

one which is intrinsically quite accurate can reduce total modeling costs

substantially by minimizing the empirical modification requirements.

During the study period, several grids were utilized to model

Apalachicola Bay. Each was generally characterized by finer grid (higher

resolution) in the East and West Bayou areas and in the East Bay region,

with moderate refinement in the Sikes Cut area and coarser grid elsewhere.

Other specialized grids could easily be developed to study other features

in the bay, although experience has shown that each new grid has special

problems in model stability and convergence to be overcome. These problems

will be discussed in more detail later.

Inputs to CAFE-1 include

1) a grid, which requires nodes and elements to be defined

2) depths at each node

3) boundary orientations

4) wave amplitude, phase and/or flux at boundaries

5) latitude

6) wind speed and direction

7) Manning roughness]
operator controlled
8) internal stress

These are relatively easy to come by, at least for the United States.

Outputs include

1) unit discharge vectors and/or velocities

2) heights above mean low water

for every node point for every timestep. In this case, the unit discharge

values (qx and qy) and a height were output for each node for a tidal

cycle of 744 minutes.

For preliminary purposes the water level was initially set to zero

and the model run until a periodic steady state was achieved (by 3 tidal

cycles). The output data for a steady tidal cycle were then written on a

disk. The data on the disk could then be manipulated and/or called for

graphical display, or used as input to the dispersion model DISPER-1. For

sinusoidal tides, repeating over each cycle, and use of a single cycle called

from memory to drive DISPER-1 result in considerable savings. The CAFE-1

program results presented later for the gridwith 281 nodes and 439 elements,

shown in Figure 3.1, typically took about 450 seconds CPU time at a cost of

about $25.00 on the University's IBM 3033 System for simulation of four tidal


The CAFE-1 model, as mentioned, was originally developed at the R.M.

Parsons Laboratory at M.I.T. It has been applied to







116.00 174.00

232.00 290.00 3118.00

Figure 3.1. Finite Element Grid System of Apalachicola Bay Used in this Study.



W. 00o


10b6. 00

64. O00


.I. I

S5 2.00 CI

1) Massachusetts Bay (Briggs and Madsen, 1973); (Christodoulou

etal.,1974); (Connor etal., 1973); and (Parker and Pearce,


2) Narragansett Bay (Connor and Wang, 1974); (Swanson and

Spaulding, undated).

3) Great Bay, N.H. (Cellikol and Reichard, 1976); (Swanson

et al., 1977).

4) Mareton Bay, Australia (Steele etal., 1977).

5) Biscayne Bay, Florida -(Sengupta et al., 1978); (Swakon

and Wang, 1977).

6) Lake Pontchartrain done at Louisiana State


All of the known successful applications appear to have been made by M.I.T.

trained personnel,with the exceptions of the Lake Pontchartrain study

and this University of Florida study.

These latter two studies,therefore,reflect the degree to which the

models are available to the public. Both LSU and UF personnel have been

able to run the programs, but it was concurred in private communications

that the package is, as yet, somewhat underdocumented for most users.

In some instances, the lack of documentation merely forces the user to

better understand the model before using it. In some cases, however,

required units or other items are unclearly specified and can cause consi-

derable confusion. In any event, current model use requires personnel

knowledge in both hydrodynamics and computer programming. Based, again,

on a comparison of respective progress at LSU and UF, it appears that a

training period of 6 to 8 months full time is necessary before satisfactory

results and progress can be obtained with this model package. This is not

inherently unreasonable, but it should underline the necessity of establishing

continuity in terms of personnel capability.

An excellent outline of model properties and capabilities at a lay

level, including potential user applications, is given in an M.I.T. Marine

Industry Advisory Services Opportunity Brief.

The models are currently available from the Parsons Laboratory,

Department of Civil Engineering at M.I.T., and from Dr. John Wang, Department

of Ocean Engineering at the University of. Miami. Users manuals are also

available from M.I.T. Acquisition costs for both models (CAFE and DISPER),

(not debugged) were about $70.00. This reasonably covers cost in 1978


The tape containing the CAFE-1 and DISPER-1 models modified for the

Apalachicola Bay is available from the Hydraulic Laboratory, University of

Florida. A copy of this tape is attached as Appendix B to this report in

its original copy submitted to the Florida Sea Grant.

3.1.2 Equations for CAFE-1

Wang and Connor (1975) present the most complete derivation of the

equations. In the process of equation development, averaging over time

(to remove turbulence terms) and space (to reduce the equations to two

dimensions) occurs. Such averaging always introduces additional coefficients

into the equations and changes the meaning of many terms. For example,

dispersion coefficients will include effects of depth averaging as well as

turbulent diffusion. Graham, Daniels, and Christensen (1979) follow the

outline of Swakon and Wang (1977) to present the basic equations. Several

variables must be introduced, including vertically integrated discharges

per unit width (qx and qy) in the two coordinate directions (x and y) and

surface displacement above mean low water, n.

The total depth, H, is then given by

H = h dz = h + n (3.1)

in which
z = vertical coordinate

h = depth at mean low water (MLW)

The basic equations can be written as follows:
Conservation of mass:

H+ + q y ql (3.2)
at x ay =I

Conservation of momentum in x-direction:

aqx 8(uqx) 8(uq )
+ + fq

1i n s aH b ah
+ [x pdz p p
9x -h x ax
s b
(F ) (Fx) + =0 (3.3)
ax xx 9y xy p x

Conservation of momentum in y-direction:

Sqv (vqx,) (vq )
y+ ( ) + +-fq
St Sx 3y x

1 f s 3H b ah
pdz p H pb h
p y -y Wy
s b
-(F ) F ) + -- -= 0 (3.4)
"y xy) -y yy y

in which
f = Coriolis parameter = 2 Wearth sin 290
corresponding to the latitude of Apalachicola Bay

p = pressure

pS = atmospheric pressure

p = bottom pressure
T = surface stress

T = bottom stress

p = local density of water in the bay

ql = volume addition rate

t = time

u = qx/H = local mean velocity in x-direction (3.5)

v = qy/H = local mean velocity in y-direction (3.6)

x,y = Cartesian coordinates (horizontal)

F ,F = internal specific stress terms (stress/density),
xy yy -sometimes termed "turbulent eddy viscosity

Mx,My = momentum addition per unit horizontal area

A Bouissinesq approximation is used in the pressure terms, viz:

n S An b 9h an H pS
pd p 7p --x" gH --- x

1 gH2 a (3.7)

in the x-direction. A similar equation is derived for the y-direction.

p(x,y,t) = po + Ap(x,y,t) (3.8)

with the usual definitions. For practical purposes the effect of

the atmospheric pressure usually induces a 1 cm change in

sea level;this may not be readily justified. The reason for doing so is

simply that these data are not available. Since most estuaries are small

in scale compared to weather systems, only a constant error results. If

the MLW is assumed to be a flat arbitrary datum (which, again, for

practical purposes is a necessary assumption in most cases), then no

detectable error results.

The term 1g H2 P is also often ignored. Based on the premise

that most passive pollutants do not change the density field. An exception,

of course is salinity.

The internal stress coefficients, can be written as

F = qi qj (3.9)
ii 26 ax ax
1.) .J 1

where Eij are the turbulent eddy viscosity coefficients, which can be

manipulated. In practice, the terms Eij serve several functions, including:

1) Expressing "true" turbulent eddy viscosity (assuming the

validity of a mixing length analogy).

2) "Internalize" Reynolds stress terms lost in averaging over

the depth.

3) "Internalize" horizontal Reynolds stress terms of such scale

that the grid cannot resolve them.

4) Expressing true molecular viscosity. This component is of

course of negligible importance at these scales.

5) An adjustment coefficient to calibrate the model.

6) Help stabilize the model.

Addition of the F.. terms essentially adds some .elegance to the model.

Leendertse's (1967) 2-D finite-difference model neglected these terms.

They were included by Wang and Connor (1975) because

"We feel that the inclusion of F.. has several attractive
properties. It allows for internal friction and thereby
energy dissipation, provided Eij is positive; it does
represent actual physical processes (although not accurately)
and it is particularly suitable for damping short wave noise
generated by numerical methods."

Note that no explicit stability criteria have yet been developed for

CAFE-1, at least as reported by Wang and Connor (1975), so inclusion of a

damping term for high-frequency noise is very useful. In comparison,

explicit stability criteria are known for a finite-difference grid of

constant element size.

The model does not seem particularly sensitive to values of Fij as

noted by Connor and Wang (1975). A comparison of Eij, values will be pre-

sented shortly. A significant test to determine whether a numerical model

application is credible and predictive lies in the values given the

"adjustment coefficients", such as Eij and n (Manning roughness coefficient)

in CAFE-1, and the dispersion coefficient Dij in DISPER-1. If these

coefficients have values which are reasonable in relation to values used

in analytic hydrodynamic approximations (taking the grid size into consi-

deration), then some confidence can be placed in the predictability of the

model. If extraordinary values of "adjustment coefficients" have to be

used to match the model output to measured data, then the application is

likely unique (or incorrect). This comparison may be a more valid "bottom

line" test than a comparison of model output to measured data.

Note that the gradient terms (aqi/xj, ... ) in Equation (3.9) are

approximated by finite values on the grid, and error results if the grid is

non-zero in size. The values of Fij and Eij are therefore functions of the

local grid scale as well as hydrodynamics.

The bottom stress terms in Equations (3.3) and (3.4) use the generally

accepted quadratic approximation:

b = Cf (q 2 q2 1/2 i- (3.10)
x =f o "x. i x H 2

in which 2
C = 1/ (3.11)
f H1/3

n = Manning's "n"

and g = acceleration due to gravity.

b 1 (3.12)
x H7/3

For shallow estuaries then, significant improvement in the output

quality is achieved by computing the velocities using instantaneous values

of total depth, since the tidal range is a significant proportion of total

depth. Note that Manning's "n" is another "adjustment coefficient".

Reasonable values are well known however, and lie in the range 0.020 to

0.040. [Most values of n have been computed for rivers. However, some study
is required to find appropriate values for flow over oyster bars and shell-

littered bottoms.]

The wind stress term T /p warrants qualitative discussion. Recent

oceanographic studies have indicated that wind is a far more important

energy input source to estuaries then had heretofore been surmised (Weisberg,

1976). Most Florida estuaries can be considered to be wind dominated. Work

at the University of Florida in coastal canals indicates wind is much more

influential in flushing than tidal action (see Morris, Walton, and

Christensen, 1978). Tidal measurements by Hydraulic Lab personnel show

that wind-set up can be up to about 3 times the tidal range.


The wind stress term in CAFE-1 is of the accepted quadratic form:

s C U2 (3.13)
Pair D U10

pair = air density
CD = a dimensionless drag coefficient

U10 = air velocity at 10 meters above water surface,
in m s-

The form given for CD is

CD = (1.1 + 0.0536 U10) x 10-3 (3.14)

This is based essentially on empirical data. Discussion of development

and validity of Equation (3.14) and (3.15) is given by Briggs and Madsen (1973).

Properly setting the boundaries may be a problem when there is a significant

wind setup. This usually requires an independent study of setup properties.

One approach is to run the model with no tide and adjust the boundaries so

a smooth setup is established and then superpose tides on the MLW level.

Alternatively, real-time wind and tidal data can be input.

Final forms of Equations (3.2) to (3.4) are developed by inserting

these approximations:

Conservation of mass:

3H + + = 0 (3.15)
at 3x ay

Conservation of momentum in x-direction:
8qx (qx2/H) 8(qx q /H)
x + + fq
Tt x dy y
s 2
+ g H-+ ,- + 1 2 + AP
ax p ax 2p 9x

+air C U U Cf (q 2 + q 2)1/2 qx
+ CD I10' 1x f x y H2

aF aF
7Txx __ yx- M = 0 (3.16)
ax ax x

Conservation of momentum in y-direction:

qy s (qxq/H) a(qy2/H)
+ + + fq
7t Dx y x

+ g H ar+ H-P a + aH2Ap
;y p Vy 2p ay

+ Pir C iU101 Uly Cf (q 2 q2)/2 2
PO D 0 10y f x y H 2
0 C
aF aF
x- 0 (3.17)

Numerical approximation procedures used in CAFE-1 will not be discussed

in detail here. For details see Wang and Connor (1975).

3.1.3 Turbulent Eddy Viscosity Coefficient Eij

The eddy viscosity coefficient Eij is the major variable parameter by

which the engineer can adjust the model results to fit a data set. For

this reason a separate section is devoted to this parameter here. Model
validity and predictability therefore hinge on the credibility of the

values assigned to F... A brief literature review was made to determine

what reasonable values might be. It has been shown that Fij is a function

both of the velocity field and characteristics of the numerical solution


A "rough" formula for estimating Exx is given by Wang and Connor

(1975), as

E ~ ag g ai (3.18)

a = 0.1 to 0.01

g = 9.81 m s-2 (3.19)

S= expected tidal range

S= expected velocity

Al = characteristic grid length

For a < 0.02 the only effect of this term is to dampen short wave noise.

Note that ui, i., Eii refer to components in direction i. A better

approximation for the velocity field contribution might be the difference

in characteristic velocities across Aki.

Few authors have provided any insight into the rational for choosing

the eddy viscosity coefficient values they have used. Some sample values

are produced in Table 3.1. Since grid sizes varied so greatly in different

applications, a nondimensionalized value, Eti, defined by

E.. At
Et ii (3.20)
ii A 2

is introduced as more appropriate for comparison. It can be seen in the

table that values for this parameter range from .0025 to .002, with credible

values being in the range .001 to .05 (depending on the value of a used).

As noted by Wang and Connor (1975), ideally Eii should be minimized with a

having a highest practical value for damping of about 0.02. High values of

Et. indicate the modeler may have forced stability on the system, although

many other factors may also be important in such cases.

Table 3.1.

Values of Turbulent Eddy Viscosity Coefficients for
Numerical Models.

Ati, t, n, ui, Eii, E.

Source m sec m m m2s-1

(Steele et al., 250-500 100 2 0.5 500 0.2-0.8

Great Bay, N.H. 250 10 1.2 1.5 36 0.0057
(Cellikol and
Reichard, 1976)

Portsmouth 150-250 5 1.25 0.5 20 0.0016-
Harbor, N.H. 0.0044
(Cellikol and
Reichard, 1976)

Massachusetts Bay 5000 100 2 0.2 10,000 0.04
(Connor and Wang,

For the Apalachicola grid some reasonable estimates are, for the

smallest element, with low velocities

a = 0.02

g = 9.81 m s"2

T 1 m
u ~ 0.1 m s

Aai ~ 350 m at least

Ei. z (0.02) (34335) = 686.7 m2s1

Eii At 686.7- 60 0 33
35o -0.33
A t 350

For the largest elements, with highest velocities,

u 0.5 m s

AP. ~ 3000 m

Ei = (0.02) (58860) = 1177 m2s-1

Eii At 1177 60
A.2 3000

In fact, values of E = E = 40 m2s1 and E = 2- m2s-1 were

selected for subsequent model runs to minimize impact of this parameter

and assure a search to define the important physical variables in the

system. Values of Eti for these values of Eii are

1) smallest of elements

Ei = .02

2) largest elements

Et. = .0002

In general the personnel on this project preferred grid manipulation

to increasing Eii as a means to achieve stability in CAFE-1. The decision

was made to use the lowest appropriate values found in the literature and

design a stable grid around these. The procedures used to obtain a stable

grid will be discussed in the next section.

3.1.4 Application and Debugging of CAFE-1 Model for Apalachicola

Application of the CAFE-1 model to Apalachicola Bay turned out to be

a trial and error process which took somewhat longer than anticipated. A

finite learning curve exists however, since it is now possible to change the

grid in a timespan of 2-3 days, whereas it took about 4 months to get the

first successful run.

Apalachicola Bay is complicated by the fact that there are several

tidal inlets and several river inflow locations. The application was further

complicated by the fact that the dispersion model required good resolution

of East Bay and East and West Bayous. This required an extreme range of

element sizes, which tends to make the model unstable.

Finding a stable grid is something which was found best learned by

trail. It becomes intuitive to some degree, and this is a drawback of the

model in its current form. An initial "ideal" grid proved to be completely

unstable. Then a grid was made in which elements were made more equilateral,

and element sizes varied slowly. This was gradually modified to the working

grid by changing several unstable areas. To accomplish this the graphics

routines were found very helpful.

Unlike DISPER-1, CAFE-1 tends to die quite rapidly when an instability

occurs. The initiationof the instabilities could be traced by reviewing

velocities at all points for time steps leading up to the instability. This

is, however, very inefficient. Therefore, graphics routines for the Gould

(essentially the same as CALCOMP) plotter were developed to enable plotting

velocity vectors at desired time steps. It was then very easy to see

where instabilities occurred. The problem was usually one of the following


1) improper node or element definition

2) total depth less than about 1 dm at low water

3) grid changed size too rapidly

4) irregular triangle shape

5) grid elements too large in zone of rapid velocity change

6) too rapid change of depth.

Since instabilities tended to occur sequentially, the procedure

usually consisted of fixing one problem at a time. As the total number of

potential instabilities was not known, the process was discouraging at


Despite the fact that an implicit time step routine is used in CAFE-1,

there appears to be a limit on the time step imposed by stability

considerations. This criterion seems to take the form of the Courant-

Friedrichs-Lewy (CFL) condition. The two-dimensional form of the Courant

number is given by Cunge (1977) as

C at gR (3.21)
r Ax

In theory, an implicit scheme should be stable for any value of Cr, even

considerably greater than 1, although the accuracy of the solution will

generally decrease as Cr approaches higher values. The theoretical

stability of the model is, however, tempered by real physical features.

For example, Cunge (1977) discusses the Leendertse (1967) scheme, which is

considered a strictly non-dissipative, second-order accuracy scheme. This


means that no numerical damping of wave amplitudes will occur, which is

desirable for real waves. However, this also means that discontinuities,

perturbations, rough geometry (rapid depth changes, etc.) will create

numerical waves which are not damped. Therefore, a theoretically stable

scheme may become unstable in some applications. It should also be noted

that the theoretical stability limits are always obtained by a linear

stability analysis of the von Neumann type which can best be considered a

guide to stability of the nonlinear equation system.

Considerable numerical experimentation led to the adoption of a

practical CFL-type criterion given by

At < 0.7 Ax /1TRW (3.22)
0 0

for this model in which x0 and H refer to the grid location giving the

critical value for At.

It was found that for all grids thus far used, Equation (3.22) yielded

values close to 60 seconds, and a time step of 60 seconds for CAFE-1 runs

was in fact found acceptable and will be used for all results presented in

this report.

3.2 The DISPER-1 Model

The selection of DISPER-1 was rather direct once CAFE-1 was

selected as the hydrodynamics model. It is 2-D, real-time, finite-element,

compatible with CAFE-1, and readily available. Because advection dominates

dispersion in tidal systems, a real-time model is required. This need is

particularly strong given the transient character of stormwater quality and

quantity. Availability of real-time capability will be extremely beneficial

when any biological, chemical, or other time-dependent constituent changes

and interactions are added to the model at a later date.

DISPER-1 is a real-time, 2-D, vertically-averaged finite element model

for solution of the convective-diffusion equation given hydrodynamic inputs

from another source (in this case from CAFE-1). The model is described by

Christodoulou et al. (1976), Leimkuhler et al. (1975), and in a user's

manual by Pagenkopf et al. (1976).

3.2.1 Equations and Boundary Conditions for DISPER-1

The model solves the classical convective-diffusion equations:

Following Leimkuhler et al. (1976)

ac + (uc) + (vc) + P (3.23)
at ax ay ax x ay y +

in which
S- apHD C aC (3.24)
ix iP H Dxx ax- P H Dy (3.24)

p H c pHD (3.25)
y H xy ax xy ay

c = p cH (3.26)

c = average concentration of a constituent

H = total depth

p = density of water

u,v = vertically averaged velocities in x and y
directions, respectively

c = vertically integrated concentrations

P = sources and sinks of mass
Dii = dispersion coefficient

Note that Equation (3.23) expresses conservation of mass of the constituent

of interest. Values of u and v came from CAFE-1.

Either a fixed concentration or a fixed mass flux can be specified as

a boundary condition. These may be specified on elements or nodes or some

combination thereof.

It is important to bear in mind that this is a finite element scheme

when loading the boundaries. Linear interpolation functions are used, so if

a smooth unit load of 1 kg s-l is to be placed across 3 points, then loading

will have to be assigned as 0.25 kg s"l on each of the outer points and 0.50

kg s-1 on the central point of the three. Also, if the grid has a depth

of 2 m, then the initial concentration of the above scheme will be 1/2 that

for a grid of 1 m depth. This is because the model solves for depth-

intergrated values, as noted in Equation (3.26).

At the present time boundary conditions may vary in DISPER, but only

if they are specified explicitly over the tidal cycle or over time. This

is inconvenient for many problems (and especially for the determination of

salinity) since the solution at the boundary must be specified before it

is known. Known applications avoid this problem by assigning a decay value,

or by having the source sufficiently weak that complete dilution occurs, so

that c is equal to zero at all times at the boundary. Another useful

approach is the floating boundary condition developed by Dailey and Harleman

(1972) which sets the boundary concentrations at ocean value during flood

tide, and specifies that the gradient of the dispersive flux during ebb

tide remains constant.

3.2.2 Stability and Convergence of DISPER-1

The best work in the areaof stability and convergence of DISPER-1 has

been done by Christodoulou et al. (1976), and experience at the University

of Florida has helped strengthen their findings for application of DISPER-1

to Apalachicola Bay. In the dispersion model, three factors combine to

change constituent concentrations. These are advection (represented by

the velocity, u), dispersion (represented by the dispersion coefficient,

Dii), and decay (represented by the appropriate decay or reaction coeffi-

cients). It has been found that typical decay terms create slower changes

than the other two terms and can generally be neglected in stability

analyses. Christodoulou et al. (1976) found that the total effect of

advection and dispersion could be represented in defining a "safe" region

given by the following inequality.

2 D. .At 2
(1.22 u t) + (8 ) < 1 (3.27)
Ax Ax2

It should be noted that Equation (3.27) is based on a theoretical analysis

assuming very regular and equilateral triangles and verification has

occurred only on moderately irregular grids. For highly irregular grids,

the allowable time step may be considerably below that given by Equation

(3.27). Equation (3.27) can be converted to
2 D.. 2
At < (1.22 -) + (8 ) ] (3.28)
Ax Ax2

2 -1
As an example, for element sizes of about 400 mm, Dii of 100 m s and a

maximum velocity of about 1.5 m s-l, Equation (3.28) yields a value of

147 s. This, too, may need to be further reduced due to other problems.

However, in general it is true that DISPER-1 is stable at longer time steps

than CAFE-1.

Once the models have been made stable for a given grid, then one

must turn to the question of accuracy of the solution or convergence to the

true solution. No specific criteria have been developed in this regard

for CAFE-1, although it is generally believed that a stable solution is

also accurate unless unexplained oscillations persist in any portion of the

flow region.

Christodoulou et al. (1976) proposed an accuracy criterion for

DISPER-1, given by inequality (3.29)

2 D..
Ax < (3.29)

Inaccuracies in DISPER-1 exhibit themselves as negative concentrations and

as ocillating values of concentrations. For Dii = 100 m2s- and a maximum

velocity of u = 1.5 m s- it can be seen that Equation (3.29) yields a

very small element size of about 130 m. If in fact such a small element

size were chosen, then values of the time increment would be drastically

lowered, as shown by Equation (3.29) for DISPER-1 and Equation (3.22) for


Fortunately, Equation (3.29) applies primarily in areas of high con-

-certration gradients and high concentrations. In regions far from the

sources, where concentrations remain lower, a larger length scale may

function satisfactorily. It should be realized, however, that the grid

size may have to be modified, and hence the time step, to eliminate oscil-

lations and non-negligible negative concentrations.

Several other features may play a role in determining model behavior

but have not been cast in any quantitative criterion. A few of interest

here will be mentioned. It has been shown that the schematization of the

sources) in DISPER-1 is very important. In one example given by

Christodoulou et al. (1976), the criteria given by Equations (3.28) and

(3.29) were both met. In two runs, all parameters were the same except

that the source was distributed over two elements in one run and over eight

elements in the other run. The run with eight elements for the source gave

completely satisfactory results, while those from the other runs exhibited

excessive oscillation. Therefore sources should be distributed over

several elements, but it should be realized that this often means that model

predictions very near the source may be less valid.

The element shapes and gradation of element sizes have been observed

by Hydraulic Laboratory personnel to be very important also. Elements too

different from an equilaterial triangle may cause problems. In addition,

problems may occur where element sizes change too rapidly. This also

includes too rapid a change in depth, even where element sizes and shapes

are fairly regular.

A factor in stability and accuracy is the initial condition chosen for

the bay. A common approach used in earlier Apalachicola Bay runs is the

so-called cold-start, with zero concentrations, e.g., specified in DISPER-1.

The model can, however, be operated in a hot-start mode, where concentrations

are specified throughout the grid at time zero. Of course, in either case,

the model is run sufficiently long (through several tidal cycles) to remove

any bias provided by the assumed initial conditions. However, instabilities

or inaccuracies could occur due to the extreme gradient of concentration

setup at the onset of model operation. If this occurs, it may be more

economical to specify a realistic initial concentration field as opposed

to severely restricting element size and time increment.


As noted in Section 3, the two models being used have been verified

by application to other waterbodies and comparing results to measured data.

Of course, as is always the case, verification may be in the eye of the

beholder. In the following,verifications of the models as used in the

Apalachicola area will be discussed.

4.1 Verification of CAFE-1

4.1.1 Field Data

Two primary types of data were obtained from the field: (1) tidal

information at various points in the bay and (2) water velocity measurements.

In addition, the wind velocity was measured several times during the data

collection process.

Most of the field data to be discussed in this report were gathered

from September 15 through September 26, 1980, and from November 7 through

November 9, 1980. Tide Recordings. Tide gages were place at each of the four

ocean boundaries: East Pass, Sikes Cut, West Pass, and Indian Pass. The

tides recorded at these gages were used as input boundary information to the


In addition, tide gages were installed at two internal points in the

bay (points A and B in Figure 4.1). The recordings from these two tide

gates were used as a comparison with the model results. Figure 4.2 shows a

tide gage and recorder in place on a private dock at Indian Pass. A sample

recording of the tide at East Pass is shown in Figure 4.3.

Additional tidal information was obtained from NOAA.

: "* i: ...i-:- :.- -. .::" D O G
: ..; .. ,. t.N .. A.. IS.L~N 5ISLAND
". .. o ". i .-.',.'" .

S"" :"giN ISLAND 0LF
:: : :.. : :

.::.P :-:ST:' ...TR. i.t O +

E3 :: AY


10 5 0 10km

Figure 4.1. Location of Verification Points in Apalachicola Bay.











-4 U,


52 rr

I i-


L re -

ri I

: I ~- i


I L~

17 K

r-t I-t

L- P r
it i I i i

Lj tv

L.. L 4. F:2' iV

JL iir

rr r

r r r

ir I.
i- t L W j hi
i r

r L r.

r 7
C I C L FV L ..~
t~I r r I

t T- t- I -



Figure 4.3. Sample Recording from Tide Gage at East Pass.





I: I

i Ei

*EF *r



r 7 .L..

r II-

1 P-I


I~--4~F~ .~-i4

-. L

.r-f +

i~T'FI-- t'-;

L i I

r : '

"" Velocity Measurements. All velocity measurements were

taken from a boat using a V "Arkansas" Ott Meter (see Figures 4.4 and 4.5).

Logarithmic velocity profiles were assumed, so the velocity readings to

be compared with the model output were taken at 36.8% of the total depth

from the bottom (which is the location of the spatial mean velocity for a

logarithmic velocity profile). Velocities were also recorded at other

depths so that velocity profiles could be plotted and the assumption of

logarithmic velocity profile checked.

At points C and D (Figure 4.1) measurements were taken at short

intervals for periods of about four hours. These values were used as a

comparison with the model results. Velocity profiles were also measured

at several other points in the bay.

Flow direction was measured with a compass after hoisting the Ott

Meter to an elevation where it could be seen (see Figure 4.6). It is

assumed that the flow direction is the same at all depths below the water


4.1.2 Model Input for Verification Runs Grid Configuration. The finite element grid shown in

Figure 3.1 was used for all of the computer runs to be discussed in this

report. It contains 281 nodes and 439 elements.

This particular grid was chosen after experimentation with several

grids because of its many desirable characteristics. First, the element

size is small enough to provide an accurate representation of the bay.

Second, most of the triangles are approximately equilateral and there is

no rapid change in element size. This is of particular importance because

rapid change in element size as well as irregular triangles, tends to produce

instabilities which cannot be predicted by quantifiable stability criteria as

discussed earlier.


Figure 4.6.

Observation of Direction of Water Flow.

-I*-CC~.~"IF*~C' IC ~--C~*r~lL~CPd

The remaining reasons for choosing this grid are related to the cost

of running the model. This cost is roughly proportional to the number of

elements and inversely proportional to the timestep.

The timestep was chosen to be 60 seconds in accordance with the earlier

discussed stability criteria. In addition, the output from CAFE-1 is used

as input to a water quality model and 60 seconds is a convenient interval

for transferring this information. Turbulent Eddy Viscosity Coefficients. As stated earlier,

internal stresses from turbulence and velocity shear are represented in

the model in the form of Equation (3.9). Since there is no way to actually

measure these stresses, the literature review provided reasonable values

for the eddy viscosity coefficients. The values currently being used are

the following:

E = 40.0 ms

E = 40.0 m3s

E = 20.0 m2s1

The sensitivity of the model to these coefficients was not known.

A single computer run was made with each of these values cut in half. The

resultant change in velocities was only about 1-2%. The turbulent eddy

viscosity coefficients are parameters that can be easily manipulated by

the operator for the purpose of calibrating the model. Therefore, a more

detailed study on the effects of making larger changes in these parameters

may easily be undertaken if necessary. Manning's n. The first set of computer runs was made with

n = 0.045. This value was calculated in accordance with the method pro-

posed by Christensen (1975) from a few scattered velocity profiles obtained

in the field. From the velocity profiles that were logarithmic (surprisingly

few of them were not), an average of the n values was taken. This average

was input as a constant value for the entire bay.

Since the first set of runs yielded velocities that were too small, a

second set of runs was made with n equal to a constant 0.030 throughout the

bay. The result of this change will be discussed in more detail in the

next section.

It is also possible to input varying values of Manning's n throughout

the bay. Although this is more time consuming than simply supplying a

constant value, it is the only way to account for differences in bottom

friction at various locations. This type of input can be used as a final

step in "fine-tuning" the model. River Flow. Although river flow is a physical characteristic

of the bay, the exact river discharges are not known for the periods during

which data was taken. However, average monthly flows for the Apalachicola

River (1961 1976) were available from the U.S.G.S. so these values were

used (see Table 4.1 and Figure 1.7). Fortunately, errors in the river

flow do not significantly affect the overall hydrodynamics of the bay, because

the river flow constitutes a very small part of the total amount of water

entering the bay. This is illustrated in the following calculation showing

the ratio of river flow (Q = 705 m s for the Apalachicola River) during a

half tidal cycle to a mean tidal prism height of 0.5 m:

Table 4.1. Average Monthly Flows for the Apalachicola River for Period

Discharge Ratio to Mean
Month m3 s1 Annual Discharge

0 387 .55

N 391 .56

D 617 .88

J 985 1.40

F 1131 1.60

M 1209 1.72

A 1082 1.54

M 687 .97

J 580 .82

J 505 .71

A 496 .70

S 392 .56

Mean Annual Discharge 705 1.00

These values correspond to the hydrograph given in Figure 1.7.

(100%) (705 m3s ) (60 s min ) (60 min h- ) (6.21 h) = 5.7
(550 km2) (0.5 m) (1,000,000 m2 km 2) Wind Velocity. The model does not provide for the input of

time-varying wind velocity, so average velocities were used. For the run

used as a comparison with the internal water surface fluctuations:

wind velocity = 2.24 m s-l (5 mph) at 2070 clockwise from North

For the two runs used as comparisons with water velocities:

wind velocity = 2.24 m s-1 (5 mph) at 0 clockwise from North Tide Information at Boundaries. The tidal curves at the

boundaries are input to the model as a set of amplitudes and times. The

model then approximates the tides as a series of sinusoidal curves. It is

important that the first half tide curve input to the model is one of in-

creasing water elevation (i.e., from low tide to high tide) since the

initial depths in the model are at MLW.

Complete tidal information at all boundaries was available for the

run used as a comparison with the internal water surface fluctuations.

However, for the two runs used as comparisons with the water velocities,

only partial boundary data was available. (Tide gages were only set-up at

East Pass and Indian River Pass when the velocity data were gathered.)

Correction factors for both amplitude and time lag were established from

the complete set of curves and used to estimate the tidal curves at West

Pass and Sikes Cut for the remaining runs from the information at East Pass. Depths at MLI. The depths input at each node were obtained

from a Navigational Chart published by NOAA. The depths input at West Pass

and Sikes Cut are slightly lower than the actual values, because the latter

causes instabilities in the model with the 60 second time step. Using a

time step small enough to permit the input of actual depths would make the

cost of running the model prohibitive. To compensate for the smaller

depths, larger than actual widths were input so that the cross-sectional

area of the outlets remained the same.

4.1.3 Comparison of Model Results with Field Data

Figures 4.7 and 4.8 show field measurements and model results for the

water surface fluctuations at two points,B and A,in the bay. Computer

results are shown for n = 0.045 and n = 0.030. At point B, both computer

runs are extremely close to the range and phase lag of the actual curve.

The change in Manning's n produced only about a 5 percent change in the

tidal range and almost no difference in the phase lag. This is important

to note, since the same change induced a much larger difference in

velocities. The difference between the range of the computer run with

n = 0.030 and the actual measured range is only about 7 percent.

At point A, there is almost exact correlation between the range of the

n = 0.030 computer run and the range of the recorded curve. However, there

appears to be a slight phase lag difference between the two (about 10


Similar correlations between measured and computed velocity and

velocity orientation are obtained with the model. Figure 4.9 shows the

comparison at point C as an example. Observed and computed values

of the vertically averaged velocities at points A through B show

similar satisfactory agreement. The predicted model velocities are

often lower than the observed velocities. However, Wang (1978) reports

that a comparison of current velocities measured by portable current

S -- MODEL, n =0.030

O /


0 T/4 T/2 3 T/4 T


crM LW \


0 T/4 T/2 3T/4 T


Figure 4.7. Water Surface Elevation vs Time at Point B in Apalachicola Bay.


0 0.2

MO \D-, 7 n.


-" 0.1
-0.1 1 I _____1
- MODEL, n=0.045
S---- MODEL, n =0.030

0 T/4 T/2 3T/4 T


Figure 4.8. Water Surface Elevation vs Time at Point A in Apalachicola Bay.

--- MODEL, n=0.045
300 _
30 MODEL, n=0.030

o //
zp 200 -

o 100 /
o-0 /


0 T/4 T/2 3T/4 T


Figure 4.9. Comparison of Observed and Computed Velocity Direction as Function of Time at Point C.

meters from boats with those measured by recording current meters mounted

on fixed structures showed that the former tend to be inflated due to wave

motion. Hence, the agreement between observed and model predicted

velocities may be even better than the observations indicate.

4.1.4 Comparison with Wang's Verification

Figures 4.10 and 4.11 show field and preliminary model results

obtained by Wang (1978) using the CAFE-1 model at Biscayne Bay. Wang

considers his results to be good for an initial computer run with no

adjustments. All of the results for Apalachicola Bay discussed in the

previous section fall well within the errors illustrated in the two

figures. Therefore, it may be assumed that CAFE-1 correctly predicts the

principal physical processes in Apalachicola Bay. However, further refine-

ment is still possible as more detailed data in specific areas are obtained.

4.1.5 Further Calibration

Because of the large change in velocities occurring when n was de-

creased from 0.045 to 0.030, it was estimated that further decreasing n

to 0.015 would cause the model to output velocities of magnitudes approxi-

mately equal to those measured in the field. However, use of n = 0.015

results stability problems, presumably because the larger velocities

produced violated the stability criterion. Since the instabilities occurred

early in the run (about halfway through the first tidal cycle), it is

possible that the velocities causing the stability problems are produced

only as a result of the "cold-start" (i.e., all initial velocities and

water surface elevations set equal to zero). Assuming that the steady-

state velocities are smaller in magnitude than these initial velocities,

it may be possible to eliminate the instabilities by inputting initial


--- I (DATA)




----- H (OATA)

Figure 4.10. Surface Displacement Verification by Wang for a Model of
Biscayne Bay Using CAFE-1 (Wang, 1978b).



Figure 4.11. Hodograph Comparison by Wang for a Model of Biscayne Bay
Using Cafe-1 (Wang, 1978b).



values for the velocities and water-surface elevations. These values can

be estimated by using the model output from past runs. Although it will

take some time to input these values initially, this change will probably

decrease computer cost because less computer time will be required for the

model to reach a steady-state. (The model is currently being run for

three to five tidal cycles before the final information is output).

4.2 Verification of DISPER-1

The verification of a transport model can take many forms. Field

data can be obtained for releases of a dye or other tracer from the river

mouth or other sources. In an estuary of this size, the problems with such

measurements can be simply too great in terms of time and manpower. A

reasonable alternative is to model some naturally occurring substance which

can be monitored by point samples through the bay, rather than following a

tracer cloud. The best choice is salinity. While large amounts of salinity

data exist on the bay, unfortunately concurrent tidal and velocity are data

rarely available. Therefore, additional data was obtained for this study.

4.2.1 Field Data for Quality Verification

The specific set of salinity data to be treated here was taken in

September, 1980, by a field crew from the Hydraulic Laboratory at the

University of Florida and analyzed in Gainesville. The data are reported

in Table 4.2.

A review of the data in Table 4.2 indicates both spatial and temporal

variation of salinity; see, for example, the values at East Pass at two

different times. Notice also some apparent small anomalies in the data

which evidently reflect measurement error, symptomatic of difficulties of

sampling in such an environment. It is also interesting that samples 22 to

26 at East Pass are inverted from that expected, indicating possible

Table 4.2. Salinities in Apalachicola Bay, September 1980. Observed by
Hydraulic Laboratory, University of Florida.

Sample Approximate Salinity Date, Time
No. Location ppt of Sample

1 Apalachicola River Mouth 1.6s 9-18-80, 1010
2 Apalachicola River Mouth 21.2 9-18-80, 1010
3 Center of Bay 30.43 9-18-80, 1100
4 Center of Bay 30.0 9-18-80, 1100
5 Center of Bay 31.6 9-18-80, 1100
6 Center of Bay 28.8 9-18-80, 1100
7 Center of Bay 31.2b 9-18-80, 1100
8 West Pass 33.2s 9-18-80, 1630
9 West Pass 34.0 9-18-80, 1630
10 West Pass 29.2 9-18-80, 1630
11 West Pass 33.8 9-18-80, 1630
12 West Pass 32.4b 9-18-80, 1630
13 Sikes Cut 29.8s 9-18-80, 1730
14 Sikes Cut 30.0 9-18-80, 1730
15 Sikes Cut 31.8 9-18-80, 1730
16 Sikes Cut 32.0 9-18-80, 1730
17 Sikes Cut 33.0b 9-18-80, 1730
18 East Pass 35.0s 9-19-80, 1140
19 East Pass 35.0 9-19-80, 1140
20 East Pass 34.4 9-19-80, 1140
21 East Pass 34.8b 9-19-80, 1140
22 East Pass 36.8b 9-26-80, 1500
23 East Pass 35.8 9-26-80, 1500
24 East Pass 38.2 9-26-80, 1500
25 East Pass 44.3 9-26-80, 1500
26 East Pass 43.5s 9-26-80, 1500
27 New River 28.7 9-26-80, 1600
28 South Span-Causeway 34.9s 9-26-80, 1820

s = surface
b = bottom

Table 4.2 continued.

Sample Approximate Salinity Date, Time
No. Location ppt of Sample

29 South Span-Causeway 36.1 9-26-80, 1820
30 South Span-Causeway 35.0 9-26-80, 1820
31 South Span-Causeway 37.8 9-26-80, 1820
32 South Span-Causeway 35.5b 9-26-80, 1820
33 South Span-Causeway 37.5s 9-26-80, 1820
34 North Span 30.4s 9-27-80, 1150
35 North Span 34.4 9-27-80, 1150
36 North Span 34.8b 9-27-80, 1150
37 St. George Sound 38.4s 9-27-80, 1300
38 St. George Sound 37.4 9-27-80, 1300
39 St. George Sound 37.8 9-27-80, 1300
40 St. George Sound 38.6 9-27-80, 1300
41 St. George Sound 41.1b 9-27-80, 1300
42 St. Vincent Sound 25.0s 9-28-80, 1200
43 St. Vincent Sound 22.6 9-28-80, 1200
44 St. Vincent Sound 24.1b 9-28-80, 1200

s = surface
b = bottom

upwelling or local disturbances. This sort of variability should be recalled

when assessing model agreement with the data. It should further be noted

that the model, being two-dimensional, yields a depth-averaged concentration


4.2.2 Model Input for Verification Runs

Due to the time span in salinity measurements, typical tides and

winds for that period were specified for a CAFE run to form hydrodynamic

input to DISPER. While this of course may lead to some inadequacies in the

simulation, it was perceived as a good test, for uncertainty in these

model parameters often exists. The 439 element, 281 node general grid shown

in Figure 3.1 was used again for the simulation. Salinities of 36 ppt were

specified at the external passes and inlets, with 10 ppt and 5 ppt at appro-
3 -1
private elements for the river mouth. River flow was taken at 337 m s

during this period,based on USGS records.

Some oscillations of concentrations were observed in the early simu-

lation. As noted earlier, both the spatial increment, Ax, and the time

increment, At, bear on this problem. Changing Ax means a change in the

grid, locating critical areas and then refining the grid there.

Oscillations in the results caused some experimentation with the dis-

persion coefficient, Dii, resulting in specifying input values which

typically yield maximum values of Dii of 150-300 m2s-1 in the bay.

These values are within the range of reported values for estuaries.

Further accuracy considerations resulted in selection of a time step, At,

of 60 seconds. These steps were chosen first in the verification process

as they represent reasonable options for the typical user.

4.2.3 Results of DISPER-1 Verification Runs

Results of one verification run are shown in Figure 4.12. Values

shown in the elements are salinities in parts per thousand. It can be seen

by comparison with Table 4.2 that there is reasonable agreement with

measured.values in the East and West Pass regions, in St. Vincent Sound,

and in the East Bay East Point region where the freshwater input from the

Apalachicola River is so important. There are, however, some regions where

values are too low, especially near the causeway island and slightly east

of that region.

There are several possible explanations for these low values in the

regions noted. First, it should be observed that the figures shows a

synoptic view of the bay, i.e., at a single time, while the observations in

Table 4.2 were taken over a range of times and therefore occur at different

parts of the tidal cycle. What is more important, however, is the fact that

measurements spanned such a long time that varying wind, tides, and other

conditions can significantly influence bay behavior. The run pictured here

is based on a constant wind speed and direction. In addition, no attempts

have been made to fine-tune CAFE results by varying Manning's n across the

bay and other such steps which might help resolve some localized errors.

It appears that verification is easier against data taken synoptically,.

such as in enhanced LANDSAT photographs. Comparison of predicted and

observed shapes and extents of plumes, i.e., pattern recognition, can provide

a good assessment of model performance. Graham, Hill and Christensen

(1978) and Hill and Graham (1980) reported on use of such data for this

In conclusion, it can be stated that an attempt at model verification

with no fine-tuning yielded acceptable results over much of the bay.


'V. a 7P.u

.6. L
... .... -t I -t t + ,,.,

ME lE ;" '. F II j 23

Figure 4.12. Salinity Run for Apalachicola Bay Using DISPER-1.
Fur 4.12.c Salinit Ru o plcioa a sn IPR1

However, some changes in element arrangement would be necessary to remove

some locally low values near the causeway island.

4.3 Satellite Verification of DISPER-1

Use of LANDSAT imagery computer enhanced for surface water color may

be utilized for verification of the DISPER-1 model, since surface water

color and water quality (e.g., acidity expressed by pH) may be related.

Since the satellite images only reveal what is going on in the upper inches

of the water column,it may be difficult to let the images relate to the

vertically averaged water quality. However, the general pattern of the

dispersion of a pollutant at the surface may be observed and compared to

the print-out of the DISPER-1 model. Assuming a well-mixed bay,the surface

water quality should indeed be indicative of the vertically average quality


Verification by such pattern recognition from LANDSAT images has given

surprisingly good results.

Figure 4.13 shows an example of the color enhanced pictures used.

Note the red colored water in East Bay and at Carrabelle. This color indi-

cates low pH water. The pattern shown in Figure 4.13 is in almost perfect

agreement with computer printouts. The plume at Sikes Cut is also note-

worthy. Also this phenomenon may be reproduced by the model.

Figure 4.13. Computer Enhanced LANDSAT Image Showing Water Quality
Patterns' in Apalachicola Bay.



To illustrate more clearly behavior in the bay, model solutions

demonstrating variation through a typical year were sought. The objective

is to provide an overview of possible behavior,with the expectation that

specific problems will still require running the model for conditions

appropriate to those problems. The material generated for the average

year is presented in the form of an atlas at the end of this report.

There are several sets of physical parameters to be selected to

provide a view of bay behavior. Major ones include the following:

1) Tides amplitude and phase lag.

2) River flows.

3) Wind speed and direction.

4) Source location and loading.

To provide a reasonable view of expected bay behavior and yet keep the

size of this report reasonable, it was decided to produce runs for each

month, or twelve in all. Therefore, parameter values selected should be

expected to represent monthly average values.

5.1 River Flows

The average monthly river flows shown in Table 4.1 were used for the

CAFE-1 simulations.

5.2 Tides

Tidal data was obtained from the National Ocean Survey for the five-

year period from 1975-1980. Tide tables and the tidal data discussed in

Section were used to estimate expected tidal height variations from

points of measured values to other boundary points and time (phase) lag

between tidal peaks.

The tide at Apalachicola is considered diurnal, in that two highs or

two lows may occur on the same day with different amplitudes. This is

clearly indicated in Figure 4.3. There are difficulties and uncertainties

associated with selecting a "typical" tide pattern for each month of the

form in Figure 4.3. In addition, the runs of interest, and many problems

of practical interest, occur over several tidal cycles. This tends to

average out the influence of tidal amplitude variation. Therefore, given

the objective of providing an overview of bay behavior, it was decided to

simply use mean tidal amplitudes to drive the model, with a respecting

sinusoidal tide specified. This should reproduce the general structure

well, although there may be small local differences at individual times

within a tidal cycle.

A review of the tidal data from 1975-1980 showed mean tidal ranges at

the Apalachicola gage rather close to one foot for all months. The lowest

value obtained was 0.86 ft and the highest 1.10 ft. No consistent trend

appeared by month, at least for the five years reviewed. It was therefore

decided to simply use a single-value of 1.0 ft (or a tidal amplitude of

0.50 ft = 0.155 m). Review of tide tables and UF measured data then led to

the tides at the bay boundaries, as shown in Table 5.1.

Table 5.1. Tidal Amplitudes and Phase Lags for Model Runs.

Tidal Phase
Location amplitude, m lag, s

East Pass 0.23 0
Sikes Cut 0.19 2520
West Pass 0.13 3780
Indian Pass 0.10 5400

Again, these values are taken as typical and should not be interpreted

as representing any specific case.

5.3 Winds

Wind is an extremely important factor in bay behavior. Wind data from

1975 through early 1981 was obtained from the Apalachicola office of the

National Weather Service in the form of monthly summaries of Local

Climatological Data, with wind speed and direction reported at three-hour

intervals, as well as daily and monthly average speeds and directions.

A review of the data led to selection of the monthly average wind values

shown in Table 5.2.

Table 5.2.

Monthly Average
Model Runs.

Values of Wind Speed and Direction Used in

a: For example, 090 is
the South, etc.

a wind from

the East, 180 from

Wind direction,
Location Speed, mi/hr deg. from N

January 8.4 030
February 8.2 073
March 8.4 186
April 7.9 192
May 7.4 182
June 7.0 197
July 6.3 217
August 6.1 130
September 6.7 110
October 7.15 063
November 7.3 078
December 7.4 045

Variation in wind speed shown in Table 5.2 seems less significant

than direction. Note the tendency of the wind to be from slightly west of

south during March-July, but more like easterly to northeasterly during

the rest of the year.

5.4 Generated Results

The model CAFE-1 is run with the grid shown in Figure 4.7 for each of

the twelve months, utilizing river flows from Table 4.1, tidal information

from Table 5.1 (the same for all months), and the wind data from Table 5.2.

The results will be presented as views of the velocity field throughout

a tidal cycle and the net velocity field over a cycle. Output from CAFE-1

runs for each of the twelve months will then be used to drive a DISPER-1

run representing a continuous discharge concentration of 100 units (ppm,

eg.) in the river water.

The results from the conservative pollutant can be scaled directly

to salinity by assuming an average salinity at the river mouth of 7.5 ppt

(typical) with 36 ppt at the ocean boundary. Then any reported concentration

can be converted to salinity by

s = 36 0.285 c (5.1)

in which c = concentration projected by model

s = salinity in ppt

Full results are presented in the attached Atlas in 204 maps depicting

hydrodynamics as well as water quality during the "average" year.


The computer prepared maps,given in the atlas showing the distribution

of conservative pollutant concentrations c in the Apalachicola Bay, are based

on the river concentration cR = 100 (e.g., ppm) and the ocean (Gulf of

Mexico) concentration co = 0.

Pollutant concentrations, cl, corresponding to other river and ocean

concentrations may be found from these maps by use of a simple conversion

formula to be established in the following.

Since the c distribution,shown in the maps and schematically in the

upper part of Figure 6.1,is a solution to the differential equations governing

the migration of conservative pollutants in the Apalachicola Bay,it is

easily proven that the concentration cI given by the linear expression

c = ac + b -(6.1)

also must be a solution. In Equation (6.1) a and b are constants.

A special case of Equation (6.1) is

c c
cI = c 00 + c (6.2)

in which cR and co are arbitrary pollutant concentrations in the river and

ocean, respectively.

This equation satisfies the boundary conditions

cl = c for cR = 100 and co = 0 corresponding

to the original c-mapping,






C eBAY P a



Figure 6.1. Schematic Maps Showing Boundary Conditions Used in Concentration
Conversion Formulas.

c = c for c = 0
1 o
c = c for cR = c
and c1 = cR for co = 100

Equation (6.2) will therefore give the pollutant concentration at any

point where c, based on cR = 100 and co = 0, is known when the river and

ocean pollutant levels are known. The maps presented in the attached atlas

will therefore, when combined with Equation (6.2), yield the vertical

average of the pollutant concentration at any location in the bay and at any

time during the tidal cycle for any values of cR and c .

In the same way the maps may be used to predict the salinity at any

point and at any time during the average tidal cycle.

The equation

s = ( c ) s (6.3)

in which s = the salinity at (x,y) at time t and s = salinity of the Gulf

of Mexico waters is of the same form as Equation (6.1). s is therefore a

solution to the equations governing the water quality of the bay. Since it

furthermore satisfies the boundary conditions

s = 0 for c = 100

and s = s for c = 0

it must represent the bay's salinity distribution. These boundary conditions

are indicated in the lower part of Figure 6.1.


The work presented in this report is a result of numerical modeling

research being performed by the University of Florida's Hydraulic Laboratory,

sponsored by the National Oceanic and Atmospheric Administration, Office of

Sea Grant, Department of Commerce, under Grant #NA80AA-D-00038,

Project #R/EM-13, and by the Engineering and Industrial Experiment Station

of the University of Florida. Support from the U.S. Army Corps of Engineers,

Mobile District, is also acknowledged.


The symbols used in this report are defined where they first appear

in the text and in the following list of symbols.

a = dimensionless constant

b = dimensionless constant

c = vertically averaged pollutant concentration at point (x,y),
at time t corresponding to cR = 100. Note in Secton 3 c
stands for the vertically integrated concentration while
the vertically averaged concentration there is denoted c

c = vertically averaged pollutant concentration at point (x,y)
at time t corresponding to cR f 100

c = vertically averaged pollutant concentration in Gulf of
S Mexico and entrances to the bay

cR = pollutant concentration of river water entering the bay

Cf = dimensionless friction coefficient

D = average estuary depth

Dii = dispersion coefficient

Dxx,Dxy,Dyy = dispersion coefficients
Et. = dimensionless turbulent eddy viscosity coefficients

Eij = turbulent eddy viscosity coefficients

f = Coriolis parameter

Fxx,Fxy,Fy = internal specific stress terms (stress/density)
g = 9.806 m s-2 = acceleration due to gravity

h = water depth below mean lowwater (MLW) at point (x,y)

H = h + n = total depth at point (x,y) at time t

H = refer to grid location giving critical At value

k = Nikuradse's equivalent sand roughness

L = estuary length

M ,M = momentum addition per unit horizontal area
x y

= Manning's n

= pressure

= pressure at

= pressure at

= tidal prism

P =

qxqy =
q =

Q =



s =


T =

u =

(u,v) =


x,y =

z =

Greek Letters

a =

AZ =

At =

Ax0 =

rl =

= kl/6/(8.25Vg) (metric)

the bottom

the water surface

sources and sinks of mass

discharge per unit width in x- and y-direction, respectively

volume addition rate

river discharge

tidal range

salinity at point (x,y) at time t

salinity in Gulf of Mexico at entrances to the bay


tidal period = 12.42 hr

expected velocity

vertically averaged velocity components in the x- and y-
direction, respectively, at point (x,y) at time t
air velocity at 10 m above the water surface in m s-

horizontal Cartesian coordinates

vertical coordinate

dimensionless coefficient

characteristic grid length

time increment

refer to grid location giving critical At value

elevation of water surface above mean low water (MLW)

expected tidal range

p = density of water

pair = density of air

Po = reference value of p

p,'p2 = characteristic extreme densities

T = bottom shear stress

T = surface shear stress

Earth = angular velocity of earth = 27r/(24 3600)



A detailed bibliography on the subject of numerical models is given

in Graham, D.S., DeCosta, K. and Christensen, B.A. (1978), "Stormwater

Runoff in the Apalachicola Estuary," Hydraulic Laboratory Report No. HY7802,

University of Florida, Gainesville, Florida.

Briggs, D.A. and Madsen, O.S. (1973), "Mathematical Models of the
Massachusetts Bay," Part II, Massachusetts Institute of Technology Report
No. MITSG 74-4, Cambridge, Massachusetts.

Cellikol, B. and Reichard, R. (1976), "Hydrodynamics Model of the Great
Bay Estuarine System," Part I, U.N.H. Mech. Res. Lab. Report No. UNH-SG-153,
Durham, New Hampshire.

Christensen, B.A. (1975), "Hydraulics and Water Resources Engineering,"
Unpublished Class Notes, University of Florida, Gainesville, Florida.

Christensen, B.A. (1979), "Rational Development and Management of Coastal
Drainage Basins by Composite Basin/Estuary Modeling," Proceedings of the
Eighteenth Congress of the International Association for Hydraulic
Research, Cagliari, Italy.

Christodoulou, G.C., Connor, J.J. and Pearce, B.R. (1976), "Modeling of
Dispersion in Stratified Waters," Report No. MITSG 76-14, Massachusetts
Institute of Technology, Cambridge, Massachusetts.

Christodoulou, G.C., Leimkuhler, W.F. and Ippen, A.T. (1974), "A Mathematical
Model for the Dispersion of Suspended Sediments in Coastal Waters,"
Massachusetts Institute of Technology, Parsons Laboratory,Technical Report
No. 179, Cambridge, Massachusetts.

Connor, J.J. and Wang, J.D. (1974), "A Numerical Evaluation of Hydrodynamic
Circulation in the Same Point Region of Narragausett Bay," Massachusetts
Institute of Technology, Parsons Laboratory, Preliminary Report, Cambridge,

Connor, J.J. and Wang, J.D. (1973), "Mathematical Models of Massachusetts
Bay," Part I, Massachusetts Institute of Technology, Report No. MITSG 74-4,
Cambridge, Massachusetts.

Cunge, J.A. (1977), "Difficulties of Open Channel Hydraulic Mathematical
Modeling as Applied to Real Life Situations," pp. 147-154, in Applied
Numerical Modeling, Ed. by C.A. Brebbia, Pentech Press, London, 1977.

Dailey, J.E. and Harleman, D.R.F. (1972), "Numerical Model for the
Prediction of Transient Water Quality in Estuary Networks," Massachusetts
Institute of Technology, Parsons Laboratory, Report No. 158, Cambridge,

Feigner, K.D. and Harris, H.S. (1970), "Documentation Report FWQA Dynamic
Estuary Model," U.S. Department of the Interior, FWQA (now EPA),
Washington, D.C.

Graham, D.S. and Christensen, B.A. (1978), "Development of Drainage Basin/
Estuary Numerical Model System for Planning in the Coastal Zone," Proceedings
of COASTAL ZONE 1978, American Society of Civil Engineers, San Francisco,

Graham, D.S., Daniels, J.P. and Christensen, B.A., (1979), "Predicting
Circulation Changes from Bathymetric Modification," Proceedings of the
Specialty Conference on CIVIL ENGINEERING IN THE OCEANS IV, American Society
of Civil Engineers, San Francisco, California.

Graham, D.S., DeCosta, K. and Christensen, B.A. (1978), "Storm Water Runoff
in the Apalachicola Estuary," Hydraulic Laboratory Report No. HY7802,
University of Florida, Gainesville, Florida.

Graham, D.S., Hill, J.M. and Christensen, B.A. (1978), "Verification of an
Estuarine Model for Apalachicola Bay, Florida," Proceedings of the 26th
Annual Hydraulics Division Specialty Conference, American Society of Civil
Engineers, College Park, Maryland.

Hill, J.M. and Graham, D.S. (1980), "Using Enhanced LANDSAT Images for
Calibrating Real-Time Estuarine Water Quality Models," Water Quality
Bulletin, Satellite Hydrology, Vol. 5, No. 1, Canada.

Huber, W.C., Heany, J.P. et al. (1975), "Storm Water Management Model
User's Manual, Version II," Department of Environmental Engineering Sciences,
University of Florida, Gainesville, Florida.

Hydroscience, Inc. (1977), "The Effects of Forest Management on the Water
Quality and Aquatic Biota of Apalachicola Bay, Florida," Report Prepared for
Buckeye Cellulose Corp., Perry, Florida.

Leimkuhler, W., Connor, J., Wang, J., Christodoulou, G. and Sungren, S.
(1975), "Two-Dimensional Finite Element Dispersion Model," Symposium on
Modeling Techniques, Vol. II, ASCE.

Leendertse, J.J. (1967), "Aspects of a Computational Model for Well-Mixed
Estuaries and Coastal Seas," RM 5294-PR, the Rand Corporation, Santa
Monica, California.

Livingston, R.L. (1981), Final Report to Florida Sea Grant, Project No.
R/EM-11, Tallahassee, Florida.

Masch, F.D. et al. (1969), "A Numerical Model for the Simulation of Tidal
Hydrodynamics in Shallow Irregular Estuaries," Technical Report HYD 12-6901,
Hydraulic Engineering Laboratories, University of Texas, Austin, Texas.

Morris, F.W., Walton, R. and Christensen, B.A. (1978), "Hydrodynamic Factors
Involved in Finger Canal and Borrow Lake Flushing in Florida's Coastal Zone,"
Vol. I, UFSG Project No. R/OE-4, Hydraulic Laboratory Report No. HY7801,
University of Florida, Gainesville, Florida.

Pagenkopf, J.R. et al. (1976), "A User's Manual for CAFE-1, A Two-
Dimensional Finite Element Circulation Model, Massachusetts Institute of
Technology, Parsons Laboratory, Report No. 217, Cambridge, Massachusetts.

Parker, B.B. and Pearce B.R. (1975), "The Response of Massachusetts Bay
to Wind Stress," Massachusetts Institute of Technology, Parsons Laboratory,
Technical Report No. 198, Cambridge, Massachusetts.

Pritchard, D.W. (1970, 1973), Coursenotes for (24)626, Lectures on
Estuarine Oceanography, of DWP by D.S. Graham, Johns Hopkins University,
Baltimore, Maryland, Unpublished.

Sengupta, S., Lee, S.S. and Miller, H.P. (1978), "Three-Dimensional
Numerical Investigation of Tide and Wind Induced Transport Processes in
Biscayne Bay," University of Miami, Sea Grant Technical Publication No. 39,
Miami, Florida.

Steele, J.G. et al. (1977), "Finite Element Modeling of Mareton Bay,
Australia,'Massachusetts Institute of Technology, Parsons Laboratory,
Technical Note No. 20, Cambridge, Massachusetts.

Swakan, E.A., Jr. and Wang, J.D. (1977), "Modeling of Tide and Wind Induced
Flow in South Biscayne Bay and Card Sound," Sea Grant Technical Bulletin
No. 36, University of Miami, Miami, Florida.

Swanson, C. and Spaulding, M. (undated), "Generation of Tidal Current and
Height Charts for Narragansett Bay Using a Numerical Model," Department
of Ocean Engineering, University of Rhode Island, Marine Technical Report- --
No. 61, Kingston, Rhode Island.

Swenson, E., Brown, W.S. and Trask, R. (1977), "Great Bay Estuarine Field
Program, 1975 Data Report, Part 1: Currents and Sea Levels," Department of
Earth Sciences, University of New Hampshire, Durham, New Hampshire.

Teal, John and Milred (1969), "Life and Death of the Salt Marsh,"
Ballantine Books, New York.

Wang, J.D. (1978), "Verification of Finite Element Hydrodynamic Model
CAFE," Proceedings, 26th Annual Hydraulics Division Specialty Conference,
University of Maryland, Maryland.

Wang, J.D. and Connor, J.J. (1975), "Mathematical Modeling of Near Coastal
Circulation," Massachusetts Institute of Technology, Parsons Laboratory,
Report No. 200, Cambridge, Massachusetts.

Weisberg, R.H. (1976), "The Nontidal Flow in the Providence River of
Narragansett Bay: Stochastic Approach to Estuarine Circulation," Journal
of Physical Oceanography, 6, (5).









Franklin County, Florida

Prepared Under
Sea Grant

Contract No. R/EM-13


Hydraulic Laboratory
Department of Civil Engineering
University of Florida
Gainesville, Florida

December 1981


This somewhat unconventional atlas is intended to give the reader

detailed information about velocity conditions, pollutant concentrations

and salinity in the Apalachicola Bay System during an average year taking

typical astronomical tides, wind intensities and river flow into consi-


It covers the entire bay from Dog Island and East Pass in the East

to St. Vincent Island and West Pass in the West and is based on the CAFE

and DISPER finite element models applied to a grid system consisting of

489 elements and 281 nodes representing this bay area. The models have

been verified by direct observations of salinities and velocities in the bay

and by pattern recognition in LANDSAT pictures computer enhanced to show

ocean and bay water quality by colors.

The Apalachicola Bay may be classified as a tide dominated well

mixed estuarine system with only sporadic and unsignificantly small areas

of stratification. Consequently, the river flow is of minor importance to

the hydrodynamics of the bay. A hydrograph representing the monthly

discharges averaged from 1971 to 1976 has therefore been used in the model

to represent rates of freshwater flow into the bay.

While the river discharge is of minor importance to the hydrodynamics

of the bay the same can not be said for pollutants brought into the bay

by the river. This influence must be and is modeled in detail and the

results are presented in such a way that the influence of any numerical

value of the river water's pollutant concentration may be evaluated at any

location-in the system. As represented in this atlas the water quality

predictions are limited to conservative pollutants, however, the numerical

model may be extended to consider nonconservative substances transported

by the water.

The atlas shows conditions during a typical tidal cycle representing

each of the twelve months of the year. Each month's events are shown on

seventeen individual sheets. The first eight show the distribution of the

vertically averaged water velocities and their circulation during that

tidal cycle in eight time increments of equal lengths beginning at low tide.

These eight sheets are marked by the name of the month and number 1 through

The ninth sheet, marked by name of the month and 9, represents the

velocities from the first eight sheets averaged over the tidal cycle. In

other words this is the net vertically averaged water velocity.

The remaining eight sheets, marked by the name of the month and

numbers 10 through 17, are reserved for water quality. They show the

distribution of the vertically averaged pollutant concentration c at the

same times during the tidal cycle the vertically averaged horizontal

velocities are given in sheets No. 1 through 8. The shown concentrations c

are based on a concentration cR of the same pollutant in the river discharge

equal to 100 (e.g., ppm) and co = 0 pollutant concentration at all inlets

to the bay from the Mexican Gulf. Simple formulas for the calculation of

c-values corresponding to other cR- and c -values and for determination

of the vertically averaged salinity s from the salinity s0 at the inlets

to the bay are given on the individual water quality sheets.
While this atlas will answer most questions concerning velocities,

their orientation, pollutant concentrations and salinities it is limited

inasmuch as it is prepared for average conditions during the year. Data

corresponding to extreme conditions such as tropical storms or periodic

excessive pollutant loads must be generated separately by use of the detailed com-
puter model. This model is stored on tape and attached to this report.


JAN 1 through JAN

JAN 9:

JAN 10 through JAN 17:

8: Vertically averaged velocities vm
given at time increments equal to
one eighth of the tidal period
beginning at low tide.

Net values of vertically averaged
velocities. Time averaged over one
tidal cycle.

Vertically averaged pollutant
concentration corresponding to
river concentration cR = 100 and
concentration co = 0 at all inlets
to bay. Concentrations given at
time increments equal to one eighth
of the tidal period beginning at low

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs