EPA/600/R-WG49
June
Thn:e Dimensional Hydrodynamic Modd
for Stratified Flows in Lakes and Eslitaries (HYDRO3D):
Theory. User Guidance, and Applications for
Superfund and Ecological Risk Assessments
by
Y. Peter Sheng, Ph,D,'? Mansour Zakikhani, Ph.D.2. Steven C, McCutcheon, Ph.D., P.E.'
with contributions by
7,. TFosseiniptmr, Ph.D,2 and Pei-Fang Wangf Ph.D,1
Donald Eliason, Ph.D,!
Douglas S- Hexui4 and Stephen F, Paiker, Ph.D.4
Earl Hayter, Ph.D.' and Phyllis KohIJ
'University of Florida
Gainesville, FJIorida 32611
3AScI Corporation
Athens, Georgia 30613
1Ecosystems Research Division
Athens, Georgia 30605
4Aeronautica! Research Associates of Princeton, Inc,
Princeton, New Jersey 08540
5Clemson University
wf
j South Carolina 29634
ECOSYSTEMS RESEARCH DIVISION
NATIONAL EXPOSURE RESEARCH LABORATORY
OFFICE OE RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ATPIENS, GEORGIA 30605
-------
DISCLAIMER.
The information in this document has been funded wholly or in part
by the U.S. Environmental Protection Agency under Cooperative Agreement
Number CR-814345-01-0 with the University of Florida. It has been subject
to the Agency's peer and administrative review, and it has been approved
for publication as an EPA document. Mention of trade names or commercial
products does not constitute endorsement or recommendation for use by the
U.S. Environmental Protection Agency.
ii
-------
FOREWORD
As environmental controls become more costly to implement and the penalties of
judgment errors become more severe, environmental quality management requires more efficient
analytical tools based on greater knowledge of the phenomena to be managed. As part of this
Division's research on the occurrence, movement, transformation, impact, and control of
environmental contaminants, the Processes and Modeling Branch develops management or
engineering tools to help pollution control officials address environmental problems.
In assessing ecological risk, models are needed to simulate the effects of complex
reversing flows in lakes, harbors, coastal areas, and estuaries and to determine where chemicals
are transported to in surface waters and where contaminated sediments accumulate. HYDR03D
is a dynamic modeling system that can be used to simulate currents in water bodies as they
respond to tides, winds, density gradients, river flows, and basin geometry and bathymetry.
Rosemarie C. Russo, Ph.D.
Director
Ecosystems Research Division
Athens, Georgia
in
-------
ABSTRACT
Increasing demands for maintaining the quality of stratified surface
waters at reasonable levels have required the development of three-dimensional
hydrodynamic models. To meet these needs, the HYDR03D program has been
documented to aid in the simulation of lakes, harbors, coastal areas, and
estuaries.
HYDR03D is a dynamic modeling system that can be used to simulate
currents in water bodies as they respond to tides, winds, density gradients,
river flows, and basin geometry and bathymetry. The code is a three-
dimensional, time-dependent, a-stretched coordinate, free surface model that
can be run in fully three-dimensional (3-D) mode, two-dimensional vertically-
averaged (x-y), and two-dimensional laterally-averaged x-z mode.
The prognostic variables are the three components (x-y-z) of the
velocity field, temperature, and salinity. The governing equations together
with their initial and boundary conditions are solved by finite difference
techniques. A horizontally and vertically staggered lattice of grid points is
used for computation. The code solves for steady-state or the time-dependent
water surface displacement, vertically-integrated velocities, 3-D velocities,
temperature, salinity, and dissolved species concentrations. The vertical
turbulence parameterization schemes include constant eddy viscosity, variable
eddy viscosities (Munk-Anderson type), and a simplified version of a second-
order closure model.
The applications provided here demonstrated that the model is capable of
realistic simulation of flow and salinity transport in complex and dynamic
water bodies. These applications include simulations of tidal circulation and
salinity transport in Suisun Bay, California and, Charlotte Harbor, Florida
and; wind-forced circulation in Green Bay, Lake Michigan. Tidal circulation
in Prince William Sound, Alaska was investigated to determine the feasibility
of applying the model under emergency conditions. Finally, the calibration of
the model for the Mississippi Sound is illustrated.
HYDR03D is a far-field model that like any other computer code, has
limitations. The present version does not contain a flooding and drying
scheme. Near field effects of cooling water discharges, diffusers, other
jets, and reservoir withdrawal cannot be adequately simulated. In addition,
short-period waves are not included in the model.
The information provided in this manual, along with the complete program
listing which will be provided separately, should be sufficient for the user
to operate the code. However, a successful model simulation of HYDR03D
requires sufficient data and familiarity with the code. The documentation
iv
-------
provides a brief review of the theory and structure of the program. Data
requirements are noted and example applications demonstrates uses of the
program.
-------
Page Intentionally Blank
-------
CONTENTS
DISCLAIMER ii
FOREWORD iii
ABSTRACT iv
LIST OF FIGURES ix
ACKNOWLEDGEMENTS .... ...... xiv
PREFACE xviii
1.0 INTRODUCTION 1
1.1 PROGRAMATIC NEEDS FOR MODELS TO SIMULATE STRATIFIED FLOWS . . 1
1.2 USE OF HYDRODYNAMIC MODELS BY TECHNICAL EXPERTS AND MANAGERS
AND HOW THIS MANUAL MAY BE OF ASSISTANCE 4
1.3 SCREENING LEVEL SIMULATIONS 6
1.4 CALIBRATION AND VALIDATION 7
1.5 DATA REQUIREMENTS 8
2.0 MODELING SYSTEM 9
2.1 OVERVIEW OF THE MODELING SYSTEM . 12
2,2 MODEL FORMULATION 13
2.2.1 Governing Equations . 13
2.2.2 Grid System 17
2.2.3 NQn-Dimensionalization of the Governing Equations . . 19
2.2.4 DimensionlessEquationsin Stretched Coordinates ... 23
2.2.5 Vertically Integrated Equations 25
2.2.6 VerticalVelocities 26
2.3 BOUNDARY AND INITIAL CONDITIONS 26
2.3.1 Vertical Boundary Conditions . . 26
2.3.2 Lateral Boundary Conditions 27
2.3.3 InitialConditions ..... 29
2.4 NUMERICAL SOLUTION ALGORITHM ....... 29
2.4.1 External Mode 29
2.4.2 Internal Mode 32
2.5 TURBULENCE CLOSURE 34
2.5.1 Constant Eddy Coefficients 34
2.5.2 Hunk-Anderson Type Eddy Coefficients 35
2.5.3 A Simplified Second-Order Closure Model 37
2.6 GRID LAYOUT 40
2.6.1 Staggered Grid 40
2.6.2 Grid Index 42
2.7 FLOW CHARTS 42
vi
-------
CONTENTS - continued
3,0 USER'S MANUAL 49
3.1 INTRODUCTION 49
3.2 DATA REQUIREMENTS OF THE PROGRAM 49
3.3 INPUT DATA DESCRIPTION 52
3.4 MODEL OUTPUT 65
4.1 COMPARISON WITH A ONE-DIMENSIONAL ANALYTICAL SOLUTION .... 69
4.2 SUISUN BAY. CALIFORNIA . 72
4.2.1 Physical Setting 72
4.2,2 CirculationPatterns 73
4.2,3 Modeling 3-D Circulation in Suisun Bay . 77
4.2.4 Results 81
4.3 CHARLOTTE HARBOR^ FLORIDA 90
4.3.1 Physical Setting ..... 90
4.3.2 Circulation in CharlotteJHarbor . 93
4.3.3 Modeling 3-D Circulationin Charlotte Harbor 93
4.3.4 Results ..... 99
4.4 GREEN BAY. LAKE MICHIGAN 115
4.4.1 Physical Setting 115
4.4.2 Two-Jimensional Simulation of Flow 116
4.4.3 3-D Simulation of Flow 124
4.5 PRINCE WILLIAM SOUND,ALASKA 130
4.5.1 Physical Setting 130
4.5.2. Modeling Parameters ... 134
4.5.3 Results 134
4.5.4. Discussion 140
4.6 CURRENTS IN MISSISSIPPISOUND 146
4.6.1 Physical Setting , 146
4.6.2 Circulation in Mississippi Sound ...... 146
4.6.3 Results 148
4.6.3.1 Tidal Simulation 148
4.6.3.2 Wind-effect on Tidal-Driven Currents ..... 153
5.0 HYDR03D PROGRAMMER'S GUIDE 148
5.1 OVERVIEW 156
5.2 HARDWARE AND SOFTWARE REQUIREMENTS 156
5.3 INSTALLATIONANDIMPLEMENTATION 156
5.4 DESCRIPTION OF THE COMPUTER PROGRAM . 156
5.5 SUBROUTINE DESCRIPTIONS 159
5.6 INPUT/OUTPUT UNITS 163
6.0 CONCLUSIONS AND RECOMMENDATIONS 165
REFERENCES 167
APPENDIX A 174
APPENDIX B 175
VII
-------
CONTENTS - continued
APPENDIX C 177
APPENDIX D ........ . 178
viii
-------
Page Intentionally Blank
-------
LIST OF FIGURES
Figure Page
1 Cartesian coordinates at the nominal water surface 16
2 Vertical stretching of the coordinates 18
3 Lateral stretching of the coordinates, ... 20
4 (a) Empirical stability functions of vertical
turbulent eddy coefficients
(b) Stability functions determined from a second
order closure model of turbulent transport . 39
5 Staggered numerical grid 41
6 Grid indices NS and MS 43
7 Flow chart of the main program EHSMML 44
8 Flow chart of the hydrodynamic subroutine EHSMHC 45
9 Flow chart of the external mode subroutine EHSMEX 46
10 Flow chart of the internal mode subroutine EHSMB3 47
11 Flow chart of the internal mode subroutine EHSMB4. 48
12 Water surface elevation and current velocity at x=5 km
(solid lines; analytical solution; dashed lines;
numerical solutions) ...... 70
13 Water surface elevation and current velocity at x=25 km
(solid lines; analytical solutions; dashed lines;
numerical solutions) .... 71
14 Map of San Francisco Bay estuarine system 74
15 Map of Suisun Bay region and the location of current-meter
moorings, tide stations, and a USGS weather station 75
16 Three-dimensional plot of the Suisun Bay bathymetry when
viewed from (a) the southwest and (b) the southeast 76
ix
-------
17 Time histories of surface elevation and depth-averaged
velocity components at (a) C26, (b) C27, (c) C28, (d) C30,
(e) C239 and (f) the left boundary during a 5-day model
simulation of tidal circulation in Suisun Bay 78-80
18 Time histories of salinity at 3 vertical levels near-bottom,
mid-depth and near-surface at C26, C27, C28, C30, C239
and the left boundary during a 5-day model simulation of
tidal circulation in Suisun Bay 83
19 Tide- and salinity-driven currents in Suisun Bay near
the bottom (a - -0.9) and near the surface (o = -0.1)
at 96 hours 84
20 Tide- and salinity-driven currents in Suisun Bay near ;
the bottom (cr - -0,1) and near the surface (cr «= -0.1) !i
at 108 hours . 85
21 Tide- and salinity-driven currents in Suisun Bay near
the bottom (a - 0,9) and near the surface (a - -0.1)
at 120 hours 86
22 Salinity distribution in Suisun Bay near the bottom
(u -0,9) and near the surface at 96 hours 87
23 Salinity distribution in Suisun Bay near the bottom
(a « -0.9) and near the surface at 108 hours 88
24 Salinity distribution in Suisun Bay near the bottom
(a = -0.9) and near the surface at 120 hours 89
25 Map of Charlotte Harbor Estuarine System ,91
26 Map of Northern Charlotte Harbor with locations of
water quality/current meter stations during the June
and July 1982 study 92
27 Tidal stage at Burnt Store Marine during July 20 to
July 22, 1982 , , 94
28 Discharge of Peace River during June 1 to July 30, 1982 .... 95
29 Initial 3-D velocity field in Charlotte Harbor for model
simulation to June 25 to June 28, 1982. ............ 97
30 Initial salinity field in Charlotte Harbor for June 25, 1982. . 98
31 Time histories of water level, surface currents, bottom
currents, surface salinity and bottom salinity at Station
10 during the 3-day model simulation period . 100
-------
32 Time histories of water level, surface currents, bottom
currents, surface salinity and bottom salinity at Station
7 during the 3-day model simulation period 101
33 Time histories of water level, surface currents, bottom
currents, surface salinity and bottom salinity at Station
22 during the 3-day model simulation period 102
34 Time histories of water level, surface currents, bottom
currents, surface salinity and bottom salinity at Station
19 during the 3-day model simulation period 103
35 Computed 3-D velocity field in Charlotte Harbor after
24 hours of model simulation 104
36 Computed salinity field in Charlotte Harbor after 24
hours of model simulation 105
37 Computed 3-D velocity field in Charlotte Harbor after
72 hours of model simulation. 106
38 Computed Salinity field in Charlotte Harbor after 72
hours of model simulation 107
39 Vertical salinity profiles at Stations 7, 22, 15 and 19
in Charlotte Harbor after 72 hours of model simulation. .... 108
40 Computed 3-D velocity field in Charlotte Harbor after 48
hours of model simulation with 1/2-km grid 109
41 Computed salinity field in Charlotte Harbor after 48 hours
of model simulation with 1/2-km grid 110
42 Time histories of water level, surface currents, bottom
currents, surface salinity and bottom salinity at Station
7 during the 2-day model simulation period Ill
43 Time histories of water level, surface currents, bottom
currents, surface salinity and bottom salinity at Station
22 during the 2-day period 112
44 Time histories of water level, surface currents, bottom
currents, surface salinity and bottom salinity at Station
19 during the 2-day model simulation period 113
45 Time histories of water level, surface currents, bottom
currents, surface salinity and bottom salinity at Station
10 during the 2-day model simulation period . . 114
46 Map of Green Bay showing relation to Lake Michigan and other
Great Lakes 117
xi
-------
47 Three-dimensional plot of Green Bay bathymetry 118
48 Grid network of Green Bay 119
49 Measured water surface elevation at the mouth of Green Bay
and at Green Bay city during September, 1969 120
50 Measured water surface elevation at the mouth of Green Bay
and at Green Bay city during October, 1969 121
51 Simulated circulation in Green Bay using October 1969
data (54 hours from at the at-rest) 122
52 Simulated circulation in Green Bay using October, 1969 dataS
(114 hours from the at-rest) , 123
53 Measured and calculated water surface elevation at Green Bay
mouth near Green Bay city during September 17-20, 1969 125
54 Measured and calculated water surface elevation at Green Bay
mouth near Green Bay city during October 8-12, 1969 . 126
55 Water surface elevation in Green Bay after 40 hours ...... 128
56 3-D vertically averaged currents in Green Bay after 40 hours. . 129
57 3-D simulation of currents in Green Bay (near the surface
layer) .131
58 3-D simulation of currents in Green Bay (near the bottom
layer) 132
59 Map of Prince William Sound, Alaska 133
60 Coarse grid of Prince William Sound, Alaska 135
61 Fine grid of Prince William Sound, Alaska 136
62 2-D vertically averaged of currents in Prince William Sound
using coarse grid (1 hour after simulation) 137
63 2-D vertically averaged of currents in Prince William Sound
using coarse grid (2 hours after simulation) 138
64 2-D vertically averaged of currents in Prince William Sound
using coarse grid (3 hours after simulation) 139
65 3-D simulation of currents in Prince William Sound using
coarse grid (1 hour after simulation) 141
xii
-------
66 3-D simulation of currents in Prince William Sound using
coarse grid (2 hours after simulation) ...... 142
67 3-D simulation of currents in Prince William Sound using
coarse grid (3 hours after simulation) 143
68 2-D vertically averaged of currents in Prince William Sound
using fine grid (2 hours after simulation) 144
69 2-D vertically averaged of currents in Prince William Sound
using fine grid (3 hours after simulation) . 145
70 Lateral Numerical Grid Used for Dynamic Simulation of Coastal
Currents Within the Mississippi Coastal Waters ........ 147
71 Transient Variation of Surface Displacements at Four Stations
Within the Mississippi Sound from 9/20/80 to 9/25/80 149
72 Surface Displacement Contours Within the Mississippi Coastal
Waters at 0 hr., 9/23/80 150
73 Transient Variation of Mid-depth Velocities at Two Stations
Within the Mississippi from 9/20/80 to 9/25/80 151
74 Horizontal Velocity Field at 0 hr and 1 m depth, 9/23/80. ... 152
75 Influence of Wind on Surface Displacement at Two Stations
from 9/20/80 to 9/24/80 154
76 Influence of wind on mid-depth horizontal velocities at
two stations from 9/20/80 to 9/24/80 155
77 Operational chart of the HYDR03D model 158
78 Steady-state wind-driven currents in an enclosed square basin
of 50 km on each side; linearly varying bottom from 3m
(South and North/to 10 m (at center) . 179
79 Three-dimensional simulation of wind-driven currents in an
enclosed basin; results are for three grid points of
(2,6), (6,6) and (10,6) . 180
xiii
-------
Page Intentionally Blank
-------
ACKNOWLEDGEMENTS
The development of a complex hydrodynamics model is an effort that
extends over a number of years and involves the cooperative efforts of many
individuals and groups. However, these models are not developed to a
practical level without the long-term commitment and intensive effort of at
least one person dedicated to the creation of a useful simulation tool. In
the case of the HYDR03D (or EHSM3D) model described in this report, that
individual is Dr. Peter Sheng. Dr. Sheng began the development of:, state-of-
the-art studies at Case Western Reserve University between 1969 and 1975. A
precursor of the HYDR03D model was developed for the simulation of lake,
estuarine and coastal circulation (Sheng et al. 1978, Sheng and Lick, 1980).
During his post-doctoral research Dr. Sheng developed a sediment transport
model (Sheng and Lick 1979, Sheng 1980) which was driven by the precursor of
EHSM3D to provide current information and a wave model to provide wave
information. These earlier works are still widely cited as the appropriate
framework for analyzing hydrodynamic and sediment transport problems in lakes.
During these landmark studies, Dr. Sheng recognized the need to develop better
hydrodynamic models which can resolve turbulent mixing and the effects of
complex geometry and bathymetry in a more rigorous fashion. During his
employment with Aeronautical Research Associates of Princeton, Dr. Sheng was
able to significantly improve the hydrodynamic and sediment transport models
through several research projects.
Among the important projects were those with the U.S. Army Engineer
Coastal Engineering Research Center (CERC) and the Hydraulics Laboratory (HL)
at the Waterways Experiment Station (WES). Dr. Sheng developed an updated
version of the original hydrodynamic and sediment transport models (Sheng,
1983; Sheng and Butler, 1982), with the cooperation of Lee Butler and others
at WES. This version was used in several studies in coastal and estuarine
waters, e.g., the Mississippi Sound and the adjacent shelf in Gulf of Mexico,
Humboldt Bay, and Los Angeles-Long Beach Harbor. The Mississippi Sound study
(Sheng, 1983) was particularly interesting because it included the modeling of
tidal and wind-driven circulation, modeling of waves, modeling of sediment
transport, and field and laboratory sediment measurements. The model
resulting from this study was notable for the use of a simplified second order
turbulence closure technique and a variable grid. Cooperative efforts between
Peter Sheng and the engineers at WES eventually produced a curvilinear grid
model (Sheng 1986a; Sheng 1987) that is similar to HYDR03D, Since 1986, Peter
Sheng and his group at the University of Florida continued to further enhance
the hydrodynamic and sediment transport models through studies on the James
River (Sheng et al. 1989a), Chesapeake Bay (Sheng et al. 1989b) and Lake
Okeechobee (Sheng et al. 1989c). The work with similar models and enhanced
versions have contributed to the improvement of the HYDRO3D code. At the same
time, the earlier modeling studies done in cooperation with WES generated
xiv
-------
several in-house projects through which the original hydrodynamics model was
further improved. Particularly noteworthy was the Chesapeake Bay project. As
a result of the Chesapeake Bay Study, a new version of the basic (Sheng 1983)
model was developed. The new version has improved the capability to maintain
channel stratification during year-long simulations and uses a z-grid or
cartesian coordinate in the vertical dimension. Additional improvements
include a new representation of the momentum-convection terms and the addition
of a spatial third-order scheme called QUICKEST, in the transport equations
for salt and heat (Johnson et al. 1989). This work has been performed
primarily by Dr. Billy Johnson of the WES Hydraulics Laboratory, and Dr. Kim
and Mr. David Marks with the Coastal Engineering Research Center. Mr. Mark
Dortch and others with the Environmental Laboratory at WES have made
significant contributions in developing linkage to larger scale water quality
models. These efforts are valuable in ensuring that this hydrodynamics model
and others are fully useful in routine water quality simulation. As noted,
the recent advances by the WES have resulted in different versions of the
model from that being documented here, but these studies have provided a
number of improvements, some of which have been incorporated into the version
of the code being documented herein. As a result, we find it very appropriate
to acknowledge the contribution of the WES to the development of this type of
hydrodynamics model for stratified flows.
The support of the U.S. Geological Survey has also been important in the
development and documentation of the EHSM3D model. Dr. Ralph Cheng of the
Water Resources Division National Research Program was the project officer of
the effort that among other things, funded the development of the salt balance
equations and testing of the model. Mr. Carl Goodwin of the Florida District
Office provided data to test the model. A draft of this manual (Sheng et al.
1986) was submitted for those studies and that manual was used as a guide for
this document. We therefore wish to acknowledge the original contributions of
Mr. Douglas Henn and Dr. Stephen Parker to elements of this document and we
have added them contributors to this document on that basis. The original
document (Sheng et al. 1986) was edited by Mr, Peter Smith of the U.S.
Geological Survey California District and Mr. Smith has made a significant
contribution to the testing and correction of an early version of the EHSM3D
model. Mr. Smith and Dr. Cheng have also recently published a satisfactory
validation of the EHSM3D model in San Pablo Bay of the San Francisco Bay
system (Smith and Cheng 1989). In addition, we wish to acknowledge the recent
review of this report by Mr. Smith. He noted, in his formal review comments
and has relayed in previous conversations with Dr. Steve McCutcheon and Ms.
Sandra Bird of the U.S. EPA Environmental Research Laboratory at Athens,
Georgia, several areas where the users manual was inadequate for general use
and several areas where the documentation does not seem to fully correspond
with or explain the equations in the code. Mr. Smith has also noted several
areas where the model could be improved, based on his tests for adapting the
code for different purposes and including getting the code to run faster on a
supercomputer. The contribution of Mr. Smith and the Geological Survey is
appreciated.
At the U.S. EPA Environmental Research Laboratory in Athens (ERL-
Athens), Georgia, Dr. Steve McCutcheon has been responsible for the
development of hydrodynamics and sediment transport models and headed the
xv
-------
competitive selection process that selected Dr. Peter Sheng and the University
of Florida to develop a more comprehensive sediment transport model. As part
of this effort, a three dimensional hydrodynamics model has been delivered by
Dr, Sheng, This model (HYDR03D or EHSM3D) provides information about fluid
velocities and shear stress for the sediment transport algorithms under
development. The hydrodynamic model has been tested in a study of Green Bay
by Ms. Sandra Bird (ERL-Athens, formerly with AScI Corp.), Dr. Mansour
Zakikhani (AScI Corp.), Dr. Pei-Fang Wang (AScI Corp.), Dr. Steve McCutcheon
(ERL-Athens), Dr. James Martin (AScI Corp.), and Dr. Edward Hosseinipour (AScI
Corp.). Ms. Bird did much of the initial testing and implementation of the
model at ERL-Athens and she developed the original linkage to larger scale
models being used in the study. Dr. Zakikhani has continued the testing and
calibration of the model and made a number of changes to achieve some
improvements in it's use (Zakikhani et al. 1989), He has developed and
written the user guidance contained in this manual with other authors listed.
More recently, Dr. Wang has continued the testing of the model arid has
developed the Green Bay Case study originally begun by Ms. Bird and Dr.
Zakikhani. Dr. Wang also developed the comparison with the analytical
solution and resoulved a number of questions about how to described the model
and input data. Drs. Wang and Zakikhani summarized the case study on the
Mississippi Sound. Dr. Hosseinipour has coordinated the efforts producing the
manual and the Green Bay Study used as an example.
The implementation and testing of the model in Green Bay has been
supported by Mr. William Richardson of the Large Lakes Research Station
attached to the U.S. EPA Environmental Research Laboratory at Duluth,
Minnesota. Data collection to complete the testing in the future has been
undertaken by the U.S. EPA Great Lakes National Program Office in cooperation
with NOAA, Sea Grant, State of Wisconsin and several other institutions.
Developed of the sediment transport and hydrodynamics models at the University
of Florida has been supported by the U.S. EPA Ecological Risk Assessment
Research Program coordinated at ERL-Athens by Dr. Harvey Holm (now with the
U.S. EPA Newport Field Station) and Dr. Craig Barber.
The development of this documentation has also been supported by the
Center for Exposure Assessment Modeling (CEAM) at ERL-Athens managed by Mr.
Robert B. Ambrose, Jr. The CEAM supports general purpose models that can be
used for Superfund site investigations, including those where complex
hydrodynamic circulation has an effect on the spread of contamination in
harbors, lakes and estuaries. Such sites include those located in Eagle Rock
Harbor of Puget Sound, New Bedford Harbor in Massachusetts, and Sheboygan
Harbor in the Great Lakes.
Support of the Civil Engineering Department of Clemson University for
Ms. Phyllis Kohl and Dr. Earl Hayter is also gratefully acknowledged. This
hydrodynamics model was selected for preliminary testing in Prince William
Sound to aid in tracking the long-term effects of the March 24, 1989 oil spill
from the EXXON Valdez. Dr. Hayter and Ms. Kohl tested the feasibility of
using the model to simulate circulation in the Sound and continued the study
as a Masters Thesis at Clemson. The continued testing and preliminary report
on the results contained herein is appreciated.
xvi
-------
Ms. Donna Hinson (AScI Corp.), Ms. Tawnya Robinson (AScI Corp.), Mr. Tim
Register (ERL-Athens), and other support staff at ERL-Athens have assisted in
in retyping, proofreading, scanning the original drafts, retyping equations,
and redoing graphics. These efforts are very much appreciated. We especially
appreciated the efforts of Ms. Robinson and Ms. Hinson in resolving last
minute software glitches in getting the document in the final format. Mr.
Robert Ryans (ERL-Athens) edited the document and assisted with revisions of
the format. As always, we appreciate his diligence and dedication to
improving the editorial quality of these manuals.
This document has been significantly improved by the reviews of a number
of colleagues. These include Dr. Raymond Walton (EBASCO Environmental,
Seattle), Dr. Alan Blumberg (HydroQual, New Jersey), Dr. Tien-Shuen Wu
(Northwest Florida Management District), Dr. Michael Amein (North Carolina
State University, Raleigh), Dr. Frederick Morris (Saint John's River Water
Management District, Palatka, Florida), and Mr. Peter Smith (U.S. Geological
Survey). The reviews of Dr. Walton, Dr. Morris, and Mr. Smith were especially
helpful in improving the usefulness of this documentation. We have mentioned
Mr, Smith's assistance above. We appreciate the various reminders from Drs.
Walton and Morris and the other reviewers about the elements that make a
documentation report useful.
Dr. McCutcheon wrote the Preface, Acknowledgements, and Introduction of
this report. He edited the report for consistency and supplemented various
sections to fully describe the work. Dr. Sheng, Mr. Henn, and Dr. Parker
wrote the original draft of Section 2. Dr. Sheng, Dr. Eliason, Dr. Zakikhani,
Dr. Hosseinipour, and Dr. McCutcheon enhanced and revised the section. Dr.
Zakikhani and Dr. Sheng wrote the User's Manual, Section 3 using Sheng et al.
(1986) as a guide. Dr. Wang performed the calculations and wrote about the
comparison of the model simulations with an analytical solution suggested by
Dr. Hosseinipour. Dr. Sheng performed the calculations and wrote about the
studies for Suisun Bay and Charlotte Harbor and Dr.s Zakikhani, McCutcheon,
and Hosseinipour revised the sections to answer reviewer comments. Dr. Wang
and Dr. Zakikhani performed the calculations and wrote about the studies for
Green Bay using simulations initially began by Ms. Bird, Dr. Hayter and Ms.
Kohl performed the calculations and wrote about the feasibility studies for
Prince William Sound. Dr. Zakikhani revised the section to address review
comments. Dr. Wang and Dr. Zakikhani excerpted a summary of the study of the
Mississippi Sound from Sheng (1983). Dr. Zakikhani and Dr. Hosseinipour wrote
the Programmers Guide for material available from the original draft by Sheng
et al. (1986). Dr. Zakikhani and Dr. McCutcheon wrote the Conclusions and
Recommendations, and the Abstract,
Dr. McCutcheon, Dr. Hosseinipour, and Dr. Zakikhani are responsible for
the editorial and technical accuracy of this document. Dr. Sheng is
responsible for the technical accuracy of the theoretical basis of the model
as described in Section 2, for the calculations and conclusions described in
the sections on Suisun Bay and Charlotte Harbor, and for the original code,
EHSM3D. Dr. McCutcheon, Dr. Wang, Dr. Zakikhani, and Dr. Hosseinipour are
responsible for changes in the code as part of the Green Bay study and for the
slightly revised code now referred to as HYDR03D.
Steve C. McCutcheon (Athens, Georgia)
xvii
-------
PREFACE
The HYDR03D computer program is one of several codes under development
at the U.S. EPA Environmental Research Laboratory in Athens, Georgia (ERL-
Athens). The development of hydrodynamic and sediment transport codes at ERL-
Athens is proceeding as follows:
Name
HYDRD3D
Dimen-
sion-
ality Type
3D Dynamic circulation model
Status
Documented in this report.
Existing and Anticipated
Level 'cf Support
as of April 1990
Level II - Code,
HYDR02D-V
20
HYDR01D-DYNHYD ID
for far-field transport
in lakes, estuaries, and
coastal areas. Employs
approximate second-order
closure scheme.
Vertically averaged finite
element hydrodynamic model
coupled with a cohesive
sediment transport code
described as SED2D-V
below.
Branched version of
Dynamic Estuary Model
involving Manning
roughness coefficient
and wind stress.
The code is expected to be
ready for release by July 1,
1990. A beta test version
is ready now for preliminary
implementation at Superfund
sites and other critical
study areas. Updates to
the code are anticipated in
September 1990 when the
final hydrodynamic and sedi-
ment transport model is
delivered by the University
of Florida.
Documentation has not been
published and is not readily
available except in draft
for beta test users. CEAM
may be able to assist select
EPA projects, especially
those involving Superfund
sites. Documentation and
code will be available by
sunnier 1990.
Documented as part of the
WASP4 code. Fully
operational and applied in a
number of studies, but the
basic equations have some
limitations that must be
understood. Case studies
include use in moderately
dynamic flows in estuaries
and rivers.
documentation, and start-up
instructions available
from CEAM. Implementation,
debugging, and
interpretation assistance
not fully available,
except on a limited basis.
Level II (anticipated) -
Code, documentation, and
start-up instructions will
be available from CEAM.
Implementation, debugging,
error correction, and
interpretation assistance
is available from Dr. Earl
Hayter, Clemson University
on a negotiated basis.
Level I - Code,
documentation, and start-up
instructions available
from CEAM. Implementation,
debugging, error
correction, and
interpretation assistance
fully available for most
studies.
xviii
-------
List of Hydrodynamic and Sediment Transport Models Available At ERL-Athens, continued
Name
Dimen-
sion-
ality
Type
Status
Existing and Anticipated
Level of Support
as of April 1990
HYDR01D-RIVMOD 10
Dynamic routing model for
single channels.
Preliminary documentation
will be available July 1990.
The code has not been for-
mally released but is avail-
able for use in select CEAM
projects and by beta testing.
Limited support is avail-
able in the early stages
of development. Users
with limited experience
are referred to models of
the Corps of Engineers and
U.S. Geological Survey if
HYDR01D-DYNHYD and kinema-
tic wave routing in HSPF
are not adequate for
the problem to be solved.
Level II support is
anticipated after the
summer of 1990.
HSPF
(kinematic
wave routing)
SED3D
SED2D-V
1D Kinematic nave
(stage-discharge) and
simple sediment routing
for dendritic branched
channels.
3D Dynamic sediment
dispersion, resuspension,
and deposition model based
on the most recent under-
standing of the important
processes. Model inte-
grated into an updated
version of the HYDR03D
code described above.
2D Finite element cohesive
sediment transport for
vertically averaged
estuaries, rivers and
other unstratified water
bodies. Linked with
HYDR02D-V to calculate
average shear stress
levels.
Documentation and code are
fully available and opera-
tional. Limitations are
well described in
documentation.
Level I - Code, and
start-up instructions
available from CEAM.
Documentation available
from NTIS. Implementation
debugging, error
correction, and
interpretation assistance
fully available for most
studies from CEAM and U.S.
Geological Survey.
Documentation and code
expected in FY91. Beta test
versions may be available by
late summer 1990 for selected
projects, especially those
involving Superfund sites.
Documentation has not been
published and is not readily
available except in draft
for beta test users. CEAM
may be able to assist select
EPA projects, especially
those involving Superfund
sites. Documentation and
code will be available by
summer 1990.
Not available but Level
support is anticipated.
II
Level II (anticipated) -
Code, documentation, and
start-up instructions will
be available from CEAM.
Implementation, debugging,
error correction, and
interpretation assistance
is available from Dr. Earl
Hayter, Clemson University
on a negotiated basis.
XIX
-------
UASP4
1D, Simple sediment mass bal-
2D, ance with advection, depo-
& sition, and resuspension
3D velocities, end eddy dif-
fusivity mixing. Precise
sediment transport
calculations require input
from other algorithms or
codes.
Documented, fully
operational and applied in a
number of studies of lakes,
estuaries, and rivers.
Level I - Code,
documentation, and start-up
instructions available
from CEAH, Implementation,
debugging, error
correction, and
interpretation assistance
fully available for most
studies.
Note: CEAM is the Center for Exposure Assessment Modeling located at U.S. EPA Environmental Research
Laboratory, College Station Road, Athens, GA 30605, (404) 546-3130, Bulletin Board Phone; (404) 546-3402.
The HYDR03D code is essentially the same code as EHSM3D (Estuarine
Hydrodynamic Software Model) developed by Peter Sheng in conjunction with the
U.S. Army Engineer Waterways Experiment Station (Sheng 1983) and the U.S.
Geological Survey (Sheng et al. 1986). The recent applications have
concentrated on investigations in estuaries as the title EHSM3D indicates.
However, the HYDR03D code is a general purpose computer program designed to
simulate complex dynamic currents in lakes, estuaries, harbors, and coastal
waters. The original code was developed in the Canadian-American Great Lakes
(Sheng et al. 1978, Sheng and Lick 1980). Prior to documenting the code with
this report, ERL-Athens investigated the feasibility of using the code in lake
settings (Zakikhani et al. 1989) and made a few minor changes to improve the
usefulness of the program. However, these recent changes by ERL-Athens are
not significant enough to warrant changing the model name except that the name
EHSM3D is misleading regarding the applicability to lakes and other waters.
-------
SECTION 1
INTRODUCTION
1-1 PROGRAMATIC NEEDS FOR MODELS TO SIMULATE STRATIFIED FLOWS
There is increasing emphasis being placed on the simulation of
stratified flows in lakes arid estuaries by a number U.S. Environmental
Protection Agency (EPA) Programs. The emphasis arises from the need to
prevent and mitigate pollution in lakes, estuaries, and other stratified
surface waters. For example, the Superfund Program is beginning to
investigate more sites where human and ecological health is affected by
contaminant transport in stratified surface water flows. In addition, the
Ecological Risk Assessment Research Program of the EPA Office of Research and
Development (ORD) is developing, in conjunction with the Office of Pesticides
and Toxic Substances, exposure assessment tools to determine the exposure of
biota in lakes and estuaries. These exposure assessments involve determining
how chemical concentrations are controlled by flows that are normally
stratified. Also, stratified flows and other complex flows control the
transport of chemicals attached to sediments. Other EPA programs that will
require some understanding of the hydrodynamic transport of sediment include,
the EPA Great Lakes National Program Office ARCS (Assessment and Remediation
of Contaminated Sediments) Program for the cleanup of contaminated sediments
from the Great Lakes toxic hotspots; the EPA Office of Water Programs and the
ORD Sediment Quality Initiative aimed at developing waste load allocation
methods for sedimentary contaminants; the ORD initiatives to investigate
eutrophication and toxic chemical fate in large lakes and marine waters; the
ORD Global Climate Program on the affects on contaminants and biogeochemical
cycles in stratified coastal waters, estuaries, and lakes; the ORD Oil Spill
Response Initiative; the development of response and cleanup plans for the ORD
Alternative Fuels Initiative; the EPA National Estuary Studies guided by the
Office of Marine and Estuarine Protection; the EPA Region IV Gulf of Mexico
Initiative, and the EPA Office of Radiation Safety programs to determine the
fate of marine sediments contaminated with radioactive elements. To support
these programs, hydrodynamics models are required to simulate the effects of
complex reversing flows at harbor entrances into lakes and estuaries, to
determine where chemicals are transported to, and where contaminated sediments
accumulate. Other programs involving sediment transport will also require
methods to determine the hydrodynamic effects on the transport and dispersion
of fine sediments enriched with nutrients, metals, radioactive elements, and
pesticides and other toxic organic chemicals. The resuspension and deposition
of fine sediment is best simulated using a hydrodynamic model to map out
levels of fluid shear stress at the bottom and throughout the water column.
-------
A list of the Programs and areas of research that this work will support
include:
Superfund site assessment and remediation,
Determining transport of nutrients and contaminants for bioremediation of
hazardous waste sites,
Ecological and human health risk assessments in the Great Lakes and other
critical stratified surface waters,
Contaminated sediment resuspension, deposition, and transport studies,
Waste load allocation for sedimentary contaminants (EPA Office of Water
Sediment Quality Criteria Programs and ORD Sediment Quality Initiative),
Waste load allocation for conventional and toxic substances (pesticides,
organic chemicals, and metals) in estuaries (see Ambrose et al. 1990) and
lakes,
Effects of global climate change on circulation in lakes, estuaries, and
coastal waters,
ORD Oil Spill Response Plan,
Development of emergency response and cleanup plans for the ORD
Alternative Fuels Initiative,
Determining Circulation and Sediment Transport for EPA National Estuaries
Studies and Regional Estuaries Programs, and
Tracing the fate of radioactively contaminated marine sediments.
The primary reasons that simulations of three-dimensional stratified
flows are of added importance is the need to describe shear stress in greater
detail to fully simulate sediment resuspension and deposition, and to predict
the effects of complex flows on contaminant transport. Scientists and
engineers have long recognized that sediment resuspension is episodic and
highly variable in spatial extent, A number of methods have been developed to
measure resuspension in the field and laboratory [sea flume devices, sediment
profile measurements (Sheng et al, 1989b), core shaking methods, and
laboratory flumes]. Experience has shown, however, that these measurements
can not be made frequently enough or at enough locations to adequately
represent a mapped history of shear stresses that cause resuspension and
control deposition. As a result, hydrodynamic models calibrated with select
measurements at a few locations are the only practical approach the
determining the flux of contaminated sediments between the water column and
benthos at the moment.
Until recently, the state of the art in contaminant transport simulation
in lakes and estuaries involved calibrating a transport model (see Ambrose et
al. 1987 for example) with measurements in the water body of interest (see
Ambrose et al. 1990 for the general procedures and other examples).
Measurements of chlorides, total salt, total dissolved solids, major ions or
cations, other conservative substances, and even some non-conservative
parameters such as water temperature are collected and used to surmise what
combination of advective circulation or flow, and dispersive mixing caused the
observed concentration distributions. Unfortunately, it has been virtually
impossible to determine if these types of model calibrations were unique and
thus representative of a wider range of conditions. In effect, modelers have
been able to describe the effects of advection (or circulation) and mixing in
-------
a black box fashion (an Input-output model calibrated to describe a system
without regard to determining Important mechanisms or processes that define
the cause and effect relationships for water quality). However, the approach
rarely leads to valid predictions.
There is now more evidence that hydrodynamic simulations are required so
that water quality can be predicted and not just described as has been done in
the past. The effects of circulation and mixing must be predicted for a
number of very dangerous chemicals proposed for wide-scale manufacture or use.
These chemicals can not be released to the environment simply to calibrate a
model. Hydrodynamic transport (and effects on sediment if the chemical sorbs
to particles) must be predicted beforehand. Likewise, spills of large amounts
of oil and other materials can not be introduced simply to learn the affect of
currents, wind, and tides on its fate and to determine its effects on the
environment.
S'
Also, it is difficult to measure circulation and mixing in all flows.
Thus it is more cost effective to calibrate a model with a few selected
measurements and use the model to extrapolate to conditions of interest.
Finally, most contamination problems in simpler river and stream flows
have been cleaned up or controlled. What remain are sedimentary contaminants
contributed by diffuse, unmeasurable sources. These sedimentary contaminants
are controlled by hydrodynamics and sediment transport from the sources, and
dispersion of in-place contaminated sediments. Effective cleanup requires
investigation of leaving the contaminated sediment in place (the "no action
alternative") and investigation of various remedial alternatives. In most
cases, hydrodynamics and sediment transport must be predicted to adequately
assess the risks to human health and ecological viability. It is therefore
clear that complex stratified flows in lakes, harbors, estuaries, and coastal
areas must be understood If;
Superfund and hazardous waste sites are to assessed and cleaned up,
If existing point and nonpoint sources are to be adequately controlled,
If future sources are to regulated on a rational basis, and
If spills are to be prevented from causing extensive damage.
To address these needs, the engineers and model developers at the U.S.
EPA. Environmental Research Laboratory at Athens (ERL-Athens), Georgia have
begun development of computer programs to simulate sediment transport and
hydrodynamics. In addition, ERL-Athens has begun the development of
hydrodynamic programs to address the complex transport of dissolved
contaminants in stratified lakes and estuaries. This involves the initial
testing and application of a hydrodynamics code originally developed by Sheng
& Lick (1979) and significantly enhanced by Sheng (1983) and Sheng et al.
(1986). This code was selected because of its use in the development of three
dimensional sediment transport and dispersion models by Dr. Peter Sheng
(University of Florida) for ERL-Athens. This computer code and documentation
(represented by this report) are being developed as an interim tool that will
be improved and expanded upon after the sediment transport model has been
developed in September 1990 if necessary.
-------
ERL-Athens intends to distribute and maintain this computer code for a
limited group of engineers and scientists who investigate the effects of
complex circulation in lakes and estuaries. As this constituency of
hydrodynamics and water quality modelers grows, ERL-Athens expects to expand
the support for these types of investigations by rigorously evaluating
alternative codes and methods, streamlining application procedures, and
publishing updated codes and supporting documentation as warranted.
1.2 USJL OFJIYDRODYNAMIC MODELS BY TECHNICAL EXPERTS AND MANAGERS AND HOff
THIS MANUAL MAY BE OF ASSISTANCE
As hydrodynamics models continue to be developed, there will be some
debate about how these computer codes are best employed to avoid misuse and
misapplication. Misuse is a more important issue for complex models like
hydrodynamic and sediment transport models, contrasted with the simpler models
that are already widely used in assessing environmental problems. It is
difficult to fully document and provide comprehensive guidance for complex
codes that will assist in preventing incorrect interpretations. Complex codes
generally require greater experience and more in-depth training that is not
readily available from many graduate study programs in environmental
engineering and science. Also, these codes are being rapidly developed and it
is difficult to maintain up-to-date documentation in such cases.
To assist in the use of hydrodynamic simulation programs, this
documentation will provide several types of information useful to technical
experts and managers of projects. First, this introductory section will
review screening level studies and the potential information available.
Second, this section will briefly review procedures for calibrating and
validating models. Third, this introductory section will briefly review data
requirements.
Section 2 of this documentation will review the theoretical basis of
this code and review the structure of the program to aid in matching the
development of the theory with sections of code implementing those equations.
This section is intended for applications experts with experience in fluid
mechanics and who need a more precise definition of the limitations and uses
of the code. The practical implications of the theoretical basis of the
computer program are discussed in the Introduction of Section 2 and in Section
2.5, Turbulence Closure. These sections should be of interest to all readers.
Most important is the discussion of modeling limitations indicating that the
model can not be properly applied to:
Waterbodies where wetting and drying occurs over larger areas (i.e.,
tidal flats and shallow reservoir embayments),
« Power plant cooling water discharges, sewage diffusers, and other jets
with excess momentum (near-field effects),
Withdrawals from reservoirs, and
Flows in which short-periods waves are important causes of mixing.
Section 2 includes several important components that establish the
theoretical validity of the code for certain hydrodynamic conditions. Section
-------
2.2 on the program formulation involves a review of the governing equations,
how the model domain is described with a grid system, how the equations are
nondimensionalized, what dimensionless variables and parameters are used, how
the governing equations are rewritten in stretched coordinates, how the
equations are vertically averaged for two-dimensional applications, and how
vertical velocities are computed. Section 2.3 describes the boundary and
initial conditions that must be specified. Section 2,4 describes the external
and internal modes for numerically solving the equations. Section 2.5
describes the three means of simulating turbulence in this program. The final
two sections give more information about the grid system used to describe the
water body of interest and provide program flow charts.
Section 3 is an abbreviated user's manual intended for the applications
expert. This section is not as comprehensive as might be desired but it does
cover data requirements to initiate and operate the code, and data
requirements for calibrating and validating the model for a specific site. In
the latter part of this valuable section, input data and their formats are
described. Unfortunately, it is not yet possible to tabulate typical ranges
of all the parameters but reference to those documents that provide some of
this useful information is cite here and elsewhere in the report. In
addition, assistance to less experienced modelers in selecting model
parameters and interpreting the results is available from the ERL-Athens
Center for Exposure Assessment Modeling (CEAM). Project managers may wish to
review the introduction to Section 3 that describes the general data
requirements. To best understand critical data requirements for a specific
study, it is recommended that the applications expert setup the program to
make preliminary simulations with estimates of the input data. This will
provide the best indication of the frequency and spatial coverage required for
a data collection study.
Section 4 lists several case studies using the model, HYDR03D, and
EHSM3D and CELC3D (codes that preceded HYDR03D and for which the simulations
are essentially the same). This section should be useful for program managers
and applications experts. Program managers may wish to note the diverse
nature of the affects on circulation simulated by the program for lakes,
estuaries, and coastal areas. The case studies show that the model has been
adequately setup for the deep waters and complex bathymetry and geometry in
Prince William Sound where tidal amplitudes are typically 4 to 5 meters (13 to
16 feet), in the extremely deep waters of the Gulf of Mexico, for the complex
bathymetry and swift currents in Suisun Bay of San Francisco Bay, for the
shallower partially mixed Charlotte Harbor estuary where the tides from the
Gulf of Mexico are on the order of 1 meter (2 to 3 feet), and in wind
dominated Green Bay attached to Lake Michigan. Finally, a detailed
calibration and testing case study from Sheng (1983) is included. This was
a study of the Mississippi Sound and the adjacent deep waters of the Gulf of
Mexico. In the final case study illustrating calibration of the model, a
second, older version of the code documented in Sheng (1983) was used to
perform the simulations. As with, EHSM3D the hydrodynamics results from this
code should be essentially the same as those that could be obtained with
HYDR03D for this site.
-------
The case studies, including a comparison to an analytical solution, are
intended to demonstrate the usefulness and validity of the program. While
minor coding inconsistencies and errors may remain in any code, the extensive
use of this code and its predecessors indicate that all major problems have
been resolved and that it is ready for application.
There are other notable applications as well and these are referenced in
the report at appropriate points. For example, Smith and Cheng (1989) have
just recently examined the use of this model in a study of San Pablo Bay in
the San Francisco Bay,
The case study of Green Bay indicates the usefulness of the model in a
large lake setting where inflows, wind driven circulation, and seiche from
Lake Michigan may be important at various times. Other studies of Lake
Okeechobee due to be published in 1990 or 1991 (also see Sheng et al. 1989c),
and older studies in other parts of the Great Lakes indicate that the model is
potentially useful throughout the Great Lakes and in large shallow lakes. A
review of the theoretical basis of the model indicates that it should be
useful in smaller lakes as well. Applications in reservoirs have not been
attempted as far as we are aware. At this time, the hydrostatic approximation
seems to preclude adequate simulation of the vertical accelerations of flow
that may be important in reservoir outflows.
Other limitations are that the model does not simulate near-field
cooling water inflows, diffuser flows, and other jet discharges. Also, the
equations do not take all short-period wave effects into account.
The section on case studies should also give a project manager some
understanding of the intensity of the calculations involved and the overall
data requirements. However, only the case study for the Mississippi Sound
involves a rigorous calibration of the model and thus Section 4 only gives
some indication of overall data requirements. This is primarily because study
objectives must be integrated into such a determination.
Section 5 is the Programmer's Guide. This is intended to give limited
assistance the applications expert to install and run the program. At this
time (May 1990), the code has a few VAX specific FORTRAN Statements that users
must modify when working on other processors. We believe this will take about
one-man week of effort and expect to resolve these problems before the final
release of the code by July 1, 1990.
1.3 SCREENING LEVEL SIMULATIOHS
One of the more important debates among hydrodynamics modeling experts
is whether or not these complex models can be applied in a screening mode with
available data. In the case studies, we illustrate how a model can be setup
for illustrative purposes by using examples from Suisun Bay and Charlotte
Harbor first used in a workshop by Sheng et al. (1986). We also investigated
the use of the model in the recent oil spill emergency in Prince William Sound
and these preliminary results are reported herein. For Prince William Sound,
we found that the code could be setup in a matter of days and made ready for
-------
follow-up studies to calibrate the model for longer-term assessments and
management of the cleanup. However, it is noted in the studies and
highlighted by the reviewers of this report that simulations based on limited
data can be misleading and highly suspect if not properly interpreted.
Nevertheless, it is clear that some useful information is to be gained and we
now include it as a useful means to obtain limited information or screen out
some alternatives. Primarily, we recommend screening level simulations using
existing data to design data collection programs for model calibration, for
extrapolation of tide and current measurements to areas not covered by
existing data, and for preliminary investigations of effects on circulation
and transport in place guesswork and approximate means even if adequate
calibration data are not available.
1.4 CALIBRATIQfl ANDVALIDATION
Calibration and validation of a computer code applied to a specific site
results in a simulation model of the site. The ability of the model to
describe or predict conditions at a site depends on how well the code is
calibrated. Calibration is the process of selecting model parameters and
configuring the computational domain to be simulated. Model simulations based
on alternative selections of parameters are compared with measured data. The
coefficients used in the simulations that best match the measurements are
chosen as the calibrated parameters. Validation is used to determine how
uncertain the results of the model are for limited ranges of conditions in the
water body of interest (lake, estuary, or coastal waters). For additional
information on calibration and validation see McCutcheon et al. (1990) for
estuary modeling, and Chapra and Rechow (1983) and Henderson-Sellers (1984)
for lake modeling.
The general procedure for calibrating and testing the adequacy of a
model is as follows:
Determine study objectives,
Define the subset of objectives to be addressed by model studies,
Collect historical data from monitoring or previous studies,
Attempt a preliminary calibration of the model,
Design a calibration data collection study based on the preliminary
calibration,
Simulate conditions during the calibration period and compare to
determine if the preliminary calibration Is sufficient (if it is not,
calibrate the model),
If the model is calibrated, collect a second independent set of data for
validation,
Validate the model for the limited range of conditions defined by the
calibration and validation data sets (if the model can not be validated
repeat the calibration step and collect more validation data for a second
attempt), and
Determine uncertainty in the calibrated model simulations by sensitivity
analysis.
See Ambrose et al. (1990) and Ambrose and Martin (1990) for guidance on these
7
-------
procedures.
1-5 DATA REQUIREMENTS
In general, data requirements for calibrating a hydrodynamics model are
extensive but not overwhelming. Typically, the following data are necessary:
Navigation charts and maps or soundings to define bathymetry and
geometry,
Measurements of current speed, salinity, and water temperature at two or
more levels at each of important boundaries of the water body (mouth of
estuary, lake outlet, fresh water inflows, rivers, etc.),
Measurements of current speed, salinity, and water temperature at two or
more levels at a number of stations throughout the water body (for
calibration and validation),
Some measurements of the initial condition of the current, salinity, and
temperature fields at the beginning of the simulation, and
Long term monitoring stations for water level to identify critical
episodes or seasonal changes.
Multiple stations may be necessary to define some open boundaries. Five to
ten stations in the domain should be sampled to calibrate and validate a
model. Long term monitoring at one or two stations in the study area is
desirable. Sampling frequencies depend on study objectives (as do sampling
location to some extent). Calibration for episodic events requires data
collection at the boundaries and internally over a period that at least
exceeds the occurrence of the event and hopefully defines conditions
beforehand and after. Simulation of seasonal changes requires multiple
deployments of current meters and water quality sampling equipment over longer
periods.
See Section 3 for more details on data requirements.
-------
SECTION 2
MODELING SYSTEM
Section 2 presents the basic model theory and describes how the
modeling system is formulated. The section is intended to assist engineers
and scientists charged with the application of the model, but some of the
introductory material may interest project managers. Unfortunately, this
document can not cover many of the basic elements of hydromechanics that are
generally needed to fully grasp the limitations a complex hydrodynamics model,
and the reader should look elsewhere for this information. See for example,
White (1977), Monin and Yaglom (1971), Reynolds (1974), Rodi (1980), Hinze
(1959), Schlichting (1979), Goldstein (1960), Turner (1973), and Tennekes and
Lumley (1972) among the few good references that can provide useful background
information. Sheng (1983) provides additional discussion ofthe theoretical
basis of the model not covered here.
HYDR03D is a FORTRAN code designed to simulate two-dimensional (2-D)
and three-dimensional (3-D) stratified (or non-stratified) flows in lakes,
estuaries, coastal waters, and harbors. In solving for the effects of density
stratification on circulation, the model also simulates the distributions of
dissolved solids (salt) and temperature. These flows and the associated mass
and heat transport are simulated dynamically. Important forces taken into
account in the simulations are those caused by tides, winds, density gradients
caused by salt (dissolved solids) and heat, and forces due to the resistance
of flow over irregular bathymetry and around irregular geometry in the water
body of interest. The grid system that users set-up for model simulations is
rectilinear in the plan view but uses a sigma stretched grid in the vertical
direction. A sigma stretched grid divides the depth into vertical layers of
equal thickness and maintains those equal thicknesses even as the total depth
of flow changes during the simulations.
The governing equations solved by this model are an approximation that
are designed to simulate some water flows but not all possible conditions that
arise in the natural environment. The important approximations in the
equations are related to the manner in which turbulence in the flow is
simulated, and the treatment of vertical velocity accelerations.
The original principles from which the governing equations for any
mechanistic hydrodynamics model are derived include:
Newton's second law that force is equal to mass times acceleration,
Conservation of water in a defined volume,
Conservation of heat, and
Conservation of dissolved solids or salt.
-------
Application of these principles results in a set of equations for the sum of
the forces acting on a fluid element in all three directions of three-
dimensional space, plus the conservation equations for water, salt, and heat.
Fluid density is calculated from dissolved solids and heat using an equation
of state.
If it is noted that Newton's law of viscosity (Streeter and Wylie 1975)
adequately relates fluid viscosity to shear stress in the fluid (viscous drag
force caused by fluid moving over the bottom or other layers of fluid), then
the original force balance equations (derived from Newton's second law) can be
expressed in an mathematically exact form known as the Navier-Stokes equations
(White 1974, Schlichting 1979, Monin and Yaglom 1971). Newton's law of
viscosity is essentially an empirical formulation but it is based on extensive
observation and can also be explained by several mechanistic and conceptual
premises. As a result of the extensive and comprehensive observations, one
can confidently treat all water flows as Newtonian, i.e., there is a linear
dependence between water viscosity and shear stress of the flow (Streeter and
Wylie 1975, Bird et al. 1977).
Unfortunately, the Navier-Stokes equations can not be solved exactly for
most turbulent flows because there is insufficient computer memory and speed
available from present day processors (including supercomputers). See Rodi
(1980) for a discussion of computing needs to solve a typical environmental
fluid mechanics problem using the Navier-Stokes equations. To derive a
practical means of solving the equations, an averaging technique dating back
to the 1894 work of Sir Osborne Reynolds (for whom the Reynolds number in
fluid mechanics is named), is typically employed (see Monin and Yaglom 1971),
This technique resolves the turbulent velocity and mass transport into two
components; a mean velocity or concentration (or temperature for the
conservation of heat equation), and a fluctuating component typically written
for the three coordinate directions, i, J, and k, as Ui -= UA + ut' , Uj <= u^ +
Uj' , and Uk = uk + uk' and for concentration or heat as C = c + c' or H <= h -»
h'. The fluctuating component about the mean velocity or constituent property
of the flow (temperature or concentration) is determine by the period of time
over which the properties are averaged. When the mean and fluctuating
components are substituted into the Navier-Stokes equations to achieve what is
known as Reynolds averaging of the equations, the resulting Reynolds equations
(White 1977) can be readily solved for many types of environmental fluid flows
with several numerical methods (Rodi 1980).
Reynolds averaging is a powerful technique that makes a number of
numerical and analytical solutions possible (Schlichting 1979).
Unfortunately, averaging introduces two severe disadvantages to overcome in
solving the resulting equations. First, the Reynolds equations are no longer
time-continuous (or, exact dynamic) expressions. Second, the substitution of
two variables (the mean and fluctuating components) for one variable, results
in a mathematical closure problem, i.e., there are no longer enough equations
to solve for all the unknown variables.
The averaged equations represent fluid motion as a series of averages
and fluctuating components that change from average values of both components
in one discrete interval to other average values in the next discrete
10
-------
Interval. This change in behavior can be handled by a number of finite
difference (used in this model) and finite element numerical schemes, as well
as by other numerical schemes such as the method of characteristics (Lai
1965), In fact, these schemes can be used to solve the equations over short
enough time intervals that the solution describes most of the important
dynamic behavior (or turbulence) of the fluid flow. As a result, the solution
is finely resolved enough in most cases to refer to the simulations as being
dynamic (as is done for this model). This is normally the cases for flows in
estuaries, lakes, harbors, and coastal waters. However, the optimum averaging
interval (model time step) can not be theoretically derived and thus must be
arbitrarily selected using judgement and experience.
In addition, the time interval for averaging effects the value of eddy
coefficients through a direct effect on the fluctuating component that will be
further described below. In effect, this introduces another correlated
parameter to the selection process that includes time intervals, spatial
segmentation (how the grid is set-up), and eddy coefficients for momentum,
mass (contaminants or salt), and heat. In certain ranges, these parameters
are not highly correlated (i.e., one is not very sensitive to values of the
others); but in general one can not select a grid to represent a water body, a
time step for the solution, and eddy coefficients independently without
understanding the mutual effects. As a result, the time averaging is
necessary to solve the equations but introduces a need for some experience and
guidance in applications.
The second drawback in the averaged equations is the addition of new
variables, making it impossible to achieve mathematical closure of the set of
equations. This becomes the problem known as turbulence closure (Rodi 1980,
Sheng 1983) that has been studied extensively (Hinze 1959, Rodi 1980, Monin
and Yaglom 1971, Reynolds 1974, Tennekes and Lumley 1972).
Turbulence closure is effectively achieved by deriving additional
equations for momentum, mass, and heat transport (see McCutcheon et al. 1990).
The equations that can be derived, range from simple expression for mixing
lengths (Prandtl 1925) and eddy coefficients (Streeter and Wylie 1975) related
to various mean flow properties. Mixing length and eddy coefficient methods
have been classified as zero-order turbeulence closure schemes (Rodi 1980),
Higher-order closure schemes can be derived from conservation of kinetic
energy, expressions for turbulence length scales, and other approaches (see
Rodi 1980 and ASCE 1988 for comprehensive reviews). Regardless of the
approach, each equation has at least one or more empirical constants that must
be determined from observations. The highest order schemes have only one or
two constants that may be widely applicable to most environmental flows (Rodi
1980, 1984, ASCE 1988) but these will never be universal constants because of
the closure problem. The more approximate schemes, i.e., eddy coefficients,
remain much more empirical (McCutcheon et al. 1990). Equation coefficients
and parameters are quite variable from one flow to another, and vary within a
flow at different locations and at different times. As a result is it is
difficult to forecast eddy coefficients. Therefore, calibration with field
measurements is necessary for precise studies. There is less uncertainty in
the parameters of the higher order schemes but these schemes are not fully
practical. Therefore, the state of the art is to use eddy coefficient methods
11
-------
(McCutcheon et al. 1990), but there are some cases where higher schemes are
useful and this model has one option involving a simplified second order
scheme that will be described later in this section. In general,
practitioners should recognize that higher order closure schemes are expected
to attain recognition as state-of-the-art in a few years and should begin
using the methods forthwith.
The basic approach used in this model to achieve turbulence closure is
to use the eddy coefficient approach. The three options available include:
Constant eddy coefficients,
Munk and Anderson type vertical eddy coefficients, and
Simplified second-order scheme expressed in terms of eddy coefficients.
The eddy coefficients are derived in the terms of the Reynolds equations where
the fluctuating components describing the turbulent momentum, mass, and heat
flux terms (Rodi 1980, ASCE 1968), These turbulent fluxes are assumed to be
proportional to the vertical gradients of mean velocity, concentration, and
temperature. The proportionality constants in these expressions, are the eddy
coefficients (see Rodi 1980) used in this model.
2,1 OVERVIEW OF THE MODELING SYSTEM
Development and implementation of a complex mathematical model such as
HYDR03D requires that the fundamental concepts be formulated in a clear and
concise fashion. To facilitate the interpretation, maintenance, and upkeep of
the computer code it is necessary to use structured and modular programming
techniques. In accordance with these criteria the HYDRO3D modeling system
consists of 65 subroutines, 2 INCLUDE files and post-processor files for
graphical presentation of results.
The governing equations of the model are incorporated in a discrete form
using a finite difference numerical method coded in VAX, FORTRAN 77. In
solving the equations the program can be run with either a fixed time step or
a variable time step.
HYDRD3D is a dynamic model that allows the specification of a variable
wind field and a variable river inflow or outflow. It also allows tidal
forcing boundary conditions at multiple sites. The vertical turbulence
closure parameterization schemes of HYDR03D include; 1) constant eddy
viscosity, 2) variable (Munk-Anderson) type eddy viscosity, 3) and a
simplified version of a second-order closure model.
Although every effort has been made to develop a general purpose
hydrodynamic modeling package, HYDR03D has its limitations like any other
computer code. In spite of the many special features contained in the code,
HYDR03D does not have the full capability to allow universal application of
the model to all water bodies under all physical conditions with arbitrarily
chosen grid patterns and time steps. To understand the limitations of this
program, those planning to apply the model are advised to read this manual and
12
-------
a previous report (Sheng 1983). Experience has clearly shown that the user
must thoroughly understand the capabilities and limitations of HYDR03D before
attempting to solve a site-specific problem,
Specific limitations of the modeling system that are obvious from the
governing equations and that have been noted from applications of the program
include:
Use of hydrostatic pressure distribution that precludes detailed
simulation of the effect of jets and other near-field mixing phenomena
from cooling water discharges, sewage diffusers and pipes, and other high
momentum discharges. It also precludes simulation of reservoir
withdrawals where vertical accelerations are significant.
Exclusion from governing equations of the effects of short period gravity
waves that describe near-shore circulation.
Lack of flooding and drying features to simulate tidal flats and deltas
in lakes where tidal amplitudes or water surface elevation changes are
large and the near shore bottom slopes are small.
Effects of grid resolution, especially near boundaries, that may slow the
speed of the calculations or cause erratic results if rapidly changing
bathymetry is not adequately resolved.
In addition, there are other less general limitations related to turbulence
modeling and other features that may cause problems in a few applications for
inexperienced users. These will be introduced in the following material as
the need becomes obvious. Since the range of experience of users is not clear
in the initial stages of development, however, these experience-related
disadvantages of the modeling system can not be fully complied in this first
edition of the documentation. These will be compiled in future documents as
feedback is received from users of the model and documentation. Therefore,
present users should have an adequate understanding of the physics of water
circulation, numerical methods, and computer programming to understand and use
the model.
2.2 MODEL FORMULATION
The governing equations used in this program and the assumptions aon
which the equations are based are described in the following subsections.
2.2.1 Governins Equations
The governing partial differential equations are based on the following
assumptions:
The hydrostatic pressure distribution adequately describes the vertical
distribution of fluid pressure,
The Boussinesq approximation Is useful (small density differences in
stratified flows are assumed to have a negligible effect on fluid
inertia), and
The eddy viscosity approach adequately describes turbulent mixing in the
flow.
13
-------
Use of the hydrostatic pressure distribution means that vertical accelerations
of the fluid are ignored. This generally limits the model to simulations of
far-field conditions. Significant deviations of simulations with actual flow
conditions may occur in near-field flows involving jets into receiving waters.
This is especially true if the high momentum effects more that 5 to 10 percent
of grid points. For more information on the hydrostatic pressure
distribution, see White (1974) and Sheng (1983).
As a result of the use of the hydrostatic pressure distribution, it is
not readily feasible to use HYDR03D to simulate the detailed behavior of
cooling water inflows (from coal-fired or nuclear plants), sewage inflows from
diffusers, outfalls, pipes, and other jet-like discharges into lakes,
estuaries, coastal waters, and harbors that involve high velocities and flow
rates. However, the existence of Intense near-field mixing in limited areas
does not preclude use of this model. For example, the existence of a diffuser
or high velocity jet from a pipe or channel into the vicinity of one or two
grid points (a few at most) of the computational domain can be simulated.
Simulations that compensate for the effects of near-field mixing usually
involve;
Specification of elevated values of the eddy coefficients at the affected
grid points, or
Design of the computational domain to exclude the high momemtum areas from
the simulation.
Near field mixing is usually computed according to the approaches in EPA
(1985), Fischer et al. (1979), Jira and Donecker (1988), Roberts (1979),
Wright (1977), and Jirka (19B2). These calculations could be used to estimate
elevated eddy coefficients or calculate the expected mixing for boundary
conditions to the hydrodynamics model when the hydrodynamics model domain
excludes the near-field effect. These procedures must be worked out on a case
by case basis, but these approaches represent to state-of-the-art at the
moment (1990).
Model users should note that neither the selection of elevated eddy
coefficients nor the selection of an ideal model domain can be accomplished
easily. Normally, users should expect to calibrate the model and collect
extra data in the vicinity of the jet(s), There are no reliable means of
relating eddy coefficients to jet mixing averaged over large scale distances.
Extra data will be required if boundary conditions are quite variable and if
the effect of the jet extends into the model domain. To reduce extra data
collection, selection of boundaries to isolate the effects of jets is
recommended. For example, moving model domain boundaries to the mouth of
embayments or arms of an estuary or lake will avoid the complications. Of
course, selection of a different model (one that solves the vertical momemtum
equation) should also be considered as well.
The use of the Boussinesq approximation does not seem to severely limit
the application of the model. The approximation is normally valid for most
water flows in the natural environment (see Monin and Yaglom 1971).
14
-------
The use of the eddy-viscosity concept indicates significant limitations
of the governing equations, as implied in the introduction of this section.
In general, a simple eddy viscosity scheme does not keep track of the
generation and dissipation of turbulence (Rodi 1980). The eddy viscosity
scheme is based on the assumption that the flow is uniform and that the
turbulence is dissipated under the same conditions under which it is
generated. Unfortunately, these are conditions that do not exist in complex
flows (multi-directional at different levels, reversing with time, etc.), and
especially in stratified flows (see McCutcheon et al. 1990 for more discussion
of these limitations). In general, turbulence is generated by complex
interactions from wind shear, fluid shear on the bottom, flow around islands
and obstructions, and internal density currents, as well as by other
mechanisms not included in HYDR03D (i.e., turbulence due to waves). This
turbulence is transported throughout the water bodies of interest and can be
dissipated under very different conditions. For example, bottom generated
turbulence from tidal flats can be swept into deeper channels where the
turbulence is of different scales and intensity. Also, salt stratified flows
are nonuniform (the vertical salinity gradient varies downstream as the
density differences decay due to mixing across the interface) and bottom
generated turbulence is dissipated under different salinity gradients
downstream. In many cases, the generation and dissipation conditions are not
radically different and the eddy coefficients are useful. But in general, the
transport of turbulence must be taken into account, especially if the model is
used in forecasting. This is one reason why the simplified second order
closure scheme employed by this code is important.
From the assumptions outlined above, the basic flow equations for an
incompressible fluid (i.e., water) can be expressed using the right hand
Cartesian coordinate system (with x,y,z) shown in Figure 1. These equations
are written as:
Continuity Equation:
cJu 3v 3w
+ + - 0 (1)
3x 3y 8z
Momentum Equations:
du du2 8uv 3uw 1 ap
B f 3u| 3
AH +
ax[ axj ay
+ - + - fv - - -- + |AH I + K | + ^IAV - I (2)
3t 3x 3y 8z p0 fix
flv duv flv2 9vw
at ax ay
i dp 9 ( avl a ( avl a f av>
_ _ fu _ _ _1 + _ AH + AH + AT <3>
PD ay a-x.{ axj ay[ ayj az[ azj
Sp_
Bz
15
-------
,- Displaced Water
Surface
*^^\
z, w
Bottom
Nominal Water Surface
Figure 1. Cartesian coordinates at the nominal water surface.
16
-------
Temperature and Salinity Equations:
9T 3uT
at 3x ay 3z
as aus avs aw
(. H +
3t 5x 8y 8z
Equation of State:
f aTl 3 { ail a f ail
Ke + KB + Kv <5)
L axj 8y[ ayj 8z[ azj
8 I as] a [ as! a f as]
= 1DB I + I Dg I + I DV I (6)
3x[ 3xJ 5y[ 3yJ flzt 3zJ
(a + 0.698P)
where
p - P/(a + 0.698P)
P = 5890 + 38T - 0.375T2 + 3S
a = 1779.5 + 11.25T - 0.0745T2 -(3.8 + 0.01T) S (7)
Equation 7 is based on the equation of state by Eckert (1958), where
temperature, T, is in degrees centigrade, salinity, S, is in ppt (part per
thousand) and density, /?, is in g/cm3.
In Equations 1 through 6, (u,v,w) are velocities in (x,y,z) directions
(see Figure 1), f is the Coriolis parameter defined as 20 sin$ (where ft is the
rotational speed of the earth, and is the latitude), p0 is a reference fluid
density, p is the variable density, p is pressure, T is temperature, S is
salinity, (AH, KH, DH) are horizontal turbulent eddy coefficients, and (Av, Kv,
Dy) are vertical turbulent eddy coefficients for momentum, mass, and heat,
respectively. In addition to the above equations, HYDR03D includes an
equation for dissolved species concentration written similar to Equation 6.
2.2.2 Grid System
HYDR03D uses a vertically stretched grid, i.e., the so-called "cr-
stretching", which leads to a smoother representation of the topography and
the same order of vertical resolution for the shallow and deeper parts of the
water body as shown in Figure 2. The transformation is done using the
following equation,
(8,
h(x,y) + f(x,y,t)
cr-stretching offers some advantages over z-grid configurations, but
there are some problems that can arise in setting up the grid for steep bottom
slopes. The a-stretching 3-D model allows the smoother resolution of
bathymetry that the stair-step configuration of z- grids (Leendertse and Liu
1975). It is necessary, however, to have sufficient resolution in the
horizontal and lateral directions of the p-grid to avoid erratic results. In
addition, Sheng et al. (1989b) and Johnson et al. (1989) report that the z-
grid is better that the cr-grid in the presence of steep bottom slopes. It was
17
-------
Figure 2. Vertical stretching of the coordinates.
18
-------
found that it is advisable to evaluate the baroclinlc gradient along the
constant z-plane (horizontal plane) and that higher order advection schemes
should be avoided.
The a-stretching introduces extra terms into the original Cartesian
equations of motion. However, most of these extra terms appear in the
horizontal diffusion terms, which are generally less significant.
In the horizontal direction, HYDR03D allows the use of either a uniform
or a non-uniform Cartesian grid. For non-uniform grids, there are two
options. The HYDR03D program will accept either an arbitrary, non-uniform
grid or a smoothly varying stretched grid (see Figure 3) which satisfies the
following equations;
x=ax + bIaC- (9)
y - ay + by i ^ (10)
where (0,7) represent the uniformly-spaced computational grid and (a^ bx, Cx,
ay, by, Cy) axe stretching constants, A uniform grid can be obtained by
setting: a^ = 0, bx <= 1, Cx = 1, ay=0, by-l, and Cy - 1. To generate a
non-uniform grid according to the transformation Equations 9 and 10, the
example in Sheng and Butler (1980) can be followed. The detailed procedure
for deriving the cr-stretched grid can be found in Sheng and Lick (1980) and
Sheng (1983).
The lateral stretching introduces stretching coefficients, ^x = dx/da
and py dy/d-y into the spatial derivative terms in the transformed equations
of motion, when a non-stretched grid is used, then nx = /iy = 1.
2.2,3 Kon-Dimensionalization of the Governing Equations
Diraensionless governing equations make it easy for a user to compare the
relative importance of various terms and to minimize numerical errors. When
properly non-dimensionalized, the part of each terra contained within a
parenthesis or bracket of Equations 13 through 17 that follow should be of
unity order and the part of the term containing the dimensionless number(s)
will indicate the order of the term.
The governing equations are nondimensionalized using the following
dimensionless variables;
f Y 1
(u*, v*, w*) = -, ] -
1 Ur Ur Zr Ur J
19
-------
y
I
*- I
Figure 3. Lateral stretching of the coordinates.
20
-------
... f x y zX
t y i ^ / i*r i s » .^^
I xt xz zptr
t* - tf
- T0) />0
gr r
Pr - Po
Tr - T0
_AH
^
A,
Avr
KH
^Hr
21
(ID
-------
Where Xr is the reference length in lateral direction, usually the maximum
dimension of the basin or water body, Zr in the reference depth, usually the
average depth of the basin, the Ur is the reference velocity, w is the
vertical velocity in the a-direction, t is time, f is the Coriolis parameter,
q is the heat flux at the surface, C is the specific heat of water, Sr
fUjJCj./g, f is the displacement of the water surface at any time as defined in
Figure 1, pD is a base fluid density usually taken as the minimum density in
the water body being simulated, p is the variable density, pr is another
reference density usually taken as the maximum density, T is water
temperature, T0 is a base water temperature usually taken as the minimum, Tr
is a reference water temperature usually taken as the maximum, AH, Kg, DH are
horizontal turbulent eddy coefficients, and Av, Kv, DV are vertical turbulent
eddy coefficients for momentum, mass, and heat, respectively. Some of these
parameters are defined in Appendix A for future reference.
The following dimensionless parameters are derived when the equations
are nondimensionalized:
Avr ti
Vertical Ekman Number: E
Lateral Ekman Number: Eg =
Vertical Prandtl Number: Prv
Lateral Prandtl Number: PrB -
Vertical Schmidt Number: Scv =
Aflr
Lateral Schmidt Number: ScH -
22
-------
Froude Number: Fr
Rossby Number: Ro -
ur
£Xr tc
gZr (Ro)
f2 X2. (Fr)2 (t,e)z
sr
Densimetric Froude Number: FrD - -
fr t0 (12)
where the various t variables appearing in the above expressions represent the
characteristic time scales associated with various physical processes. The
inertia oscillation time scale is t±, vertical turbulent diffusion time scales
are tv tvdh ^vte» the lateral turbulent diffusion time scales are tidm, ttdh,
tids, the convection time scale is te, the gravity wave time scale is tge, and
the internal gravity wave is tgi. These are better defined in Appendix C.
2.2.4 Dimensionless Equations in Stretched Coordinates
Using the dimensionless variables and parameters in Sections 2.2.3 and
dropping the asterisks for clarity, the following dimensionless equations can
be derived:
8j da
23
-------
1 9Hu 1 Sf E., 3 f 5u > Rof 1 9Huu 1 aHuv 3Huw
J i ^V 1 A i i i ^ a_ i ^ _s
Rof 1
-rl^
H 9t nx ax W- da [ da ) H J_ /ix 5x >jx 3y
+ E r-£-L-^-J +-^ fA H +HOT!
^ fiH| .. a,, I^H .. 3i + ,. a., IAH ,. a.J * n.U.1.1
f" (*fl - X1
-------
2.2.5 Vertically Integrated Equations
For vertically mixed estuaries, the governing equations can be
integrated over the depth. The resulting a non-dimensional form is given as:
au
(19)
au H ar
+ r T + v - Ro
at u, 9x
au
Wsx
au
Ro H2 3/5
Frg 2
H aj-
- -f D
(20)
av
at
H
- U - Ro
[ 9 (
\ h
L^x [
Ro H2
Frn 2
H
(21)
where the variables are defined in Appendices A and C.
The nonlinear inertia, lateral diffusion, and baroclinic pressure
gradient terms in Equations 20 and 21 are obtained by vertically integrating
the corresponding terms in Equations 14 and 15, respectively. However, these
terms are obtained by assuming that horizontal velocity and density are
uniform in the vertical direction, an assumption which is not always valid.
In addition, the vertically integrated equations ignore the baroclinic terms.
Thus, the above forms of vertically integrated equations are not recommended
for all flows, especially those involving baroclinic circulation. When
barolclinic circulation is important, the fully 3-D version of HYDR03D should
be used. Details of the vertically integrated equations are briefly discussed
later in this section.
25
-------
2.2.6 Vertical Velocities
Equations describing the vertical velocity in the transformed coordinates
are:
ar Iff
at Hj.it
1-Hr ar Iff SEu 3Hv 1
= ----- - + - dcr (22)
(i +
-------
where qs is the heat flux at the surface, rbx is bottom shear stress in the x
direction, rby is the bottom shear stress in the y-direction, and 9u/3a and
Sv/9a are the partial derivatives of longitudinal and lateral velocities with
respect to the a coordinate defined earlier.
Idu 3v 1
, are in tensor form notation. The other
da 80 J
parameters were defined earlier.
The boundary conditions at the bottom (CT = -1) are:
f3u 5v 1 H ur , ,
A I /, , \ _ U 7 T fn2 -1- -rr21!1/2 fn V 'I
^v i <^bx rby; ~ « iri>d (.% 1- v^; {Ui , v1(>
[da da J EL f^^
^ ^ __ rt
5S
- 0 (25)
So
where Cd is the drag coefficient applied to the bottom surface, and (ulf vx)
represents the velocity vector components at the first grid point above the
bottom. The drag coefficient is related to the Manning roughness coefficient
for bottom roughness as shown by McCutcheon et al, (1990). Also see Sheng
(1983), McCutcheon et al, (1990) compiles representative values of the
Manning for estuaries,
2.3.2 Lateral Boundary Conditions
Along the shoreline where river inflow or outflow may occur, the
boundary conditions are;
u - u(x,y,cr,t)
v -= v(x,y,<7,t)
w ~ 0
T - T(x,y,ff,t)
S = S(x,y,P,t) (26)
where u, v, T, and S are velocity in the x-direction, velocity int he y-
direction, water temperature, and salinity, respectively varying dynamically
with time t and spatially in the (x, y, a) coordinate system, u> is the
vertical velocity, assumed to be zero.
27
-------
At the solid boundary, both the normal and tangential velocity
components are equal to zero. In addition, the normal derivatives of
temperature and salinity are also set equal to zero.
Along an open boundary, either f or the velocity can be specified. For
the water surface elevation f, there are currently three options:
(27)
or
- 0 and = 0 (28)
dx 3y
or
ft t c |£ - 0 «d J£ ± c | = 0 (29,
where AQ, Tn, $n are the amplitude, period, and phase angle of the tabular
tidal constituents, respectively; ^^ is the maximum number of these tidal
constituents; and c is the dimensionless phase speed of surface gravity wave
at the open boundary. See the Case Studies for Suisun Bay and Charlotte
Harboer (Equation 62) in Section 4 for an illustration of the application of
Equation 27.
When open boundary conditions are given in terms of f , the normal
velocity component is assumed to be of zero slope. The tangential velocity
component may be either zero, or of zero slope, or computed from the momentum
equations .
The salinity or total solids along an open boundary or river entrance is
computed from a 1-D advection equation during the outflow. For example, along
an open boundary perpendicular to the x-direction.
3HS 3HuS
- + RO - = 0 (30)
St pxdx
where the spatial salinity flux is evaluated from the salinity values at the
boundary and the interior grid point via a one-sided differencing scheme.
During the inflow, however, the salinity value at the open boundary can
either take on a prescribed value or be determined from the 1-D advection
equation while using the boundary salinity value and the prescribed salinity
value to evaluate the spatial flux term. Section 4.3 on the Charlotte Harbor
case study regarding the application of Equation 30 and other options.
28
-------
2.3.3 Initial Conditions
To start a simulation, the initial spatial distributions of f, u, v, w,
T, and S must be specified. However, when initial data are completely
unknown, there is usually little choice but to start with zero initial fields.
This is a process referred to as spin-up. It invilved starting at resta nd
running the model until reason conditions are attained. These attained
conditions become the initial conditions of the next simulation if the
simulation is stopped and restarted. If the simulation is continued, this
point in time defines the beginning of the results that will be used to
investigate circulation and mass transport. See Case Study 4.3 on Charlotte
Harbor for a brief review of the procedure. When initial data are known at a
limited number of locations an initial field can be generated by an
appropriate interpolation scheme. In principle, the interpolated field must
satisfy the conservation law governing that field variable.
For practical simulations of barotropic flow in the absence of salinity
and temperature variation, the HYDR03D code usually assumes zero initial flow
if few initial data are known. This is reasonable because the spin-up time of
a barotropic flow field is relatively short due to the use of a variable
time-stepping scheme. In case of a baroclinic simulation where salinity and
temperature varies with space and time, the spin-up time is longer and an
interpolation routine is provided to produce a reasonable initial field from
limited data points.
2.4 NUMERICAL SOLUTION ALGORITHM
2.4.1 External Mode
In the external solution mode the model solves for the surface
displacement f and the vertically integrated velocities U and V in Equations
19 through 21. To speed up the model simulation, all the terms in Equations
19 through 21 related to the propagation of surface gravity wave are treated
implicitly. The time derivatives and surface slopes in the momentum equations
are generally treated implicitly, whereas the bottom stresses are computed
explicitly from the latest vertical profiles of horizontal velocities.
The dimensionless finite-difference equations needed to obtain the
external mode solution given in matrix notation as:
[A] {F]n+1 - [I] {F]n + At ID}" (31)
29
-------
where:
[A]
HAt
HAt
(D)
[I] -
1 0 0
010
0 0 1
[F) -
f
u
V
(32)
30
-------
where Ax, Az, At, Sx, and Sy are space intervals in the x-and y-directions,
the time interval, and delta notations in x, y directions, respectively.
It is convenient to express the matrix [A] as the sum of three matrices,
These are the identity matrix [I], a matrix [A^J and a matrix [Ay). The
matrices [\] and [Ay] are written as:
[A,]
HAt
sx
0
[Ay] =
HAt
My
0
/9At
(33)
The factorization of Equation 31 and neglecting of terms of order At2 yields
the following x-sweep and y-sweep expressions:
x-sweep:
([I] + [AJ] IF}* = ([B] - [Ay]) (F}H + At {D)N
y-sweep:
[I] + [Ay])
IF}* -t- [Ay]tF)N
(34)
(35)
31
-------
where N, and K+l are the successive time step counters; [F]N+1, [F]N and
[F]* are the solution vectors at time steps N+l, N, and the intermediate
between the x-sweep and y-sweep, respectively. {D)N is the residual( or
forcing) matrix at time step N.
Alternatively, the x sweep and y sweep are also listed here can be
written as:
x-sweep:
and
f* +
U*
(36)
U* +
HAt
° At
(37)
y-sweep:
S v
,
n+1
SAt
Vn
(38)
and
HAt
Vn
At
(39)
where IF1, Vn are velocity matrices at time step N; and D£, and D£ are the
forcing or residual matrix in the x and y-sweep at time step n, respectively.
2.4.2 Internal-Mode.
The internal mode solution is obtained by defining deficit velocities as
u = u - u and v =* v - v, where u and v are vertically averaged by subtracting
the vertically averaged momentum equations from the three-dimensional momentum
equations multiplied by H. The resulting differential equations are:
1 3Hu Dx K. 8
= B_ - +
H 3t H H2 8(7
a
(u + u)
(40)
H at
H
H
A, (v + v)
8z
32
-------
which can be written in a finite difference form as:
(Hu)
n+1
(Hu)n + At (HBX
Hn At 2
(u + u)1
8z
(42)
(Hv)n+1 - (Hv)n + At
- Dv)n + HD At -^
y (Hn)2 5r
A, (v +
8z
(43)
The vertical diffusion terms in the momentum equations are treated implicitly
to ensure numerical stability in shallow water. It also is important to
ensure that the vertically integrated deficit velocities always equal zero or:
Ua
(44)
(45)
k-l
where Kj^ is the maximum number of the vertical layers. To ensure that
Equations 44 and 45 are satisfied, the nonlinear inertia, baroclinic, and
horizontal turbulent diffusion terms in the vertically integrated Equations 20
and 21 must be evaluated by summing the corresponding terms in the 3-D
equations at all vertical levels.
Once un+1 and vn+1 are obtained, u and v can be obtained from:
- un+1
(46)
H"
(47)
33
-------
Following these, the vertical velocities wn+1 and w""1"1 can be computed
from:
n+1 nn Aa f ar)n+1 Aa f f 3Hu )n
t,> u? - I I - I I
k *-* 0H [ fltj H L I M* Jk
(48)
av, ln-t-1
n+l .. ^4.1 4-TU I "S I I on on
r - tf*1 «£ i + - + .7 u-+v- (49)
k k-l a I *- I I ,. a,, ., a,.
where the vertical velocity w should be almost zero at the free surface,
2.5 TURBULENCE CLOSURE
Turbulence parameterization is necessary because of the averaged nature
of the governing equations, as was discussed in the introductory part of this
section. As noted, there is some art to selecting turbulence parameters.
This selection process is therefore aided if there are several methods
available to give the user some flexibility to chose a method tailored to his
experience and the conditions at the site of interest. This code offers
relatively good flexibility ranging from a simple constant eddy coefficient to
a simplified second-order closure approach. In all, three methods are
available to determine the vertical eddy coefficients. These are:
Specification of constant eddy coefficients,
Calculation by the Munk-Anderson formulation, and
Calculation using a simplified second-order closure scheme.
There are a number of other closure schemes, but all have various theoretical
and practical disadvantages, as do these options. See Sheng (1983), Blumberg
(1986), Rodi (1980) and ASCE (1988) for more detailed information on these
methods and other alternative methods. At the moment, we can not offer the
users comprehensive guidance on turbulence modeling in this document.
However, when simulations are sensitive to turbulence parameters and
assistance is need, users should contact the Center for Exposure Assessment
Modeling (CEAM), U.S. EPA ERL-Athens (see Preface). The CEAM can assist in
calibrating a model or refer users to experts in the field to handle difficult
problems. In addition, users may wish to consult any readme files associated
with the distribution of this code or any information that may be available on
the CEAM bulletin board to learn of supplements to this manual on turbulence
parameterization. Several reviewers have noted the need to supplement the
information available in this report and the CEAM will attempt to do this as
time permits,
2.5.1 Constant Eddy Coefficients
For this option, constant eddy coefficients are specified in the
vertical direction. These constant values are equal for momentum, mass, and
34
-------
heat transfer (Peter Smith, in review) despite the knowledge that this is not
precisely correct. (Monin and Yaglom 1971, Turner 1973). Given the
approximate nature of the eddy viscosity approach however, this can be a
practical means of making reasonable calculations.
The use of constant eddy coefficients is not normally recommended for
precise calibrations of a hydrodynamic model (e.g., see McCutcheon et al. 1990
for guidance). Such an application, will only crudely approximate
circulation, transport, and mixing in a. stratified flow or a complex flow. In
some screening analyses, however, approximate descriptions may be adequate for
a preliminary investigation of circulation patterns. In addition, it may be
useful to begin calibration of a model with this option to determine if a
reasonable and stable solution is possible. As a result, the ability to
specify a constant eddy coefficient is occasionally useful, but results can be
inaccurate and misleading.
Because the use of constant eddy coefficients results in very
approximate results, it is difficult to estimate the coefficients for various
types of flows in different water bodies. Furthermore, published values can
not be adequately assessed without knowledge of the sensitivity of model
results to eddy coefficients. Unfortunately, this is rarely investigated or
discussed in the available published reports. As a result all guidance on
ranges of values must be used with care.
To assist in selecting eddy coefficients, see McCutcheon et al. (1990)
for values obtained in estuaries, harbors, and coastal areas. Bowie et al.
(1985) offers guidance for the relative order of magnitude of eddy
coefficients that may be useful for selecting tentative values. Chapra and
Reckhow (1983) and Henderson-Sellers (1984) discuss eddy coefficients for
lakes and reservoirs. Typical values for North American large lakes are
(Lynch 1986, Cheng et al. 1976, Csanady 1975):
AH 103 to 105 cm2 s"1 (horizontal eddy viscosity)
Av = 1 to 100 cm2 s"1 (vertical eddy viscosity)
U = 10 cm s"1 (horizontal velocity)
L - 106 to 107 cm (horizontal length scale)
H <= 103 to 10A cm (vertical length scale)
Ap/p 10~3 (relative density difference
over the depth)
Also see Sheng (1986b) for additional guidance on values of eddy coefficients
that may be appropriate when studying lakes.
2.5.2 Munk-Anderson Type Eddy Coefficients
There are a number of different formulas for the vertical eddy
coefficients for momentum, mass, and heat transfer (see McCutcheon et al. 1990
Henderson-Sellers 1982, Blumberg 1986, and Sheng 1983 for a review).
Unfortunately, there is no clear guidance on which of these divergent
formulations are best adapted to certain water bodies. McCutcheon et al.
(1990) indicates that a Munk-Anderson type formulation among others may be
useful for studies of estuaries when calibration may not be possible. Similar
35
-------
guidance for lakes, coastal areas, and harbors is not readily available. When
calibration is possible, HcCutcheon et al. notes that there may be several
alternative forms that are useful, but there has not been sufficient study to
show that any form is significantly better that the Munk- Anderson form used in
this and a number of other hydrodynamics codes (McCutcheon, 1983), As a
result, the Munk-Anderson formulation offers useful flexibility and
consistency with other models.
The Munk-Anderson formulation is based on the observation that eddy
coefficients for stratified flows are a fraction of the eddy coefficients for
non- stratified flows under the same conditions. The ratio of the stratified
to non-stratified coefficients are equal to stability functions [^(Ri) and
<£2(Ri)] of the gradient Richardson number (Richardson 1921, Turner 1973),
expressed as :
Av - A i (Ri)l K* - K^, *Z(RI>; and Dv -= Dvo 02(R1) (50)
where :
Uzr
dp I f 3u }z ( 8v
+ -
9a [ [ da J I da
-1
(51)
e is the dimensionless density and the other terms in the gradient Richardson
number are also as have been described beforehand. The variables Avo, K^,, and
Dvo are the eddy coefficients for momentum, mass, and heat transport under
non-stratified conditions, respectively. For this code, it is assumed that
the eddy coefficients for mass and heat transport are equivalent. This is not
exactly the case and some small discrepancies may arise in the calculations of
salinity and temperature. Also the usual practice is to assume that A^0 = K^
= Dvo (e.g., see McCutcheon 1983) as is done in this code. These variables
are computed from the vertical transport of momentum as:
^2o
H
(52)
where A0 is the length scale assumed to be a linear function of a that
increases with distance above the bottom and below the water surface, with a
peak value at raid-depth, but not exceeding a certain defined fraction of the
local depth.
Equation (52) is limited to describing the generation and dissipation of
turbulence in boundary-layer like shear flows over the bottom. Among other
effects, it does not include the transport of turbulence from different flow
conditions, nor does it include the effects of surface waves. Waves can cause
more intense mixing in the upper layers of deeper flows. In these cases where
the flow is not a boundary-layer type, Equation (52) may be of limited
validity. When flows are not boundary-layer types, it is typically assumed
that the eddy coefficients are constant in the upper layers of flow,
especially down to the thermocline in lakes (see Sheng et al. 1986a,
Henderson-Sellers 1984) or constant below the thermocline (McCormick and
Scavia 1981, McCutcheon 1983) depending on relative depths. For limited
36
-------
additional information about the effects of wind mixing on the mixing
coefficients see Kent and Pritchard (1959) .
As noted above, the stability functions are of diverse forms but the
Munk- Anderson formulations seems to be the best available for use. The
general form of the Munk-Anderson (1948) formula is generally written as
(McCutcheon et al, 1990).
^ - (l+^iRi)"1 and #2 - (l+a-jRi)*2 (53)
where oi, al, az> a^ are empirical coefficients that vary from one water body
to the next (McCutcheon et al . 1990), and seem to be spatially and temporally
variable in a given body of water. See McCutcheon et al. (1990) for a
tabulation of the coefficients that have been used in past studies of selected
estuaries and coastal waters. Henderson-Sellers (1982) and Blumberg (1986)
tabulate a few more values related to studies in lakes. Munk and Anderson
(1948) found for the thermoclrne in the ocean that Equation 53 was best
written as :
^ - (l+10Ri)'1/2 and ^2 - (1+3 . 33Ri)"3/z (54)
Not only the form of the stability function may vary from site to site,,
but different Richardson numbers may also be used for different types of flow
conditions. For example, the formation and deepening of the thermocline in a
relatively shallow basin depends strongly on the relative importance of wind
stress and heat flux at the free surface. In such a case, the following
Richardson number could be used:
KZHZ£72£ S
(55)
U2ru2.(l+e/0 da
where ic is the von Karman constant and u* is the dimensionless friction (or
shear) velocity at the free surface. McCutcheon et al (1990) review a number
of other stability functions as well. These include gross Richardson numbers,
Froude number, the Monin-Obukhov scaling length, the Nyquist Buoyancy
frequency, and Richardson numbers based on shear velocity and average velocity
along with various definitions of linear and nonlinear density gradients.
2.5.3 A Simplified Second-Order Closure Model
The simplified second- order closure model is derived from a complete
Reynolds stress turbulent transport model by assuming a local equilibrium
condition (Sheng 1982, 1983, 1984, 1986a) . In addition to the mean flow
equations, a set of algebraic equations are solved for the second- order
correlations are solved to obtain the stability functions ^ and 2 in terms
of the mean flow variables. These Cartesiaon coordinate equations, when
written in dimensional and tensor forms, are:
37
-------
- Si
- gj
- 2
q
A
(56)
dp
- u
3u,
--1
Sx
2
- 0.75 q -
(57)
2u3'p'
0.45 q p'p'
(58)
The subscripts i, j, and k correspond to the coordinate directions (xi( Xj,
xk) , usually numbered 1, 2, and 3. Therefore, UL' , Uj', and uk' are velocity
fluctuations about the mean velocities uis u j , and uk, respectively. Hence
Equation 56 represents six separate formulas and Equation 57 represents three
separate formulas. Also, p and p' are the mean density and density
fluctuation, respectively; p0 is a reference density, usually taken as the
minimum density; g is the gravitational vector; fi is the vector for the
rotational speed of the earth; £ is the pernutation tensor; q is the total
turbulent micro -length velocity; and A is a turbulent length scale,
The detailed derivation of Equations 56 through 58 can be found in Sheng
et al. (1989b) . Also see Sheng (1983). A graphical comparison of this
formulation versus some of the semi-empirical forms discussed in the last
section is shown in Figure 4.
The length scale, A0, is assumed to be a linear function of the vertical
distance above the bottom or below the free surface. In addition,
stratification is assumed to modify the length scale through the following
empirical relationship:
A = A0 ( 1 + Sx Ri)S2
(59)
38
-------
(a)
Blumberg
Kent & Prilchord
Bowden B Hamilton
Munk B Anderson
74.0
2.5 5.0 7.5 10
.6 -4 -.2 0 .2 .4 j6
R!
Figure 4. (a) Empirical stability functions of vertical turbulent eddy
coefficients from Slumber (1975. (b) Stability functions determined from a
second-order closure model of turbulent transport.
39
-------
where Sx and S2 are arbitrary coefficients.
The length scale A is then adjusted such that the following relationships
are satisfied (Sheng and Chin 1986):
||fl * 0.65 (60)
A < CA H (61)
A < CA Hp
(62)
where CA is usually on the order of 0.1 to 0,2, H is the total depth, and Hp
is the depth of the pycnocline.
The simplified second-order closure model presented above is strictly
valid when the turbulence time scale (A/q) is much less than the mean flow
time scale and when turbulence does not change rapidly over A. It has been
found, however, to be quite useful in simulating vertical flow structures in
estuarine and coastal waters.
Figure 4 illustrates the behavior of the resulting stability functions
as a function of the gradient Richardson number (see Equation 51). This
response is similar to that obtained from the Munk-Anderson type vertical eddy
coefficient functions of the gradient Richardson shown in Figure 4a. Figure 4
was not designed for an exact comparsion, but clearly the trend for the Munk-
Anderson formulation [and similar forms by Kent and Pritchard 1959 and Bowden
Hamilton (1975)] and the simplified second order closure formulation are the
same for Ri>o (stable stratfication; if water is stratified, it is almost
always stably stratified). The curve in Figure 4a by Blumberg, (1975) is from
an alternative closure secheme.
2.6 GRID LAYOUT
The grid system used to describe any computational domain is explained
below.
2.6.1 Staggered Grid
A staggered grid is used in both the horizontal and vertical directions
of the computational domain (see Figure 5), This grid is often referred to as
the "C-grid". In the horizontal directions, a unit cell consists of a f-point
at the center (fj,i)i a U-point to its right (UJ(i) and a V-point at its top
(Vj(i). In the vertical direction, the vertical velocities are computed at
the "full" grid points including the free surface (k=kmax).
40
-------
o U,u
D V, V
A t , W, T, p
(
-B-
-B-
A O A O A O
G G n
)AoAoAoAo
u a rj n
) A o A O A O A O
u G U a
)AOAOAOAO
-B-
n a
) A O A O
-B B
-* cx
u
A W
o T, P
A
O
A
O
A
A
O
A
O
.A.
oc
7777//////7T/
Figure 5, Staggered numerical grid.
-------
Horizontal velocities, temperature, salinity and density are computed at the
"half grid points (half grid spacing below the full points).
2.6.2 Grid Index
Two arrays, each of dimension (J^,^, I,nax) , are used to index the grid
cells. The array NS indicates the condition of the left and right cell
boundaries, whereas the array MS denotes the condition of the top and bottom
cell boundaries (see Figure 6). JUl(I) and JU2(I) indicate the first and the
last water points for computing U along the I-th column. JV1(I) and JV2(I)
denote the second and the second to last water points for V. IU1(I) and
IU2(I) indicate the second and the second to last water points for U along the
I-th row. IV1(J) and IV2(J) denote the first and the last water points for V.
2.7 FLOWCHARTS
Flow charts of the major programs EHSMML, EHSMHC, EHSMEX, EH8MB3, and
EHSMB4 are shown in Figures 7 through 11. The names of the major variables
are listed in Appendix A,
-------
NS
MS
NS
\
V
N.
MS
\ \ \ \
..J-l-..
6
Figure 6. Grid indices NS and MS
43
-------
EHSM3D Mainline Routine - EHSMML
CHSMIN
Input Paramal«rt
Inftlatlze Variable*
i
f
No
EHSMHC
Advanc* Hydtodynamtc
Variable*
EHSMCN
Advance
Dissolved Species
6HSMDT
Output Variable*
Set (STOP
EHSMHR
R»ad Hydredynamlc
Variable* from Disk
Figure 7. Flow chart of the Main Program EHSMML
-------
EHSMHC Subroutine
EHSMDE
Compute Denolty
Barocllnlc Terms
EHSMB3
Advance
Veloeltle*
EHSMED
Compute Vertical
Eddy Coefficients
IHSMEX
Advance
External Variables
Set Velocities
u-U/H
V-V/H
EHSMB4
Compute Inertia! and
Diffusion T»rnll
EHSMSA
Advance
Sallnlly
EHSMTE
Advance
Temperature
Return
EHSMD4
Compute (nerllal and
Dllfuelon Terms
Figure 8, Flow chart of the Hydrodynamic Subroutine EHSMHC
45
-------
EHSMEX Subroulifie
No
EHSMWS
Set
Wind Stresses
T
EHSMm
Set RiverFlow
Vetocily, Salinity
EHSMTD
Set Tidal
Surface Elevations
EHSMFF
CompulB X-Sweep
Forcing Terms
EHSMZS
Advance
External Variables
EHSMFF
Compute Y-Sweep
Forcing Terms
EHSMZS
Advance
External Variables
X - Sweep
Y - Sweep
, N * t
Yes
Figure 9. Flow chart of the External Mode Subroutine EHSMEX
46
-------
X - Sweep of EHSMB3 Subroutine
Comput*
Corlolli Term
1
i
EHSMTB
Compute
Bottom Stre»»
i
Comput* Matrix
Coefficient* (ot
u' (K):K 1,KM
EHSMMt
Invert Matrix tor
u'(KJ:K 1,KM
i
Compute u Velocity
u(K»fromu'(K), U
K 1,KM
Sot Boundary
Conditions on
Open Boundaries
Figure 10. Flow chart of the Internal Mode Subroutine EHSMB3
-------
X - Sweep of EHSMB4 Subroutine
Begin DO
I- 1,IM
J- 1, JM
Compute Local and
Integrated Inertia!
Terms: Ft, Fll
Compute Local and
Integrated Dlf lualon
Termt:FOABC,FHD
Ye»
Figure 11. Flow chart of the Internal Mode Subroutine EHSMB4
48
-------
SECTION 3
USER'S MANUAL
3.1 INTRODUCTION
This section briefly describes the operation of the program, data
requirements, and the form of the output. There are several operational modes
for 3-D and 2-D simulations discussed below. The input data requirements are
briefly reviewed in general terms in Section 3.2 for project managers and
applications experts. Specific data formats are reviewed in detail for
applications experts in Section 3.3. Output information can be printed or
recorded in dimensionless or dimensional forms and a number of different
alternative output data files can be produced that are reviewed in Section 3.4
for applications experts.
As discussed in previous sections, the HYDR03D program can simulate
time-dependent currents in coastal, estuarine, harbor, and lake waters.
Parameters simulated by the program include surface displacement (£")>
vertically integrated velocities (U,V), 3-D velocities (u,v,w), temperature
(T) , salinity (S), density (/?) , and dissolved species concentration (C). The
code can be run as a fully 3-D model, or as a 2-D vertically integrated (x-y)
model. In addition, the code has been designed to simulate 2-D laterally
integrated flows as an x-z model. However, this option has not been fully
tested and it is not presently recommended for use (if such calculations are
necessary, users should contact the CEAM for guidance).
Changing from 3-D mode to 2-D mode or vice versa, requires changing
three parameter statements in the include file, HYDR03D.INC, and a few input
parameters in the input file (*.INP, where * is an arbitrary file name
assigned by the user). The file, HYDR03D.INC, is included with the source
code and has sufficient comments to explain to the user what parameters must
be changed. See Section 3.3 for the parameters that should be specified in
the input file (*.INP).
3.2 Data .Requirements of the Program
Data required to initiate the simulation:
1) l|j|l§:|i||l!|;|fflal|iii§S that defines the horizontal boundaries and bottom
topography of the computational domain. The spatial scales of
physical processes that the model can properly resolve depend on the
49
-------
grid information as well as the governing equations.
2) . The temporal scales of physical processes that
the model can properly resolve depending on the time step information
at the beginning of the simulation.
These include the flow variables as well as the water quality
parameters (salinity, temperature, and concentration of a dissolved
species) ,
Data required to operate the program:
1) ^^^ttlil^l^HlSl^iiiiSi- These include the specification of
fluxes of momentum, heat, and dissolved species at the air-sea
interface as well as the bottom. Alternatively, these conditions
could be given in terms of the state variables instead of their
fluxes .
/* \ *f-y₯^^'W'₯^'₯S*%f»y'*-f.*t-^ mt » 11 .1 * .c " .«_» j- T * j
2) m6gs§|?||!SJ2n|ga^|||0|^|.:|g3p||Si. These include the specification of solid
boundaries, river flows, and open boundary conditions. To run a
successful simulation, valid boundary conditions must be provided at
all times throughout during the simulation.
With these data, the model can be used in a screening mode to develop an idea
about the important processes that control circulation at a site. This is
useful to aid in preliminary investigations and for designing calibration data
collection studies. The simulations, however, must be interpreted with care
until the model is tested with calibration and validation data,
The general procedure for calibrating and testing the adequacy of a
model is as follows:
Determine study objectives,
Define the subset of objectives to be addressed by model studies,
Collect historical data from monitoring or previous studies,
Attempt a preliminary calibration of the model,
Design a calibration data collection study based on the preliminary
calibration,
Simulate conditions during the calibration period and compare to
determine if the preliminary calibration is sufficient (if it is not,
calibrate the model) ,
If the model is calibrated, collect a second independent set of data for
validation,
Validate the model for the limited range of conditions defined by the
calibration and validation data sets (if the model can not be validated
repeat the calibration step and collect more validation data for a second
attempt) , and
Determine uncertainty in the calibrated model simulations by sensitivity
analysis.
See Ambrose et al. (1990) and Ambrose and Martin (1990) for guidance on these
procedures ,
Calibration is the process of changing model parameters until the
50
-------
simulations match measured data. Calibration is necessary because several
critical model parameters such as the Manning roughness coefficient and eddy
coefficients can not be adequately related to conditions at a site or forecast
without simulating the site and changing the parameters as need to match
selected conditions in the modeled domain. This is especially necessary if
precise calculations are required, or if cause and effect relationships must
be defined with some care (and these relationships are sensitive to
hydrodynamics and transport).
Calibration and validation consists of collecting two independent data
sets that define the distribution of currents, salinity, and temperature in
the study domain to be modeled. Data must be collected at enough important
locations and with sufficient frequency to adequately define the phenomena of
interest. Data collection procedures are essentially the same for collecting
calibration and validation data but the data sets must be collected
independently. Occasionally, collection of only one data set for validation
may be possible if there is sufficient data available for a preliminary
calibration of the model.
In practice, some studies only collect one set of data for model
calibration because of resource limitations. These studies are useful but
care must be exercised if model results directly influence resource decisions.
Without validation testing, it is not possible to accurately report
uncertainty in the model simulations (see Chapra and Reckhow 1983) .
When calibrating and validating a model, some care is necessary to
compare model results and data. Both the simulations and data collected
should describe the flow phenomena of interest at the same temporal and
spatial scales. Field data should be collected over long enough periods at a
number of stations and properly averaged. If necessary, simulations should be
averaged for consistent comparisons. Field data collection sites should be
representative of selected parts of the water body being simulated or data
from several sites averaged to provide representative data. In addition to
comparing averages, variances should be compared as well to evaluate the
dynamic response of the model. See McCutcheon et al. (1990) for limited
guidance on statistical testing methods and criteria for simpler models.
At this time, there is only limited guidance on the data necessary to
calibrate a hydrodynamics model. Ambrose and Martin (1990) and Ambrose et al.
(1990) provide guidance on data collection for hydrodynamic model calibration
and validation for estuaries and that information is useful for lakes as well.
For this specific model, Sheng (1983) provides a detailed study of the
Mississippi Sound that should be reviewed to determine the frequency and
spatial locations for calibration and validation sampling. Each site will be
different, however. As a result, it is not possible to foresee all
contingencies and recommend comprehensive data collection procedures.
In practice, the amount of data that should be collected for a specific
study depends on many factors. The primary factors include the resources
available, objectives of the study, flow conditions under study, uncertainties
in the data, and limitations of the model. Generally, data should be
collected over sufficiently long periods of time to define the phenomena that
51
-------
control the water body hydrodynamics. If wind driven episodic events are of
interest, these should be documented with measurements from start to finish.
If spring tides, or other tidal conditions are important, sampling should take
place during these occurrences. If long-term simulations are of interest,
several cruises or data collection studies over important seasons or periods
may be required.
The best practice, is to spend one to two weeks recording and
reviewingstudy objectives, developing a subset of objectives to be addressed
by modeling, and then design a modeling plan. From the modeling plan, it
becomes clearly how much and what kinds of data are necessary to calibrate and
validate the model. Also, the planning clarifies how much time will be
necessary. Tentative sampling plans can be formulated and costs estimated at
this point. The first component of the modeling plan should include a
preliminary calibration with historic data. Following this the sampling plans
and cost estimates should be revised. This is the optimine time to estimate
resource requirements. After some experience is gained, this process can be
streamlined somewhat. However, it does not seem possible to fully estimate
data requirements and costs using a manual like this without knowledge of
specific study objectives and a modeling plan,
3,3 Input Data Description
Most of the data are in a free format and this is denote by (*). Other
data formats are as noted.
At this time, this manual does not provide sufficient guidance on the
ranges of data that parameters can be selected from. For guidance, users are
referred to McCutcheon et al. (1990), Grey (1986), Sheng (1983) and a number
of studies using the model and earlier version of the code (Sheng et al. 1978,
1986, 1989a, 1989b, 1989c; Sheng and Lick 1979, 1980; Sheng 1975, 1980, 1982,
1984, 1986, 1987; Johnson et al. 1989; Smith and Cheng 1989). In addition,
other 2-D and 3-D modeling studies should be consulted as well. Finally,
users may wish to consult with the Center for Exposure Assessment Modeling,
U.S EPA Environmental Research Laboratory, Athens, GA 30506 [(404) 546-3130,
bulletin board (404) 546-3402].
Following are the data formats for each record line:
#1 TITLE CARD: ISTART(I4), TITLE(A64)
START; Start up flag.
- 0 New start, initial flow variables read from input device (file).
1 Restart, initial flow variables and salinity values read from
discfile IR.
-= 2 Special restart, initial flow variables read from discfile IR but
initial salinity values evaluated by interpolation from data at
several stations.
TITLE: A brief description of the run, e.g., Circulation in Green Bay.
#2 PHYSICAL CONSTANTS: XREF, ZREF, UREF, COR, GR, ROO, ROR, TO, TR (*)
XREF: Reference length in lateral direction, usually the maximum
52
-------
ZREF:
UREF:
COR:
GR;
ROD;
ROR:
TO:
TR:
dimension of the basin (cm).
Reference depth, usually the average depth of the basin (cm).
Reference velocity (usually 10 cm/sec for estuaries).
Coriolis acceleration (f 20 sin ^; Q - angular velocity of the
earth; is the latitude).
Gravitational acceleration (usually 981 cm/sec2) ,
Base water density or minimal density in the model domain (g/cm3) .
Reference water density, e.g., density at 20°C and 30 ppt or
maximum density in the model domain (g/cm3) .
Base water temperature, e.g., 1°C or minimum temperature (°C) in
the model domain.
Reference water temperature, e.g., maximum temperature (°C) in the
model domain.
#3 EQUATION FLAGS: IVLCY, ITEMP, ISALT, ICC, IFI, IFA, IFB, IFC, IFD (*)
IVLCY: Velocity flag.
0 Does not compute velocities and other hydrodynamic-variables,
- 1 Computes velocities and other hydrodynamic variables.
ITEMP: Temperature flag.
- 0 Does not compute temperature distribution.
- 1 Computes temperature distribution.
ISALT: Salinity flag.
0 Does not compute salinity.
1 Computes salinity.
ICC: Concentration flag
0 Does not compute dissolved species concentration.
= 1 Computes dissolved species concentration.
IFI: Nonlinear inertia flag for the momentum equations.
= 0 Does not compute nonlinear inertia terms in the equations.
-= 1 Computes nonlinear inertia terms in conservative form with central
differencing scheme.
- 3 Computes nonlinear inertia terms in conservative form with second
upwind differencing scheme.
4 Computes nonlinear inertia terms in conservative form with
combined central and upwind scheme.
IFA: Coefficient for group A of the higher-order lateral
diffusion terms.
= 0 Does not include one group A of the higher-order lateral diffusion
terms.
1 Includes one group A of the higher-order lateral diffusion terms.
IFB : Coefficient for group B of the higher-order lateral
diffusion terms.
IFC : Coefficient for group B of the higher-order lateral
diffusion terms.
IFD : Coefficient for the leading-order lateral diffusion terms
0 Does not include lateral turbulent diffusion.
1 Include the leading-order lateral diffusion terms.
#4 TEMPERATURE PARAMETERS: BVR, SI, S2, PR, PRV, TWE, TWH, FKB, TQO (*)
BVR: Reference turbulent thermal eddy diffusivity (cm2/sec).
S1,S2: Empirical constants used in the simple variable vertical eddy
coefficients.
53
-------
PR:
PRV:
TWE:
TWH:
FKB:
TQO:
#5
CMAX,
IVER:
- 1
- 2
ICON:
- 0
= 1
Turbulent Prandtl Number (typically assigned a value of 1),
Vertical turbulent Prandtl Number (typically assigned a value of
1).
Initial temperature in the epilimnion or upper layer(°C),
Initial temperature in the hypolimnion or upper layer(°C),
Vertical grid index of the initial thermocline location (not in
use at this time)
Initial surface heat flux (cal/cm/cm/sec).
IUBO
IBL
IBR
JBM
JBP
CREF
CMAX
CO:
IC1:
IC2:
JC1:
JC2:
ID1:
ID2:
JD1:
JD2:
CONCENTRATION PARAMETERS: IVER, ICON, IUBO, IBL, IBR, JBM, JBP, CREF,
CO, IC1, IC2, JC1, JC2, ID1, ID2, JD1, JD2 (*)
Vertical diffusion flag.
Explicit vertical diffusion term for water quality equations.
Implicit vertical diffusion term for water quality equations.
Advection flag for water quality equations (similar to IFI).
Does not compute advection terms in the equations.
Computes advection terms in conservative form with central
differencing,
- 3 Computes advection terms in conservative form with second upwind
differencing scheme.
= 4 Computes advection terms in conservative form with combined
central and upwind differencing scheme.
Bottom orbital velocity flag (enter value of 0; not used at this
time 0).
Concentration computation does not have to he performed
for the entire computational domain. Instead, it can be
done for a window that covers an area from I~IBL to
I-IBR and from J=JBM to J=JBT, initially.
Reference species concentration (units determined by user).
Maximum concentration allowed by the code (The run stops if Cmax
is exceeded)
Initial concentration (in units determined by user).
Initial concentration field may be specified to
be zero everywhere in the computational domain except
within two windows: the first one covers an area
from I=IC1 to I-IC2 and from J-JCl to J-JC2.
the second one from I=ID1 to I=ID2 and J=JDl to J=JD2.
#6 TURBULENCE PARAMETERS: IEXP, IAV, AVR, AVI, AV2, AVM, AMR(*)
IEXP: Vertical eddy coefficient flag (see EHSMED.FOR and EHSMEZ.FOR for
details).
= 0 Constant eddy coefficient. Must also set ISPAC(9) = 0,
The following options are used with variable eddy coefficients, i.e., when
ISPAC(9) is nonzero. See record #12 following.
-1 Richardson-number dependent on eddy coefficients with length scale
linearly increasing from the bottom and surface.
= 2 Richardson-number dependent eddy coefficients with length scale
linearly increasing from the bottom.
= -3 Eddy coefficients determined from simplified second-order
54
-------
turbulence closure model,
IAV; Reference vertical eddy viscosity flag.
0 Input parameter AVR is used as reference eddy viscosity.
- 1 Reference eddy viscosity is computed form AV1+TXY*AV2, where TXY
is the total wind stress and AVI and AV2 are input parameters,
AVR: Reference vertical eddy viscosity (cmz/sec),
AVI: Background vertical eddy viscosity when wind is zero (cma/sec).
AV2: If IAV 1, urxstratifled vertical eddy viscosity is computed from
AV1+TXY*AV2.
AVM; Minimum allowable vertical eddy coefficient (cm2/sec)
AHR: Reference lateral turbulent eddy viscosity (cmz/sec).
#6A MORE TURBULENCE PARAMETERS: FM1, FM2, ZTOP, SLMIN, QQMIN (*)
FM1: Empirical constant used in Richardson-number dependent eddy-
coefficient formula.
FM2: Empirical constant used in Richardson-number-dependent eddy-
coefficient formula.
ZTOP: Distance between the top of the computational domain and the free
surface (cm).
SLMIN: Minimum scale length (cm).
QQMIN: Minimum total turbulent velocity (cm/sec)).
#6B MORE TURBULENCE PARAMETERS: QCUT, ICUT, GAMAX, GBMAX, FZS, KSMALL (*)
QCUT: (Not used for now).
ICUT: Eddy coefficient parameter.
= 0 Cutoff not operating.
* 1 Eddy coefficients below a sharp density gradient are not allowed
to exceed that at the sharp gradient.
GAMAX: Maximum vertical eddy viscosity (for u, v variables),(cm2/sec).
GBMAX: Maximum vertical eddy diffusivity (for B, t. variables), (cm2/s
ec) .
FZS: The turbulence length scale is not allowed to exceed the product
of FZS and the depth.
KSMALL: A non-zero value of KSMALL implements a routine to the vertical
eddy coefficients (a value of 1 is recommended).
#1 WIND PARAMETERS: IWIND, TAUX, TAUY (*)
IWIND : Wind stress flag.
- 0 Uniform wind stress specified by TAUX and TUAY.
= 1 Variable wind stress read from file unit IR4.
TAUX: Wind stress at the air-sea interface in the x-direction.
TAUY: Wind stress at the air-sea interface in the y-direction.
#8 VERTICAL BOUNDARY CONDITION PARAMETERS: ISMALL, ISF, ISIE, IBTM, ITB,
HADD, HMIN, ZREFBN, CTB, BZ1, HI, H2 (*)
ISMALL: Small amplitude flag.
- 0 Small amplitude assumption is invoked. Surface elevation is not
added to the depth to obtain the total depth.
- 1 Small amplitude assumption is not invoked. Surface elevation is
included in the total depth.
IBTM: Bottom topography flag.
= 0 Depth changes linearly from HI along the western boundary to H2
55
-------
along the eastern boundary.
1 Depth changes linearly from HI along the southern boundary to H2
along the northern boundary.
- 2 Depth deck (series of records) is read in at the end of this input
stream.
- 3 Components of depth (HU,HV,HS) are read from discfile (unit 11).
ITB : Bottom friction flag.
1 Linear stress law with no slip condition is employed.
> 3 Quadratic stress law is employed.
HADD: Constant datum added to the depth at all locations.
HMIN: Minimum depth,
0. No adjustment on the depth data.
> 0. Depth cannot be less than HMIN.
ZREFBN: Reference height above the bottom (ZREFB - ZREFN * BZ1/ZREF).
CTB: Constant bottom friction coefficient (.004 to 0.4). Also see
JSPAC(2) on record #12. CTB is only used if a constant friction
or drag coefficient is requested by setting JSPAC(2) = 0.
BZl : Constant bottom roughness height (0.1 cm to 0.5 cm).
HI: Depth along one boundary (cm).
H2: Depth along the opposing boundary (cm).
#8A ZREFTN, TZl, SSSO (*)
ZREFTN; Reference height at the top (cm).
TZl : Constant surface roughness height (cm).
SSSO: Initial uniform surface elevation (dimensionless divided by Sr =
fu). sr = fUrXr/g.
#9 LATERAL BOUNDARY CONDITION FLAGS: ITIDE, IOPEN, JWIND, IJLINE (*)
ITIDE : Tidal forcing flag.
=0 No tidal forcing.
1 With tidal forcing and constituent tide boundary condition (This
option is currently inactive).
2 With tidal forcing and tabular tide boundary condition.
3 Surface water elevations are prescribed.
IOFEN : Open boundary flag. A 4-digit number that indicates whether there
are open boundaries along the west-south-east-north sides of the
computational domain. Zero indicates no open boundary and 1
indicates open boundary. For example, 1010 means open boundaries
on west and east.
JWIND : Open boundary flag in case of wind forcing only.
IJLINE: Number of open boundary lines along which tidal forcing
information is to be specified,
#10 FOR EACH IJLINE, IF IJLINE.GT.O, READ:IJGAGE, IJDIR, IJROW, IJSTRT,
IJEND (*)
IJGAGE: Gage number. Not important.
IJDIR: Direction of the line segments along which tidal data are
prescribed. In case of constituent tides, IJDIR = 1 indicates the
x-direction while IJDIR 2 indicates the y-direction. In case of
tabular tides, IJDIR = (1,2,3,4) indicates
(west,south,east,north),
IJROW : Row/column index of the line segment,
IJSTRT: Grid index of the starting point on IJROW.
56
-------
IJEND : Grid index of the ending point on IJROW.
#11 (ISPAC(I),1-1,10) (*)
(1) : Disk output flag,
- 0 No output to diskfile (unit IW).
> 0 Every ISPAC(l)-th time-step, output are written to diskfile.
(2) : Smoothing flag,
- 0 No smoothing is applied.
> 0 Every ISPAC(2)-th time-step, apply smoothing to salinity.
(3) : Open Boundary flag for elevation.
- 0 Prescribed surface elevation on open boundary.
> 0 Surface elevation has zero slope on open boundary.
< 0 Surface elevation satisfies radiation condition on open boundary.
(4) ; Open Boundary flag for mass flux.
- 0 Compute tangential mass flux from equations of motion.
< 0 Tangential mass flux - 0.
> 0 Tangential mass flux has zero slope in normal direction.
(5) : Open Boundary flag for advection/diffusion terms.
«= 0 Computes advection/diffusion terms on open boundary.
> 0 Set advection/diffusion terms to zero on open boundary.
(6) : 2-D flag.
= 0 Performs the x-sweep and the y-sweep.
> 0 Performs the x-sweep only.
< 0 Performs the y-sweep only.
(7) : Basin geometry flag.
- 0 Grid indices NS, HS, JU1, JU2, JV1, JV2, IU1, IU2, IV1, IV2, are
determined from the depth arrays in EHSMI*.FOR routines.
= 1 Previously determined grid indices are read from discfile (unit
12).
(8) : Smoother flag.
=0 No smoother is applied.
> 0 Every ISPAC(8)-th step, smoothing is applied to the velocity
field.
(9) : Vertical eddy coefficient flag.
= 0 Constant vertical eddy coefficients,
- 1 Variable vertical eddy coefficients computed from EHSMED and
EHSMEZ routines.
(10): Residual current flag.
0 Does not compute Eulerian residual currents,
1 Computes Eulerian residual currents. See subroutine EHSMRS.
#12 (JSPAC(I), 1-1, 10) (*)
(1) : Dimensionality flag.
= 0 All output in dimensionless units,
- 1 All output in c.g.s. unit.(Does not work at this time),
(2) : Bottom friction coefficient flag,
- 0 Constant bottom friction coefficient is specified by CTB, See
Record #8.
= 1 Variable bottom friction coefficient is computed in EHSMEX and
EHSMBC routines based on the law of the wall.
(3) : Coriolis acceleration terms.
1 No Coriolis flag,
57
-------
- 0 Coriolis acceleration terms are evaluated.
(4) : Dummy flag used to check steady state termination.
(5) : Open boundary salinity flag.
- 0 Prescribed salinity is used along open boundary during inflow.
- 1 Prescribed salinity is used in a 1-D advection equation to obtain
the salinity along open boundary during inflow.
(6) Dummy flag.
(7) Dummy flag.
(8) Bottom friction flag for the 2-D vertically-integrated version
(KM-1).
- 0 Explicit bottom friction.
1 Implicit bottom friction.
(9): Open boundary flag for salinity and temperature (To be used later
for thermally stratified flow. > 1 will include salinity time
varying data; - 2 Freeze open boundary salinity to initial
values).
(10): Include initial values for salinity at interior nodes (stations).
#13 (RSPAC(I), 1=1, 10) (*)
(1): Manning's n in association with other parameters in c.g.s. units.
(2): Dummy parameter.
(3): An small number used in checking the convergence to steady-state
(0.0001 is recommended).
(4) : An small number used in checking the convergence to steady-state
(0.0001 is recommended).
(5) ; Dummy parameter.
(6) ; Dummy parameter.
(7): Depth below which the bottom friction coefficient follows a ramp
function (see EHSMTB.FOR).
(8); When depth falls below RSPAC(7), the bottom friction coefficient
is
linearly interpolated between the one computed in EMSKTB and
RSPAC(B).
(9): Coefficient for the spatial smoother (0.25 is recommended).
(10): Coefficient for the curvature check of the spatial smoother (4 is
recommended, See Sheng (1983)), p. 258.
#14 TIME-STEPPING PARAMETERS: ISTEP, 111, IT2, ITS, DELT, DELTMIN, DELTMAX,
EPSILON, BUFAC, WTS, WTU, WTV (*)
ISTEP: Time-stepping flag.
0 Constant time-step is used in the time-integration of the finite-
difference equations.
- 1 Dynamic time-stepping is used. At the end of each time-step, the
total weighted maximum rate of change of the major variables is
compared with EPSILON. If the rate of change is less than
EPSILON, the time-step is allowed to increase by 10% to 20%.
Otherwise, the time-step is cut back proportional to the ratio of
EPSILON and the rate of change.
IT1: Initial time index. NOTE IT1 - 1 when beginning a simulation
(ISTART -0). If a simulation is being restarted (ISTART - 1 or
2), IT1 should be set equal to the value IT2 in the previous run
plus 1.
58
-------
IT2: Final time index,
ITS: Ratio of the internal time-step to the external time- step.
BELT: Initial time-step (seconds). For constant time stepping this is
the time step used.
DELTMIN; Minimum allowable time-step when dynamic time-stepping is used.
DELTMAX: Maximum allowable time-step when dynamic time-stepping i s used.
EPSILON: Maximum allowable rate of change of major variables (5% or 0.05 is
recommended).
BUFAC; When rate of change exceeds BUFAC*EPSILON, the run stops.
WTS: Weighting factor for surface elevation when computing EPSILON (0.1
is recommended)
WTU: Weighting factor for surface u-velocity when computing EPSILON (1.
is recommended).
WTV: Weighting factor for surface v-velocity when computing EPSILON (1.
is recommended).
#15 PRINTOUT PARAMETERS: ITEST, IP1, IP2, IP3, IPU, IPW, IPA, IPB, ID, JPA,
JPB, JD, KPA, KPB, KD (*)
ITEST : Testing/debugging flag.
= 0 Operational run with minimal output information to printer.
= 1 Test run with extra output to the printer.
- 3 Creates time history file (*.SUV) that contains major variables at
selecting stations and vertical levels. See Cards #20 and #21.
Time index interval for brief printout.
Time index interval for total printout.
Time index interval for printout within each internal step
IP1
IP2
IF3
(unused]
IPU
IPW
IPA
IPB
ID:
JPA
JPB:
JD:
KPA
KPB
KD:
IPU
0 turns off printing and
0 turns off printing and
Horizontal velocity printout flag.
IPU = 1 activates printing.
Vertical velocity printout flag. IPW
IPU = 1 activates printing.
First index for x-direction printout.
First index for y-direction printout.
The printout does not have to cover the entire computational
domain.
Instead, the printout goes from I=IPA to I=IPB, every ID-th
spacing in the x-direction, print information in y-direction
from J-=JPA to J-JPB, every JD-th spacing in the y-direction.
First index for z-direction printout.
Last index for z-direction printout.
Time index interval for z-direction printout.
#16 PRINTOUT FORMAT FLAGS: IGI, IGH, IGT, IGS, IGU, IGW, IGC, IGQ, IGL,
IGR, IGRI, IGTB (*)
IGI : Printout format flag for initial data.
- 0 Procedures digital printout (via EHSMWR routine) of initial data.
= 1 Procedures simple contour plot (via EHSMGR routine) of initial
data.
IGH : Printout flag for depth arrays.
=-1 Does not print depth arrays.
= 0 Procedures digital printout of depth arrays.
59
-------
- 1
IGT
IGS
IGU
IGW
IGC
IGQ
IGL
IGR
IGRI
IGTB
Procedures graphical printout of depth arrays.
Printout format flag for temperature variables (Same as IGI).
Printout format flag for surface elevation (Same as IGI),
Printout format flag for mass flux and velocity and integrated
velocity in the horizontal direction (Same as IGI).
Printout format flag for and velocity and integrated velocity in
the vertical direction (Same as IGI).
Printout format flag for concentration variables (Same as IGI).
Printout format flag for turbulent velocity. Set to <0 or >1 to
turn off printing.
Printout format flag for turbulent scale. Same as IGQ.
Printout format flag for density field.
Printout format flag for Richardson number.
Printout format flag for bottom stress.
ICI
IWC
ICONC).
#17 DISKFILE INFO: IRD, IV, IWR, ICI, IWC, ICONC, IWS, IREAD, IR4 (*)
IRD; Unit number of input diskfile for storing arrays of major flow
variables.
IW; Unit number of output diskfile for storing arrays of major flow
variables.
IWR : Output flag for arrays of minor flow variables.
- 0 Does not writes to diskfile unit IW.
= 1 Writes to diskfile unit IW every ISPAC(l)-th steps. In case of
dynamic time-stepping, writes to diskfile every time break
specified by TBRK on data set (record) #19.
Concentration input flag.
0 Does not read concentration field from diskfile unit ICONC.
1 Reads concentration field from diskfile (unit ICONC).
Concentration output flag.
0 Does not write concentration field to the output diskfile (unit
- 1 Writes concentration field to output diskfile unit ICONC.
ICONC : Unit number of diskfile for storing concentration variables.
IWS ; Unit number of diskfile for storing Eulerian residual variables.
IREAD : If IVLCY = 0 and ICC > 0, reads flow variables from unit IR
diskfile every IREAD-th steps.
IR4 ; Unit number of the wind stress file.
#18 3 MAJOR I/O FILENAME: FNAME (6A4)
FNAME : A six-element vector specifying the names of (1) the input
diskfile for flow variables, (2) the output diskfile for flow
variables, and (3) the diskfile for concentration variables,
#19 TIME BREAKS FOR STORING MAJOR OUTPUT ARRAYS: (TBRK(I),1=1,10) (*)
TBRK: Time breaks (in hours) at which major flow output will be stored
on file unit IW. Only used when ISTEP = 1.
> 0. Hours at which flow data are written to file unit IW.
< 0. Run is stopped after flow data at ABS(TBRK) hour are dumped to
file unit IW.
#20 TIMEFILE (UNIT 18) GAGE STATIONS: NSTA, NRANGE, NFREQ (*)
NSTA: Number of stations (not to exceed NSTATS in HYDR03D.INC) where
60
-------
major flow variables are stored on a timefile (*.SUV). No data
stored on timefile if NSTA - 0.
NRANGE: Dummy flag.
NFREQ : Flow data are stored on timefile every NFREQ steps.
#21 IF(NSTA.GT.O) READ : IST(K), JST(K) , KST(K), STATID(K)(314, A48)
1ST : x-Grid index of the timefile gage stations. Proceeds to read river
info cards (records) when a zero 1ST is detected.
JST : y-Grid index of the timefile stations.
KST ; z-Grid index of the timefile stations.
STATID: A 48-character title for each of the timefile stations.
#22 RIVER INFO: NRIVER (*)
NRIVER: Total number of river stations in the computational domain, not to
exceed NRIVRS specified in HYDR03D.INC.
#22A FOR EACH NRIVER.GT.O, READ: IRIVER, JRIVER, LRIVER, URIVER, VRIVER (*)
IRIVER: Grid index (of a river station) in the x-direction.
JRIVER: Grid index (of a river station) in the y-direction.
LRIVER: Alignment index of a river station.
- 1 River flows in the x-direction.
- 2 River flows in the y-direction.
< 0 Read time-varying river flow rate from a disc file.
URIVER: Volumetric flow rates in ft3/sec for rivers with LRIVER=1.
VRIVER: Volumetric flow rates in ft3/sec for rivers with LRIVER=2.
#23 INITIAL VERTICAL PROFILES OF SALINITY, TEMPERATURE, AND CONCENTRATION
ALONG THE OPEN (WEST-SOUTH-EAST-NORTH) BOUNDARIES:
IF (ISALT.NE.O) READ (SABW(K),SABS(K),SABE(K),SABN(K),K=1,KM) (*)
(Note that ISALT Is defined in Record # previously)
#23A,B,C IF ISALT.NE.O, READ NUMBER OF INITIAL SALINITY STATIONS: NISS
IF (NISS.GT.O) READ THE FOLLOWING CARDS FOR N=1,NISS
ISS(N), JSS(N), NDEPTH(N), TDEPTH(N)
(ZDEPTH,(K,N), K-l, NDEPTH(N))
(ZSAL(K,N), K-l, NDEPTH(N))
ZSEDI: Concentration data at each depth level of a station
ISS : I-grid index of the salinity station.
JSS : J-grid index of the salinity station.
NDEPTH: Number of depth levels at which salinity data are to be specified.
TDEPTH: Total water depth at the salinity station.
ZDEPTH: Depth measured from the water surface at each depth level of a
station. The units of TDEPTH and ZDEPTH must be the same (c.g.s.
or f.p.s.).
ZSAL: Salinity data at each depth level of a station.#24INITIAL VERTICAL
PROFILES OF TEMPERATURE ALONG THE OPEN (WEST-SOUTH-EAST-NORTH) BOUNDARIES:
IF (ITEMP.NE.O) READ (TBtf(K), TBS(K), TBE(K), TBN(K), K-l,KM
#24A,B,C IF ITEMP.NE.O, READ NUMBER OF INITIAL TEMPERATURE: NITT
IF (NITT.GT.O) READ THE FOLLOWING CARDS FOR N=l, NITT
ISST(N), JSST(N), NDEPTT(N), TDEPTT(N)
(ZDEPTT(K,N), K-l, NDEPTT(N))
61
-------
(ZTEM(K.N), K-l, NDEPTT(N))
ISST: I-Grid index of the temperature station.
JSST: J-grid index of the temperature station.
NDEPTT: Number of depth levels at which temperature data are to be
specified.
TDEPTT; Total water depth at the temperature station.
ZDEPTT: Depth measured from the water surface at each depth level of a
station. The units of TDEPTT and ZDEPTT must be the same (c.g.s.
or f.p.s.)
ZTEM: Temperature data at each depth level of a station.
#25 INITIAL VERTICAL PROFILES OF CONCENTRATION ALONG THE OPEN (WEST-SOUTH-
EAST-NORTH) BOUNDARIES
IF (ICONC.NE.O) READ (CBW(k), CBS(k), CBE(k), CBN(k), K-l, KM)
#25A,B,C IF ICONC.NE.O, READ NUMBER OF INITIAL CONCENTRATIONS; NISSS
IF (NISSS.GT.O) READ THE FOLLOWING CARDS FOR N=l, NISSS
ISSS(N), JSSS(N), NDEPTS(N), TDEPTS(N),
(ZDEPTHS(K,N), K-l, NDEPTHS(N))
(ZSEDI(K,N), K-l, NDEPTHS(N))
ISSS; I-grid index of the concentration station
JSSS: J-grid index of the concentration station
NDEPTHS: Number of depth levels at which concentration data are to be
specified
TDEPTHS: Total water depth at the concentration stations
ZDEPTHS: Depth measured from the water surface at each depth level of a
station. The units of TDEPTHS and ZDEPTHS must be the same
(c.g.s, or f,p.s.)
ZSAL: Salinity data at each depth level of a station.
#26 TIDAL BOUNDARY CONDITION PARAMETERS: NCG, NCONST, XYEAR, XMONTH, XDAY,
XHOUR (*)
NCG : When ITIDE = 1, total number of tidal gages (not to exceed NTIDE
specified in HYDR03D.INC) where constituent tide information are
to be specified. Since the constituent tide option does not work,
there is no need to specify the following 4 groups of cards (#25
through #28). If ITIDE = 2, tabular tide information are to be
read (#29 through #31).
NCONST: Total number of tidal constituents to be used (not to exceed NCNST
in HYDRO3D.INC).
XYEAR : The year at the beginning of the constituent tidal record.
XMONTH: The month at the beginning of the constituent tidal record.
XDAY: The day at the beginning of the constituent tidal record.
XHOUR : The hour at the beginning of the constituent tidal record.
#26A INDEX NUMBER OF TIDAL CONSTITUENTS: (NCST(NCI),1=1.NCONST) (*)
NCST: Index number of tidal constituents.
FOR EACH (J=1,NCG), READ #26B,27,28:
#26B KNGAGE(J),HO(J),XLONG(J) (*)
#27 TIDAL AMPLITUDES: (AMP(I,KNGAGE(J)),1-1,NCONST) (*)
62
-------
#28 TIDAL PHASES: (XKAPPA(I,KNGAGE(J)),1-1,NCONST) (*)
KNGAGE: Tidal gage number.
HO: Gage datum or mean sea level relative to depth datum,
XLONG : Longitude of tidal gage.
AMP : Constituent tidal amplitudes, in the same order as chosen by NCST,
XKAPPA: Constituent local epochs. Set XLONG
used.
0 if Greenwich epochs are
#29
1
#30
#31
(*)
TP:
J:
TABULAR TIDE DATA (PERIOD): (TP(I),1-1,NCONST) (*)
Tabular tidal periods along four boundaries (in seconds).
(AMP, PHASE)ON(WEST AND EAST): J, NC, AMPW, PHW, CAW, AMPE, PHE, CAE
(AMP, PHASE)ON(SOUTH AND NORTH): I, NC, AMPS, PHS, CAS, AMPN, PHN, CAN
If J < 0, exit
y-grid along the western and eastern boundaries.
the read loop.
NC: Tidal constituent index. Must be less or equal to NCONST.
AMPW Tidal amplitude at certain J along the western boundary,
PHW Tidal phase at certain J along the western boundary.
CAW Constant tidal amplitude added to AMPW.
AMPE Tidal amplitude at certain J along the eastern boundary.
PHE Tidal phase at certain J along the eastern boundary.
CAE Constant tidal amplitude added to AMPE.
I: X-grid index along the southern and northern boundaries. If I<0,
exit the read loop.
NC: Tidal constituent index. Must be less than or equal to NCONST,
AMPS: Tidal amplitude at certain I along the southern boundary.
PHS ; Tidal phase at certain I along the southern boundary.
CAS : Constant tidal amplitude added to AMPS.
AMPN: Tidal amplitude at certain I along the northern boundary.
PHN : Tidal phase at certain I along the northern boundary,
CAN : Constant tidal amplitude added to AMPN.
#32 NUMBER OF THIN-WALL BARRIER: NEAR (*)
NEAR: Total number of thin-wall barriers.
#32A FOR EACH NBAR.GT.O. READ: IJBDIR(I), IJBROW(I), IJBSTR(I), IJBEND(I)
(*)
IJBDIR; Direction of the thin-wall barrier
= 1 The thin-wall barrier is along the x-direction
- 2 The thin-wall barrier is along the y-direction
Grid index of the row/column where the thin wall barrier is
IJBROW
located.
IJBSTR
IJBEND
Starting grid index of the barrier along IJBROW
Ending grid index of the barrier along IJBROW
#33 LATERAL GRID MAPPING: IGRID, XMAP, ALREF, ALYREF (*).
IGRID : Horizontal grid index,
- 0 Uniform horizontal grid. Skips cards (records) #34 through #36.
1 User-specific non-uniform horizontal grid. Grid information read
from unit 14 file.
= 2 Exponentially stretched horizontal grid.
XMAP: Mapping ratio of the physical domain and the computational domain.
63
-------
ALREF : Reference length in the x-direction of the computational domain.
ALYREF: Reference length in the y-direction of the computational domain.
#34 VARIABLE GRID MAPPING IN X DIRECTION: NCR, ALPHA1 (*)
NRG : Total number of mapping regions in the x-direction.
ALPHAl: Counter index of the first cell (usually 1).
#34A FOR EACH NRG, READ VARIABLE GRID MAPPING IN X DIRECTION: LPR, A, B, C
(*)
LPR : Total number of cells in a grid region.
A: Coefficient of the coordinate-stretching equation:
X - A + B (ALPHA) ** C
B: Coefficient of the coordinate-stretching equation.
C: Exponent of the coordinate-stretching equation.
#35 VARIABLE GRID MAPPING Y DIRECTION: NRG, ALPHAl (*)
NRG ; Total number of mapping regions in the y-direction.
ALPHAl: Counter index of the first cell (usually 1).
#36 FOR EACH NRG,READ VARIABLE GRID MAPPING IN Y DIRECTION: LPR, A, B, C
(*)
LPR : Total number of cells in a grid region.
A: Coefficient of the coordinate-stretching equation :
Y - A + B* (ALPHA) ** C
B: Coefficient of the coordinate-stretching equation.
C: Exponent of the coordinate-stretching equation.
#37 IF IBTM=2, READ BATHYMETRY DECK: ((HS(J, I), 1=2, IM), J=2, JM) (*)
64
-------
3.4 MODEL OUTPUT
This section discusses the creation of output files and the form of
information available to the user. During the simulation, certain files will
be used as input, some files will be created as temporary scratch files, and
there will be some output files created. Manipulation of unit numbers and the
names of these files are controlled in the input file (*.INP, reviewed in the
last section) and/or the run file (RUN.COM). There are options to control the
dimensions of the output data and these are also covered in the previous
section. For some additional detail, the reader should also refer to Section
3.3 and Section 5.
In the HYDR03D code, there are two options to control the dimensions of
the output data. One option records the output in dimensionless form and the
other records the information in a dimensional form (i.e, with physical units
of measurement). The option to record the output in dimensional units
involves a separate program that reads the standard output in dimensionless
units, and converts to the final dimensional form. The conversion is
controlled by flags in the input data file.
The HYDR03D program reads the input files and produce several output
files. These output files are controlled and produced by flags defined in the
input file (*.INP) and by the proper files assigned in the run file (RUN.COM).
To aid in bookkeeping and cataloging of results from different simulations,
the output files created by the model use the name of the input data set with
a unique extension as follows:
65
-------
Created for
File name all Runs
Purpose
*.OUT
yes
*.suv
*.DAT
yes
GRID.PAR yes
CONG.OUT
WIND.INP
*.IWX yes
*.TEMP
*.SAL
*,VEL
BOUND.DAT
*,RES.
yes
Contains an "echo" of the input data including Initial
and boundary conditions, some computed results, and any
error messages. Values are either in graphic or
numeric form depending on flags assigned in the input
file (*.INP, see Section 3.3),
Contains time-dependent information calculated at
assigned stations where the outputs are desired.
Needed when ISTART - 1 or 2 and contains initial
computations that can be used in subsequent simulation
runs.
Contains bathymetry information needed for graphic
presentation.
Contains output of dissolved species concentration and
will be created if applicable.
Contains information regarding wind shear stress
calculated in advance using WINDSHEAR.FOR and wind data
(from file WIND.DAT).
Stores the total number of runs applied to a specific
problem.
Created in during the simulation of thermally
stratified flow and contains temperature data required
for graphic presentations.
Contains output of salinity computations.
Contains general output information to be used In
vector and contour plotting programs.
Contains water surface elevation data prescribed at the
open boundaries.
Stores Eulerian residual velocity arrays.
Note that * denotes a user defined name that will be selected by the program
to match the input data file name.
Files not created in all HYDR03D runs are created only for specialized
simulations. For a detailed explanation of these files the reader is referred
to Section 5.6 of this manual.
66
-------
SECTION 4
CASE STUDIES
In this section, several applications of the model are presented for the
purpose of illustrating the utility of the program and the feasibility of its
use in several diverse settings. Although limited field data are available
for direct comparison with model results, much can be learned about the model
performance by analyzing the general and specific results from site specific
simulations. The sites chosen for these case studies include, Suisun Bay of
San Francisco Bay, Charlotte Harbor in Florida, Green Bay of Lake Michigan,
Prince William Sound, and Mississippi Sound. In addition, the case studies
begin with a comparison with a simple analytical solution that illustrates the
general validity of the program. Appendix D is a related type of case study
that is used to illustrates the structure and format of the input and output
data sets. The study described in Appendix D is a simulation of wind driven
currents in an hypothetical enclosed basin.
The specific objectives of this Section are to:
Illustrate the feasibility of using the program in several important types
of water bodies,
Describe the results that may be obtained from screening-level studies,
Demonstrate the options available,
Document the validity of the program, and
Show how the model has been calibrated in at least one study.
To illustrate the feasibility of using this model in diverse bodies of water
we have selected studies from Prince William Sound, Alaska (site of the March
24, 1989 Oil Spill), Suisun Bay of the San Francisco Bay, Mississippi Sound
and adjacent deep waters of the Gulf of Mexico, the partially mixed Charlotte
Harbor on the west coast of Florida, and Green Bay off of Lake Michigan.
Prince William Sound is relatively deep, though not as deep as the waters of
the Gulf of Mexico off the Mississippi Sound. Tide ranges in Prince William
Sound are typically 3 to 5 m (12 to 15 feet). Use of the model in a screening
mode for the emergency response to the EXXON Valdez spill was a primary
consideration in the selection of Prince William Sound as a case study. The
complex bathymetry and the unique influence of San Francisco Bay on many
aspects of the American economy and culture are the appealing attributes of
the study in Suisun Bay. Charlotte Harbor is a partially mixed, shallow
estuary of a classical type (see Ambrose and Martin 1990). The Harbor is
significantly influenced by freshwater flow and has moderate tides of a few
feet. Both San Francisco Bay and eastern estuaries like Charlotte Harbor
(i.e., Tampa Bay and Sarasota Bay) are expected to be important in the U.S.
EPA National Estuaries Studies. The Green Bay study explores the effects of
67
-------
wind and river flow on circulation in a seriously contaminated part of the
U.S.-Canadian Great Lakes where a number of studies are underway.
In addition to the case studies selected, there are a number of other
studies that demonstrate the adaptability of the model to many different
sites, including other critically important sites. These include recent work
by Smith and Cheng (1989) in San Pablo Bay adjacent to Suisun Bay. Johnson et
al. (1989) used a significantly modified z-grid model in Chesapeake Bay.
Sheng et al. (1978) made a realistic application of an early version the model
to Lake Erie that produced very satisfactory agreement with measured data.
Curvilinear versions of the model have also been applied to the James
River Estuary (Sheng et al. 1989a), in Lake Okeechobee (1989c) where a version
of this o--stretched model is also being implemented for a comprehensive test,
and other sites by the Corps of Engineers (Los Angeles-Long Beach Harbor and
Humboltd Bay in California).
The validity of the code is supported by the comparison with a simple
analytical solution in Section 4.1, as well as the other case studies where
various options have been explored. In addition, Appendix D illustrates that
reasonable results are possible for shallow water conditions. However, these
are not the primary studies by which the validity of the code has been
assessed. Other comparisons with analytical solutions and laboratory data are
presented in Sheng and Lick (1980) and Sheng (1983).
Studies by Smith and Cheng (1989), including others to be reported in
the summer of 1990; studies reported by Zakikhani et al. (1989); and studies
by the Corps of Engineers (e.g., see Johnson et al. 1989), are also
noteworthy. These studies, although influenced by the primary architect of
the model, indicate that the code is sound and useful enough for other
investigators to implement. It is hoped that this documentation will
accelerate the use of the program and others like it by other investigators
who need to understand effects of circulation in water-quality studies.
There are several case studies included in this section that are not
complete calibrations of the program. The study in Prince William Sound was
purposely designed as a feasibility investigation and it is reported as such.
The study of Green Bay is in an early stage of calibration where the model has
been calibrated and checked out with historic data, as recommended earlier in
this manual. The studies of Suisun Bay and Charlotte Harbor were selected to
demonstrate the use of several option and the general implementation of the
program. These studies were examples from a workshop found in Sheng et al.
(1986) and the reviewers point out, as do the written sections, that these
studies do not provide definitive conclusions about the circulation in Suisun
Bay and Charlotte Harbor. The same can be said of the feasibility study in
Prince William Sound and the preliminary calibration in Green Bay. However,
the case study for Mississippi Sound does represent an adequate calibration of
the model and is useful for that reason. Both project mangers and
applications experts should benefit for the brief review given for the
Mississippi Sound calibration and applications experts may wish to examine the
68
-------
details of the application. Review of the case study in this section should
provide an adequate idea of the data requirements and the intensity of the
calculations.
4.1 COMPARISON WITH A ONE-DIMENSIONAL ANALYTICAL SOLUTION
In this section a comparison was made between the applications of
HYDR03D and a one -dimensional analytical solution for standing waves in
rectangular basins. A square basin with sides of 50 km in length and a depth
of 10 meters deep, and with one open boundary along y-axis was used. The open
boundary was subjected to a sinusoidally oscillating wave forcing with an
amplitude of f0 - 10 cm and a wave frequency w = 2 w/12.
The linearized governing equation for a one dimensional standing wave in
a rectangular basin, assuming that the flow is incompressible and inviscid,
and the water depth is shallow compared to wave amplitude, is given as (Lynch
and Gray, 1978) :
i£ + h - o
5t ** (63)
flu ar n
-at + & a x °
The boundary conditions are :
f (x=L, t) = ro
U (x=0, t) - 0 (65)
The variables were defined in Section 2 and Appendix A.
A steady-state solution these equations may be obtained using the method
of separation of variables as:
Cos(kx)
f(x,t) - To - sinwt (66)
Cos(kl)
gk sin(kx)
u(x,t) - - j- - Cos(wt) (67)
« Cos(kl)
w
where h is the water depth, and the wave number k «= - is the wave
(gh)1'2)
number.
To obtain a numerical solution, the flow domain was divided with 5 grid
lines in each lateral directions resulting in 25 grid cell s of size of 10 x
10 km. Other parameters are assigned as h - 10 m, g = 980 cm2/sec, and L = 50
km. The initial conditions for the numerical model simulation were selected
to be the at-rest conditions. The results of water surface elevation and
velocities for two locations are plotted in Figures 12 and 13.
69
-------
O
i l
M
CL
in
I
12
24 36 48
TIME (HOURS)
60
72
CO
X.
21 O
CJ -
" 0
12
24 36 48
TIME (HOURS)
60
72
Figure 12, Water Surface elevation and current velocity at x - 5 km (solid
lines represent the analytical solution and dashed lines represent
the numerical solution).
70
-------
24 36
TIME (HOURS
a*
en
o
Qi
" 0
12
24 36
TIME (HOURS
48
60
72
Figure 13. Water surface elevation and current velocity at x - 25 km (solid
lines represent the analytical solution and dashed lines represent
the numerical solution).
71
-------
The analytical solutions are given as solid lines and the numerical
solutions as dashed lines. As can be noted, there is good agreement between
the numerical (HYDR03D) and analytical solutions. The initial differences
between the results are due to the selection of the at-rest initial condition
for the numerical solution. After the initial transient response lasting
approximately 4 to 5 cycles, the model results closely mimics the steady state
analytical solution. The relative maximum differences between the model
results and the analytical solution, after 5 cycles, is less than 5 percent.
Other factors which contribute to the difference between the numerical and
analytical results may be due to the assumptions used to linearize the one-
dimensional wave equations used to derived the analytical solution.
4.2 SUISUN BAY. CALIFORNIA
A previous study of Suisun Bay (Sheng et al. 1986) was chosen to
illustrate the feasibility of applying this model in areas of San Francisco
Bay. The model has not been calibrated or validated, but enough work has been
done to establish that the model is potentially useful. The follow-up study
by Smith and Cheng (1989) in the adjoining San Pablo Bay, is an additional
indication of the feasibility of this model for estuary studies in San
Francisco Bay.
This case study is based on simulations done prior to 1986 with the
EHSM3D code. That code is essentially the same as the code being documented
in this report (see Preface), but there are a few changes and modifications
that have been made since that time. However, none of the changes invalidate
the use of these results to show feasibility of the present code.
One advantage of this study is that it provides a brief review of the
process involved in initially setting up a study. In fact, it has been used
for that purpose in a training workshop (Sheng et al. 1986). Briefly reviewed
in this illustrative example are:
Investigation of pertinent data and previous studies,
Investigation of the processes that may influence circulation,
Initial model setup, including selection of boundary conditions, initial
condition options, and model parameters, and
Interpretation of preliminary results.
Although the study was not carried to the calibration stage, the process
of implementing the model in a feasibility study is well illustrated.
Another advantage in selection this study for illustrative purposes, is
that it may be possible to do additional work to calibrate and validate the
model. Studies have continued in San Francisco Bay that may provide data to
evaluate these initial results.
4.2.1 Physical Setting
Suisun Bay is part of the San Francisco Bay and Delta system in
72
-------
California, which is one of the world's largest and most complex estuarine
systems (Figure 14). The central bay (mean depth 10,7m) is connected to the
Pacific Ocean at Golden Gate (depth of 110m). To the north and northeast,
the system extends to the extremely shallow San Pablo Bay (more than 50% of
the Bay has a depth less than 2m), through the Carquinez Strait (mean depth
8.8m) to Suisun Bay (mean depth 4.3m) and finally into the Delta (See Figures
15 and 16). Each of the embayments of the system usually consist of a deep
navigation channel (depth > 9m) surrounded by shallow shoals. In spite of
their similarity, the various embayments exhibit very distinctive features.
Suisun Bay (Figure 15) is quite complex. It consists of several deep
navigation channels surrounding numerous shoals and islands (see Simmons,
Chips, Van Sickle, Sherman, and other minor islands), and includes two very
shallow sub-embayments, Grizzly and Honker Bays (mean depth <2m). Suisun Bay
has an area of 94 km2 and a mean depth of 4.3m. The main navigation channel
depth is between 9 and 14m and it connects Carquinez Strait and the Delta. The
Delta, which provides 90% of the freshwater in the San Francisco Bay system
has a volumetric outflow rate between 50 and 150 m3/sec i-n summer and 8,000
and 12,000 m3/sec in winter.
4.2.2 Circulation Patterns
Observations and analysis indicate that circulation in Suisun Bay is
affected by four major factors: (1) tides, (2) salinity gradients, (3)
meteorological forcing and (4) bathymetry and geometry. These factors are
explored in this section to indicate what phenomena the model should simulate.
Ocean tides enter the Bay System at Golden Gate and travel a significant
distance through Suisun Bay into the north by northeast end of the Bay system.
Extensive field studies on tidal circulation in the Bay system have been
performed by the U.S. Geological Survey and others (e.g., Conomos and others
1978; Patchen and Cheng 1979; Cheng and Conomos 1980; Smith, 1980; Cheng and
Gartner 1984). The recent report by Cheng and Gartner (1984) provides the
most comprehensive database on tides, tidal currents, and residual currents in
the San Francisco Bay system. Water levels were measured at several stations,
and currents were measured at several current meter located within Suisun Bay
as shown in Figure 15. Salinity intrusion within the navigation channels of
San Francisco Bay and in particular, the Delta system, was investigated by a
U.S. Army Corps of Engineers contractor Kinnetics Laboratories, Inc., (1981),
A tidal harmonic analysis of Suisun Bay data indicate that the major
constituents are the M2 and the Kj tides. Figure 17 graphically illustrates
the spatial distribution of properties of M2 and Kj^ tides, at Stations C26,
C27, C28, C30, C239, and at the east boundary of Suisun Bay. At Station 5103,
the M2 and Kj amplitudes are 52.4 cm and 30 cm, respectively. Values of the
same parameters are 43 cm and 25,4 cm at Station 5112. There is a net phase
shift of approximately 33 degrees between the two stations.
73
-------
San Francisco Bay System
South
San Francisco Bay
37"30* -
122°30
122°
T2T°30'
Figure 14. Map of San Francisco Bay estuarine system.
-------
Current Meter Station
Stage Station
U.S.G.S. Weather Station
>C22i
Martinez
Figure 15. Map of the Suisun Bay region and the location of current-meter
moorings, tide stations, and a USGS weather station.
-------
(a)
(b)
Figure 16, Three-dimensional plot of the Suisun Bay bathymetry when viewed
from (a) the southwest and (b) the southeast.
76
-------
Tidal currents within Suisun Bay were analyzed to determine the Olt Ka,
N2, M2, S2 and M^ tides. Harmonic constants were given along the major and
minor axes of the tidal current ellipses, A strong bi-directional tendency
was observed at most stations and the basin bathymetry was found to
significantly affect the principal current direction also see Cheng and
Gartner 1984), In addition, it was noted that the tidal current speed can
vary up to a factor of two between spring and neap tides.
Salinity varies from the ocean value of approximately 30 ppt at Golden
Gate to the freshwater value of approximately 0 ppt upestuary at the Delta,
The location of the salinity front and the detailed salinity distribution
within the northern reach of the Bay system is significantly influenced by the
Delta outflow. During the low flow summer months, the salinity front may reach
into the confluence of the Sacramento and San Joaquin Rivers (see Figure 14).
During winter, the salinity front may retreat into San Pablo Bay. Salinity
within Suisun Bay varies between 10 to 15 ppt at the western end to between 0
to 10 ppt at the eastern end. Salinity data also exhibit significant daily and
temporal variations. From this it seems clear that the salinity distribution
may have a significant effect on the currents within Suisun Bay and must be
taken into account to properly simulate circulation.
The bathymetry of Suisun Bay is another element that influences
circulation, and it is relatively complex as described above. To illustrate
the complexity, Figure 16 shows 3-D plots of the bathymetry viewed from the
southwest and southeast. In addition and perhaps most important is the
findings of a tidal analysis by Cheng and Gartner (1984) that indicates a
strong influence of bathymetry.
In the initial report, Sheng et al. (1986) does not review any direct
influence of meteorological forcing. It is clear from the initial discussion
above, however, that indirect affects on fresh water inflows are at least an
important seasonal effect on the location of the salt water wedge. In
addition, it is likely that the open shallow Grizzly Bay and Honkers Bay are
subject to some wind-driven circulation. Since a decision was made to ignore
short-term episodic events in the initial study (Sheng et al, 1986),
consideration of meteorological forcing is of lessened importance for this
case study,
4.2.3 Modeling 3-D Circulation in Suisun Bay
For illustrative modeling in Suisun Bay, a grid was defined, initial
conditions were specified, selected data were used to provide a reasonable
representation of the tide and salinity boundary conditions, model parameters
were initially selected, and the resulting simulations were investigated.
Major features and trends of the circulation became the focus of the initial
study. Short-term events episodic were ignored.
The model grid was setup to simulate five vertical layers. Each layer
was divided by grid line spacings of 1/2 km. The resulting network had a
total of 46 x 26 x 5 grid points. The bathymetry arrays were smoothed by the
option provided to eliminate sharp bottom slopes,
77
-------
(n)
SURFF1CE
FIT C25 ( ID. 8 1
SURfnCE ELEVHTI0M HT C27 t 12,11 )
-BO
TIHE ( HOURS I
TIHE
H0URS )
DErTH-nVERnGED CURRENT RT CZ5 I ID, 8 I
DEFTH-HVERRGED CURRENT HT CZ7 ( 12,11 I
TIME ( H0UR5 I
TIME I HCURS 5
OEPFH-flVERnGED CURRENT RT C2B ( ID, 8
DEPTH-PVERRGFD CURRENT flT CZ7 ( 12. 1
TIHE I H0URS !
TIME ! HBURS )
Figure 17. Time histories of simulated surface elevation and
depth-averaged velocity components at (a) C26( (b) C27, (c)
C28, (d) C30, (e) C239 and (f) the west boundary during a
5-day model simulation of tidal circulation in Suisun Bay,
78
-------
(d)
ELEVflTIHH HI C2B I 18.22 I
SURFflCE ELEVnTIHN P? C30 ! 26.
TINE
H0UR5
TIME
H0URS )
CURRENT RT c2B i IB.22 )
DEPTH-nVERRGED CURRENT RT C30 ! Z5-17
LJ
UJ
U
TIME f H0URS )
TIME
tTEPTH-nVERPGED CURRENT flT C28 ( 18.22
FlFPTH-qVERRGEO CURRENT PT C30
u 2
LJ
-2
-3
IM
r-
TIME
03
I H0URS )
to
01
OD
O
TIME t H0URS )
Figure 17. Time histories of simulated surface elevation and
drpth-averaged velocity components at (a) C76, (b) C27, (c)
C2R, (rl) C30, (e) C239 and (f) the west houndary during a
5-day model simulation of tidal circulation In Sulsun Bay.
79
-------
(p)
(O
sutrncE ELEVFUIBH nt C239 i 29. 9
ELEVRTIBN FH LEFT B0UW]PRY I
I H0URS 1
TIME ( H0URS )
UErTH-nvEKnGEO CURRENT RT CZ39 I 29, 9 1
OEFTH-RVERFIGEO CURRENT IT LEFT B0UNDRRT
UJ
5
-10
u
m
_j
UJ
>
TINE (
TIME I H0URS }
UJ
10
-10
CURRFNT nT C239 ( 29. 9
r\i w
rsi
CD a
m i
TIME
I H0URS )
a
CNJ
OEPTH-flVERflGED CURI^NT flT LETT BBUNRflRY f ?,<\ 1
LJ
UJ
UJ
>
TIME ( H0UR5 1
Figure 17. Time histories of simulated surface elevation and
dfpth-averaged velocity components nt (a) C26, (b) C27, (c)
C7B, (d) C30, (e) C739 and (f) the west boundary during a
5-dny model simulation of tidal circulation In Sulsun Bay.
80
-------
Turbulence and friction losses were simulated in straight forward
manner. The simplified second order turbulence closure option was select to
represent vertical mixing. A horizontal eddy viscosity of 100 mz/sec was
selected. The related bottom roughness height of 0.4 cm was chosen.
Tidal boundary conditions were developed from a synthetic tide based on
the previous analysis. A synthetic tide composed of the M2 an(* Kx
constituents only, was applied to the west and east boundaries. Along the
west boundary, the Mz and K^ tides were assigned amplitudes 55 cm and 31 cm,
respectively, and a phase angle of 90 degrees. Along the east boundary,
amplitudes were assigned as 43 cm and 23 cm, respectively. A phase shift of
38 degrees between the west and east boundaries for both the M2 and Kx tide
constituents, was selected. This was necessary because the west boundary of
the computational domain was established westward of the Benicia tide station
(see Figure 15).
Salinity boundary conditions were selected to provide reasonable
representations of observed vertical salinity gradients at the west and east
boundaries, and of the horizontal gradient across the computational domain.
At the west boundary, top layer values of 18 ppt and bottom layer values of 20
ppt were specified. At the east boundary, a top to bottom variation of 14 ppt
was specified.
Boundary fluxes were described with an advection calculation for outflow
and a constant concentration inflow. This was necessary because the EHSM3D
model only allows for outflow-inflow open boundary conditions. When outflows
occurred, the flux was calculated from the one-dimensional advection equation.
Inflows were computed assuming that the flow originated from a. constant
concentration "reservoir" (i.e., inflows were assumed to have a constant
concentration regardless of the history of the outflows).
Simplified initial conditions were selected for circulation and
salinity. The at-rest option as selected to represent a quiescent flow
conditions in the beginning of the simulations. A linear salinity gradient
was assumed to describe salinity across the computational domain.
4.2.4 Results
Simulation of 120 hours (5 days) provided a number of notable results.
Time history results for simulated water level and depth averaged velocity,
illustrated in Figures 18 to 24 for stations C25, C26, C28, C30, C239, and at
the west boundary, show some of the same features noted in the tidal
measurements. These include:
That simulated ebb currents are stronger that simulated flood currents at
many stations, and
That bi-directionality is also apparent in the depth-averaged simulations
of currents.
Simulated current speeds are comparable to observed neap tide
81
-------
measurements In Cheng and Gartner (1984), but spring tide simulations of
current speed are generally smaller than observations. This is probably due
to the idealized synthetic tidal boundary condition, and the simplified
salinity boundary conditions. In addition, there are concerns that the
salinity boundary conditions selected for this simulation do not necessarily
allow the complete internal propagation of baroclinic perturbations at the
boundaries.
Although the water level and currents appeared to reach a dynamic steady-
state in relatively short time periods, salinity values at various stations
were still very slowly changing at the end of 120 hours of simulation. This
is illustrated in Figure 18 that shows the salinity at three vertical levels
(near-bottom, mid-depth, and near-surface) gradually increasing at each of the
stations except at the west boundary that is constrained as shown. The
vertical structure of the velocity fields at 96, 108 and 120 hours shown in
Figures 19, 20 and 21, respectively, indicated that a dynamic steady-state may
have been achieved in the simulations. Other notable occurrences include flow
reversals at some locations and near-surface currents well in excess of 1
m/sec at 108 hours. The lest satisfactory simulations involved the vertical
salinity structure. The simulated salinity fields at 96, 108, and 120 hours
are shown in Figures 22, 23, and 24, respectively. Figure 18 shows the
simulated vertical structure of the salinity field at six stations over the
full course of the simulation. The simulated vertical stratification is not
very pronounced and there are indications that greater degrees of
stratification actually exist (Peter Smith, in review).
The simulation of limited stratification is cause for further investigation
and indicates that present results can only be used for illustrative purposes
until additional calibration is possible. In this case, it does not seem
possible to reach into preliminary conclusions about the flow and salinity
distributions. When calibration is undertaken, the specifications for the
open boundaries, initial conditions, and the coarse five-layer grid spacing
should be examined. Ten or more layers are likely to improve vertical
resolution and a nonuniform Cartesian grid could be used to better resolve the
steep topography between the navigation channel and the shallow areas.
Comparison with additional calibration data is also likely to point out other
input data should be investigate, however, the results for simulated
circulation in Suisun Bay and the test of the model in the adjacent San Pablo
(Smith and Cheng 1989) indicate there seem to be no insurmountable problems
that would prevent calibration.
82
-------
m C2$ t to.e )
20
SRLINITY RT C30 ( 26,17 1
tn
tn
t
z
TIME ( H0URS 1
TIME T H0URS )
II
- 18
16
SRLiNirr Rr C27 r 12,11 )
16
SRLINITY BT C239 I 29,9 1
5 16
TIME ( H0URS 1
TIME ( H0URS )
18
m
m
16
ff
in
SRLINITr m C28 ( 18,22 1
WtOeOOCM^-CDCOO
rsiorw
TIME t HBURS 1
20
SRLINITY RT THE LEFT BOMDRRY ( 2.4 J
18
IB
v\f\f\fv\f
TIME ( H0URS 1
Figure 18. Time histories of simulated salinity at 3 vertical
levels near-bottom, mid-depth and near-surface during
a 5-day model simulation of tidal circulation in
Suisun Bay.
83
-------
SUISUN BOY ; SflLINITY RUN : 96 H0UR5
VEL0CITY F0R 5IGMR = -0.9
2.506*01
mxinun
VE1.0CJTY F0R SIGMR - -0. !
7.S1E»01
HRXIHUM VECT0R
Figure 19. Simulated Tide- and salinity-driven currents in
Suisun Bay near the bottom (a - -0.9) and near the
surface (a = -0.1) at 96 hours.
84
-------
SU1SUN BHY : SRLINJTY RUN : 10P H0UR5
VELBHTY F0R SIGHO = -0,9
VETTDR
VEL0CJTY F0R SIGMR = -D. !
i.26E«02
MRXJHUM VECT8R
Figure 20. Simulated Tide- and salinity-driven currents in
Suisun Bay near the bottom (a = -0.1) and near the
surface {a -0.1) at 108 hours.
85
-------
5UISUN BflV : SBLIN1TY RUN : 120 H0URS
VEL0CITY F0R SIGHfl r -0.9
Z.SOE'0!
VECTUR
VEL0CJTY F0R SJGHR = -0.J
tiff// SS~
, I I /
i i i f r /.
f f y x- y f r ii
₯tf%'',':\\\
il./S/SSs,..'.1^. N.v V ,S, X-
f {/,',',' :d^^;-s^Q^7
,1 u//, !)>-:::T3;::;::-.n-^^f f:
' '. *it,f?ff ... iv ...... .CI. vH>.vN\\\*\ 1 .
/ / l/S/s ) 1. . . .V\ \ If
WXlHUf VECT8R
Figure 21,
Simulated Tide- and salinity-driven currents in the
Suisun Bay near the bottom (o = 0.9) and near the
surface (a - -0.1) at 120 hours.
86
-------
SUISUN BRY : SPLJNITY RUN : 96 H0URS
SPLIN1TY F0R S1GMP = -0.9
IE CW78LRE
LEVEL? FRBtt 13,D 1C 20.C
rNTERVflL PT .500
SflLINITY F0R SIGMR = -D.]
19 CBKTBJRS
CBNTBUT LTVtLS FTWH 12.0 Tfl 1B.D
CWTBJR INTWVfl W .SOD
Figure 22. Simulated salinity distribution in Suisun Bay near
the bottom (a = -0.9) and near the surface (a = -0.1)
at 96 hours.
87
-------
SUISUN BflY ! SRLINJTY RUN : 106 H0UR5
SRL1N1TY F0R S1GMR = -0.9
19 CW7BJRS
UVELS rra« is.5 Te 20.0
C«NTRJR fNTERVflL BT .500
SRLINIJTY F0R SIGMR r -0.1
LFVtLS nWH 1Z.3 T« JB.O
CIN1BUR INIERVfl F .500
Figure 23. Simulated salinity distribution in Suisun Bay near
the bottom (a = -0.9) and near the surface (a = -0.1)
at 108 hours.
88
-------
BUI SUN env : SRLINJTV RUN : izo HBURS
SflLlNJTY F0R SIGHP = -0,9
1C
WM 13.0 ie 20.0
CfNTRW IHTEKVfl BT .500
SRLINITY F0R BIGMR = -0.1
J9
LtVtLS MfflM 12,0 Tf 16.0
CBKTBIH INTERVRL IT .500
Figure 24. Simulated salinity distribution in Suisun Bay near the
bottom (a = -0.9) and near the surface (a = -0.1) at
120 hours.
89
-------
4.3 Charlotte Harbor. Florida
A previous study of Charlotte Harbor (Sheng et al. 1986) was also chosen
to Illustrate the feasibility of applying this model, as was done with the case
study for Suisun Bay. The model has also not been calibrated or validated for
Charlotte Harbor, but enough work has been done to establish that the model is
potentially useful in partially mixed U.S. Gulf Coast estuaries.
The simulations for this case study were also performed prior to 1986 with
the EHSM3D code. As noted before (see Preface), that code is essentially the
same as the code being documented in this report. There are a few changes and
modifications that have been made since that time, but none invalidate the use
of these results to demonstrate the feasibility of using the present code.
This case study also provides a brief review of the process involved in
initially setting up a study. It has also been used in a training workshop
(Sheng et al. 19S6) for that purpose. In addition, this study includes other
illustrative investigations of interest. Briefly reviewed in this illustrative
example are:
Investigation of pertinent data and previous studies,
Processes that may influence circulation,
Initial model setup, including selection of boundary conditions, initial
condition options, and model parameters,
Investigation of pronounced freshwater effects,
Analysis of affects of the grid scale,
A study of the affects of initialization, and
Interpretation of preliminary results.
Although the study was not carried to the calibration stage, the process
of implementing the model in a feasibility study is well illustrated. The study
focussed on investigation of general trends and behavior, and the sensitivity
of the results to various model options and approaches. Short-term episodic
events, including the affects of tropical storms, have not been considered.
4.3.1 Physical Setting
The Charlotte Harbor area of southwest Florida is a shallow water body
with complex boundaries and flow patterns. The estuary (Figure 25) receives
discharges from the drainage of 16 percent of the State of Florida through the
Peace, Myakka and Caloosahatchee Rivers. The estuarine system is connected with
the Gulf of Mexico through various inlets between the barrier islands on the west
of the system. The northern area of Charlotte Harbor (Figure 26), which has a
maximum water depth of approximately 7m, is of particular importance because of
the rapid development of adjacent land areas. The Pine Island Sound to the south
of Charlotte Harbor is extremely shallow (maximum depth 2 m). The hydrodynamics
of the system are complicated by islands, shoals, and multiple openings to the
Gulf of Mexico.
90
-------
CHARLOTTE HARBOR
DE SOTO
SARASOTA
r KOI
CHARLOTTE
CHARLOTTE
FORT MEYERS
GULF OF MEXICO
Figure 25. Map of Charlotte Harbor Estuarine System.
91
-------
Figure 26, Map of northern Charlotte Harbor with locations of
water quality/current meter stations during the June
and July 1982 study.
92
-------
4.3.2 Circulation in Charlotte Harbor
Like that in most partially mixed estuaries, circulation in Charlotte Harbor
is affected by ocean tides propagating through the Harbor entrances, salinity
gradients caused by freshwater inflows, meteorological forcing, and Harbor
bathymetry and geometry. Tidal records indicate a primary diurnal tide with some
semi-diurnal influence. Measurements of water surface elevations on which the
analysis of the tides were based, were made at the Harbor entrance, Brunt Store
Marina, and north and south of Pine Island. Figure 27 shows the tidal elevations
measured during July 20 and 21, 1982 at Brunt Store Marina.
Water quality and current meter data were collected at a number of stations
shown in Figure 26. The data of interest were collected during June and July
1982.
The freshwater inflow from the Peace River strongly influences the
circulation within Charlotte Harbor. The volumetric flow rate can vary from less
than 2,000 cubic feet per second (cfs) to 18,000 cfs within a week (Figure 28).
This kind of flow variation can significantly affect the location of the salinity
intrusion front. Potentially the tides can propagate upstream into the Peace
River.
Due to the shallowness of the estuarine system, tropical storms from the
Gulf of Mexico can significantly affect the tides and circulation within
Charlotte Harbor, In this study, however, these extreme events have not been
investigated. Other meteorological effects and the effects of bathymetry and
geometry were also not investigated in detail in the initial phases of this
study (Sheng et al. 1986).
4.3.3 Modeling 3-D Circulation in Charlotte Harbor
For illustrative modeling, a grid was defined, initial conditions were
specified, selected data were used to provide a reasonable representation of the
tide and salinity boundary conditions, model parameters were initially selected,
and the resulting simulations were investigated. Select sensitivity analyses
of grid resolution and other options were performed.
The model grid consisted of nine vertical layers and horizontal grid line
spacings of 1 km. The computational domain was extended from the northern part
of the Harbor to between stations 17 and 19 shown in Figure 26. The domain
extended about 1 km into the Peace River beyond station 10 shown in Figure 26.
This resulted in a domain of 11 x 11 x 9 grid points.
The tidal boundary condition at the southern end of the computational
domain was based on water elevation data collected during June 25 to 27, 1982
at the Harbor entrance. From these data, the water elevations (f) for the
boundary condition were determined to be;
93
-------
TIDAL STAGE AT BURNT STORE MARINA
CHARLOTTE HARBOR, FLORIDA
E
3
2.0
2 1-°
5
t
S o
LU
O
H
-J-1.0
g
H
July 20,1982
July 21,1982
M
N 18 M 6
TIME (hrs.)
N 18
M
Figure 27. Tidal stage at Burnt Store Marina during July 20 and
July 21, 1982.
-------
DISCHARGE OF PEACE RIVER AT ARCADIA, FLORIDA
I I I
16,000
12,000
LU
a
CC 8,000
3:
o
4,000
10
Recurrence Interval (years)' '
20 30 10
TIME (days)
20
30
Figure 28. Discharge of Peace River during June 1 to July 30,
1982,
95
-------
f - A sin (2n-t/T) + C (62)
where A is 33 cm, T is 24 hours, and C is 11 cm.
It was suspected a priori that the inflow from the Peace River should be
influenced by tides and stratification. As a result, the boundary condition
was selected to modify the boundary velocities accordingly. It was assumed for
illustrative purposes that the inflow velocity at the surface was 2.5 times the
average velocity known from measurements upstream. The bottom velocity was
assumed to be -1.5 times the average velocity, where the negative sign denotes
that the bottom waters are assumed to move upstream and out of the domain at the
boundary, partially in compensation for the increased surface inflow rates.
Inflow rates of freshwater were specified as follows:
Peace River: 15,000 cfs, and
Myakka River; 1000 cfs.
The salinity boundary conditions were approximately represented for the
southern boundary and assumed to be constant at the Peace River boundary.
Salinity along the open southern boundary is computed from the one-dimensional
advection equation with a prescribed valued obtained from measurements at station
17 (see Figure 26) located outside the model domain. No tidal variation in
salinity was specified for the river inflow and outflows.
Specification of initial conditions involved several steps, including
preliminary simulations to set up the final simulations. Initially the
simulations were begun with a quiescent flow field (at-rest conditions) and no
salinity stratification. A simulation was conducted for 24 hours and this new
condition used to establish the initial velocity field for the next series of
simulations. The simulated currents after 24 hours are shown in Figure 29.
These are idealized presentations that are difficult to interpret, but the strong
surface currents from the Peace River are simulated as expected.
The salinity initial conditions for the next series of simulations was
obtained by quadratic interpolation from the measurements at seven stations in
northern Charlotte Harbor over the 2-day period between June 25 and 27, 1982.
These interpolated salinity fields are shown in Figure 30, where bottom
salinities of up to 21 ppt were derived in the bottom layer of the southeastern
section of the domain. The interpolations indicated that the surface waters were
relatively fresh.
The effects of finer grid resolution were explored by switching to a 1/2
km grid line spacing in the lateral and horizontal directions. This led to
differences in the model domain and the open boundary condition at the southern
end of the Harbor. In testing the finer grid, the domain was extended to station
17 (Figure 26) in the southern part of the Harbor and extended further past
station 10 into the Peace River.
A better representation of the southern tidal boundary condition was
employed. Two tidal constituents were used to approximate tidal forcing.
Vertical salinity profiles were estimated from the data collected at stations
96
-------
UV VELOCITY AT SIGMA OF -0.944
0
B
7.1BE«00
MAXIMUM VECTOR
UV VELOCITY AT SIGMA OF -0.056
B
3.3«E»01
MAXIMUM VECTOR
Figure 29. Initial 3-D velocity field in Charlotte Harbor for a
model simulation from June 25 to June 27, 1982.
97
-------
SALINITY AT SIGMA OF -0.944
o
in
o
7 CONTOUnS
CONTOim LEVELS FROM 3 OH TO 21.0
CONTOUR INTERVAL OF 3.00
SALINITY AT SIGMA OF -0.056
o
J>
B
COHTOUH LEVELS mOM .600 TO 4.20
CONTOUR INTERVAL Or .800
Figure 30, Initial salinity field in Charlotte Harbor for June
25, 1982.
98
-------
Vertical salinity profiles were estimated from the data collected at stations
15, 19, 20, and 17 and employed in the simulations. Consequently, stratification
is weaker along the southern boundary in the fine grid simulation. Initial
conditions were selected in the same manner as for the coarser grid.
4.3.4 Results
A series of records for velocity and salinity distributions were produced
during the model simulation of 72 hours. Surface and bottom, velocity and
salinity distributions are presented in Figures 31 through 34 for stations 10,
7, 22, and 19, From these results there are several observations worth noting.
First, a dynamic steady state was obtained. Second, the Peace River flow is a
dominate influence. The relatively strong currents at stations 10, 7, and 22,
and the initial reduction in the high salinities at the bottom locations of
stations 7 and 22, are related to the freshwater inflow.
The velocity and salinity fields are shown in Figures 35 and 36,
respectively. In these illustrations, the surface currents are distinctly
different from the bottom currents and significant stratification exists, as
would be expected in a partially mixed estuary. Conditions are similar at the
end of the 72-hour simulation period as shown in Figures 37 and 38.
Dimensionless depth profiles of salinity shown in Figure 39 illustrate the
stratification simulated by the model. These results are for stations 7, 22,
15, and 19 at the end of the 72-hour simulation.
The 48-hour simulations were performed after the initial 24-hour
simulation, using the finer grid and the boundary conditions described above.
These results were not interpreted (Sheng et al. 1986), but are presented for
the reader to investigate here. The resulting circulation pattern at the end
of the 48-hour simulation is shown by the surface and bottom currents in Figure
40. The salinity distributions at the same time are shown in Figure 41, The
time variations of surface elevation, surface currents, bottom current, surface
salinity, and bottom salinity at a number of stations are shown in Figures 42
through 45.
99
-------
9JnrncE FLEvnnoM ni ( LO. 9 i : srnio
sn ,,.,.,,.,,^
ac
u
10
£ n
UJ
"'-in
0 5 12 15 24 30 36 12 18 54 60 66 72
TIME ( M0URS I
stmrncE CURRENT RT
STRIO
SURfRCE CURRENT RT ! 10. 9 ) : STF110
0 6 12 19 24 30 36 42 49 54 60 66 72
TIHE < H0URS 1
0 6 12 18 24 30 36 42 48 54 60 66 72
TIME ( H0UR3 )
B0H0H CURRENT RT S 10. 9
STRIO
B0TT0H CURRENT RT ( 10. 9 ) : STRIO
ui
-4.0
-3.6
§3.2
^z.n
=2.1
0 6 12 IB 24 30 3B 12 48 54 60 66 72
TIHE ( HWJRS I
s- _
0 6 12 18 24 30 36 42 48 54 60 66 72
TIHE f H0URS )
.15
-.10
£-35
- .30
-.20
5.I5
j . ID
£.05
3URFHCE SflLINITY RT ( 1C. 9 1 : 3TRIO
0 G 12 18 24 30 36 12 40 54 60 66 72
TIHE I H0URS I
t
0.
10
9
e
7
6
5
4
3
2
1
B0TTHH SflLINITY RT [ 10. 9 1 : STRIO
0 5 12 18 24 30 3B 42 48 54 60 66 17
TIME ( H0URS I
Figure 31. Time Histories of water level, surface currents,
bottom currents, surface salinity and bottom
salinity at Station 10 during the 3-day model
simulation period.
100
-------
SUP.mCE ELfcVMUHH HI
5
tr
i.
UJ
-in
0 6
18 21 30 36 12 18 51 50 55 72
TIME « H0URS )
-5
-ID
-15
-20
-25
-30
-35
-10
-15
-50
suRrncE CURRENT nr t 7. 9 ) : sin?
0 6 12 18 21 30 36 12 18 51 60 66 72
TINE ( H0UR5 1
o 9
uj g
3* 7
" G
,_ 4
t 3
H 2
CD 1
^ o
^-i
CURRENT OT t 7, 9 1
0 6 12 18 21 30 36 12 18 51 60 66 72
TIME t H0URS )
B0TT0H CURRENT fll { 7. 9 ) ! STfl?
B0TT0H CURRE^fT flT ! 7, 9 I : STF17
0 6 12 IB 21 30 36 12 18 51 GO 66 72
TIME ( H0LRS J
0 6 12 18 21 30 36 12 18 51 60 66 72
TIME < H0URS !
SURFflCE SflUNITY flT I 7. 9 ) : 3TP7
B0TT0H SRLINITY RT ( 7, 9 ) ; STR7
-2.4
£2.0
" 1.6
>- 1.2
z .9
d . T
0
0 6 12 18 21 30 36 12 18 51 60 66 72
TIME ( H0URS )
0 6 12 18 24 30 36 12 18 51 60 66 12
TIME I H0URS )
Figure 32. Time histories of water level, surface currents,
bottom currents, surface salinity and bottom
salinity at Staion 7 during the 3-day model
simulation period.
101
-------
,URrncE ELEvnnaM tu < 5, 7 ) : 57022
0 6 12 18 21 30 36 12 18 51 60 66 12
TIME ( HflURS 1
10
5
0
-5
-10
-15
-20
25
-30
-35
-10
SURFF1LE CURRENT mi 5, 7 1 : STR22 . SURFflCE CURRENT RT t 5, 7 I : STR2Z
0 6 [2 JB 24 30 36 12 49 54 60 66 72
TIHE I H0UR3 )
8
sc
cj
^ 0
a-8
d~'2
r-te
0 6 J2 10 24 30 3G 42 48 54 60 G6 72
TIME I H0URS )
B0II0H CURRENT RT ( 5, 7 J : 5TR2Z
0 6 12 18 24 30 3B 12 48 54 SO 66 72
TINE ( H0URS I
S0TT0H CURRENT RT 1 5, 7 I : STR22
5 2
0
0 6 12 18 24 30 36 42 48 54 60 66 72
TIME I H0URS 1
7.0
6'0
3URFRCE SflLINITY RT ( 5. 7 ) : 3TR22
K
Ul
1.0
0 B 12 18 24 30 36 12 48 54 60 66 72
TIHE I H0URS I
90TTBM SRL1NITY flT t 5. 7 I : 3TH22
0 6 12 18 24 3D 36 42 48 54 60 SB 77
TIHE ( HOURS !
Figure 33. Time histories of water level, surface currents,
bottom currents, surface salinity and bottom
salinity at Station 22 during the 3-day model
simulation period.
102
-------
EiE'vnngM nr i 5. 2 i *. smia
0 6 12 18 21 30 36 12 18 51 60 66 72
TIME f HOURS )
SURFnCE CURRENT FIT [ 5. 2 1 : STfllS
fe -6
'J .1
>- 0
tr. -,2
H -.1
d -.6
> -.9
3-1.0
0 6 12 18 21 30 36 12 16 51 60 66 72
TIME I H0UR5 1
25
20
15
10
5
0
-5
-10
-15
5URFRCE CURRENT Rt ( 5. 2 1 : SIR 19
0 6 12 IB 21 30 3G 12 IB 51 60 66 72
TIME ( H0URS )
1,0
"-
-------
CHRRLSTTF HRRB0R : FILE - CH3001D
UV VEL0CITY RT TIME 0F 24.0 H0UR5 HMO SIGHR 0F -0.91-1
1 t
\ 1
T \
< \
\
/
t
N 1 /
N \ \
X \ \
CHRRL0TTE HRRB0R I FILE = CH3D01D
UV VEL0CITY flT TIME 0F 21.0 H0URS RND S1GMR 0F -0.056
) t
_| I
MflXIMUH VECT0R
Figure 35. Computed 3-D velocity field in Charlotte Harbor after
24 hours of simulation.
104
-------
CHARLOTTE HHRB0R : FILE = CH3DOID
BfiLiNITY RT TIME 0F 24.0 H0URS RND SIGMR 0F -0.9*4
9 C0NT0LJRS
LEVELS FRflH 6.00 T0 22.0
CflNTBUR tmCRVaL BF 2.00
CHRRL0TTE HRRB0R : FILE = CH3D01D
SRL1NTTY RT TIME 0F 24.0 HBURS fiNO SIGHfl 0F -0.056
6 CftWTBURS
LEVEL3 FR0N J.DO 10 fl.OO
INTERWL BF 1.00
Figure 36. Computed salinity field in Charlotte Harbor after 24
hours of simulation.
105
-------
CHRRL0TTE HRRB0R t FILE - CH3D01D
UV VEL0CITY RT TIME 0F 72.0 H0URS RND SIGHfl 0F -0.914
I I
\ \ /
\ \ 1
\ s \
V «-. N
N \
i '
' X V
i '
B
VECTBR
CHPRL0TTE HRRB0R : FILE = CH3001D
UV VEL0C1TY PT TtHE 0F 72.0 H0UR5 flND SIGMfl 0F -0.056
\ t
t
HBX1HUH VECT0R
Figure 37. Computed 3-D Velocity field in Charlotte Harbor after
72 hours of simulation.
106
-------
CHRRL0TTE HRRB0R : FILE = CH3D01D
SflLINITY flT TIME 0F 72.0 H0URS RNO SIGMfi 0F -0.944
9 CBMT0URS
CflNTflUR LEVELS FKOH 6.00 T0 22.0
CflNTBUR IHTERVOt flf 2.00
CHRRL0TTE HRRB0R : FILE = CH3D01D
SRLINITY FIT TIHE 0F 72.0 H0URS flND SIGMfl 0F -0.056
\
C*)TOJR LEVELS flWH 1.00 T0 9,00
tMTERWL 0r !.00
Figure 38, Computed salinity field in Charlotte Harbor after 72
hours of simulation.
107
-------
FlLE:CH3DOtD
0
RUN: DODO
X =72.0
a
I
-.1
,z
,3
-.1
-.5
-.9
-1.0
2 3 1 5 6 7 6 9 10 tl
BT STR7 (PPTI
FILE:CH3D01D
0
RUN: 0000
X =72.0
H
-.7
-.9
9
-1.0
S 6 7 8 9 10 11 12 13 14 15
SFLINITY flT STR22 tPPT)
FILC;CH3D010
0
RUN: ODOO
=7Z.O
inoinoiftotrtomai/joinoin
(Dr~r-03eooiO)aa-«»*f>iN»nin
AT STP1S (Wf)
P1LE.-CH3DOID
0
ftUNlOODO
X =7Z.O
_.!
-.2
-.3
uj ,4
-.5
£ -.6
-.7
-.9
-t.D
4 6 Q 10 12 H IB 18 20 22
SfLINITY RT STR1S IPPTJ
Figure 39. Vertical salinity profiles at Stations 7, 22, 15 and
19 in Charlotte Harbor after 72 hours.
108
-------
Charlotte Harbor : File = CH3DHIO
UV Velocity al Time of 18.0 Hours and Slgmo or -0.9 11 t/Sfff" I
%a
VtClOT
Figure 40. Computed 3-D velocity field in Charlotte Harbor
after 48 hours of simulation with a 1/2-km grid.
109
-------
Chorlotte Horbor : File = CH3DHID
Salinity al Tine of 48.0 Hours ond Sigma or -0.941
1 ContoiH
Contour l8v*U fro* 3.00 to 21.0
ConLou- Interval of 3.00
Charlolle Harbor : File = CH30H1D
Sollnlly al Tine of 48.0 Hours and Sigma or -0.056
II Contort
Contour lev«U tra* 1,00 to It.O
Cantor Inltrval of I.M
Figure 41. Computed salinity field in Charlotte Harbor after 48
hours of simulation with a 1/2-km grid.
110
-------
so
40
30
20
ID
SURFfiCE ELEVOTJ0N PT ( 18.27 J : STJJ7
-20
-30
0 6 1Z 18 Z4 30 35 42
H0UR5
SURFFCE CURRENT flr ( 18.27 » ; STP7
SURFnCE CURRENT RT t 10.27 ) : STR7
w
~ -9
M5
M -1Q
g-zo
g-2<
= -28
0 6 12 18 24 3D X 42 48
H0UR3
12 18 24 30 3S 42 48
B0TT0H CURRENT RT ( 10,27 J : STFI7
12 18 24 30 36 42 48
B0TT0H CURRENT flT I 18.27 \ : STR7
0 6 12 18 24 30 36 42 48
3,6
2.0
I!:?
i:!
o
SURfflCE 3flLlHITT RT ( 18.27 } : 3Tfl7
0 6 12 18 24 30 36 42 48
H0URS
22
20
16
B0TTBM SflLtNtTY flT I 18,27 > : 9TR7
0 6 12 IB 24 30 36 42 46
H0URS
Figure 42. Simulated time histories of water level, surface
currents, bottom currents, surface salinity and bottom salinity
at Station 7 during the 2-day simulation period.
Ill
-------
30
25
- 20
6 '5
_ to
5
o: 0
fc -5
^-10
-15
-20
SURFRCE ELEVRTI0N HT I 13.23 ) ! STR22
0 6 12 IB 21 30 36 12 4fl
TINE ( H0UR5 )
SURFRCE CURRENT RT ( 13.23 I : STR22
SURFRCE CURRENT RT ( 13,23 1 : STR2Z
0 6 12 18 24 30 36 42 48
TtHE I H0UR3 )
.-1
0 S 12 16 24 30 36 42 45
TIMt ( H0UR3 I
u
-30
.20
B0TTBH CURRENT RT I 13,23 1 I STRZ2
-.io
-.30
0 6 12 18 24 30 3E 42 48
T1HE ( HOURS 1
B0TT0M CWWENT RT I 13,23 ) : STR22
0 E 12 19 24 30 36 42 48
TIME ( H0URS J
5.
it
SnttMlTT Hr f 13,23 J : 3TR22
B0TTBM 3RL1N1TY flT ( 13.23 J : 3TR22
0
5
0
3.S
3.0
2.5
2.0
t.S
1,0
0 6 12 IB 24 30 36 42 46
TIME ( H0URS J
0 6 12 18 24 30 35 42 48
TIME t H0URS )
Figure 43. Simulated time histories of water level, surface
currents, bottom currents, surface salinity and
bottom Salinity at Station 22 during the 2-day period.
112
-------
30
25
- 20
5 is
_ 10
5
-------
30
25
- 20
5S
5
tt 0
fc -5
^-10
-15
-20
SURFRCE EIEVPT1BN FIT t 30.28 ): STfUO
0 6 1Z 18 Z4 30 36 42 48
H0URS
SURFHCE CURRENT RT ( 30.28 ) : STR10
5URFRCE CURRENT RT I 30.28 } ! STfllO
-
y
d
=3
-I.B
-1.8
-2.0
-2,2
-2.4
-2.6
-2.8
-3.0
-3.2
-3.4
0 6 12 IB 24 30 36 42 48
HWJR3
E-
H
d.
12 IB 24 30 36 42 48
B0TT0H CURRENT RT I 30,28 I : STfllO
6 12 18 24 30 3E 42 48
B0TT0H CURRENT flT t 30,28 ) : STfllO
6 12 10 24 30 36 42 48
.18
.10
£.08
2 -06
3-M
S.02
SURFflCE 9flLtNlTY RT I 30.29 J : 3TRIO
BflTTBM 3RLINITY flT ( 30,28 I : STfllO
0 6 12 18 21 30 36 42 48
HBURS
0 6 12 18 24 30 36 42 48
HBURS
Figure 45. Simulated time histories of water level, surface
currents, bottort currents, surface salinity and
bottom salinity at Station 10 during the 2-day model
simulation period.
114
-------
** GREEN BAY. LAKE MICHIGAN
In this case study HYDR03D is applied as a part of a comprehensive study
of the effects of PCB's in Green Bay sediments. Historically, the Fox River
in Wisconsin has contributed a significant amount of PCB,s to the environment
and it is suspected that much of this contaminant has migrated into Green Bay.
The Fox River has one of the largest concentrations of pulp and paper mills in
the world, from which PCB's are suspected to have been discharged. Most of
the fisheries are presently closed because of the PCB levels in fish.
The study involves a preliminary calibration of the hydrodynamics model
and a sediment model with the historic data available. Following the
preliminary calibration, both hydrodynamics and sediment transport models will
be calibrated with data being processed from the 1989 field season. If the
study is fully successful, it will involve the linkage of hydrodynamics,
sediment transport, and large scale water-quality models of the box type
(Ambrose et al. 1987). Presently, large scale box modeling suffers from an
inability to simulate and predict the effects of complex stratified flows on
transport. An important test will be calibrating the model with historic data
and checking the calibrated parameters with data collected in the summer of
1989 to see if the simulations are predictively valid.
In this case study, results from the preliminary calibration are
reported to illustrate the use of the model in a large lake setting with
complex wind driven circulation. The dynamic nature of Green Bay (Miller and
Saylor, 1985) indicates the need for unsteady state two- and three-dimensional
simulations. HYDRO3D, which is capable of treating stratification and lake-
bay interactions, is used to model the flow and transport processes in the
Bay. In this case study the model is applied to simulate 2-D and 3-D
circulation patterns and these results are found to be similar to general and
specific observations of the Bay.
4.4.1 Physical Setting
Green Bay is a long and relatively shallow water body in northern Lake
Michigan. The Bay is separated from Lake Michigan by the Door Peninsula and
connected to the lake by four main channels near its northern end. These
channels are Martin Island Passage, Rock Island Passage, Porte des Morts
Passage, and Poverty Island Passage (see Figure 46). The Bay is approximately
40 km wide and 190 km long and its main axis is oriented from the north by 38
degrees to the east. More than a dozen streams drain the area around and in
the vicinity, and discharge into the Bay. Major tributaries that contribute
water and sediment to the Bay include the Fox, Oconto, Peshtigo, Menominee,
and Escanaba rivers. The upper part of Green Bay is generally deeper than 20
m, with a maximum depth of 48 m west of Washington Island. The lower half of
the Bay, south of Chambers Island, is 30 m deep near the island, but very
shallow (few meters deep) at the southern end (see Figure 47). Several small
islands exist in the Bay. However for these simulations, only the effects of
Chambers Island will be simulated. Chambers Island has an area of 12 km2 and
is located midway between the mouth of the Bay and Green Bay city.
The flow and circulation are controlled by wind, Lake Michigan water
115
-------
levels (also wind dominated), and river inflows. As the historic data to be
presented in the following sections indicate, seiche is an important phenomenon
in the Bay. In addition, winds significantly control local circulation. There
is a counterclockwise gyre northeast and a clockwise gyre southwest of Chambers
Island in Green Bay that typically describe the general circulation patterns.
These general trends must be simulated for the hydrodynamics model to achieve
general usefulness in this study.
4.4.2 Two-Dimensional Simulation of Flow
The HYDRD3D model is tested in the 2-D mode using historical data from
Green Bay given by Heaps et al, (1982). The main objective of these tests are
to show the response of the Bay to wind forcing and Lake Michigan water level
changes. The limited data available for a 2-D simulation consist of:
a) Water level observations at the mouth of the Bay (St. Martin Island, and
Plum Island) and at two other stations: Menominee and Green Bay cities.
b) Hourly wind observations at the airport in Green Bay city. No current data
are available due to the malfunction of current meters.
These data are used to specify the boundary conditions at the passages into Lake
Michigan and the wind shear on the Bay. The simulations of water movements in
the Bay are performed for two periods when data were collected. These periods
are September 17-20 and October 8-12, 1969. The simulation for each period was
started from the at-rest condition to define the initial velocity field for the
simulations.
Hourly wind data taken at the airport (a short distance from the south
end of the Bay), are used to calculate the time-varying wind stresses acting
over the entire water surface of the Bay. Winds were variable during the two
simulation periods. During the September period, the wind directions were mostly
northward for the first day of simulation (September 17), northeasterly during
September 18, easterly during September 19, and southwesterly during September
20. The maximum wind speed during this period was about 8.5 m/sec and occurred
on September 8. Wind speeds during the October period were stronger than those
for the September period. The maximum wind speed during the October period was
about 12 m/sec and occurred on October 9. The wind directions during this period
were westerly during October 8, southerly during October 9, southwesterly during
October 10, northerly during October 11, and northeasterly during October 12.
Two dimensional simulations of water surface elevations and depth-averaged
velocities are performed on a 2 km x 2 km horizontal grid network. This network
has a total of 21 x 96 grid cells, as shown in Figure 48.
The water surface elevation data measured at the mouth of the Bay and near
Green Bay city for the September and October periods are given in Figures 49 and
50, respectively. These data are used as input for the model. As shown, the
variation of the water surface at the mouth ranges from 5 to 10 cm for each
period of the data set. The major force that distinguished the October results
from the September results is the wind force. The wind direction also plays a
major role on the general circulation in the Bay. Figures 51 and 52 show these
116
-------
GREEN BAY
rjti.lt
Hi S Q 10 M It".
l « J in ii a Mi
14
10
"IBT
Figure 46, Map of Green Bay showing relation to Lake Michigan
and other Great Lakes
117
-------
118
bathymetry
-------
x
E
Figure 48. Grid network of Green Bay
s
119
-------
GREEN BAY
SEPT 17-20 1969
UJ
111
UJ
O
to
cc
LEGEND
GBay City
GB'av"Mbu(fi
o.o
12,0
24.0
86.0
48.0
TIME (HR)
60.0
72.6
V»6.C
Figure 49 . Measured vater surface elevation at the mouth of
Fox River near Green Bay city during September,
1969
120
-------
GREEN BAY
Oct 8-12 1969
o
to
O
" 'o
Z o
O"
I
UJ o
uJ?
UJ
O
1=
05
CC
UJ
2: °
5?
LEGEND
Boundary
GTeehBiy"
o
1C
0.0 12.0 24,0 36,0 46.0 60.0 72.0 84.0 86.0
TIME (HR)
108.0
120.
Figure 50. Measured water surface elevation at the mouth of
Green Bay and at Green Bay city during October, 1969
121
-------
CIRCULATION IN GREEN BAY
October 1969
TIME = 54,00 HOURS
TIME STEP =324
200 n
1800
1800
HO 0
120 0
s
3
1
3, 100 o
a
a
BO 0
90 0
40 0 -
20 0
0 0
50
Velocity Scale
00 100 200 300 400 SO (t lo 0 TOO 800 BIO 1000
X in Kilometers
Figure 51. Simulated circulation in Green Bay using October
1969 data (54 hours from at the at-rest)
122
-------
CIRCULATION IN GREEN BAY
Oclober 1969
TIME = 114.00HOURS
TIME STEP =684
zoo.o
leoo
1800
140.0
130.0
c
I
.3
a
60.0
flOO
400
800
50
Velocily Scale
0.0 100 200 300 400 500 *0 0
X in Kilometers
TOO BOO 900 1000
Figure 52. Simulated circulation in Green Bay using October
1969 data (114 hours from the at-rest)
123
-------
effects for southwesterly and northeasterly winds during the October period.
During certain times in October, a counterclockwise circulation was calculated
due to southwesterly wind. These two-dimensional results, however, should not
be compared with the general pattern of circulation in the Bay which was reported
by Miller and Saylor (1985) and Modlin and Beeten (1970). This counterclockwise
circulation is shown in Figure 51 where, along the western shore of the Bay flow
is southward, along eastern side flow is northward in respect to long axis of
the Bay, and near Green Bay city the flow is from left to right parallel to the
x-axis of the grid. With the forcing of northeasterly wind, the currents run
along both shores, producing counterclockwise flow along Chambers Island in the
northern part of the Bay (Figure 52).
Figure 53 shows the computed and measured water surface elevations
versus time at a station near the City of Green Bay for September 17-20, 1969.
During the first day of the simulation a lag was observed between the measured
and computed water surface elevations. Because zero initial values for dependent
variables such as water surface elevation and velocities were used, it was
concluded that this lag is due to stabilization time, or the time that a column
of water will take to absorb the inertia of a suddenly applied wind stress.
Differences between the measured and computed water surface elevations during
the last day of simulation may be attributed to the use of the over-land wind
data measured at the Green Bay airport for over-water wind data in the Bay,
For the October simulation shown in (Figure 54) substantial
differences are apparent between the measured and computed water surface
elevation during the first two days and last day of the simulation.
Nevertheless, there is good agreement between the observed and computed
oscillatory patterns. Again as mentioned by Heaps, et al. (1982), the large
differences between the simulated and measured water level elevations may be due
to some forces that affected the observed values and were not included in the
computations. But in general, the model responded fairly well considering the
inadequacy of available data used in the simulation.
4.4.3 3-D Simulation of Flow
Since circulation in Green Bay seems to be three dimensional, the
application of the 3-D mode is expected to provide improved simulations. This
is investigated in this case study by applying the model in the 3-D mode.
Heaps et al. (1982) studied water motion in Green Bay by analyzing
the measured field data for September and October 1969. These investigators
pointed out that the main external forcing mechanisms to the Bay water included
the wind, the semidiurnal tide, and the first free longitudinal mode of
oscillation of the Lake Michigan, with the latter two forcing components acting
at the Bay mouth. Using a vertically averaged 2-D numerical model, Heaps et al.
was able to simulate the water surface and the vertically integrated currents
due to specific external forces. With the currents and temperature measured at
different depths and locations within the Bay, Miller and Saylor (1985) analyzed
the data and found strong variations of water motion and temperature in both the
horizontal directions and in the water column.
124
-------
GREEN BAY
SEPT 17-20 1969
ALEGEND\
/ ISIMULATEC
o.o
12.0
24.0
36.0
48.0
TIME (HR)
Figure 53. Measured and calculated water surface elevation at
Green Bay mouth near Green Bay city during
September 17-20, 1969
125
-------
GREEN BAY
Oct 8-12 1969
o
o-
Og
Z
O
I"
f^ CD -
01 *
-J
UJ
LU
O
03
tr
ILI
LEGEND
SIMULATED
MEASURED
0.0 12.0 24.0 36.0 48.0 60.0 72.0
TIME (HR)
84.0
96.0
108.0
120.1
Figure 54. Measured and calculated water surface elevation at
Green Bay mouth near Green Bay city during October
8-12, 1969
126
-------
Due to the 3-D characteristics of water motions within the Bay, both Heaps
et al. (1982) and Miller and Saylor (1985) pointed out that a 3-D numerical
hydrodynamic model is essential to accurately simulate the circulations in the
Bay. The HYDR03D model was used to simulate the 3-D currents within the Bay;
in which the domain was divided into 21 grid points in x-direction and 96 grid
points in y-direction, and 5 vertical layers resulting 21 x 96 x 5 square cells.
An experiment using the model, was conducted to determine the
responses of water motion in Green Bay under the action of a uniform wind field.
Starting from a zero initial condition, we applied a uniform wind of 8 m/sec to
the Bay and held it constant throughout the simulation period. The wind was
primarily directed along the main axis of the Bay (negative x-axis) , 38°
clockwise from the north toward the City of Green Bay. For the sake of simplicity
and the unavailability of boundary conditions at the mouth, the Bay was assumed
to be an enclosed domain.
With the above assumptions and grid configurations, a 3-D simulation
was performed for a duration of 40 hours. Figure 55 shows the water surface
elevations in Green Bay after 40 hours. Due to the direction of the wind, a
positive water surface profile is maintained toward the City of Green Bay. The
profile decreases to zero somewhere in the middle of the main axis and then to
negative values in the two northern gulfs. Figure 56 shows the 3-D, vertically
averaged currents in the Bay at the same time. The currents along the shallow
shore regions are driven by the wind. Currents against the wind in deeper
central regions are driven by the pressure gradient associated with the positive
surface setup as was shown in Figure 55. This phenomenon is often seen in the
studies of estuarine and lake hydrodynamics.
Several small gyres are distributed in the Bay. These gyres are associated
with the bathymetry and geometry of the Bay. Comparing Figures 55 and 56 with
the Heaps's 2-D model results, we find that the agreement between surface
elevations from both models is remarkably good. The (maximum) surface setup at
the City of Green Bay is 11.7 cm from HYDR03D and 11 cm from Heaps' model. The
general patterns of 2-D circulation in Figure 56 are very similar to those of
Heaps' model except in regions near the mouth. In Figure 56, from the mouth to
the northern, shore, there exists two types of circulation. One is
counterclockwise near the mouth and the other clockwise near the northern shore.
In Heaps' results, these two gyres are merged into one large counterclockwise
gyre extending from the mouth to the northern shore. This difference of local
circulation may be attributed to the fact that in Heaps' model, the mouth is not
a boundary; rather, it is continuously connected to the lake. An open boundary
is assumed to exist in the central region of the lake, which is far away from
the mouth and hence the Bay. In the present model, however, we assume rigid,
closed boundaries along the mouth. Figure 57 shows the currents in the near-
surface layer and Figure 58 shows those in the near-bottom layer.
127
-------
Figure 55. Water surface elevation in Green Bay after 40 hours
128
-------
GREEN BflY
CONSTRNT WIND (RFTER 40 HOURS)
30.6 CM/SEC
3D - VERTICALLY flVERflGED
WIND (3 M/SEC1
Figure 56. 3-D vertically averaged currents in Green Bay after
40 hours
129
-------
Similar to the wind-driven currents along the shallow shore regions,
the near-surface currents also are driven by the wind and hence follow the
direction of the wind. These unidirectional surface currents cause a gradient
of water surface elevation along the wind direction. To balance the pressure
gradients caused by the water surface setup, the currents return in the lower
layers (near the bottom), in particular, the deeper regions, as shown in Figures
57 and 58.
4.5 Prince William Sound. Alaska
To further test the capability of HYDR03D under a different situation it
was applied to Prince William Sound in Alaska, to simulate water circulation
during the recent oil spill from the EXXON Valdez that began on March 24, 1989.
4.5.1 Physical Setting
Prince William Sound (Figure 59) lies on the southern coast of Alaska. The
sound covers an area of approximately 8000 square kilometers (3090 square miles)
and includes many islands of various sizes. Most of the islands are concentrated
in the western half of the sound, leaving a large open area of approximately 1800
square kilometers (700 square miles) in the eastern half. The terrain in the
area is very rough, creating numerous bays and causing the shoreline of the sound
and islands to be quite irregular.
The sound is separated from the Gulf of Alaska by Montague and Hinchinbrook
Islands, which form the south-eastern boundary, and has two major connections
to the Gulf of Alaska. The Hinchinbrook Entrance is 11.4 kilometers (7.1 miles)
wide and opens directly to the gulf between Hinchinbrook and Montague Islands,
at about the middle of the eastern side of the sound. At the southern end of
Montague Island and of the sound, Montague Strait forms an 8.4 kilometer (5.2
miles) wide passage parallel to the main shoreline. The average depth in the
Hinchinbrook Entrance and Montague Strait is 300 meters (980 feet) and 195 meters
(630 feet), respectively.
A navigation channel extends from the port of Valdez in a bay at the
northern end of the sound, across the previously mentioned open stretch of water,
and out through the Hinchinbrook Entrance. Depths along this channel are
primarily in the range of 275 to 460 meters (900 to 1500 feet). These depths
are typical of the more open, western half of the sound. In the eastern half,
a scattering of islands separates the sound into a network of passages of widely
varying widths and depths. The two major passages lie on either side of the
largest interior island, Knight Island. The passage between Montague and Knight
Islands averages about 6.3 kilometers (3.9 miles) wide and 180 meters (600 feet)
deep; the narrower passage between Knights Island and the main shoreline averages
about 10.1 kilometers (6.3 miles) wide and 400 meters (1300 feet) deep. The
maximum depth in the sound is approximately 870 meters (2860 feet) and occurs
in the western half of the sound, off the northern end of Knight Island.
130
-------
GREEN BRY
CONSTRNT NINO (RFTER 40 HOURS!
59.o OH/SEC
3D VELOCITY
LflYER NO. - 5
HfNO 18 M/SECt
Figure 57. 3-D simulation of currents in Green Bay (near the
surface layer).
131
-------
GREEN BRY
CON5TRNT WIND (flPTER 40 HOURS)
M,2 CM/SEC
3D VELOCITY
LflYER NO. - i
WIND (8 M/SEC)
^-^*-:*1 \/ ^±sr?zrr:r.3A1 ..,^zry-TTl?^rgr=^y?^T*-r °
Figure 58. 3-D simulation of currents in Green Bay (near the
bottom layer).
132
-------
Figure 59. Map of Prince William Sound, Alaska,
133
-------
4.5.2. Modeling Parameters
Two finite difference grid networks were used to describe the Prince William
Sound area. Initially a coarse uniform grid network (Figure 60) of 35 x 28 square
blocks was used to minimize data processing and computation time. In this grid
system the area of each grid block was 25.8 square kilometers (10.0 square
miles). This relatively coarse grid failed to represent the highly irregular
nature of the shoreline, omitting many small islands, passages, and bays. The
second grid was (Figure 61) four times as fine as the first, consisting of 70
x 56 square blocks where the area of each grid block was 6.5 square kilometers
(2.5 square miles). This finer grid was much more successful in representing
the features missed by the coarser grid.
Both circumscribing grids had an open boundary on their western and southern
sides (approximately corresponding to the south and east of the map) . On the
western side, the open boundary extended from blocks 2 to 10 for thecoarse grid
and blocks 2 to 21 for the fine grid. On the eastern side, the open boundary
extended from blocks 2 to 29 for the coarse grid and blocks 2 to 63 for the fine
grid.
Three simulation runs were performed to calculate the flow field in the
sound. One run for each grid was completed in the 2-D mode in addition to a 3-
D run (with three layers) using the coarse grid. For simplicity in these
comparative test runs, a tidal amplitude of 3.0 meters (9.8 feet) with no phase
angle was assumed along the open boundaries. All other factors, such as wind
stresses and river inflows, were neglected. A time step of 1.0 and 1.5 minutes
was used for the coarse grid and fine grid runs, respectively. These runs were
used to compare the results from a 2-D and 3-D analysis and from a coarse and
fine grid analysis (2-D only).
4.5.3 Results
Scale vector plots of the calculated velocity field were obtained at hourly
intervals for each run. These plots are shown in Figures 62 to 64. The scale
of the velocity vector is given in terms of its horizontal, and vertical
components one inch equals 0.77 meters per second (2.5 feet per second). An
arrow with no stem indicates that the velocity is too small to be revealed at
this scale. The map scale is 1 inch equals 23,4 kilometers (14.6 miles).
Several velocity vectors in the lower right corner were unexpectedly large.
At one point in particular, the velocity was so great that, for several plots,
this vector was truncated at the boundary of the plot. Extensive mud flats exist
here, causing some of this area to be declared as land (the blank area) and
others to be very shallow. The isolated large velocities in this corner may be
attributed to this shallowness combined with the boundary effects.
134
-------
Figure 60, Coarse grid of Prince William Sound, Alaska.
135
-------
nniiijiEirr icwi
itii
J4J4+H-HHf»H!-H4tt*f4i4+-H-H-
,JL,l_Lt_i_Lji_j_JLjLiJlwti._i_i.i.J.j.J,ji,j.Xa.~LLLl~4_l,J-.j_-t,
4- -J-H-J-H-4J- i-Hfir
:L:i±J:fm4::±|:t:i:i:j4
Figure 61. Fine grid of Prince William Sound, Alaska.
136
-------
anoq -[) pjaS SSJEOO 3ui:sn punog
UT s^uaa^no jo
-------
SET
uf
sanoq g) P"pj9 BSJBOD 9ujsn punog
jo pa3eaaA.B
-£9
-------
Figure 64. 2-D vertically averaged of currents In Prince William Sound using
coarse grid (3 hours after simulation).
139
-------
Some velocity vectors are plotted over land. This is due to the
discretization process. Attempting to represent the shoreline with a fixed grid
resulted in some blocks containing both land and water. Velocity vectors are
plotted at the center of the block, occasionally resulting in the vector being
plotted on land. This phenomenon occurs most often in the coarse grid due to
its poorer representation of the actual shoreline. (This grid has also been
slightly rotated from its original position in order to match the alignment of
the fine grid.)
4.5.4. Discussion
In comparing the results from the 2-D and 3-D runs utilizing the coarse
grid (Figures 64-66), we note two factors. First, the flow direction varied
with time and space in a similar manner for each run. Thus, at any given point
in time, the current pattern for the sound was the same for the 2-D and 3-D runs.
Second, the magnitude of the flow velocities at any given point in time and space
varied little for these runs. Only at the third hour (Figures 64 and 67) can
any detectable difference be discerned. At this time, the velocities at the two
main entrances to the sound are slightly greater for the 2-D run. Thus, it
appears that the 2-D and 3-D options differ most in the calculated magnitudes
of the flow field with little or no impact on the directionality.
A comparison of the 2-D coarse grid and fine grid results (Figures 64-69)
yields a similar conclusion. Due to the longer time step and better resolution
of the fine grid, it took slightly longer for the flow field to stabilize so no
comparison could be made in the first hour. At subsequent hours, however, it was
again seen that the flow direction varied with time and space in a similar manner
for each run, resulting in similar current patterns. Also, some difference was
noted in the magnitude of the flow velocities at any given point in time and
space. One must compare the velocities with care as the greater density of the
vectors in the fine grid plots tends to exaggerate any differences in magnitude.
By examining individual, corresponding velocity vectors in each run, it may be
seen that within the two main entrances to the sound, only a small difference
occurs whereas, on either side of the entrances, the differences are greater.
At a given time for the fine grid run, the flow gains speed more rapidly
approaching the entrances, reaches a slightly greater maximum velocity within
the entrance, and loses speed more slowly on the other side. The fine grid, as
would be expected, provided a much more detailed visualization of the flow field,
including the representation of flows through several smaller entrances that
parallel Montague Strait. Thus, it appears that the velocity magnitude is
affected by the degree of resolution of the finite difference grid while the
directionality remains unchanged.
140
-------
Figure 65. 3-D simulation of currents in Prince William Sound
using coarse grid (1 hour after simulation).
141
-------
1 I I
Figure 66, 3-D simulation of currents In Prince William Sound
using coarse grid {2 hours after simulation),
-------
Figure 67. 3-D simulation of currents in Prince William Sound
using coarse grid (3 hours after simulation).
143
-------
» t
t«
»»
»tti»***
i,»<» »i» 11
» *'it * i i 11
I I M t t t I l
4, 44 *i t 4 t 1
« -i4 »,
4 *
* 4 1J4*-
Figure 68. 2-D vertically averaged of currents in Prince William
Sound using fine grid (2 hours after simulation).
144
-------
...I
Figure 69. 2-D vertically averaged of currents in Prince William
Sound using fine grid (3 hours after simulation).
145
-------
4. 6 Currents in Mississippi Sound
In this application, an earlier version of the model has been applied to
simulate the tide- and wind-driven currents in Mississippi Sound and adjacent
continental shelf waters of the Gulf of Mexico (Sheng, 1983).
4.6.1 Physical Setting
Mississippi Sound and adjacent areas (Figure 70) is a region that has
received greater attention due to increasing utilization of its resources,
including the dredging of shipping channels and the disposal of dredged
materials. The Mississippi River is located at the Western end of the sound
and it dominates flow and sediment transport in area. Other major tributaries
which discharge into the sound include the Pearl, Pascagoula, and Mobile Rivers.
4.6.2 Circulation in Mississippi Sound
The circulation in Mississippi Sound is affected by (1) ocean tides
propagating from the Gulf of Mexico through the sound entrance, (2)
meteorological forcing, and (3) bathymetry and geometry.
Gulf tides in the area consist of the diurnal components Kl, 01, and PI
collectively over the semi-diurnal components M2 and S2, except along the Western
Florida Coast. Platzman (1972) and Hansen (1974) found that the period of the
lowest mode of long gravity waves in the Gulf might be quite close to the diurnal
tide period, hence suggesting a quasi -resonant condition. Reid and Whitaker
(1981) developed a numerical tide model for the Gulf based on the vertically-
integrated linearized, Laplace tidal equations in spherical coordinates to
portray the barotropic response of the Gulf to tidal forcing. Their study on
the Gulf tides may provide a useful option to supply seaward boundary conditions
for this application.
The water level response for a given tidal constituent is usually expressed
in the following form (Shureman, 1941) in terms of the surface displacement f :
f - F(t) A(A,^) cos [ wt + x - G (X,^)] (68)
where X is the longitude, is the latitude, A is the mean amplitude over 18.6
years and G the Greenwich phase or epoch at given position (A,^), is tidal
frequency, x is the astronomical argument, while F is the nodal factor (a slowly
varying function of time). Tides at particular stations are characterized by
A and G for individual constituents. In Sheng's study (1983), A's and G's for
5 constituents (01, Kl, PI, S2 and M2) along the open boundaries of our grid are
supplied from Reid and Whitaker 's model. Surface displacements at the open
boundary stations are determined from a linear combination of those due to the
five tidal constituents.
146
-------
Figure 70. Lateral Numerical Grid Used for Dynamic Simulation of
Coastal Currents within the Mississippi Coastal Waters,
147
-------
4.6.3 Results
Simulations were performed for tide and wind-driven current during September
20 to September 25, 1980 and briefly are described below.
4.6.3.1 Tidal Simulation
In this example, the model was run using water surface displacement
(Equation 68) as the boundary conditions. The model-simulated water surface
displacements at four stations (see Figure 70 for locations) within the
Mississippi Sound are compared with measured data in Figure 71. The measured
data shown on those figures have been filtered such that variations due to short
period oscillations on the order of a few hours or less are not included. Over
the simulation period, diurnal tides dominate over the semi-diurnal tides.
Towards the end of the five-day period, the diurnal tides become somewhat less
predominant while the semi-diurnal tides became gradually more apparent. Good
agreement is found at all stations.
Surface displacement over the coastal area at the end of the third day of
simulation is shown in Figure 72, The results exhibit variation in surface
displacement from nearly zero along the open boundary to -7 cm within the
Mississippi Sound, indicating the phase difference in tide.
In this simulation, a relatively large time step of 12 minutes was used
for both the external and the internal modes. Seven grid points were used in
the vertical direction. A relatively smooth bottom with a roughness length,
Z^, of 0.1 cm was assumed. A parabolic length scale, A, no more than 25% of
the local depth, was assumed in the vertical direction. River inflows from six
rivers were considered: Pearl, Jourdon-Wolf, Biloxi, W. Pascagoula, Pascagoula,
and Mobile,
The tide-driven horizontal currents at mid-depth are shown in Figure 73
for two stations in the Mississippi Sound. Currents on the order of 1 ft/sec
(30 cm/sec) exists at both stations. Again, reasonable agreement is found
between data and model results.
The horizontal velocity field at 1 m depth, after 3 days of simulation, is
shown in Figure 74. Relatively large currents exist at the various tidal inlets
and In the area between Ship Island and Chandeleur Island. Except in these
areas, at this instant of time, bottom shear stress generated by the tidal
currents are generally less than 0.3 dyne/cm2.
148
-------
ftuftvStiftN Bt IB. 44 - STATION I
Figure 71. Transient Variation of Surface Displacements at Four Stations
Within the Mississippi Sound from 9/20/80 to 9/25/80.
149
-------
X 9
H1591SSIPP1 SflUND t FILE s
SURFACE DISPUCEMENt fit TIME = 72.0 H0URS
-7.0
I OWTIUM i tOU
-5. a -3.0 -1.0 O.SO
Figure 72, Surface Displacement Contours Within the Mississippi
Coastal Waters at 0 Hr., 9/23/80.
150
-------
Model Retultt
Dolo
Tmo.
STATION 5
a. e
-o. to
-l . 6
I .
I .»
a.
e.o
STATIONS
Mo
, 4~ 1 » ,-1 * _Jl
Model
. oafa
Reiulft_
_
i
**a.
ze.es STATION 6
i.a
a.*
a.a
-o.fc
-l, o
.(d)
'"«
too.
FIGURE 73, Transient Variation of Mid-Depth Velocities at Two
Stations Within the Mississippi from 9/20/80 to 9/25/80,
151
-------
MISSISSIPPI SWNO : FILE = 11300X4
UV VEL0C1TY HT TIME = 72.0 H0URS RND DEPTH =
l.OM
x 3
2.66E+01 ICM/SECJ
VECTIR
FIGURE 74, Horizontal Velocity Field at 0 hour and 1 in depth, 9/23/80.
152
-------
4.6.3.2 Wind-effect on Tidal-Driven Currents
The results presented above did not contain any wind-driven effect. During
this study, wind data were collected at several meteorological stations
surrounding the Mississippi Sound. The wind during the 5-day period between
September 20 and September 25, 1980, was generally quite mild ( 10 mph) blowing
from the southeast. To examine the effect of wind on the currents, Sheng (1983)
carried out a three-day simulation from September 20 using a uniform wind stress
of 1 dyne/cmz from the southeast. As shown in Figure 75, the southeasterly wind
caused water to pile up within the Mississippi Sound at (I,J)=(22,62), outside
Pascagoula Harbor along the northern shore. The wind resulted in a set-up of
0.4 ft. The wind set-up at (I,J)=(30,56), however, is only 0.2 ft. due to the
shielding effect of the Horn Island.
The influence of wind on the current also depends on the location. Figure
76 shows the along-shore velocity at 2 locations over the 3-day period. At (I,J)
= (33,28), off Cat Island, the presence of the wind did not have an appreciable
effect on the tidal current. At (I,J) - (26,88), within the pass between the
Mississippi Sound and the Mobile Bay, the wind caused significant flow from the
Mobile Bay into the Sound. This resulted in a significantly larger bottom shear
stress, which leads to the reduction in the amplitude of the tide-driven
currents. For detailed information on this application, the reader is referred
to Sheng (1983).
153
-------
ELEVRTICJN RT 22.62
120.
EIl-ElVHTieiN RT 3O.56
(HBunai
1X0.
FIGURE 75. Influence of Wind on Surface Displacements at Two Stations from
9/20/80 to 9/24/80.
154
-------
MID-DEPTH u AT 31
Ha)
- CT . H
1.0
Mb)
.Wind B Tide
-Tide Only
T I M
MfO-DEPTK
MID-DEPTH u
He)
Winrf B Tide i
Tide Only F
MID-DEPTH » AT
CC
(d)
Mo
T I MR I MtBUHS J
Figure 76. Influence of Wind on Mid-depth Horizontal Velocities at Two
Stations from 9/20/80 to 9/24/80.
155
-------
SECTION 5
HYDR03D PROGRAMMER'S GUIDE
5.1 OVERVIEW
This section of the manual provides information for the operation
of the program on a computer system such as the Digitial Equipment Corporation
VAX. This section will also explain the various subroutines in the model
which should facilitate modification of the program for specific application
and design of a specialized input/output by adding new modules. A description
of to the programming aspects of the code will also help users in linking the
hydrodynamic program to water quality modeling packages.
5.2 HARDWARE AND SOFTWARE REQUIREMENTS
At this time, the model is operational only on the DEC VAX
computer systems. The program modifications and test runs have been done on
the VAX and therefore the model operations on the VAX system are described
here. The program code is written in VAX FORTRAN 77 and requires about 3000
blocks of hard disk storage, which increases proportionately with the 2-D or
3-D mode of operation and the length of simulation time. For output in the
graphic forms the CA-DISSPLA graphic software package is used.
5.3 INSTALLATION AND IMPLEMENTATION
Although the program is designed for operation on a VAX computer
system, It can be run with some modifications on other computer systems that
support the VAX FORTRAN programming language. For VAX operation, the supplied
program on tape must be installed on the computer system according to the
instructions in the README file accompanying the program codes. The executable
code should then be tested with the sample input file supplied with the model
and the output compared with the sample output file to ensure that the program
is installed properly on the computer system. If it is desired to modify the
program or add extra subroutines to perform specialized calculations, then the
source code must be re-compiled after the modification and linked before it
can be used in performing hydrodynamic simulations.
5.4 DESCRIPTION OF THE COMPUTER PROGRAM
The model consists of 64 subroutines which enable the code to
perform various tasks in a structured fashion. These subroutines facilitate
156
-------
the input of data to the program, perform the mathematical calculations, and
output the simulation results in either numerical or graphical form. The main
routine supervises the overall model operation. It opens input files, calls
subroutines, closes input files and opens output files. For graphical
presentation of the simulation results the software package DISSPLA is used.
The graphical outputs are the basin topography, temporal and spatial variation
of velocities, elevations, temperature, and salinity. Figure 77 illustrates
the functional relationships among the different modules of the program.
157
-------
EHSM3D
Main Program
X-Y MODE
2-D 3-D X-Y-Z MODE
1 1
X-Z MODE
INPUT
PARAMETERS
INPUT
PARAMETERS
INPUT
PARAMETERS
OUTPUT
PARAMETERS
GRAPHIC
NUMERICAL
Figure 77, Operational Chart of the EHSM3D model.
158
-------
5.5 SUBROUTINE DESCRIPTIONS
This section describes the characteristics of each individual
subroutine of the HYDR03D code.
EHSM3D : Main program that supervises the overall model simulation, as
shown in the flowchart in the previous section.
EHSMAI : Sets the lateral turbulent eddy viscosities on the computational
star used to compute the lateral diffusion terms of the horizontal
(u) velocity in subroutine EHSMD4 (2-D runs only).
EHSMAJ : Sets the lateral turbulent eddy viscosities on the computational
star used to compute the lateral diffusion terms of the horizontal
(v) velocity in subroutine EHSMD4 (2-D runs only).
EHSMAS ; Sets the lateral turbulent eddy diffusivities on the computational
star used to compute the lateral diffusion terms of the water
quality parameters in subroutine EHSMC4.
EHSMAU ; Sets the lateral turbulent eddy viscosities on the computational
star used to compute the lateral diffusion terms of the horizontal
u velocity in subroutine EHSMB4(3-D runs only).
EHSMAV ; Sets the lateral turbulent eddy viscosities on the computational
star used to compute the lateral diffusion terms of the horizontal
(v) velocity in subroutine EHSMB4(3-D runs on]y).
EHSMB3 ; Advances the 3-D velocity fields. Using a vertically implicit
scheme, the horizontal perturbation velocities (u',vf) and
computes. These are then combined with the horizontal vertically
integrated velocities (U,V) to obtain the horizontal velocities
(u,v). The continuity equation then is used to compute the
vertical velocity on both the vertically stretched grid and the
original grid.
EHSMB4 : Computes the explicit advection and horizontal diffusion terms of
the momentum and vertically integrated momentum equations. These
are then saved for use by EHSMB3 and EHSMFF for advancing the
velocity fields (3-D runs only).
EHSMC4 : Computes the explicit advection and horizontal diffusion terms of
the concentration, salinity or temperature equation. These are
then saved for use by EHSMCN, EHSMSA or EHSMTE for advancing the
fields (3-D runs only).
EHSMCN : Advances the concentration field using a vertically implicit
scheme and the explicit terms computed by EHSMC4.
EHSMCS : Sets the field values on the computational star used to compute
the explicit terms of the water quality parameters in subroutine
EHSMC4.
EHSMCU : Computes the coefficients and inverts the matrix for advancing the
water quality parameters,
EHSMD4 : Computes the explicit advection and horizontal diffusion terms of
the vertically integrated momentum equations. These are then
saved for use by EHSMFF for advancing the vertically averaged
velocity fields (2-D runs only).
EHSMDE : Computes the water density field and the baroclinic pressure
gradient terms for the horizontal momentum equations.
EHSMDP : Dumps step number information to a disk file (DUMP.TMP)when
159
-------
switch-1 is sent while the program is running,
EHSMDT : Computes the individual parts of horizontal diffusion terms for
EHSMB4, EHSMC4 and EHSMD4. Entry point EHSMDO computes the basic
diffusion term. Entry point EHSMDF computes the higher order
terms.
EHSMED : Computes the lateral turbulent eddy viscosity and diffusivity
fields. This routine also computes the Richardson number,
square-root of the turbulent energy and turbulent scale fields.
EHSMEX : Advances the external variables (surface elevation and vertically
integrated velocities). The river flows, tidal conditions and
wind stresses are first set. Then using a horizontally implicit
scheme (implicit in the x direction only) the surface elevation is
partially advanced. From this the vertically integrated u
velocity is advanced. Then using a similar horizontally implicit
scheme (implicit in the y direction) the advance of the surface
elevation is completed. From this the vertically integrated v
velocity is advanced.
EHSMEZ : Computes the lateral turbulent eddy viscosity or diffusivity for a
given water column (called by EHSMED).
EHSSFF : Computes the explicit terms for the x- and y-sweeps of the surface
elevation equation (called by EHSMEX).
EHSMPN : A subroutine composed of functions CONCEN and DIFFUS, which
provide values of concentration and diffusion coefficient at grid
points according to grid indices MS and NS.
EHSMGA Tridiagonal matrix inversion routine called by EHSMCU and EHSMZS.
EHSMGR Generates a printer plot of 2-D or 3-D field.
EHSMHC Supervises the computation of the hydrodynamic variables(surface
elevation, vertically integrated velocities, velocities, salinity
and temperature).
EHSMHR : Reads the hydrodynamic variables (surface elevation,vertically
integrated velocities, velocities, salinity and temperature) from
disk. Used to compute concentration fields from a previously made
run,
EHSMIF Provides the initial 2-D and 3-D fields at the beginning of a run.
EHSMIH Initializes the bottom topography fields at the
beginning of a run.
EHSMI1 Initializes the index fields at the beginning of a run.
EHSMIN Supervises the input and initialization of the fields.
EHSMIR Reads the input parameters.
EHSMIS Initialize the salinity field based on quadratic interpolation of
salinity data at up to 10 stations.
EHSMPT To initialized temperature field based on quadratic interpolation.
EHSMIT Computes individual advection terms for subroutines EHSMB4, EHSMC4
and EHSMD4.
EHSMIW Outputs input parameters and initial fields,
EHSMMI Matrix inversion routine for momentum equations (called by
EHSMB3).
EHSMND Computes nondimensional parameters and normalizes initial fields.
EHSMOT Output routine. This routine supervises the output and checks for
160
-------
run termination. If variable time steps are used this routine
computes the maximum change of the controlling variables and the
maximum Courant number and adjusts the time step accordingly.
EHSMRF ; Reads the 2-D and 3-D fields from disk then restarting a run.
EHSMRI : Computes the river flows and advances velocity and salinity fields
at river points (called by EHSMEX).
EHSMRS : Computes velocity and bottom stress residuals in the output
routine.
EHSMSA : Advances the salinity field using a vertically implicit scheme and
the explicit terms computed by EHSMC4.
EHSMSB : Computes salinity value at the open boundaries using a linear time
interpolation.
EHSMSC : Controls the smoothing of fields. This routine may be called by
EHSMED to smooth the lateral turbulent eddy viscosities and
diffusivities (KSMALL. NE. 0) or by ENSMOT to smooth the velocity
fields (ISPAC(8), NE. 0) and/or the water quality parameters
(ISPAC(2). NE. 0).
EHSMBT: Computes temperature values at the open boundaries using a linear
time interpolation.
EHSMSE : Sets the surface elevations and depths for all computational
stars. This routine is called by most routines that need the
surface elevation or depth when ISMALL - 1.
EHSMSM 1-D spatial smoother routine (called by EHSMSC).
EHSMSS Routine called by EHSMOT to check for steady state.
EHSMTB Computes the bottom stress (called by EHSMFF).
EHSMTD Computes the tidal surface elevations and advances the surface
elevation field on tidal points (called by EHSMEX).
EHSMTE Advances the temperature field using a vertically implicit scheme
and the explicit terms computed by EHSMC4.
EHSMTP Generates test output (ITEST flag).
EHSMU4 Computes u velocity at v points (contains EHSMV4 to compute v
velocities at u points, and EHSMW4 to compute v velocities at w
points).
EHSMVI : Sets the velocities on the computational star used to compute the
advection and lateral diffusion terms of the horizontal u velocity
in subroutine EHSMD4 (2-D runs only).
EHSMVJ : Sets the velocities on the computational star used to compute the
advection and lateral diffusion terms of the horizontal v velocity
in subroutine EHSMD4 (2-D runs on]y).
EHSMVS : Sets the velocities on the computational star used to compute the
advection terms of the water quality parameters in subroutine
EHSNC4.
EHSMVU : Sets the velocities on the computational star used to compute the
advection and lateral diffusion terms or the horizontal u velocity
in subroutine EHSMB4 (3-D runs only).
EHSMW : Sets the velocities on the computational star used to compute the
advection and lateral diffusion terms of the horizontal v velocity
in subroutine EHSMB4 (3-D runs only).
EHSMW3 : Generates numerical printout of 3-D fields.
EHSMWR : Generates numerical printout of 2-D fields.
EHSMWS ; Reads surface wind stress from disk and interpolates the
161
-------
surface wind stress field (called "by EHSMEX).
EHSMWW
EHSMXX
EHSMXY
EHSMZE
EHSMZS
Computes the vertical velocity field from the continuity equation
(called by EHSMB3).
This is the external routine name called to set surface elevations
and depths. If ISMALL - 0 then EHSMXX is EHSMZE else it is
EHSMSE.
Computes x and y grids if not read from disk.
Sets the depths for all computational stars. This routine is
called by most routines that need the depths when ISMALL 0.
Computes the matrix coefficients for the surface elevation
equation, inverts the matrix, sets new surface elevation and
computes new vertically integrated velocity field (called by
EHSMEX). In addition to the above subroutines of the HYDR03D
code, the following programs are used to generate the depth arrays
and the various grid indices:
DEPTH_FILE_CREATE ; This program reads the depths at the corner points of
the grid lines and creates the 3 depth arrays for the
U, V and f points.
This program reads the depth file created by
DEPTH_FILE_CREATE and produces the grid indices NS,
MS, JU1, JU2, JV1, JV2, IU1, IU2, IV1, and IV2.
INDEX FILE CREATE
162
-------
5.6 INPUT/OUTPUT UNITS
Unit 4 - This is the main Input file providing the essential input
information via formatted card images that are described in detail
in Section 2.
Unit 6 - This is the file containing the major printouts.
Unit 6 - This is a formatted sequential input/output file which stores the
surface displacements, vertically- integrated velocities and
three-dimensional velocities at selected stations and time
intervals. It is created in EHSMOT by:
WRITE( 8,911) TIME, IT
911 FORMAT (1PE13.6, OPI13)
WRITE(8,912) (S(JST(I),IST(I)), I-l.NSTA)
WRITE(8,912) (UI(JST(I),IST(I), I=1,NSTA)
WRITE(8,912) (VI(JST(I),IST(I), I=1,NSTA)
912 FORMAT (1P10E13.6)
DO 10120 KZ=1,KM
WRITE(8,912) (U(KZ,JST(I),IST(I)) ,I-1,NSTA)
10120 CONTINUE
DO 10130 KZ=1,KM
WRITE(8,912) (V(KZ,JST(I),IST(I)),I-1,NSTA)
10130 CONTINUE
Unit 11 - This input file contains the variable bottom topography provided by
the user. It is an unformatted sequential file containing HU, HV,
and HS each dimensioned as (JM, IM) . It is read in EHSMIR by:
READ (11) HU.HV.HS.
This unit is required if IBTM=3 (IBTM is used in input file).
Unit 12 - This is an unformatted sequential output file containing the grid
parameters NS , MS, JU1, JU2, JV1, JV2 , IU1, IU2, IV1, and IV2. It
is read in EHSMII by:
READ(12) NS,MS
READ(12) JU1,JU2,JV1,JV2
READ (12) IU1,IU2,IV1,IV2
Unit 13 - This is an unformatted sequential file that stores the major
species concentration data at desired time instants. It is created
in EHSMDT by:
WRITE ( ICONC ) TIME , IT , FNAME ( 5 ) , FNAME ( 6 ) , IM , JM , KM , XREF ,
ZREF.UREF.COR
WRITE(ICONC) XS,XU,YS,YV,HU!HVJHS
WRITE (I CONG) C
Unit 14 - This file contains user -generated non-unifora grid when IGRID=1.
It is created by
WRITE (14) XU
WRITE(14)YV
Unit 16 - This is a formatted sequential file that contains the run number
and two indices needed for restarting a run:
Unit 18 - This is a formatted sequential file that stores the salinity at
selected stations and time intervals. It is created in EHSMOT by:
WRITE(18,911) TIME, IT
DO 10140 KZ - l.KM
WRITE(18,912) (SA(KZ,JST(I),IST(I)), I-l.NSTA)
10140 CONTINUE
163
-------
Unit 20 - This is the sequential file that contains the river inflow data at
selected time instants. It is read in EHSMRI by:
DO N - l.NRIVER
READ(20,*) IDAY,IHOUR,URIVER(N),VRIVER(N)
END DO
Unit IRQ - This is an unformatted sequential file that contains all the
necessary information when initiating or restarting the run. The
structure of this file can be found in subroutine EHSMRF.
|Jnit IW - This is an unformatted sequential file that stores the major flow
output data at desired time instants. The structure of this file
is similar to that of unit IRD file and is created in subroutine
EHSMOT by:
WRITE(IW)TIME1IT,FNAME(3),FNAME(4))IM,JM1KM>
XREF,ZREF,UREF,COR,AVO
WRITE(IW) XS, XU, YS, YV, HU, HS, FMU, FWV, FMS, FMSV
WRITE(IW) U, V, W, WW
WRITE(IW) UI, VI
WRITE(IW) S
WRITE(IW) T
WRITE(IW) SA
WRITE(IW) GA,GB
WRITE(IW) TBX.TBY
WRITE(IW) QQQ, SL
Unit IWS - This is the output file for storing residual flow data and contains
the same variable groups as the unit IRD and IRW files.
Unit IR4 - This is the sequential file that contains the wind stress field at
selected time instants.
164
-------
SECTION 6
CONCLUSIONS AND RECOMMENDATIONS
This report documents existing and developing programmatic needs in the
U.S. EPA for a stratified flow model to simulate complex flows in lakes,
estuaries, harbors, and coastal waters. This model is needed to assist in
ecological assessments, risk assessments, and exposure predictions for
dissolved and sediment-bound contaminants. The model is needed for estuary
studies to support the National Estuary Studies and wasteload allocations.
There is a need to assist in the clean up of contaminated sediments in the
Great Lakes and other lakes, A number of research programs ranging from oil
spill initiatives to investigation of the effect of global climate change
could benefit from a simulation tool designed to determine circulation and
transport in lakes, estuaries, and near coastal marine waters.
To meet these well defined needs, the U.S. EPA Environmental Research
Laboratory located at Athens, Georgia, has worked cooperatively with others to
document the hydrodynamics model, HYDR03D. This model is a dynamic three
dimensional circulation model. The model simulates water circulation,
dissolved solids or salinity, water temperature, and a dissolved species
concentration. The model will also be used with the SED3D sediment
resuspension and dispersion model due to be completed in FY 1991 (see
Preface).
HYDR03D has been tested in a number of estuaries and lakes. In this
documentation, the model is used to simulate diverse water bodies that include
Prince William Sound in Alaska, Suisun Bay of the San Francisco Bay, Charlotte
Harbor in Florida, Green Bay of Lake Michigan, and the Mississippi Sound in
the Gulf of Mexico. In addition, these illustrative examples and other
hypothetical cases are reviewed to demonstrate the validity of the code and
the flexibility of the program to simulate different conditions.
This documentation provides other important elements to aid the user as
well as establishing the validity and flexibility of the program. This report
reviews the data required, and the form that the data must be transformed
into. It reviews the structure of the program and provides information about
the derivation of the governing equations that form the basis of the model.
From all of this, one can conclude that a useful and necessary tool is
available to support U.S. EPA studies and other environmental investigations.
It should be noted, however, that this is a complex model that may
require assistance beyond that available in the manual. It is recommended
that potential users, including program managers and applications experts,
consult the Introduction (Section 1) and introduction to the major sections
for guidance on how best interpret and use this document. For further
assistance, contact the Center for Exposure Assessment Modeling (CEAM). CEAM
165
-------
can assist in the design of studies involving stratified flows, aid in the
development of data collection programs, and provide expert assistance in
implementing and interpreting the results for Superfund investigations and a
number of other different types of studies.
166
-------
REFERENCES
Ambrose, R.A., T.A. Wool, J.P. Connolly, R. W. Schanz, 1988. WASP4, A
Hydrodynamic and Water Quality Model--Model Theory, User's Manual,
and Programmer's Guide, EPA/600/3-87/039, U.S. Environmental Protection
Agency, Athens, Georgia.
Ambrose, R.A., Jr., and J.L., Martin, 1990. Technical Guidance Manual for
Performing Waste Load Allocations, Book III: Estuaries, Part 1, Estuaries
and Waste Load Allocation Models, U.S. EPA, Office of Water, Washington,
D.C.
Ambrose, R.A., Jr., Martin, J.L., and McGutcheon, S.C., 1990, Technical
Guidance Manual for Performing Waste Load Allocations, Book III:
Estuaries, Part 2, Application of Estuarine Waste Load Allocation Models,
U.S. EPA, Office of Water, Washington, D.C.
ASCE Task Committee on Turbulence Models in Hydraulics Computations, 1988.
Turbulence modeling of surface water flow and transport: Parts I to V,
Journal of Hydraulics Engineering, American Society of Civil Engineers,
114(9), 970-1073.
Bird, R.B., R.C. Armstrong, and 0. Hassager, 1977. Dynamics of Polymeric
Liquids, Vol. 1, Fluid Mechanics, Wiley, New York.
Blumberg, A. F. 1975. A Numerical Investigation into the Dynamics of Estuarine
Circulation. Chesapeake Bay Institute Report No. 91, Baltimore MD.
Blumberg, A.F., 1986. Turbulent Mixing Processes in Lakes, Reservoirs, and
Impoundments, in Physics-Based Modeling of Lakes, Reservoirs, and
Impoundments, W.G. Grey, ed., American Society of Civil Engineers, New
York.
Bowden, K. F. and P. Hamilton. 1975. Some Experiments with a Numerical Model
of Circulation and Mixing in a Tidal Estuary. Est, Coastal Mar. Sci., 3,
281.
Bowie, G.L, W.B, Mills, D.B, Porcella, C.L. Campbell, J,R. Pagenkopf, G.L.
Rupp, K.M. Johnson, P.W.H. Chan, and S.A. Gherini, 1985. Rates,
Constants, and Kinetics Formulations in Surface Water Quality Modeling,
2nd Edition, EPA/600/3-85-040, U.S. Environmental Protection Agency,
Athens, Georgia,
167
-------
Butler, H. L. 1978, Coastal Flood Simulation In Stretched Coordinates. Proc.
16th Coastal Engineering Conference, ASCE, Hamburg, Germany.
Chapra, S.C., and K.H. Reckhow, 1983. Engineering Approaches for Lake
Management, Volume 2: Mechanistic Modeling, Butterworth Publishers,
Boston.
Cheng, R. T., and T. J. Conomos, 1980. Studies of San Francisco Bay by the
U.S. Geological Survey, Institute of Environmental Sciences Proceedings,
299-303.
Cheng, R. T., and J. W. Gartner, 1984, Tides, Tidal and Residual Currents in
San Francisco Bay, California - Results of Measurements, 1979-1980. 5-
Part Report, U. S. Geological Survey Water Resources Investigation Report
84-4339.
Cheng, R.T., T.M. Powell and T.M, Dillon, 1976. Numerical Models of Wind-
Driven Circulation in Lakes, Appl. Math. Modeling 1, 141-159.
Conomos, T. J. and others. 1978. Field and Modeling Studies of San Francisco
Bay. Coastal Zone '78 Proceedings, ASCE.
Cox, R. A. and others. 1967. The Electrical Conductivity/Chlorinity
Relationship in Natural Seal Water. Deep-Sea Research, 14, 203-220.
Csanady, G.T., 1975. Hydrodynamics of Large Lakes, Ann, Rev, Fluid Mech,, 7,
357-386.
Doneker, R.L,, and G.H. Jirka, 1968, CORMIX1: An Expert System for Mixing Zone
Analysis of Convential and Toxic Single Port Discharges, U.S. EPA Report
EPA/600/3-88/013, Athens, Georgia.
Ekeart, C., 1958. Properties of Water, Part II. The Equation of State of
Water and Sea Water at Low Temperature and Pressure, American Journal of
Science, 256, 225-240.
Fischer, H.B., List, E.J., Koh, R.C.Y., Imberger, J., and Brooks, N.H., 1979.
Mixing In Inland and Coastal Waters. Academic Press, New York.
Goldstein, S., 1960. Lectures on Fluid Mechanics, Interscience London.
Grey, W.G., ed., 1986. Physics-Based Modeling of Lakes, Reservoirs, and
Impoundments, American Society of Civil Engineers, New York.
Hansen, K.W., 1974. Calculation of Normal Modes for the American
Mediterranean Seas, Technical Report No. 26, University of Chicago.
Heaps, N.S., C.H. Mortimer, and E.J. Fee, 1982. Numerical Models and
Observations of Water Motion in Green Bay, Lake Michigan. Phil. Trans.
Royal Society of London, Series A, 306, 371-398.
168
-------
Henderson-Sellers, B., 1982. A simple formula for vertical eddy diffusion
coefficients under conditions of non-neutral stability, Journal of
Geophysical Research, American Geophysical Union, 87(C8), 5860-5864.
Henderson-Sellers, B., 1984. Engineering Limnology, Pitman, London.
Hinze, J.O., 1959. Turbulence, 2nd Ed., McGraw-Hill , New York.
Jirka, G.H., 1982. Multiport Diffusers for Heat Disposal - A Summary, American
Society of Civil Engineers Journal of the Hydraulics Division, Vol. 108.
Johnson, B., Y.P. Sheng, K. Kim, and R, Heath, 1989. Development of a Three
Dimensional Hydrodynamic Model of Chesapeake Bay, Proceedings of the
American Society of Civil Engineers Conference on Estuarine and Coastal
Modeling, H.L. Spaudling and J.C. Swanson, Eds., American Society of
Civil Engineers, In Press.
Kent, R. E., and D. W. Pritchard, 1959. A Test of Mixing Length Theories in a
Coastal Plain Estuary, J. Mar, Res,,18, 62-72.
Kinnetics Laboratories, Inc., 1981. In-Site Field Data Gathering Stations San
Francisco-Delta Salinity Intrusion With Navigation Channels. KLI-81-1,
Final Report to San Francisco District, U. S. Army Corps of Engineers.
Lai, C., Flows of Homogeneous Density in Tidal Reaches; Solution by Method of
Characteristics, U.S. Geological Survey Open File Report, 1965.
Leendertse, J.J., and S.-K. Liu, 1975. A Three-Dimensional Model for Estuaries
and Coastal Seas: Volume II, Aspects of Computation, Rand Corp., Report
R-1764-OWRT, prepared for the Office of Water Research and Technology,
U.S. Department of Interior, Santa Monica, California.
Lynch, D. R. and W.G., Gray, 1978. "Analytic Solutions for Computer Flow Model
Testing" Journal of the Hydraulics Division, American Society of Civil
Engineers, 1409-1428.
Lynch, D.R., 1986, Basic Hydrodynamic Equations for Lakes, in Physics-Based
Modeling of Lakes, Reservoirs, and Impoundments, W.G. Grey, ed., American
Society of Civil Engineers, New York.
McCormick, M.J., and Scavia, D., Calculation of vertical profiles of
lake-averaged temperature and diffusivity in Lakes Ontario and
Washington, Water Resources Research, 17(2), 305-310.
McCutcheon, S.C., 1983. Vertical mixing in models of stratified flow, in
Frontiers in Hydraulic Engineering, Proceedings of the Hydraulics
Specialty Conf., American Society of Civil Engineers, Cambridge, Mass.
15-20.
McCutcheon, S.C., 1989. Water Quality Modeling, Volume I Transport and Surface
Exchange in Rivers, CRC Press, Inc., Boca Raton, Florida,
169
-------
McCutcheon, S.C., D.-W. Zhu, Bird, S., 1990. Chapter 5 Model Calibration,
Validation and Use, Technical Guidance Manual for Performing Waste
Load Allocations, Book III: Estuaries, Part 2: Application of Estuarine
Waste Load Allocation Models, Ambrose, R.A,, Jr., Martin, J.L. and
McCutcheon, S.C., Eds., U.S. Environmental Protection Agency. Office of
Water, Washington, D.C,
Miller, G.S., and J.H. Saylor, 1985. Currents and Temperature in Green Bay,
Lake Michigan, J. of Great Lakes Res,, 11(2), 97-109
Moldin, R., and A.M. Beeton, 1970. Dispersal of Fox River Water in Green Bay,
Lake Michigan, Proceedings of Thirteenth Conference on Great Lakes
Resources, International Assoc. Great Lakes Res., 468-476
Monin, A.S., and A.M. Yaglom, 1971. Statistical Fluid Mechanics; Mechanics of
Turbulence, Vol. 1, The MIT Press, Cambridge.
Munk, W. H. and E. P. Anderson. 1948. Notes on the Theory of the Thermocline,
J. Mar. Res,, I, 276-295.
Prandtl, L. 1925. Bericht uber untersuchungen zuraus-gebildete turbulenz, Zs,
Angew. Math. Mech., 5(2), 136-139.
Patchen, R. C., and R. T. Cheng. 1979. Current Survey of San Francisco Bay by
NOS/USGS, EOS, American Geophysical Union, 60, p. 843.
Platzman, G.W., 1972. Two-Dimensional Free Oscillations in Natural Basins,
J. Phys. Oceanography, 2(2), pp. 117-130.
Reid, R.O., and R.E. Whitaker, 1981. Numerical Model for Astronomical Tides
in the Gulf of Mexico, Vol. I: Theory and Application, Dept.
Oceanography, Texas A & M University, College Station, TX.
Reynolds, A.J., 1974. Turbulent Flows in Engineering, A Wiley, London, New
York.
Richardson, L.F., The supply of energy to and from atmospheric eddies,
Proceedings of the Royal Society, Series A, 97, 354-373.
Roberts, P.J.W., 1979. Line Plume and Ocean Outfall Dispersion, American
Society of Civil Engineers Journal of the Hydraulics Division, Vol. 105.
Rodi, W., 1980. Turbulence Models and Their Application in Hydraulics,
International Association for Hydraulic Research, Delft, The Netherlands.
Rodi, W. , 1984. Examples of turbulence model applications, in Simulation of
Turbulence Models and Their Applications, Vol. 2, Collection de la
Direction des Estudes et Recherches, Electricite de France, editions
Eyrolles, Paris, France.
Schlichting, H., 1979, Boundary-Layer Theory, McGraw-Hill, New York.
170
-------
Sheng, Y.P., 1975. Wind-driven Currents in Contaminant Dispersion in The
Nearshore of Large Lakes, Ph.D. Dissertation, CASE Western Reserve
University, Cleveland, OH.
Sheng, Y.P., 1980. Modeling Sediment Transport in Shallow Waters, Estuarine
and Wetland Processes with Emphasis on Modeling, P. Hamilton and K.B.
Mocdonald, Ed.s, Plenum Press, 299-337.
Sheng, Y.P., 1981. Modeling the Hydrodynamics and Dispersion of Sediments in
the Mississippi Sound, Technical Report No. 455, Aeronautical Research
Associates of Princeton, Princeton, NJ.
Sheng, Y.P., 1982. Hydraulic Applications of a Second-Order Closure Model
of Turbulent Transport, in Applying Research to Hydraulic Practice, P.
Smith, Ed., American Society of Civil Engineers, New York, 106-119.
Sheng, Y.P., 1982. A Review on LA-LB Harbor Study, Technical Memo No. 82-11,
Aeronautical Research Associates of Princeton, Princeton, NJ.
Sheng, Y.P., 1983. Mathematical Modeling of Three-Dimensional Coastal Currents
and Sediment Dispersion: Model Development and Application, Technical
Report CERC-83-2, U.S. Army Corps of Engineers Waterways Experiment
Station, Vicksburg, MS.
Sheng, Y.P,, 1984. A Turbulent Transport Model of Coastal Processes,
Proceedings 19th International Conference on Coastal Engineering,
American Society of Civil Engineers, 2380-2396.
Sheng, Y.P., 1986a. A Three-Dimensional Hydrodynamic Model in Generalized
Curvilinear Coordinates, Technical Report 587, Aeronautical Research
Associates of Princeton, Princeton, NJ.
Sheng, Y.P., 1986b, Finite Difference Models for Hydrodynamic Lakes and
Shallow Seas, in Physics-Based Modeling of Lakes, Reservoirs, and
Impoundments, W.G. Grey, ed., American Society of Civil Engineers, New
York.
Sheng, Y. P., 1987. On Modeling Three-Dimensional Estuarine and Marine
Hydrodynamics, in Three-Dimensional Models of Marine and Estuarine
Dynamics, J. C. J. Nihoul and B. M. Jamart, Eds., Elsevier, 35-54.
Sheng, Y.P. and W. Lick, 1979. The Transport and Resuspension of Sediment in
a Shallow Lake, Journal of Geophysical Research, Vol. 84, pp. 1809-1826.
Sheng, Y. P. and W. Lick. 1980. A Two-Node Free-Surface Numerical Model for
the Three-Dimensional Time Dependent Current in Large Lakes, U.S.
Environmental Protection Agency, Duluth MN., EPA-600/3-80-047, 62.
171
-------
Sheng, Y. P., W. Lick, R. T. Gedney, and F, B. Molls. 1978. Numerical
Computation of Three-Dimensional Circulation in Lake Erie; A comparison
of a Free-Surface Model and a Rigid-Lid Model, J. Phys. Oceanography, 8,
713-27.
Sheng, Y. P., and H. Butler, 1982. Modeling Coastal Currents and Sediment
Transport, Proc. IBth Coastal Engineering Conference, American Society of
Civil Engineers Capetown, South Africa, pp. 1127-1148.
Sheng, Y.P. and S.S. Chin, 1986. Tropical Cyclone Generated Currents, Proc. of
the 20th International Conference on Coastal Engineering, American
Society of Civil Engineers, 737-751.
Sheng, Y. P., S.F. Parker and D.S. Henn. 1986. A Three-Dimensional Estuarine
Hydrodynamic Software Model (EHSM3D), Aeronautical Research Associates of
Princeton, Inc. Princeton, NJ 08543-2229.
Sheng, Y.P., J.K, Chou and A.Y. Kuo, 1989a, Three Dimensional Numerical
Modeling of Tidal Circulation and Salinity Transport in James River
Estuary, Proceedings of the American Society of Civil Engineers
Conference on Estuarine and Coastal Modeling, M.L. Spaudling and J.C.
Swanson, Eds., American Society of Civil Engineers, In Press.
Sheng, Y.P., H.K. Lee and K.H. Wang, 1989b. On Numerical Strategies of
Estuarine and Coastal Modeling, Proceedings of the American Society of
Civil Engineers Conference on Estuarine and Coastal Modeling (M.L.
Spaudling and J.C. Swanson, Eds.), American Society of Civil Engineers,
In Press.
Sheng, Y.P., V. Cook, S. Peene, D. Eliason, P.F. Wang and S. Schofield, 19B9c,
A Field and Modeling Study of Fine Sediment Transport in Shallow Waters,
Proceedings of the American Society of Civil Engineers Conference on
Estuarine and Coastal Modeling, M.L. Spaudling and J.C. Swanson, Eds.,
American Society of Civil Engineers, In Press.
Shureman, P., 1941. Manual of Harmonic Analysis and Prediction of Tides.
U.S. Department of Commerce, Coast and Geodetic Survey, Special Pub. No,
98.
Smith, P.E. and R.T. Cheng, 1989. Recent Progress on Hydrodynamic Modeling of
San Francisco Bay, Proceedings of the American Society of Civil Engineers
Conference on Estuarine and Coastal Modeling, M.L. Spaudling and J.C.
Swanson, Eds., American Society of Civil Engineers, In Press.
Smith, R. A., 1980. Golden Gate Tidal Measurements, 1854-1978. ASCE J. of
Waterways, Port, Coast and Ocean Division, ,106, WW3, pp. 401-410.
Tennekes, H. and J.L. Lumley, 1972. A First Course in Turbulence, The MIT
Press, Cambridge. Mass.
Turner, J.S., 1973. Buoyancy Effects in Fluids, Cambridge University Press.
172
-------
U.S. EPA, 1985. Technical Support Document for Water Quality-based Toxics
Control, Office of Water, Washington, D.C.
White, P.M., 1974. Viscous Fluid Flow, McGraw-Hill, New York.
Wright, S.J., 1977. Mean Behavior of Buoyant Jets in a Crossflow, American
Society of Civil Engineers Journal of the Hydraulics Division, Vol. 103.
Zakikhani, M., S. L. Bird, and S. C. McCutcheon, 1989. Hydrodynamic and
Sediment Transport Modeling of Green Bay, Lake Michigan.
Presented at 32nd Conference of International Association for Great
Lakes Research, Madison, Wisconsin.
173
-------
APPENDIX A
LIST OF SYMBOLS AND MAJOR VARIABLES
Symbol
T
U
V
P
V » KV ' Dv
AH s Kg , DH
U
V
u
w
T
S
C
h
^sx i rsy
rbx'rsy
X
y
a
Mx
Pi
FORTRAN
label
S
UI
VI
R, RU
GA, GB
AH
U
V
w
ww
T
SA
C
HU,HV,HS
TX, TY
TBX, TBY
XS, XU
YS, YV
Z, SG
FMS, FMU
FMSV, FMV
Array Size
JM, IM
JM, IM
JM, IM
KM, JM, IM
KM, JM, IM
JM, IM
KM, JM, IM
KM, JM, IM
KM, JM, IM
KM, JM, IM
KM, JM, IM
KM, JM, IM
KM, JM, IM
JM, IM
JM, IM
JM, IM
IM
JM
KM
IM
JM
Definition
Surface Displacement
Vertically- Integra ted velocity
Vertically- integrated velocity
Density
Vertical eddy coefficients
Lateral eddy coefficients
Velocity in x direction
Velocity in y direction
Vertical velocity in f direction
Vertical velocity in Z direction
Temperature
Salinity
Species Concentration
Depths
Wind Stresses
Bottom Stresses
X Position of f and u points
Y Position of f and v points
Z Position of u and w points
X-stretching coefficients at f and u
points
Y- stretching coefficients at f and v
points
174
-------
APPENDIX B
DERIVATION OF THE GOVERNING EQUATIONS
FOR A a-STRETCHED COORDINATE SYSTEM
The governing equations for hydrodynamics can be expressed in terms of
the a-stretched coordinates. The following is an illustration of the
derivation for the continuity equation.
The continuity equation in a stretched coordinates can be derived from
the continuity equation in Cartesian coordinates (x,y,z), a definition of the
a coordinate system, and the chain rule of differentiation. The continuity
equation in Cartesian coordinates (x,y,z) is:
au dv 9w n , ,.
+ + - 0 (A-l)
3x 3y 8z
The <7-stretched coordinate system is defined as follows;
z-f(x,y;t)
a(x,y,x;t) = i (A-2)
H(x,y;t)
where H(x,y;t) -= h(x,y) + T(x,y;t) is the total instantaneous water depth,
The chain rule is:
8_ d_ da 3_
fiy" dy dy do
_ _ 4. _ fA ^
9y " ay 3y flf l }
dz
w =
dt
w - (f (x,y;t) -I- CT(x,y,z;t) H(x,y;t>
at
dr ah ch
w -= Hu + (1+a) -^ + CT(U 4 a ) (A- 5)
dt crx oy
175
-------
where
da
«-dr (A-6)
Substituting Equations 2 through 6 into Equation 1, we obtain the
continuity equation in the new coordinate system (x,y,f)
8jL + i. (HU) + *- (Hv) + - 0 (A-7)
at ox 3y da
Non-dimensionalization is based on the definition of the following
nondimensional variables:
(x*,y*) - (x,y)/xr)
(u*,v*) - (u.v)/^) (A-8)
w* - wxr/ur
f* -
t* - tf
Substituting Equation (A-8) into (A-7) yeilds the continuity equation in the
non-dimensional form:
9Hv 3w y nx
+ jg H - - 0 (A-9)
at ^ ax ^ ay
The Equation A-9 is the same as Equation 13 in Section 2.2.6 with ^x = py
= 1. Following a similar approach, may obtain the non-dimensional forms of
the momentum equations in the (x,y,cr) coordinates as given in Equations 14 to
17 in Section 2.2.6.
176
-------
APPENDIX C
CHARACTERISTIC TIME SCALES OF VARIOUS
PHYSICAL PROCESSES IN ESTUARIES
Physical Time Order of
Process Scale Magnitude
Periodic Forcing tf 1/w
Convection tc ^r/^r
Inertia Oscillation ti 1/f
Vertical Turbulent
Diffusion t^, t^, tvds Z,2, A^, Zrz/K,f, Zt2
Lateral Turbulent
Diffusion tedml ttdh, ttds Xj-VAHr.V/Kflr.XrVDHr
Gravity Wave tge Xr/gZr
Internal
Gravity Wave tgi Xr/A/5/p0
Ekman Layer
Diffusion te Zr/2f
177
-------
APPENDIX D
SAMPLE INPUT AND OUTPUT FOR THE SIMULATION
OF WIND-DRIVEN CURRENTS IN AN ENCLOSED BASIN
In this appendix a sample input/output is described to be used to examine
the code. This test example solves a wind-driven currents in an enclosed
basin of 50 km x 50 km. The depths at the north and south are 3 meters and
vary linearly to 10 meters at the central region of the basin. A uniform wind
with speed of 9m/sec (1 dyne/cm2) is blowing from east to west (toward the
negative x-axis) 1 dyne/cm2 starts blowing from a zero initial condition. The
simulation stops after 24 hours. The domain is divided into 10 x 10 uniform
grid cells in the horizontal plane. The local depth is divided into 5 layers
with equal length. The vertical eddy viscosity is assumed to be constant and
5Cn2/sec. The boundary conditions at the lateral boundaries and bottom are
assumed to be 'no slip' condition. Starting with the action of the east wind,
the surface elevation reaches steady state within 24 hours (Figure 12), In
this run the Coriolis effect is not considered and since the basin geometry is
symmetrical, the responses of currents and surface elevation also exhibit
symmetrical characteristics.
178
-------
, , * *
4 ^ -* -
* -* -* -
i 4 * ^ > *
3-D : Slgmo - -0.9
-4~- «--
«,
/ *<*" ^ 1
v»4::: *''
> .,.^_<Ii-
- - ' -^-*-~»-«^4^l
* i
t
3-0 » Sigma - -0.5
<~^_^
3-D : Sigma - -0.I
4""*
M*
V I
\«~,«-
«~*
1 -*»*»*-*-» »
> -»-» 1_*-»-* -. i
3-D ! Hcon
42.6 CM/SEC
Fieure 78. Steady-state Wide-driven currents in an. enclosed square basin of
50 km on each side; linearly varying bottom from 3 m (South and
North) to 10 m (at center)
179
-------
o
LJ
ill
LJ
O
ft
L.,
IE
LJ
ct
(t
o
LJ
LJ UJ
O V
(T ^
u. n
cc o
io ~
o;
cr.
₯
u
CC
IX
R r,
u
in
C
O
o
m
LC
(C
n
LJ
WIND-DRIVEN CURRENTS IN RN ENCLOSED BflSIN
V
a 12 is
TIME (HOURS)
20
24
\.
-A4-
.-4
a 12 is
TIME (HOURS!
20
24
4
t 12, B)
< (6,6)
i 110,6)
12 16
TIME (HOURS)
20
24
Figure 79. Three-dimensional simulation of wind-driven currents in an
enclosed basin; results are for three grid points of (2,6),
(6,6) and (10,6)
180
-------
EXAMPLE INPUT DATA FILE
FOR AN ENCLOSED BASIN
-------
TR
#1 ISTART(I4),TITLE(A64) ***TITLE CARD***
0 3-D RUN FOR SAMPLE RUN WITH WIND=-1 DYNE/CM**2 (1 DAY)
#2 XREF ZREF UREF COR GR ROO ROR TO
*PHYSICAL CONSTANTS*
1000000. 500. 10. .00009 981. 1. 1.001 0. 20.
#3 IVLCY IFI IFA IFB IFC IFD ICCl(icon) ICC2(isal) ICC3(itemp)
ICC4(isedi)
1100010 0 0
0
#4: BVR SI S2 PR PRV TWE TWH FKB TQO **TEMPERATURE
PARAMETERS
1. 0. 0. 1. 1. 0. 0. 5. 0.
#5: IVER ICON IUBO IBL IBR JBM JBP
IC1IC2JC1JC2ID1ID2JD1JD2*CONCEN. PARAMETERS*
2 4 0 1 26 1 29 1. 10000. 1
0 0
#6: IEXP IAV AVR AVI AV2 AVM AHR
PARAMETERS***
CREF CMAX CO
000000
***TURBULENCE
#6As
f6B:
#7:
0
FM1
-.5
QCUT
.15
IWIND
0 5.
FM2
-1.5
ICUT
1
TAUX
0.
ZTOP
-.05
GAMAX
300.
TAUY
0 . 1 . :
SLMIN QQ
1.
GBMAX FZS
300. .10
***WIND P
10<
MI:
01
AR
KSMALL
0
-1.
0.
#B:ISMALL IBTM ITB HADD HMIN ZREFBN CTB
B C
1 250. 1. 5. .004 .4
BZ1
HI
,25
H2 **VERT
.75
#8A:ZREFTN TZ1 SSSO ** MORE VERT B.C.**
5. .4 .0
#9: ITIDE IOPEN JWIND IJLINE ***LATERAL B.C. FLAGS***
0 0000 0 0
#10: IJGAGE IJDIR IJROW IJSTRT IJEND ***IF(IJLINE.GT.0), FOR EACH
IJLINE***
#11 (1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
***ISPAC(I),1=1,10***
000-100000 0
#12 (1) (2) (3) (4) (5) (6) (7)
***JSPAC(I),1=1,10***
01-100
#13 (1) (2) (3) (4)
***RSPAC(I),1=1,10***
0. .0 100. 1. .25 4,
0010
(5) (6) (7)
(8) (9) (10)
0
(8) (9) (10)
.020 .1 .00001 -.0001
#14: ISTEP IT1 IT2 ITS DELT DELTMIN DELTMAX EPSILON BUFAC WTS WTU
WTV*TIMESTEP*
-------
0 1 144 1 600. 1. 900. .075 10. 1. 1.
1.
#15: ITEST IP1 IP2 IP3 IPU IPW IPA IPB ID JPA JPB JD KPA KPB
KD**PRINTOUT INFO**
3 72 72 72 1 1 1 11 1 1 11 1 1 5 4
#16: IGI IGH IGT IGS IGU IGW IGC IGQ IGL IGR IGRI IGTB**PRINTOUT
FORMAT FLAGS***
1111111-1-1-1-1-1
#17: IRD IW IWR ICI IWC ICO ISED IWS IREAD IR4 ***DISCFILE
INFO***
99100 0141 0 15
#18: FNAME(6A4) ***3 DISCFILE
NAMES(FLOWIN f FLOWOUT,CONCENTRATION)* * *
6 6 sal 6 6 sal 6 6 sal
#197(1) "(2) (3) (4) (5) (6) (7) (8) (9) (10)
***TBRK(I),TIMEBREAKS***
240. -480. -480. -480. -480. -480. -480. -480. -480. -480.
#20: NSTA NRANGE NFREQ ***TIMEFILE GAGE STATIONS***
50 2
#21+:IST(K) JST(K) KST(K) STAT ID(K) ( 3 I 4 ,A48 )
***IF(NSTA.GT.O),STATION INFO***
261
461
661
861
10 6 1
#22: NRIVER ***NUMBER OF RIVERS***
0
#22A: IRIVER(K) JRIVER(K) LRIVER(K) URIVER(K) VRIVER(K)
***IF(NRIVER.GT.O)***
#23: (SAB(K),K=1,KM)**VERTICAL SALINITY PROFILE ALONG W,S,E,N (IF
ISALT.NE.O)**
#23A: NISS **NUMBER OF STATIONS WITH INITIAL SALINITY DATA
#23B:(ISS(N),JSS(N),NDEPTH(N),TDEPTH(N),N=1,NISS)**I,J,NO. OF
PTSfTOTAL DEPTH
#23C: (CONT'D) / DEPTH / SALINITY /
#24: (TB(K),K=1,KM) ***VERTICAL TEMPERATURE PROFILE (IF
ITEMP.NE.O)***
#24A: NITT **NUMBER OF STATIONS WITH INITIAL TEMPERATURE DATA
#24B:(ISST(N),JSST(N),NDEPTT(N),TDEPTT(N),N=1,NITT)*I, J,NO.OF
PTS,TOTAL DEPTH
#24C: (CONT'D) / DEPTH / TEMPERATURE /
#25: (CB(K),K=1,KM) ***VERTICAL CONCENTRATION PROFILE (IF
ICC.NE.O)***
#25A: NISSS **NUMBER OF STATIONS WITH INITIAL CONC. DATA
#25B: (ISSS(N) , JSSS(N) ,NDEPTHS(N) ,TDEPTHS(N) ,N=1,NISSS)*I,J,NO.OF
PTS,TOTAL DEPTH
#25C: (CONT'D) / DEPTH / Dissolved CONC. /
#26: NCG NCONST XYR XMONTH XDAY XHR XMIN*TIDAL
PARAMETERS(ITIDE.NE.l SKIP 27THRU30)
01 82 7 20 14 5
-------
#26A: (NCST(I),I=1,NCONST) ***INDEX NUMBER OF TIDAL
CONSTITUENTS***
#26B: KNGAGE(J) HO(J) XLONG(J)(*) ***IF NCG>0,READ 28,29,30
FOR J=1,NCG***
#27: (AMP(I,KNGAGE(J)),I=1,NCONST) ***TIDAL AMPLITUDES***
#28: (XKAPPA(I/KNGAGE(J))/I=1/NCONST) ***TIDAL PHASES***
#29 : (TP(I),I=1,NCONST) ***TABULAR TIDE DATA***
#30: J NC AMPW(J,NC) PHW(J,NC) CAW(J,NC) AMPE(J,NC) PHE(J,NC)
CAE(J,NC)
#31: I NC AMPS(I,NC) PHS(I,NC) CAS(I,NC) AMPN(I,NC) PHN(I,NC)
CAN(I,NC)
#32: NBAR ***NUMBER OF THIN-WALL BARRIERS***
#32A: IJBDIR(I),IJBROW(I),IJBSTR(I),IJBEND(I) ***IF NBAR.GT.0, FOR
EACH NBAR**
#33: IGRID XMAP ALREF ALYREF ***LATERAL GRID MAPPING***
0 1.
#34: NRG ALPHA1
#34A: LPR A B
IN X DIR***
#35: NRG ALPHA1
#36: LPR A B
IN Y DIR***
#37: IF(IBTM.EQ.2) READ ((HS(J,I),1=2,IM),J=2,JM)
* * *BATHYMETRY* * *
300.
300.
450.
450.
700.
700.
850.
850.
1000.
1000.
1000.
1000.
850.
850.
700.
700.
450.
450.
300.
300.
300.
450.
700.
850.
1000.
1000.
850.
700.
450.
300.
300.
450.
700.
850.
1000.
1000.
850.
700.
450.
300.
300.
450.
700.
850.
1000.
1000.
850.
700.
450.
300.
5000000. 5000000.
***VARIABLE GRID MAPPING IN X DIR***
C ***FOR EACH NRG, READ VARIABLE GRID MAPPING
***VARIABLE GRID MAPPING IN Y DIR***
C ***FOR EACH NRG, READ VARIABLE GRID MAPPING
(12F6.1)
300. 300. 300. 300. 300.
450. 450. 450. 450. 450.
700. 700. 700. 700. 700.
850. 850. 850. 850. 850.
1000. 1000. 1000. 1000. 1000.
1000. 1000. 1000. 1000. 1000.
850. 850. 850. 850. 850.
700. 700. 700. 700. 700.
450. 450. 450. 450. 450.
300. 300. 300. 300. 300.
-------
EXAMPLE OUTPUT DATA FILE
FOR AN ENCLOSED BASIN
-------
***THREE-DIMENSIONAL MODEL OF: 3-D RUN FOR SAMPLE RUN WITH
WIND=-1 DYNE/CM**2 (1 DAY) RUN: 135 DATE: 12-APR-90
*THERMALLY HOMOGENEOUS *NO SALINITY *NO SEDIMENT *NO
RIVER *NO TIDE *NO WIND * OPEN BDRY*(IM,JM,KM)= 11 11 5
***** PHYSICAL CONSTANTS AND REFERENCE LENGTHS(IN CGS
UNITS):
XREF ZREF UREF COR GR
ROO ROR TO TR
l.OOOOE+06 5.0000E+02 l.OOOOE+01 9.0000E-05 9.8100E+02
l.OOOOE+00 1.0010E+00 O.OOOOE+00 2.0000E+01
***** FLAGS GOVERNING HYDRODYNAMIC EQUATIONS :
IVLCY ITEMP ISALT IFI IFA
IFB IFC IFD
10010
001
ICC(l) ICC(2) ICC(3) ICC(4)
0000
***** TEMPERATURE PARAMETERS :
BVR SI S2 PR PRV
TWE TWH FKB TQO
l.OOOOE+00 O.OOOOE+00 O.OOOOE+00 l.OOOOE+00 l.OOOOE+00
O.OOOOE+00 O.OOOOE+00 5.0000E+00 O.OOOOE+00
***** CONCENTRATION PARAMETERS :
IVER ICON IUBO IBL IBR
JBM JBP CREF CMAX CO
2 4 0 1 26
1 29 l.OOOOE+00 l.OOOOE+04 l.OOOOE+00
IC1 IC2 JC1 JC2 ID1
ID2 JD1 JD2
00000
000
***** TURBULENCE PARAMETERS :
IEXP IAV AVR AVI AV2
AVM AHR
0 0 5.0000E+00 O.OOOOE+00 O.OOOOE+00
l.OOOOE+00 l.OOOOE+04
FM1 FM2 ZTOP SLMIN QQMIN
-5.0000E-01 -1.5000E+00 -5.0000E-02 2.0000E-03 l.OOOOE-03
QCUT ICUT GAMAX GBMAX FZS
KSMALL
1.5000E-01 1 3.0000E+02 3.0000E+02 l.OOOOE-01
0
***** WIND PARAMETERS :
-------
IWIND TAUX TAUY
0 -l.OOOOE+00 O.OOOOE+00
***** VERTICAL BOUNDARY CONDITION PARAMETERS l
ISMALL ISF ISIE IBTM ITB
10025
HADD HMIN ZREFBN BZ1 HI
H2 ZREFTN TZ1 SSSO
O.OOOOE+00 2.0000E-03 5.0000E+00 4.0000E-01 2.5000E-01
7.5000E-01 5.0000E+00 4.0000E-01 O.OOOOE+00
***** LATERAL BOUNDARY CONDITION PARAMETERS :
IOPEN ITIDE JWIND IJLINE
0000
***** LATERAL BOUNDARY INFO :
J IJGAGE IJDIR IJROW IJSTRT
IJEND
***** ISPAC(I),1=1,10 :
000-10
00000
***** JSPAC(I),1=1,10 :
01-100
00100
***** KSPAC(I),I=1,10 ;
00000
00000
***** RSPAC(I),1=1,10 :
2.0000E-02 l.OOOOE-01 l.OOOOE-05 -l.OOOOE-04 O.OOOOE+00
O.OOOOE+00 l.OOOOE+02 l.OOOOE+00 2.5000E-01 4.0000E+00
***** DERIVED D-LESS PARAMETERS:
RB EV EH FR FRD
DX DY DZ DT DTI
1.1111E-01 2.2222E-01 1.1111E-04 1.4278E-02 4.5151E-01
5.0000E-01 5.0000E-01 2.0000E-01 5.4000E-02 5.4000E-02
***** DERIVED REFERENCE QUANTITIES:
WREF SREF TAUR
5.0000E-03 9.1743E-01 4.5000E-01
***** TIMESTEP INFORMATION :
ISTEP IT1 IT2 ITS DELT
DELTMIN DELTMAX
0 1 144 1 6.0000E+02
l.OOOOE+00 9.0000E+02
EPSILON BUFAC WTS WTU WTV
7.5000E-02 l.OOOOE+01 l.OOOOE+00 l.OOOOE+00 l.OOOOE+00
-------
***** PRINTOUT INFORMATION :
m'c* c*m
TEST
3
JD
1
IGW
1
IGTB
-1
IP1
72
I PA
1
IGI
1
IGC
1
IP2
72
IPB
KPA KPB
11
1 5
IGH
1
IGQ
-1
IP3
72
ID
KD
1
4
IGT
1
IGL
-1
IPU
1
JPA
1
IGS
1
IGR
-1
IPW
1
JPB
11
IGU
1
IGRI
-1
***** DISCFILE INFORMATION :
IRD IW IWR ICI IWC
ICO ISED IWS IREAD IR4
99100
0 14 1 0 15
***** MAJOR DISCFILE NAME :
UVINPUT UVOUTPUT SEDIMENT
6_6_sal 6_6_sal 6_6_sal
***** TIMEBREAKS FOR MAJOR OUTPUT TO DISC :
TBRK1 TBRK2 TBRK3 TBRK4 TBRK5
TBRK6 TBRK7 TBRK8 TBRK9 TBRK10
2.4000E+02 -4.8000E+02 -4.8000E+02 -4.8000E+02 -4.8000E+02
-4.8000E+02 -4.8000E+02 -4.8000E+02 -4.8000E+02 -4.8000E+02
***** TIMEFILE GAGE STATIONS : (NSTA,NRANGE,NFREQ) «= 5
0 2
STATIONID
STATION
1
2
3
4
5
1ST
2
4
6
8
10
JST
6
6
6
6
6
KST
1
1
1
1
1
***** HORIZONTAL DISTANCES AND SLOPES
XMAP = 1. ALREF =
5000000.
***** x-DIRECTION
CELL
5000000. ALYREF
-------
NUMBER
XU FMS
1
0.0000000
2
1.00000
0.5000000
3
1.00000
1,0000000
4
1.00000
1.5000000
5
1.00000
2.0000000
6
1.00000
2.5000000
7
1.00000
3.0000000
8
1.00000
3.5000000
9
1.00000
4.0000000
10
1.00000
4.5000000
11
1.00000
FMU
FACE
CENTER
FACE
1
CENTER
FACE
1
CENTER
FACE
1
CENTER
FACE
1
CENTER
FACE
1
CENTER
FACE
1
CENTER
FACE
1
CENTER
FACE
1
CENTER
FACE
1
CENTER
FACE
1
CENTER
FACE
NOT USED
NOT USED
1
.00000
.00000
.00000
.00000
.00000
.00000
.00000
.00000
.00000
.00000
xs
IN COMPUTATION
IN COMPUTATION
5.0000000 1.00000
***** y-DIRECTION:
CELL
NUMBER
YV FMSV FMV
FACE NOT USED
1 CENTER NOT USED
FACE 1
0.0000000 1.00000
0.2500000
0.7500000
1.2500000
1.7500000
2.2500000
2.7500000
3.2500000
3.7500000
4.2500000
4.7500000
YS
IN COMPUTATION
IN COMPUTATION
-------
0
1
1
2
2
3
3
4
4
2
1.00000
.5000000
3
1.00000
.0000000
4
1.00000
.5000000
5
1.00000
.0000000
6
1.00000
.5000000
7
1.00000
.0000000
8
1.00000
.5000000
9
1.00000
.0000000
10
1.00000
.5000000
11
1.00000
CENTER
FACE
CENTER
FACE
CENTER
FACE
CENTER
FACE
CENTER
FACE
CENTER
FACE
CENTER
FACE
CENTER
FACE
CENTER
FACE
CENTER
FACE
5.0000000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
***** z-DIRECTION:
K Z
5 O.OOOOE+00
4 -2.0000E-01
3 -4.0000E-01
2 -6.0000E-01
1 -B.OOOOE-01
DZZ
2.0000E-01
2.0000E-01
2.0000E-01
2.0000E-01
2.0000E-01
0.2500000
0.7500000
1.2500000
1.7500000
2.2500000
2.7500000
3.2500000
3.7500000
4.2500000
4.7500000
SG
-l.OOOOE-01
-3.0000E-01
-5.0000E-01
-7.0000E-01
-9.0000E-01
***** DEPTHS AT CELL CENTERS (HS)
MAXIMUM MODULUS: 2.0000E+00 PLOT INCREMENT
1.3333E-01
************************
-------
*44444444444*
*66666666666*
*AAAAAAAAAAA*
*CCCCCCCCCCC*
*EEEEEEEEEEE*
*EEEEEEEEEEE*
*CCCCCCCCCCC*
*AAAAAAAAAAA*
*66666666666*
*44444444444*
*44 -4 4444444 4*
************************
***** DEPTHS AT D CELL FACES (HU)
MAXIMUM MODULUS: 2.0000E+00 PLOT INCREMENT
1.3333E-01
************************
*44444444444*
*66666666666*
* AAAAAAAAAAA*
*CCCCCCCCCCC*
*EEEEEEEEEEE*
*EEEEEEEEEEE*
*CCCCCCCCCCC*
*AAAAAAAAAAA*
*66666666666*
*44444444444*
*44444444444*
************************
***** DEPTHS AT V CELL FACES (HV)
MAXIMUM MODULUS: 2.0000E+00 PLOT INCREMENT
1.3333E-01
*44444444444*
*55555555555*
*88888888888*
*BBBBBBBBBBB*
*DDDDDDDDDDD*
*EEEEEEEEEEE*
*DDDDDDDDDDD*
*BBBBBBBBBBB*
*88888888888*
*55555555555*
*44444444444*
************************
***** SLOPES IN THE X-DIRECTION (AT U PTS)
MAXIMUM MODULUS: O.OOOOE+00 PLOT INCREMENT
O.OOOOE+00
-------
EMPTY FIELD - NO PLOT GENERATED
***** CURVATURES IN THE X-DIRECTION (U)
MAXIMUM MODULUS: O.OOOOE+00 PLOT INCREMENT
O.OOOOE+00
EMPTY FIELD - NO PLOT GENERATED
***** SLOPES IN THE Y-DIRECTION (AT V PTS)
MAXIMUM MODULUS: l.OOOOE+00 PLOT INCREMENT :
6.6667E-02
************************
* *
* -8-8-8-8-8-8-8-8-8-8*
* EEEEEEEEEE*
* -9_9_9_9_9_9_9_9_9-9*
* _8_8-8-8-8-8-8-8-8-8*
* *
* 8888888888*
* 9999999999*
* EEEEEEEEE E*
* 8888888888*
* *
************************
***** CURVATURES IN THE Y-DIRECTION (V)
MAXIMUM MODULUS: 1.2000E+00 PLOT INCREMENT :
8.0000E-02
* CCCCCCCCCC*
* *
* -5-5-5-5-5-5-5-5-5-5*
* 777*7*7'7*7*7*7 *7 *
~~,i~~,i~~,j~~,j~~jm~»j~~*j~~j~~jm~,j «
* -E-E-E-E-E-E-E-E-E-E*
* _7_7_7_7_7_7_7_7_7_7*
* _5_5_5_5_5-5_5-5_5_5*
* , *
* CCCCCCCCCC*
* *
************************
***** MOD DEPTHS AT CELL CENTERS (HS)
MAXIMUM MODULUS: 2.0000E+00 PLOT INCREMENT
1.3333E-01
************************
*44444444444*
*66666666666*
-------
*AAAAAAAAAAA*
*CCCCCCCCCCC*
*EEEEEEEEEEE*
*EEEEEEEEEEE*
*CCCCCCCCCCC*
*AAAAAAAAAAA*
*66666666666*
*44444444444*
*44444444444*
************************
***** MOD DEPTHS AT U CELL FACES (HU)
MAXIMUM MODULUS! 2.0000E+00 PLOT INCREMENT
1.3333E-01
************************
*44444444444*
*66666666666*
*AAAAAAAAAAA*
*CCCCCCCCCCC*
*EEEEEEEEEEE*
*EEEEEEEEEEE*
*CCCCCCCCCCC*
*AAAAAAAAAAA*
*66666666666*
*44444444444*
*44444444444*
************************
***** MOD DEPTHS AT V CELL FACES (HV)
MAXIMUM MODULUS: 2.0000E+00 PLOT INCREMENT
1.3333E-01
************************
*44444444444*
*55555555555*
*88888888888*
*BBBBBBBBBBB*
*DDDDDDDDDDD*
*EEEEEEEEEEE*
*DDDDDDDDDDD*
*BBBBBBBBBBB*
*88888888888*
*55555555555*
*44444444444*
************************
***** NS
12345678901234567890123456789012345678901234567890123456789012345
678901234567890
11 02111111115
10 02111111115
-------
9 02111111115
8 02111111115
7 02111111115
6 02111111115
5 02111111115
4 02111111115
3 02111111115
2 02111111115
1 00000000000
***** MS ARRAY
12345678901234567890123456789012345678901234567890123456789012345
678901234567890
11 05555555555
10 01111111111
9 01111111111
8 01111111111
7 01111111111
6 01111111111
5 01111111111
4 01111111111
3 01111111111
2 02222222222
1 00000000000
JU1 JU2 JV1 JV2 I
1
1
1
1
1
1
1
1
1
1
1
11
11
11
11
11
11
11
11
11
11
11
2
2
2
2
2
2
2
2
2
2
2
10
10
10
10
10
10
10
10
10
10
10
1
2
3
4
5
6
7
8
9
10
11
IU1 IU2 IV! IV2
2
2
2
2
2
2
2
2
2
2
2
10
10
10
10
10
10
10
10
10
10
10
1
1
1
1
1
1
1
1
1
1
1
11
11
11
11
11
11
11
11
11
11
11
11
10
9
8
7
6
5
4
3
2
1
-------
TX
MAXIMUM MODULUS: 2.2222E+00 PLOT INCREMENT
1.4815E-01
************************
*_F-F-F-F-F-F-F-F-F-F-F*
*_F-F-F-F-F-F-F-F-F-F-F*
*_F-F-F-F-F-F-F-F-F-F-F*
*_FFFFFFFFFFF*
*_F-F-F-F-F-F-F-F-F-F-F*
*_FFFFFFFFFFF*
*_F-F-F-F-F-F-F-F-F-F-F*
*_F-F-F-F-F-F-F-F-F-F-F*
*_F-F-F-F-F-F-F-F-F-F-F*
*_FFFFFFFFFFF*
*xxxxxxxxxxxxxxxxxxxxxx*
************************
TY
MAXIMUM MODULUS; 1.2340E-07 PLOT INCREMENT
8.2267E-09
************************
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
************************
HU AT (1,2):0.6000E+00
BZO
MAXIMUM MODULUS: 4.0000E-01 PLOT INCREMENT
2.6667E-02
************************
*XX FFFFFFFFFF*
*XX FFFFFFFFFF*
*XX FFFFFFFFFF*
*XX FFFFFFFFFF*
*XX FFFFFFFFFF*
*XX FFFFFFFFFF*
*XX FFFFFFFFFF*
*XX FFFFFFFFFF*
*XX FFFFFFFFFF*
-------
*XX FFFFFFFFFF*
*XX FFFFFFFFFF*
B20U
MAXIMUM MODULUS: 4.0000E-01 PLOT INCREMENT
2.6667E-02
*FFFFFFFFFFF*
*FFFFFFFFFFF*
*FFFFFFFFFFF*
*FFFFFFFFFFF*
*FFFFFFFFFFF*
*FFFFFFFFFFF*
*FFFFFFFFFFF*
*FFFFFFFFFFF*
*FFFFFFFFFFF*
*FFFFFFFFFFF*
*xxxxxxxxxxxxxxxxxxxxxx *
************************
** TIME INDEX = 72 ** TIME - 12.0000 HOURS
SSUM = 2.58684158E-05 : SASUM = 0.OOOOOOOOE+00
SURFACE ELEVATION
MAXIMUM MODULUS; 5.8933E+00 PLOT INCREMENT ;
3.9289E-01
************************
*XX F A 6 3 1-1-4-6-9-D*
*XX D 9 6 3 1-1-4-6-9-C*
*XX C 9 6 3 1-1-4-6-9-B*
*XX C 8 6 3 1-1-4-6-9-A*
*XX B 8 6 3 1-1-4-6-8-A*
*XX B 8 6 3 1-1-4-6-8-A*
*XX C 8 6 3 1-1-4-6-9-A*
*XX C 9 6 3 1-1-4-6-9-B*
*XX D 9 6 3 1-1-4-6-9-C*
*XX E A 6 3 1-1-4-6-9-D*
*xxxxxxxxxxxxxxxxxxxxxx *
X MASS FLUX (UI)
MAXIMUM MODULUS: 2.2921E+00 PLOT INCREMENT
1.5280E-01
************************
* _5_7_8-8-8-8-8-7-5 *
* -5-8-9-A-A-9-8-7-4 *
-------
* -.-.-1-1-1-1-1-.-. *
* 355555542 *
* 8CEEEEDB7 *
* 8CEEEEDB7 *
* 355555542 *
* -.-.-1-1-1-1-1-.-. *
* -5-8-9-A-A-9-8-7-4 *
* -5-7-8-8-8-8-8-7-5 *
*xxxxxxxxxxxxxxxxxxxxxx*
************************
Y MASS FLUX (VI)
MAXIMUM MODULUS: 1.8115E+00 PLOT INCREMENT
1.2077E-01
************************
*xx *
*XX-7-3-.-.-. . .126*
*XX-E-6-2-.-. . 1 3 6 C*
*XX-F-6-3-l-. . 1 3 6 D*
*XX-A-5-2-.-. . 1249*
*xx . .-.-.-.-.-.-. .-.*
*XX A 5 2 . .-.-1-2-4-9*
*XX E 6 3 1 .-.-1-3-6-D*
*XX E 6 2 . .-.-1-3-6-C*
*XX 73.. .-.-.-1-2-6*
*XX *
************************
U-VELOCITY IS
K = 1
MAXIMUM MODULUS: 1.1710E+00 PLOT INCREMENT
7.8066E-02
************************
* -3-7-8-9-9-9-8-7-4 *
* 3 .-1-2-3-2-1 1 3 *
* BAAA99AA9 *
* DCCCCCCCA *
* FEEEEEEEB *
* EEEEEEEEB *
* DCCCCCCCA *
* BAAA99AA9 *
* 3 .-1-2-3-2-1 1 3 *
* -3-7-8-9-9-9-8-7-4 *
*XXXXXXXXXXXXXXXXXXXXXX*
K = 5
MAXIMUM MODULUS: 4.1082E+00 PLOT INCREMENT
2.7388E-01
************************
-------
* -A-D-E-E-E-E-E-D-A
* -B-D-D-E-E-D-D-B-9
* _9_9_9_g_9_9_9_9_8
* _g_7_7_7_7_7_7_8-8
* -9-6-5-5-5-5-6-7-7
* _9_6-5-5-5-5-6-7-7
* _g_7_7_7_7_7_7_8-8
* -B-D-D-E-E-D-D-B-9
* -A-D-E-E-F-E-E-D-A
* xxxxxxxxxxxxxxxxxxxxxx *
************************
V-VELOCITY
K = 1
MAXIMUM MODULUS:
4.6765E-02
7.0147E-01
PLOT INCREMENT
K = 5
MAXIMUM MODULUS;
1.3188E-01
************************
*xx
*XX-D-7-2-.~. .
*XX-E-7-3-l-. .
*XX-B-4-2-.-. .
*XX-6-2-l-.-. .
*XX-. .-.-.-.-.
*
. 2 6 C*
1 3 5 D*
. 1 3 A*
.126*
_ _ *
*XX 621.
*XX B 4 2 .
*XX E 7 3 1
*XX D 7 2 .
*XX
,-.-.-1-2-6*
.-.-.-1-3-A*
.-.-1-3-5-D*
.-.-.-2-6-C*
*
************************
1.9781E+00
PLOT INCREMENT
************************
*xx
*XX-A-4-l-
*XX-E-5-2-
*XX-B-4-2-
*XX-6-2-l-
*xx - -
*XX 621
*XX B 4 2
*XX F 5 2
*XX A 4 1
*XX
*
.-. . . 1 4 B*
. . . 1 3 7 E*
. . . 1 3 6 A*
. . .1246*
*
._._. -1-2-4-6*
.-.-.-1-3-6-A*
.-.-.-1-3-7-E*
. .-.-.-1-4-B*
*
W-VELOCITY IS
K = 1
MAXIMUM MODULUS:
3.6162E-01
PLOT INCREMENT
-------
2.4108E-02
************************
*XX-A
*XX-E
*XX-E
*XX-E
*XX-E
*XX-F
*XX-E
*XX-A
1 1 .
1 . .
2 . .
3 . .
3 . .
2 . .
1 . .
1 1 .
.-.-.-1-1
.-.-.-. 1
.-.-.-. 1
.-.-.-. 1
._._._. i
.-.-.-1-1
. . . . 1
9*
r*
r*
A*
A*
r*
r*
9*
6*
*xxxxxxxxxxxxxxxxxxxxxx *
************************
K = 5
MAXIMUM MODULUS:
8.2267E-09
1.2340E-07
PLOT INCREMENT
************************
*XX *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xxxxxxxxxxxxxxxxxxxxxx *
************************
** TIME INDEX - 144 ** TIME = 24.0000 HOURS
SSUM = 7.93933868E-05 i SASUM = O.OOOOOOOOE+00
SURFACE ELEVATION
MAXIMUM MODULUS:
3.9861E-01
5.9791E+00
PLOT INCREMENT
************************
*XX E A 6 3
*XX D 9 6 3
*XX C 8 6 3
*XX B 8 6 3
*XX B 8 6 3
*XX B 8 6 3
*XX B 8 6 3
*XX C 8 6 3
*XX D 9 6 3
*XX F A 6 3
-1-4-6-9-C*
-1-4-6-9-C*
1-4-6-9-B*
-l_4_6-8-A*
-1-4-6-8-A*
1-4-6-8-A*
.1-4-6-8-A*
.1-4-6-9-B*
1-4-6-9-C*
1-4-6-9-C*
-------
*xxxxxxxxxxxxxxxxxxxxxx *
************************
X MASS FLUX (UI)
MAXIMUM MODULUS! 3.0327E+00 PLOT INCREMENT :
2.0218E-01
************************
* -4_5_6_6-6-6-6-6-4 *
* -6-8-9-9-9-9-8-6-4 *
* -1-2-2-2-2-2-2-1-1 *
* 344444331 *
* 9DEEEEDB7 *
* 9DEFEEDB7 *
* 344444331 *
* -1-2-2-2-2-2-2-1-1 *
* -g-8-9-9-9-9-8-6-4 *
* _4_s_6-6-6-6-6-6-4 *
*xxxxxxxxxxxxxxxxxxxxxx*
************************
Y MASS FLUX (VI)
MAXIMUM MODULUS: 2.5646E+00 PLOT INCREMENT :
1.7098E-01
************************
*xx *
*XX-4-2-.-.-. ... 2 5*
*XX-C-4-l-. . .1259*
*XX-F-5-2-. . . 1 3 6 B*
*XX-B-4-l-. . .1248*
*xx-. .-.-. .-. .-. . .*
*XX B 4 1 .-.-.-1-2-4-8*
*XX E 5 2 .-.-.-1-3-6-B*
*XX C 4 1 .-.-.-1-2-5-9*
*XX 42.. .-.-.-.-2-5*
*XX *
************************
U-VELOCITY IS
K*~. -I
~ JL
MAXIMUM MODULUS: 1.2088E+00 PLOT INCREMENT :
8.0586E-02
************************
* -2-7-8-8-9-9-8-8-4 *
* 1-3-5-5-5-5-4-1 2 *
* A99999998 *
* CCCCCCCBA *
* EEEEEEEEB *
* EEEEEEEEB *
-------
K = 5
MAXIMUM MODULUS:
3.0329E-01
* CCCCCCCBA *
* A99999998 *
* 1-3-5-5-5-5-4-1 2 *
* -2-7-8-8-9-9-8-8-4 *
*XXXXXXXXXXXXXXXXXXXXXX*
************************
4.5494E+00
PLOT INCREMENT
************************
* -9-C-D-D-D-D-D-C-9 *
* -C-E-E-E-E-E-D-B-9 *
* _g_9_9_A-9-9-9-9-8 *
* _B-fi-6-6-6-6-7-7-7 *
* _g_3_3_3_3_3_4-5-6 *
* -6-3-3-3-3-3-4-5-6 *
* -8-6-6-6-6-6-7-7-7 *
* -9-9-9-A-9-9-9-9-8 *
* -C-E-E-E-E-E-D-B-9 *
* -9-C-D-D-D-D-D-C-9 *
*xxxxxxxxxxxxxxxxxxxxxx *
************************
V-VELOCITY
K - 1
MAXIMUM MODULUS:
5.4502E-02
K = 5
MAXIMUM MODULUS:
1.6339E-01
8.1753E-01
PLOT INCREMENT
************************
*xx
*XX-B-6-2-.-.
*XX-E-7-3-l-.
*XX-C-4-l-.-.
*XX-7-2-.-.-.
*
. . 1 5 B*
. 1 3 6 B*
. .139*
. ,127*
*XX- -. .-.-. .*
*XX 7 2
*XX C 4 1 .
*XX F 7 3 1
*XX B 6 2 .
*XX
._._, -1-2-7*
..... -1-3-9*
*
************************
2.4508E+00
PLOT INCREMENT t
************************
*XX
*XX-8-3-l-.-. .
*XX-E-4-l-. . .
*XX-C-3-l-. . .
*XX-7-2-.-. . .
*xx . .-.-. .-.
*
.149*
1 3 7 C*
1 3 6 A*
1247*
_ _ *
-------
*XX 7 2 .
*XX C 3 1
*XX E 4 1
*XX 831
*XX
,-.-.-1-2-4-7*
.-.-.-1-3-6-A*
.-.-.-1-3-7-C*
, .-.-.-1-4-9*
*
************************
W-VELOCITY IS
K = 1
MAXIMUM MODULUS:
2.3619E-02
3.5428E-01
PLOT INCREMENT
************************
K = 5
MAXIMUM MODULUS:
8.2267E-09
*XX-9-l-,
*XX-A 1- ,
*XX-E 2 ,
*XX-E 2 ,
*XX-D 3 ,
*XX-D 3 ,
*XX-F 2 ,
*XX-E 2 ,
*XX-A 1-,
*XX-9-l-,
1
.-. . . .-.-2
, . . .-.-. 1
. . .-.-.-. 1
. . .-.-.-. 1
, . . .-.-. 1
.-. . . .-.-2
,-.-. . . . 1
6*
8*
C*
r*
A*
A*
C*
r*
8*
6*
*xxxxxxxxxxxxxxxxxxxxxx*
************************
1.2340E-07
PLOT INCREMENT
*XX *
*XX *
*XX *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xx *
*xxxxxxxxxxxxxxxxxxxxxx*
************************
***** MAXIMUM STEP REACHED : IT
144
-------
NEEL-ATH-1Z85
TECHNICAL REPORT DATA
(Plate rectf Inttntcliont on ih? reverse be fate completing)
\, REPORT NO,
600/R-99/Q49
3, BtClH fcW 1 'S ACCtSSi ON NO.
4, TITLE AND SUBTITLE
Three Dimensional Hydro dynamic Model for Stratified
Flows itt Lakes and Estuaries (HYDR03D)
OATt
June
6-PERFORMING ORGANIZATION CODE
7.AUTHOWS) y. Peter sheng, Mansour
Steven C. McCutcheon, E.Z- Hosseinipour, Pei-Fang
B. CEHFQHMING ORGANIZATION Ri-FQ-HT NO
Donalj. Eliasqn,. D.^j
P. PERFORMING ORGANi£ATic»N
University of FLorida
Gainesville, FL
ANS
Parker.E.
ihL.
i u. PROGRAM
no.
11. CONTKACT^GrtANT (VO.
12. SFONSOfllNO AGENCY NAME AND AODH6SS
Ecoaystems Hesearch Division - Athens, GA
Office of Research and Development
U.S, Environmental Protection Agency
Athens, GA 30605-2700
13. t YP6 OF HfcPORT ANOPERlOO tOVE RE D
14, SPONSORING AGEFVCY CODE
iiFA/600/01
15. SUPPLEM£h|TAflY NOTE?
16. ABSTHACT
Increasing dtmatida tar maintaining the quality of stratified surface waters
at reasonable levels ha^e required the development of threo-dicietisior^al hydTodymsfiiit
models. To meet these needs, the HYUR03D prograji lias been documented to aid in the
simulation of lakes, harbors, coastal areas, stt<1 eptuar.ies.
HnCDROSl} is a dynamic modeling system that can be used to simulate currents in
water bodies as they respond to tides, winos, density gradients, river flows, and
basin geometry arid bathymetry. The code is a three-dimensional, time-dependent,
o-stretehed coordinate > free surface tcodel that can be run in irully three-d ioiensiuiia
(3-B) mode, two-dimensional vertically^averaged (x-y), snd tvo-dlT^ensional laterally
averaged X-Z made.
Thf, applications provided here demonstrated that tfre model is depable of
realistic simulation of flow and salinity transport in complex and dynamic water
bodies. These applications include simulations of tidal circulation and salinity
transport in Sulsiin Bay> California and Charlotte Harbor, Florida and wind-Toreed
circulation in Green Bay, Lake Hichigaii. Tidal circulation in Prince William Sound,
Alaska was investigated to determine the feasibility of applying the model under
emergency conditions. Finally, the calibration ot Che model for the Mississippi
Sound is illustrated.
KFV WORM AN& DOCUMENT ANALYSIS
Water Quality
HYDR03TD
Mode1ing
Estuaries
Lakes
S/OPEN ENC1ED TERM'S
C. CO SAT I
STATEMENT
RELEASE TO PUBLIC
18, SECURITY CLASS {ThisReport,
^CLASSIFIED
ZO. SECUHIT V CLASS /
UNCLASSIFIED
31. NO. OF CAGES
223
SB. PRICE
Fwm 3120-
T. 4-77) PREVIOUS
------- |