|
MODELING
THE APALACHliCOLA SYSTEM
.. S 0 SOUTH
S\ CAROLINA
S\ ^Atlanta
ALABAMA GEORGIA .::'
S. *'. Savannah
Chattahoochee Flint Rive .
.River :::::. ATLANTIC
I- *.:-- OCEAN
-.. :. -. ._. . '.
Mobile :. .. Jacksonville
Mobile Ny^^
: *.*.: .
Apalachicola Bay Apalachicola River -.' '-.*
GULF. .
0 F
MEXICO :. .. .
Tampa L*.* ..
St. Pete. -
S":.*:. ;:: Miami
'.. :.. ...:
By Carol Conner
Ann Conway
Barry A. Benedict V*
B. A. Christensen
TECHNICAL PAPER NO. 23 MAY 1982
I
MODELING THE APALACHICOLA SYSTEM*
A Hydrodynamic and Water Quality Model
With a Hydrodynamic and Water
Quality Atlas of Apalachicola Bay
Project No. R/EM-13
By
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.
TABLE OF CONTENT
Page
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. SELECTION OF MODELS
2.1 Criteria for Selection . .... .. . . 22
2.2 Models Considered. . . .... . 23
2.3 Models Selected . . .... . . . 24
3. SPECIFIC FEATURES OF MODELS USED . . .... . 26
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
4.1.1.1 Tide Recordings . . . 50
4.1.1.2 Velocity Measurements . . .... 54
4.1.2 Model Input for Verification Runs . . .. 54
4.1.2.1 Grid Configuration . . .. 54
4.1.2.2 Turbulent Eddy Viscosity Coefficients .58
4.1.2.3 Manning's n . . . . 59
4.1.2.4 River Flow . . . . 59
4.1.2.5 Wind Velocity . . . . 61
4.1.2.6 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
Page
5. THE ATLAS. SELECTION OF CONDITIONS FOR A TYPICAL YEAR .... .76
5.1 River Flows . . . ... .... .. 76
5.2 Tides . . . . .. .. . 76
5.3 Winds . . . . ....... 78
5.4 Generated Results . . . . ... 79
6. DETERMINATION OF POLLUTANT CONCENTRATION (AND SALINITY) AT AN
ARBITRARY LOCATION FOR OTHER RIVER CONCENTRATIONS AND OCEAN
CONCENTRATIONS . . . .. . . 80
7. ACKNOWLEDGEMENTS . . . . ... ...... 83
8. LIST OF SYMBOLS . . . . ... ...... 84
9. REFERENCES . . . . ... .. .. .. 87
APPENDIX A: A HYDRODYNAMIC AND WATER QUALITY ATLAS OF THE APALACHICOLA
BAY SYSTEM, FRANKLIN COUNTY, FLORIDA
APPENDIX B: USER'S MANUAL AND THE COMPLETE VERIFIED MODEL ON COMPUTER
TAPE. (1 reel)
1. INTRODUCTION
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
- DESERT
L'-COASTAL
ESTUARINE
WET AGRICULTURE
-DRY AGRICULTURE
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
resources.
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.
LANDSAT INPUT
LAND USE RIVER BASIN PREDICTIVE MODEL
LANDSAT INPUT 2 ESTUARINE PREDICTIVE MODEL
VERIFICATION
E 3 CORRELATION FROM BIOLOGICAL
DATA COLLECTION
4 ECONOMICAL ANALYSIS
5 PLANNING
6 POLITICAL SYSTEM
Figure 1.2. Logistics of Model Approach.
COMBINED
RIVER/ESTUARY
PREDICTIVE MODEL
(calibrated by existing data)
RIVER BASIN MODEL ESTUARINE
I MODEL
IMPRpVEMENT
BY
RETENTION/DETENTION
CORRELATION
FROM
BIOLOGICAL
DATA COLLECTION
r-
Detailed Flow Chart of the Model Approach to Rational Utilization of Coastal Drainage Basins.
ECONOMICAL
ANALYSIS
r
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
Bay.
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
*Atlar
/ ALABAMA GEC
SOUTH
CAROLINA
S". S ":'. Savannah
Flint Rive .
SATLANTIC
OCEAN
.: Jacksonville
la River :.. .':
GULF
OF
MEXICO
Tampa
St. Pete
[Miami
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 ,\ -
1200
1100
1000
T 77/ ^ ------- --
"E /- QAVE = 705 m's '
800
S700 ------
300
U 600 //
X 500
400 -
Q: !
200
100
100 -O -- -- -- -
JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC
MONTH
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
County.
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.
APALACHICOLA
BAY CA
GULF
OF
MEXICO
I I
-1
10km
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-
rrr;'
~: 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.
Gulf.
VIEW E. Sikes Cut Seen from the Mexican
3r. gIPtc'w4-Ci~ *c 53\.~j-
--- c~~r~--- .- -
3,C~t~t-.-.
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.
~if5~
i"
''
li- ~ L -:'--.
U-
r~rs~***~~r: ';l'-==~,iriiL
-~L ~-iT~
i;
.~L-n
---~K
~C*~Llc
%~f~
Figure 1.19. Apalachicola Bay.
VIEW K. East Bay Seen from Southeast.
...-.2
Figure 1.20. Apalachicola Bay.
from North.
VIEW L. St. George Island Causeway Seen
y---j
-1
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 (
TQ L D
where
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.
DIMENSION
SALT WATER WEDGE
PARTIALLY MIXED
VERTICALLY HOMOGENEOUS
SECTIONALLY HOMOGENEOUS
Figure 1.22.
Pritchard's Estuarine
(1970, 1973).
Classification Scheme.
TYPE
TYPE
TYPE
TYPE
Pritchard
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.
2. SELECTION OF MODELS
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
proficiency
7) model must be able to be run at a reasonable cost
8) model must be flexible, not site-specific, in its
applicability
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
selected.
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
net,
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.
3. SPECIFIC FEATURES OF MODELS USED
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.
26
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
cycles.
The CAFE-1 model, as mentioned, was originally developed at the R.M.
Parsons Laboratory at M.I.T. It has been applied to
439 ELEMENTS
281 NODES
RIVER MOU
APALACHICOLA
EAST PASS
1 DOG
L ISLAND
H-.
In
C
LLJ
C-
LL.J L
116.00 174.00
232.00 290.00 3118.00
METERS EAST(X102 1
Figure 3.1. Finite Element Grid System of Apalachicola Bay Used in this Study.
GULF
SIKES CUT OF
MEXICO
WEST PASS
1z:1
W. 00o
S
56.00
10b6. 00
64. O00
d
.I. I
S5 2.00 CI
1) Massachusetts Bay (Briggs and Madsen, 1973); (Christodoulou
etal.,1974); (Connor etal., 1973); and (Parker and Pearce,
1975).
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
University.
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
dollars.
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)
-h
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
-g
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
coefficients"
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.
Here
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.
Therefore,
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.
36
The wind stress term in CAFE-1 is of the accepted quadratic form:
s C U2 (3.13)
Pair D U10
where
pair = air density
CD = a dimensionless drag coefficient
U10 = air velocity at 10 meters above water surface,
-1
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
grid.
A "rough" formula for estimating Exx is given by Wang and Connor
(1975), as
E ~ ag g ai (3.18)
u
where
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
1977)
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,
1977)
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
-1
u ~ 0.1 m s
Aai ~ 350 m at least
Ei. z (0.02) (34335) = 686.7 m2s1
and
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
and
Eii At 1177 60
J0.0078
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
11
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
Bay
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
items:
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
times.
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
43
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)
where
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
o
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
-0.5
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
CAFE-1.
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.
4. MODEL VERIFICATION
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.
4.1.1.1 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
model.
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 .-.',.'" .
IDIAN -<>AL GULF
S"" :"giN ISLAND 0LF
ST. VINCENT ISLAND WEST E
PASS SKES ::cu PASS
:: : :.. : :
.::.P :-:ST:' ...TR. i.t O +
10. SAND LLAN
ST. VINCENT ISLAND WEST A LA
E3 :: AY
GAPE ST. EORGESLAND
10 5 0 10km
Figure 4.1. Location of Verification Points in Apalachicola Bay.
1
U,
U,
Co
-o
C,
la
-o
L
Co
-or
~I
a)
1..
-4 U,
SiQ
52 rr
I i-
I-
L re -
ri I
: I ~- i
:i.:I/.l-
'HtjL>
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
rr r
r r r
LI
ir I.
i- t L W j hi
rr.
i r
r L r.
r 7
C I C L FV L ..~
t~I r r I
t T- t- I -
t-
I
Figure 4.3. Sample Recording from Tide Gage at East Pass.
53
Ii-
:TE
li-K
I: I
i Ei
*EF *r
-TT
i-"S
41-
4.
r 7 .L..
r II-
1 P-I
il.LI
I~--4~F~ .~-i4
-. L
.r-f +
i~T'FI-- t'-;
L i I
r : '
""
4.1.1.2 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
surface.
4.1.2 Model Input for Verification Runs
4.1.2.1 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.
4.1.2.2 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
xy
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.
4.1.2.3 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.
4.1.2.4 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
1961-1976.
Mean
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)
4.1.2.5 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
4.1.2.6 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.
4.1.2.7 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
percent).
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
E
2
O /
0
S0.2-
-.1
UL
0 T/4 T/2 3 T/4 T
PORTION OF TIDAL CYCLE
;L
crM LW \
I-
-0.1
0 T/4 T/2 3T/4 T
PORTION OF TIDAL CYCLE
Figure 4.7. Water Surface Elevation vs Time at Point B in Apalachicola Bay.
E
U)
z
0 0.2
MO \D-, 7 n.
I.
-" 0.1
-0.1 1 I _____1
SML /--- FIELD DATA
- MODEL, n=0.045
S---- MODEL, n =0.030
0 T/4 T/2 3T/4 T
PORTION OF TIDAL CYCLE
Figure 4.8. Water Surface Elevation vs Time at Point A in Apalachicola Bay.
o FIELD DATA
--- MODEL, n=0.045
300 _
30 MODEL, n=0.030
04-
o //
zp 200 -
WE
U)
o 100 /
o-0 /
LL
0
0 T/4 T/2 3T/4 T
PORTION OF TIDAL CYCLE
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
10000
--- I (DATA)
0 I (MODEL)
30000
50000
(SECONDS)
70000
----- H (OATA)
A H (MODEL)
Figure 4.10. Surface Displacement Verification by Wang for a Model of
Biscayne Bay Using CAFE-1 (Wang, 1978b).
0
---MODEL
---DATA
Figure 4.11. Hodograph Comparison by Wang for a Model of Biscayne Bay
Using Cafe-1 (Wang, 1978b).
N
A
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
value.
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
verification.
In conclusion, it can be stated that an attempt at model verification
with no fine-tuning yielded acceptable results over much of the bay.
-.7
'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
parameter.
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.
75
5. THE ATLAS. SELECTION OF CONDITIONS FOR A TYPICAL YEAR
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 4.1.1.1 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
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.
6. DETERMINATION OF POLLUTANT CONCENTRATION (AND SALINITY) AT AN ARBITRARY
LOCATION FOR OTHER RIVER CONCENTRATIONS AND OCEAN CONCENTRATIONS
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,
POLLUTANT
SITUATION
INDIAN PASS
C=0
WEST PASS
000
-A
SALINITY RIVER
SITUATION Y 8 /
C eBAY P a
INDIAN PASS
s=s
WEST PASS tN
CA
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.
7. ACKNOWLEDGEMENTS
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.
8. LIST OF SYMBOLS
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 =
R=
S
s =
t=
T =
u =
(u,v) =
U10
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
time
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
-I
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)
86
9. REFERENCES
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,
Massachusetts.
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,
Massachusetts.
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,
California.
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).
r
APPENDIX A
A
HYDRODYNAMIC
AND
WATER QUALITY
ATLAS
of
THE APALACHICOLA BAY SYSTEM
Franklin County, Florida
Prepared Under
Sea Grant
Contract No. R/EM-13
by
Hydraulic Laboratory
Department of Civil Engineering
University of Florida
Gainesville, Florida
December 1981
FOREWORD
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-
deration.
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
8.
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.
JANUARY
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
tide.
|
