|
UFL/COEL-TR/103
WIND-DRIVEN CIRCULATION IN LAKE
OKEECHOBEE, FLORIDA: THE EFFECTS OF
THERMAL STRATIFICATION AND AQUATIC
VEGETATION
by
Hye Keun Lee
Dissertation
1993
WIND-DRIVEN CIRCULATION IN LAKE OKEECHOBEE. FLORIDA: THE
EFFECTS OF THERMAL STRATIFICATION AND AQUATIC VEGETATION
By
HYE KEUN LEE
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1993
ACKNOWLEDGEMENTS
I would like to express my sincere gratitude to my advisor Dr. Y. Peter Sheng for
his continuous guidance, encouragement and financial support throughout my study.
I would also like to extend my thanks and appreciation to my doctoral committee
members, Dr. Robert G. Dean, Dr. Donald M. Sheppard and Dr. Ulrich H. Kurzweg,
for their patience in reviewing this dissertation. My gratitude also extends to Dr.
Robert J. Thieke who reviewed my dissertation.
I must thank Dr. Paul W. Chun for reviewing my dissertation and the
great guidance during my stay in Gainesville while he served as a faculty advisor
of the Korean Student Association.
Financial support provided by the South Florida Water Management District,
West Palm Beach, Florida, through the Lake Okeechobee Phosphorus Dynamics
Project is appreciated.
I would like to dedicate this dissertation to my late father and my mother.
Finally, I would like to thank my loving wife, Aesook, for her support and patience,
and my beautiful daughter, Mireong, and my smart son, David.
TABLE OF CONTENTS
ACKNOWLEDGEMENTS ............................ ii
LIST OF FIGURES ...............
LIST OF TABLES ................
. . .. vi
. . . xiii
ABSTRACT .. .. .. .. ... .. .... .. . .. .. .
CHAPTERS
1 INTRODUCTION ...............................
2 LITERATURE REVIEW ...........................
2.1 Numerical Models of Lake Circulation ..................
2.1.1 One-Dimensional Model .....................
2.1.2 Two-Dimensional Model .....................
2.1.3 Steady-State 3-D Models ......................
2.1.4 Time-Dependent 3-D Models ..................
2.2 Vegetation M odels ............................
2.3 Thermal M odels ..............................
2.4 Turbulence Model ..................... ..........
2.4.1 Eddy Viscosity/Diffusivity Concept ..............
2.4.2 Constant Eddy Viscosity/Diffusivity Model ...........
2.4.3 Munk-Anderson Type Model ...................
2.4.4 Reynolds Stress Model ......................
2.4.5 A Simplified Second-Order Closure Model: Equilibrium Closure
M odel . . . .. . . . .
2.4.6 A Turbulent Kinetic Energy (TKE) Closure Model . .
2.4.7 One-Equation Model (k Model) .................
2.4.8 Two-Equation Model (k e Model) . . . .
2.5 Previous Lake Okeechobee Studies .................. .
2.6 Present Study .. .. .. ... .. .... . .. .. .. .
3 GOVERNING EQUATIONS . . . .
3.1 Introduction .....................
3.2 Dimensional Equations and Boundary Conditions
ordinate System . . . . ..
3.2.1 Equation of Motion . . . .
3.2.2 Free-Surface Boundary Condition (z = r) .
3.2.3 Bottom Boundary Condition (z = -h) .
3.2.4 Lateral Boundary Condition . .
3.3 Vertical Grid ....................
3.4 Non-Dimensionalization of Equations .. .
in a Cartesian Co-
3.5 Dimensionless Equations in cr-Stretched Cartesian Grid . ... 27
3.5.1 Vertically-Integrated Equations . . . ... 28
3.5.2 Vertical Velocities ......................... 29
3.6 Generation of Numerical Grid. . . . . 29
3.6.1 Cartesian Grid .......................... 29
3.6.2 Curvilinear Grid ......................... 29
3.6.3 Numerical Grid Generation . . . . ... 30
3.7 Transformation Rules........................... 31
3.8 Tensor-Invariant Governing Equations . . . ... 34
3.9 Dimensionless Equations in Boundary-Fitted Grids . ... 36
3.10 Boundary Conditions and Initial Conditions . . ... 37
3.10.1 Vertical Boundary Conditions . . . ... 37
3.10.2 Lateral Boundary Conditions . . . ... 37
3.10.3 Initial Conditions ... .................... ..... .. 37
4 VEGETATION MODEL ........................ .... .. 39
4.1 Introduction ... . .. .. .. .. ..... .. .. .. 39
4.2 Governing Equations ........................... 41
4.2.1 Equations for the Vegetation Layer (Layer I) . ... 41
4.2.2 Equations for the Vegetation-Free Layer (Layer II) .. .43
4.2.3 Equations for the Entire Water Column . . ... 44
4.2.4 Dimensionless Equations in Curvilinear Grids . ... 45
5 HEAT FLUX MODEL ............. ... ....... .... .. 48
5.1 Introduction .. . .. .. .. ... .... .. .. 48
5.2 The "Equilibrium Temperature" Method . . . ... 48
5.2.1 Short-Wave Solar Radiation . . . . ... 49
5.2.2 Long-Wave Solar Radiation . . . . ... 49
5.2.3 Reflected Solar and Atmospheric Radiation . . ... 49
5.2.4 Back Radiation .............. ... ... .. .... .. 51
5.2.5 Evaporation ............................ 51
5.2.6 Conduction ............. ... .... ....... 51
5.2.7 Equilibrium Temperature . . .... . 52
5.2.8 Linear Assumption ........................ 52
5.2.9 Procedure for an Estimation of K and T . . ... 53
5.2.10 Modification of the Equilibrium Temperature Method . 54
5.3 The "Inverse" Method ......................... 54
5.3.1 Governing Equations . . . . . ... 55
5.3.2 Boundary Conditions . . . . . ... 55
5.3.3 Finite-Difference Equation . . . . ... 56
5.3.4 Procedure for an Estimation of Total Heat Flux . ... 57
6 FINITE-DIFFERENCE EQUATIONS . . . . ... 58
6.1 Grid System .. . . .......... .. .... 58
6.2 External Mode ............................. 58
6.3 Internal M ode ............... ................ 61
6.4 Temperature Scheme .......................... 62
6.4.1 Advection Terms ........................ 63
6.4.2 Horizontal Diffusion Term . . . . ... 66
7 MODEL ANALYTICAL TEST .............
7.1 Seiche Test .. .. .. .. .. . ..
7.2 Steady State Test ..................
7.3 Effect of Vegetation .................
7.4 Thermal Model Test ................
8 MODEL APPLICATION TO LAKE OKEECHOBEE
8.1 Introduction ....... ..... .........
8.1.1 Geometry ..................
8.1.2 Temperature. . . . .
8.2 Some Recent Hydrodynamic Data . . .
8.2.1 W ind Data .................
8.2.2 Current Data ................
8.2.3 Temperature Data .............
8.2.4 Vegetation Data. ...............
8.3 Model Setup ..... ..... ...... ...
8.3.1 Grid Generation ....... ........
8.3.2 Generation of Bathymetry Array . .
8.4 Model Parameters .. . . . ...
8.4.1 Reference Values ..............
8.4.2 Turbulence Model and Parameters .
8.4.3 Bottom Friction Parameters . .
8.4.4 Vegetation Parameters . . .
8.4.5 W ind Stress .................
8.5 Steady State Wind-Driven Circulation . .
8.6 Wind-Driven Circulation without Thermal Stratification . .
8.6.1 Tests of Model Performance . . . . .
8.6.2 Model Results ............... ..............
8.7 Wind-Driven Circulation with Thermal Stratification: Te Method .
8.8 Simulation of Currents with Thermal Stratification: Inverse Method .
8.8.1 The Diurnal Thermal Cycle . . . . .
8.9 Sensitivity Tests ..... .. .. .. .. .. .. .. .... .
8.9.1 Effect of Bottom Stress . . . . . .
8.9.2 Effect of Horizontal Diffusion Coefficent . . .
8.9.3 Effect of Different Turbulence Model . . . .
8.9.4 Effect of Advection Term . . . ..... ...
8.10 Spectral Analysis .......................... ...
9 CONCLUSION .................................
APPENDIX
A SIMULATED CURRENTS BY INVERSE METHOD ...... .. ..
BIBLIOGRAPHY .................................
BIOGRAPHICAL SKETCH ...........................
...........
...........
...........
...........
...........
. . . .
...........
...........
...........
...........
. . . .
........ ...
........ ...
...........
. . . .
. . . .
76
76
76
76
76
77
81
81
83
83
83
84
88
88
89
91
92
96
98
99
99
99
123
150
150
156
167
168
168
168
169
174
LIST OF FIGURES
3.1 A computational domain and a transformed coordinate system.
4.1 Schematics of flow in vegetation zone . . . .
5.1 Meteorological data at Station L006 . . . .
6.1 Horizontal and vertical grid system . . . .
7.1 Model results of a seiche test . . . . .
7.2 Surface elevation contour when the lake is steady state with uni-
form wind stress of -1 dyne/cm2 . . . . .
7.3 Effect of vegetation on surface elevation evolution in a wind-driven
rectangular lake. Solid line is without vegetation, broken line is
with low vegetation density, and dotted line is with high vegetation
density. . . . . . . . .
7.4 Time history of wind stress and currents at the center
all five levels. Thermal stratification is not considered.
7.5 Time history of wind stress and currents at the center
all five levels. Thermal stratification is considered. ..
8.1 Map of Lake Okeechobee. . . . .
8.2 Wind rose at Station C ..................
8.3 Computation domain of Lake Okeechobee . .
8.4 Curvilinear grid of Lake Okeechobee. . . .
8.5 Depth contour of Lake Okeechobee when the lake stage
Unit in cm ... .. ...... . . .. .
8.6 Distribution of vegetation height in Lake Okeechobee.
8.7 Distribution of vegetation density in Lake Okeechobee.
of lake at
of lake at
. . 75
. . 78
. . 82
. . 85
. . 86
is 15.5 ft.
. . 87
. . 95
. . 97
Steady-state depth-integrated currents (cm2s-1) in Lake Okee-
chobee forced by an easterly wind of 1 dyne/cm2 . . .
8.9 Steady-state surface elevation contour (cm) in Lake Okeechobee
forced by an easterly wind of 1 dyne/cm2. . . ... 101
8.10 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 1: East-West direction). Thermal stratification
was not considered in model simulation . . . ... 104
8.11 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 1: North-South direction). Thermal stratification
was not considered in model simulation . . . ... 105
8.12 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 2: East-West direction). Thermal stratification
was not considered in model simulation . . . ... 106
8.13 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 2: North-South direction). Thermal stratification
was not considered in model simulation . . ..... 107
8.14 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 3: East-West direction). Thermal stratification
was not considered in model simulation . . ..... 108
8.15 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 3: North-South direction) . . .. 109
8.16 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 1: East-West direction). Thermal stratification
was not considered in model simulation. . . ... 110
8.17 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 1: North-South direction). Thermal stratification
was not considered in model simulation. . . .... 111
8.18 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 2: East-West direction). Thermal stratification
was not considered in model simulation. . . .. 112
8.19 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 2: North-South direction). Thermal stratification
was not considered in model simulation. . . ... 113
8.20 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 1: East-West direction). Thermal stratification
was not considered in model simulation. . . ... 116
8.21 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 1: North-South direction). Thermal stratification
was not considered in model simulation . . .. 117
8.22 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 2: East-West direction). Thermal stratification
was not considered in model simulation. . . ... 118
8.23 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 2: North-South direction). Thermal stratification
was not considered in model simulation . . .. 119
8.24 Stick Diagram of wind stress, measured currents, and simualted
currents at Station E ...................... 120
8.25 Simulated (solid lines) and measured (dotted lines) currents at
Station A (Arm 1: East-West direction). Thermal stratification
was not considered in model simulation. . . ... 121
8.26 Simulated (solid lines) and measured (dotted lines) currents at
Station A (Arm 1: North-South direction). Thermal stratification
was not considered in model simulation . . .. 122
8.27 Simulated (solid lines) and measured (d6tted lines) currents at
Station D (Arm 1: East-West direction) . . .... 124
8.28 Simulated (solid lines) and measured (dotted lines) currents at
Station D (Arm 1: North-South direction) . . .. 125
8.29 Simulated (solid lines) and measured (dotted lines) currents at
Station D (Arm 2: East-West direction) . . .... 126
8.30 Simulated (solid lines) and measured (dotted lines) currents at
Station D (Arm 2: North-South direction) . . .. 127
8.31 Simulated (solid lines) and measured (dotted lines) currents at
Station A (Arm 1: East-West direction) when thermal effect is
considered ... .. .. .. .. .. .. .. .. .. 129
8.32 Simulated (solid lines) and measured (dotted lines) currents at
Station A (Arm 1: North-South direction) when thermal effect is
considered . . . . .. . .. 130
8.33 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 1: East-West direction) when thermal effect is
considered ... .. .. .. ... .... ... . .. 131
8.34 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 1: North-South direction) when thermal effect is
considered ... .. .. .. ............. .. ... 132
8.35 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 2: East-West direction) when thermal effect is
considered ... . .. ... .. .. ... . .. 133
8.36 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 2: North-South direction) when thermal effect is
considered .... .. ............... ... .. 134
8.37 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 1: East-West direction) when thermal effect is
considered .. .. . . . .. .. .. 135
8.38 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 1: North-South direction) when thermal effect is
considered. ................... ........... 136
8.39 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 2: East-West direction) when thermal effect is
considered. ....... ........... .. ......... 137
8.40 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 2: North-South direction) when thermal effect is
considered ..... .. .. .. .. .. .. ...... 138
8.41 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 3: East-West direction) when thermal effect is
considered. ... .. . . . .. .. .. .. 139
8.42 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 3: North-South direction) when thermal effect is
considered. ... .................... ....... 140
8.43 Simulated (solid lines) and measured (dotted lines) currents at
Station D (Arm 1: East-West direction) when thermal effect is
considered ............ ..... ............... 141
8.44 Simulated (solid lines) and measured (dotted lines) currents at'
Station D (Arm 1: North-South direction) when thermal effect is
considered. ... .. .. ... . .. ... .... .. 142
8.45 Simulated (solid lines) and measured (dotted lines) currents at
Station D (Arm 2: East-West direction) when thermal effect is
considered. .. .. . . .. .. .. .... .. 143
8.46 Simulated (solid lines) and measured (dotted lines) currents at
Station D (Arm 2: North-South direction) when thermal effect is
considered ............................ 144
8.47 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 1: East-West direction) when thermal effect is
considered. ..... .. .. .. .. .. .. .. ... .. . 145
8.48 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 1: North-South direction) when thermal effect is
considered. ..... .. .. . .. .. .. ..... 146
8.49 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 2: East-West direction) when thermal effect is
considered. ............................. 147
8.50 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 2: North-South direction) when thermal effect is
considered . .. .. . . . . .. 148
8.51 Time history of eddy viscosity at Station C between Julian days
147 and 161. .... .. .... . . .. .. .. .. .. 151
8.52 Time history of wind stress and measured currents between Julian
days 150 and 152 ............................ 152
8.53 Time history of simulated currents at Station C between Julian
days 150 and 152 ............................ 153
8.54 Time history of eddy viscosity at Station C between Julian days
150 and 152 .. .. .. .... .. .. . .. ... .. 154
8.55 Time history of heat fluxes at Station C between Julian days 147
and 161 . . . . . . . . 155
8.56 Temperature contours of data and model at Station C between
Julian days 152 and 155.......................... 157
8.57 Simulated and measured temperatures at Station A. ...... ..158
8.58 Simulated and measured temperatures at Station B. ...... ..159
8.59 Simulated and measured temperatures at Station C. ...... ..160
8.60 Simulated and measured temperatures at Station D. ...... ..161
8.61 Simulated and measured temperatures at Station E. ...... ..162
8.62 Vertical profiles of currents and temperature during a typical day.
Base temperature is 25 C. .................... 163
8.63 Vertical profiles of currents and temperature during a typical day.
Base temperature is 25 C. .................... 164
8.64 Vertical profiles of currents and temperature during a typical day.
Base temperature is 25 C. .................... 165
8.65 Vertical profiles of currents and temperature during a typical day.
Base temperature is 25 C. .................... 166
8.66 Spectrum of wind stress and surface elevation . . ... 171
8.67 Spectrum of measured and simulated currents (east-west direc-
tion) at Station C .......................... 172
8.68 Spectrum of measured and simulated currents (north-south direc-
tion) at Station C .......................... 173
A.1 Simulated (solid lines) and measured (dotted lines) currents at
Station A (Arm 1: East-West direction). Inverse method was
used for the estimation of heat flux . . . .. 177
A.2 Simulated (solid lines) and measured (dotted lines) currents at
Station A (Arm 1: North-South direction). Inverse method was
used for the estimation of heat flux . . . .. 178
A.3 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 1: East-West direction). Inverse method was
used for the estimation of heat flux . . . .. 179
A.4 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 1: North-South direction). Inverse method was
used for the estimation of heat flux . . . .. 180
A.5 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 2: East-West direction). Inverse method was
used for the estimation of heat flux . . . .. 181
A.6 Simulated (solid lines) and measured (dotted lines) currents at
Station B (Arm 2: North-South direction). Inverse method was
used for the estimation of heat flux. . . . . 182
A.7 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 1: East-West direction). Inverse method was
used for the estimation of heat flux. . . . ... 183
A.8 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 1: North-South direction). Inverse method was
used for the estimation of heat flux . . . .. 184
A.9 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 2: East-West direction). Inverse method was
used for the estimation of heat flux. . . . ... 185
A.10 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 2: North-South direction). Inverse method was
used for the estimation of heat flux . . . .. 186
A.11 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 3: East-West direction). Inverse method was
used for the estimation of heat flux . . . .. 187
A.12 Simulated (solid lines) and measured (dotted lines) currents at
Station C (Arm 3: North-South direction). Inverse method was
used for the estimation of heat flux. . . . . 188
A.13 Simulated (solid lines) and measured (dotted lines) currents at
Station D (Arm 1: East-West direction). Inverse method was
used for the estimation of heat flux. . . . ... 189
A.14 Simulated (solid lines) and measured (dotted lines) currents at
Station D (Arm 1: North-South direction). Inverse method was
used for the estimation of heat flux . . . .. 190
A.15 Simulated (solid lines) and measured (dotted lines) currents at
Station D (Arm 2: East-West direction). Inverse method was
used for the estimation of heat flux . . . .. 191
A.16 Simulated (solid lines) and measured (dotted lines) currents at
Station D (Arm 2: North-South direction). Inverse method was
used for the estimation of heat flux . . . .. 192
A.17 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 1: East-West direction). Inverse method was
used for the estimation of heat flux . . . .. 193
A.18 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 1: North-South direction). Inverse method was
used for the estimation of heat flux . . . .. 194
A.19 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 2: East-West direction). Inverse method was
used for the estimation of heat flux . . . ... 195
A.20 Simulated (solid lines) and measured (dotted lines) currents at
Station E (Arm 2: North-South direction). Inverse method was
used for the estimation of heat flux . . . .. 196
LIST OF TABLES
Selected features of lake models. . . . ... 21
Application features of lake models. . . . ... 22
Installation dates and locations of platforms during 1988 and 1989. 79
Instrument mounting, spring deployment. . . ... 80
Reference values used in the Lake Okeechobee spring 1989 circu-
lation simulation ............................ 89
Vertical turbulence parameters used in the Lake Okeechobee spring
1989 circulation simulation. . . . ..... 92
Vegetations in Lake Okeechobee (From Richardson, 1991) . 94
Index of agreement and RMS error at Station C. . ... 103
Index of agreement and RMS error at Station B.
Index of agreement and RMS error at Station E.
Index of agreement and RMS error at Station A.
Index of agreement and RMS error at Station D.
Index of agreement and RMS error at Station A
effect is considered ..................
Index of agreement and RMS error at Station B
effect is considered ..................
. . 114
. . 115
. . 123
. . 123
when thermal
when thermal
.o.......
Index of agreement and RMS error at Station C when thermal
effect is considered......................... ..
Index of agreement and RMS error at Station D when thermal
effect is considered...........................
Index of agreement and RMS error at Station E when thermal
effect is considered...........................
Parameters used in sensitivity tests . . . .
149
149
149
149
149
167
8.7
8.8
8.9
8.10
8.11
8.12
8.13
8.14
8.15
8.16
8.17 Index of agreement and RMS error. . . . ... 170
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
WIND-DRIVEN CIRCULATION IN LAKE OKEECHOBEE, FLORIDA: THE
EFFECTS OF THERMAL STRATIFICATION AND AQUATIC VEGETATION
By
HYE KEUN LEE
August 1993
Chairman: Dr. Y.P. Sheng
Major Department: Coastal and Oceanographic Engineering
Wind-driven circulation in Lake Okeechobee, Florida, is simulated by using a
three-dimensional curvilinear-grid hydrodynamic model and measured field data. Field
data show that significant thermal stratification often develops in the vertical water
column during daytime in the large and shallow lake. Significant wind mixing due to
the lake breeze, however, generally leads to destratification of the water column in the
late afternoon and throughout the night. Thus, thermal effects must be considered
in the numerical simulation of circulation in shallow lakes.
During daytime the lake is thermally stratified and wind is relatively weak, mo-
mentum transfer is generally limited to the upper layer and hence the bottom currents
are much weaker than the surface currents. During the initial phase of significant
lake breeze, strong surface currents and opposing bottom currents are developed, fol-
lowed by oscillatory motions associated with seiche and internal seiche, until they are
damped by bottom friction.
Lake Okeechobee is covered with submerged and emergent aquatic vegetation over
much of the bottom on the western portion of the lake (20 % of the surface area). The
presence of the vegetation causes damping of the wind, wave and current fields. To
provide realistic simulation of wind-driven circulation in the presence of vegetation,
this study developed a simplified vegetation model which parameterizes the effect
of vegetation in terms of added "form drag" terms in the momentum equations.
Simulated currents in the open water region in the vicinity of vegetation compare
quite well with data. This physical process is successfully modeled by parameterizing
the vertical turbulence with a simplified second-order closure model. Model simulation
which assumes homogeneous density structure fails to represent the stratification and
destratification cycle. On the other hand, simulation which includes thermal effects
faithfully reproduced field data.
CHAPTER 1
INTRODUCTION
Lakes are valuable resource for a variety of human needs: drinking water, agri-
cultural use, navigation, waste water disposal, recreation sites, cooling reservoirs for
power plants, etc. Ninety-nine percent of Americans live within 50 miles of one of
37,000 lakes (Hanmer, 1984). Lakes can be dangerous to people under certain cir-
cumstances such as during flooding and the deterioration of water quality due to the
excessive loading of contaminants into the lake.
Water movement in lakes is driven by wind, density gradient, waves, and tributary
flow, but primarily is influenced by wind action. During periods of strong wind, severe
flooding can be caused by the storm surge. Some examples are the flooding in the
Lake Okeechobee area during the 1926 and 1928 hurricanes. The 1926 hurricane
caused a storm surge of 7 ft at Moorehaven on the western side of Lake Okeechobee,
and about 150 people lost their lives (Hellstrom, 1941).
Water quality of lakes is of utmost importance. So long as human activities are
limited to a small part of a lake, it may appear that the lake has an unlimited capacity
of self-purification. However, as population and human development increase, a lake
may not be able to endure the excessive stresses caused by human actions, and water
quality may become deteriorated. Typical evidence of poor water quality includes
sudden algal bloom, colored water, fish kill, taste and odor in drinking water, and
floating debris of plants. Eutrophication is the process in which excessive loading of
nutrients, organic matter, and sediments into lakes results in an increase of primary
production. Sources of eutrophication are increased use of fertilizer, waste water
discharge, and precipitation of polluted air.
2
Earlier studies on physical processes in lakes concentrated on the observations of
periodic up and down motion of water level, i.e., seiche motion. When the wind blows
over a certain period, water builds up near the shoreline. After wind ceases, water
starts to oscillate as a free long wave. Subsequently, this wave is damped out due to
bottom friction. Seiche can be initiated by sudden change of wind speed or direction.
the passage of a squall line, an earthquake, or resonance of air and water columns.
When the velocity of a squall line is close to the speed of gravity wave, resonance
occurs and damages can be more severe. It was reported that severe storm damages
in Chicago were caused by a squall line over Lake Michigan on June 26, 1954 (Harris,
1957).
Another important feature in deep temperate lakes is the temperature variation
over the depth, which is called thermal stratification. Starting in the early spring, the
lake attains a temperature of 4 OC and is more or less isothermal. During the summer
season, the surface water starts to become warmer because of increased solar radiation,
so that, gradually, a sharp temperature gradient, i.e., thermocline, is formed. The lake
remains thermally stratified during summer with a warm surface layer (epilimnion),
a thermocline, and a cold layer (hypolimnion). Though strong wind action tends to
lower or break the thermocline, the lake generally remains stratified during the entire
summer season. During the fall, as the air temperature drops, the net daily heat flux
at the water surface becomes negative, i.e., the lake loses heat daily. Hence, water
density in the epilimnion often becomes heavier than that in the hypolimnion and
causes convective mixing which, in combination with strong wind action, causes the
lake to become isothermal again in the winter. This process repeats itself annually. It
is important to know the location of thermocline at different times of the year so that
water can be withdrawn to a desirable height in deep lakes or reservoirs for various
agricultural and municipal uses.
Many processes are influenced by currents and temperature in lakes. For example,
3
the growth rate of all organic matter in lakes is governed by temperature (Goldman,
1979). The growth rate generally increases between some minimum temperature
and an optimum temperature, and decreases until it reaches maximum temperature.
Cooling water from power plants is mixed with surrounding water by the currents and
turbulent mixing which depend on the temperature field as well. Thus, predicting
currents and temperature are essential to understanding the transport of various
matters and their effects on the ecology.
In the early days, simple analytical models were used to study physical processes
in simplified conditions. For example, a set-up equation was used to predict the
storm surge height (Hellstrom, 1941). Since analytical models could not realistically
consider such effects as advection, complex geometry, and topography, they had been
applied to limited problems to understand certain basic processes.
Numerical models are valuable tools for simulating and understanding water
movement in lakes. Once a rigorously developed model is calibrated with measured
data, it can be used to estimate the flow near a man-made structure or to predict
the movement of contaminants including oil spill, sediments, etc. During the 1970s
and 1980s, vertically averaged two-dimensional numerical models, which can compute
only the depth average currents and surface elevations, were widely used because they
were simple and needed little computer time. However, since they could not give ac-
curate results for cases where the vertical distributions of currents and temperature
are required, three-dimensional models are needed.
Numerical modeling requires the discretization of the computation domain. Past
numerical models which were developed during 1970s generally used a rectangular grid
(for example, Sheng, 1975). However, to represent the complex geometry such as the
shoreline and the boundary between the vegetation zone and the open water in Lake
Okeechobee, a very fine rectangular grid is required. On the other hand, boundary-
fitted grids can be and have been used in recent models (for example, Sheng, 1987)
4
to represent the complex geometry with a relatively smaller number of grid points.
In some shallow lakes, aquatic vegetation can grow over large areas. The vegeta-
tion can affect the circulation significantly because it introduces additional friction on
the flowing water. For example, Lake Okeechobee has vegetation over an area which
covers 25% of the total lake surface. Because vegetation consists of stalks with differ-
ent heights and diameters, a representative diameter and height over each discretized
grid cell must be introduced in the model. Additional drag terms must be introduced
in the momentum equations to represent the form drag introduced by the vegetation.
The consideration of vegetation is necessary to compute the flow and transport of
phosphorus between the vegetation area and the open water.
Previously developed numerical models which were applied to deep lakes, e.g.,
the Great Lakes, cannot be readily applied to shallow lakes such as Lake Okeechobee,
since many shallow water processes are not included in these models.
Since 1988, with funding from the South Florida Water Management District and
U.S. Environmental Protection Agency, the Coastal and Oceanographic Engineering
Department of the University of Florida (under the supervision of Dr. Y. Peter
Sheng) has conducted a major study on the hydrodynamics and sediment dynamics
and their effects on phosphorus dynamics in Lake Okeechobee. The primary purpose
of the study was to quantify the role of hydrodynamics and sediments on the internal
loading of phosphorus and the exchange of phosphorus between vegetation zone and
open water. As part of the study, field data (wind, air temperature, wave, water
current, water temperature, and suspended sediment concentration) were collected
over two one-month periods in 1988 and 1989. Ahn and Sheng (1989) studied the wind
waves of Lake Okeechobee. Cook and Sheng (1989) studied the sediment dynamics
in Lake Okeechobee. This study focuses on the influences of vegetation and thermal
stratification on lake circulation. The objectives of this study are
(1) to obtain a general insight in the wind-driven circulation in Lake Okeechobee,
5
(2) to develop a numerical model which can simulate the effect of vegetation on
Lake Okeechobee circulation,
(3) to develop a numerical model which can simulate the thermal stratification and
its effect on circulation, and
(4) to determine the important factors for producing successful simulation of circu-
lation in large shallow lakes.
The literature review will be presented in the Chapter 2, after which the for-
mulation of the three-dimensional model will be given in Chapter 3. A vegetation
model will be explained in Chapter 4, and a thermal model follows in Chapter 5.
Finite-difference formulation will be presented in Chapter 6. After the model test in
Chapter 7, application to Lake Okeechobee will be described in Chapter 8. Finally,
a conclusion will be given in Chapter 9.
CHAPTER 2
LITERATURE REVIEW
2.1 Numerical Models of Lake Circulation
Sheng (1986) reviewed numerous models and suggested that models could be
classified according to numerical features (dimensionality, horizontal grid, vertical
grid, numerical scheme, etc.) and physical features (forcing function, free surface
dynamics, spatial scale, turbulence models, etc.). As an example, some models are
described in the following.
2.1.1 One-Dimensional Model
One-dimensional (1-D) models include a single spatial coordinate (longitudinal or
vertical). When the lake is elongated and well-mixed in the directions perpendicular
to the longitudinal axis, a longitudinal 1-D model can be used. A longitudinal 1-D
system of equations is derived by integrating the continuity and momentum equations
over the cross section. Sheng et al. (1990) developed a longitudinal one-dimensional
model of Indian River Lagoon. Sheng and Chiu (1986) developed a vertical one-
dimensional model for a location in Atlantic Ocean.
2.1.2 Two-Dimensional Model
Two-dimensional (2-D) models include horizontal 2-D models, which assume ver-
tical homogeneity, and vertical 2-D models, which assume transverse homogeneity.
The equations of motion are obtained by performing integration or averaging in the
vertical or transverse directions.
Hsueh and Peng (1973) studied the steady-state vertically-averaged circulation
in a rectangular bay by solving a Poisson equation with the method of successive
7
over relaxation (SOR). Their model included the terms of bottom friction, advection,
bottom topography, and lateral diffusion, while assuming steady state and homoge-
niety in density. The specification of the vertical eddy viscosity is not required in the
two-dimensional model but can give only depth-averaged velocities.
Shanahan and Harleman (1982) developed a transient 2-D model which assumed
vertical homogeneity. When the lakes are long, deep but relatively narrow, laterally
averaged 2-D model can be applied (for example, Edinger and Buchak, 1979).
2.1.3 Steady-State 3-D Models
An early study on wind-driven circulation was conducted by Ekman (1923) who
solved momentum equations analytically while neglecting the nonlinear terms. We-
lander (1957) developed a theory on wind-driven currents based on an extension
of Ekman's theory. After neglecting inertia terms and horizontal diffusion terms,
steady-state momentum equations were combined with the introduction of complex
variables. After applying boundary conditions, a solution was obtained in terms of
the imposed wind stress and unknown pressure gradient term. By introducing the
stream functions for vertically-integrated flow, the continuity equation could be sat-
isfied unconditionally. The final equation to be solved was reduced to a second-order
partial differential equation for stream function, 0, as follows:
V2 = a + bo + c (2.1)
Ox 9y
Once is found, the currents can be found by taking the derivatives, 2 and 0.
Gedney and Lick (1972) and Sheng and Lick (1972) applied Welander's theory
to Lake Erie. The equation for stream function was solved by the successive over
relaxation method. The agreement between the field data and model results was
good. Eddy viscosity was assumed to be constant but varies with wind speed.
Thomas (1975) used a depth-varying form of vertical eddy viscosity as follows:
v = o(1 + ) = ku.(h + z) (2.2)
8
where v is eddy viscosity, v0 is eddy viscosity near surface, and k is a constant (0.4).
2.1.4 Time-Dependent 3-D Models
Mode-Splitting
In order to solve the dependent variables with the unsteady three-dimensional
model, Simons (1974) used a so-called "mode splitting" method for Lake Ontario
while Sheng et al. (1978) used a somewhat different method for Lake Erie. Defining
the perturbation velocity fi = u U, i = v 7 where ui, are depth-averaged veloci-
ties, and u,v are instantaneous velocities, Sheng and Butler (1982) derived governing
equations for ii, v by subtracting the vertically-averaged equations from the momen-
tum equations. Therefore the solution procedure consists of an external mode, which
includes the surface elevation and U and F, and an internal mode, which includes ii,v
and temperature.
Time Integration of 2-D equations
Time integration is important for improving the efficiency of numerical models.
When the explicit method is used, the time step is limited by the Courant condition,
which is C t < 1. Therefore, explicit method is not desirable for long-term simu-
lations. Leendertse (1967) used the ADI (Alternate Direction Implicit) method to
simulate tidal currents in the southern North Sea. All terms in the continuity equa-
tion and pressure terms in the momentum equation were treated implicitly, while
the other terms were expressed explicitly. After factorization of the finite-difference
equations, the resulting unknowns are solved by inversion of tridiagonal matrices in
the x sweep and y sweep.
Vertical Grid
Various types of vertical grids are used in numerical models of lake circulation.
The earlier models generally used multiple vertical layers of constant fixed thickness
(z-grid) which do not change with time ("Eulerian grid") as used by Leendertse (1975).
This type of model needs a large number of vertical grid points in order to accurately
9
represent the shallow regions. More recent models use the so-called o--grid which
was originally applied in the simulation of atmospheric flow by Phillips (1957). This
vertical c-stretching uses the same number of vertical grid points in both the shallow
and deep regions, with the vertical coordinate defined as follows:
z ((x, y, t)
6 = (2.3)
h(x,y) + C(x,y,t)
where h(x, y) is the water depth, and ( is the water surface elevation. Governing
equations are transformed from the (x, y, z, t) coordinates to (x, y, r, t) coordinates
by use of chain rule and become somewhat more complex because of the extra terms
introduced by the stretching.
Other types of models (e.g., Simons, 1974) use a "Lagrangian grid" which consists
of layers of constant physical property (e.g., density) but time-varying thickness.
These models could resolve vertical flow structure with relatively few vertical layers.
However, parameterization of the interfacial dynamics is often difficult.
Horizontal Grid
One of the challenges in numerical models is the accurate representation of com-
plex geometry. Most models (e.g., Leendertse, 1967) use a rectangular uniform grid
to represent the shoreline of a lake or estuary. Thus, a large number of grid points
are needed to achieve a fine resolution near the shoreline or islands. Because com-
putational effort is directly related to the number of grid points, grid size should be
as small as possible to maintain required resolution near the interest area, so long
as the computational effort is not excessive. Therefore, to achieve a balance between
resolution and computational efficiency, a nonuniform grid method could be used.
Sheng (1975) used smaller grid size used near areas of importance but coarse grid
elsewhere.
Use of a boundary-fitted grid is another viable alternative. Johnson (1982) used
a boundary-fitted grid to solve depth-integrated equations of motion for rivers. Using
chain rules, he transformed the governing equations for a boundary-fitted grid which
10
was generated by using the WESCORA code developed by Thompson (1983). John-
son (1982), however, transformed the equations in terms of the Cartesian velocity
components.
The boundary-fitted grid has recently been adapted to three-dimensional nu-
merical models. Sheng (1986) applied tensor transformation to derive the three-
dimensional horizontal equations of motion in boundary-fitted grid in terms of the
"contravariant" velocity vectors (a "contravariant" vector consists of components
which are perpendicular to the grid line) and the water level. Sigma grid is used
in the vertical direction. The resulting equations in the boundary-fitted and sigma-
stretched grid are rather complex. However, numerous analytical tests were conducted
to ensure the accuracy of the model (Sheng, 1986 and Sheng, 1987). The model has
been applied to Chesapeake Bay (Sheng et al., 1989a), James River (Sheng et al.,
1989b), Lake Okeechobee (Sheng and Lee, 1991a, 1991b), and Tampa/Sarasota Bay
(Sheng and Peene, 1992). However, the earlier study on Lake Okeechobee (Sheng and
Lee, 1991a) did not consider thermal stratification in their model.
2.2 Vegetation Models
Vegetation can affect the aquatic life and also the water motion in the marsh area.
Early studies on the effect of vegetation on flow were conducted in the open channels.
Ree (1949) conducted laboratory experiments to produce a set of design curves for
vegetated channels. Kouwen et al. (1969) studied the flow retardance in a vegetated
channel in the laboratory and proposed the following equation:
U A
S= C, + C2n( ) (2.4)
where U is average velocity, u* is shear velocity, and C1 and C2 are coefficients. A
is a cross-sectional area of the channel, and A, is the cross-sectional area blocked by
the vegetation.
Reid and Whitaker (1976) considered the vegetation effect on flow as an ad-
ditional term, which is proportional to the quadratic power of the velocity, in the
11
depth-integrated momentum equation. Details of their vegetation models are given
in Chapter 4. Their model, however, considered only the linearized equations of
motion. In the present study, fully non-linear equations are considered.
Sheng (1982) developed a comprehensive vegetation model by including the effect
of vegetation on mean flow and second-order correlations in a Reynolds stress model.
Although the model was able to faithfully simulate the mean flow and turbulence in
the presence of vegetation, it was not used for the present study due to the extra
computational effort required when it is combined with a 3-D circulation model.
Roig and King (1992) formulated an equivalent continuum model for tidal marsh
flows. Neglecting leafiness, flexibility, and vegetation surface roughness, the net re-
sistance force due to vegetation is thought to be related to the following parameters:
V = f (p,, g u, 1,, d, s) (2.5)
where i is the viscosity of water, u is depth-averaged velocity, d is the average diameter
of vegetation, 1, is the vegetation height, and s is the spacing between vegetations.
Through a dimensional analysis,
T= pu f (F, R, ) (2.6)
where F is the Froude Number and R is the Reynold's Number.
To determine the function f, they conducted a simple flume experiment. For each
value of s/d, the dimensionless shear parameter I~ was plotted as a function of R
and F.
2.3 Thermal Models
Sundaram et al. (1969) used a one-dimensional vertical model to demonstrate
the formation and maintenance of thermocline in a deep stratified lake. The surface
boundary condition was given as follows:
aT
q, = -pcpKh & = K(Te T,) (2.7)
12
where q, is the heat flux, K is heat-exchange coefficient,Kh is vertical eddy diffusivity,
c, is specific heat of water, p is density of water, T, is an equilibrium temperature,
and T, is a surface water temperature. They assumed that the annual variation in
heat flux can be approximated by the cyclic form of the equilibrium temperature:
T, = Te + asin(wt + ) (2.8)
Following Munk and Anderson (1948), the eddy diffusivity Kh was expressed as the
product of the eddy diffusivity under neutral condition and a stability function, which
is one under neutral condition but becomes less than one under stable stratification
(positive Richardson number).
Price et al. (1986) studied the diurnal thermal cycle in the upper ocean using
field data and a vertical 1-D thermal model. Their measured data include currents,
temperature, and salinity, as well as meteorological data. Field data were collected
between April 28, 1980, and May 24, 1980, at about 400 km west of San Diego,
California.
Their major findings are the trapping depth of the thermal and velocity response
is proportional to r Q1/2, the thermal response is proportional to Q3/2, and the diurnal
jet amplitude is proportional to Q'/2, where Q is the heat flux and r is the wind stress.
They also simulated the diurnal thermal cycle using the vertical one-dimensional heat
equation coupled with the momentum equations.
Gaspar et al. (1990) determined the latent and sensible heat fluxes at the air-
sea interface using the inverse method. They stated that the total heat flux can
be divided into a solar part and a nonsolar part. While the solar radiation data
is usually available from direct measurement, the nonsolar part is usually indirectly
estimated from the meteorological data. However, this estimation of the nonsolar part
involves many empirical formulas and may contain large errors. Gaspar et al. (1990)
found that, by using the measured temperature data and solving the vertical one-
dimensional momentum equation and temperature equation, it is possible to estimate
13
the nonsolar heat flux more accurately. They called this method the "inverse method."
The major advantage of this "inverse" method is that it allows one to use the
usual satellite data (wind stress, surface insolation, and sea surface temperature) for
the estimation of heat flux at the ocean surface. The disadvantage of the method is
that horizontal advection effect is neglected in the analysis.
2.4 Turbulence Model
Flows in natural water bodies are often turbulent, although they can relaminar-
ize during periods of low wind and for tide. Although direct numerical simulation of
turbulence can now be performed for simple flow conditions, it is still computation-
ally prohibitive for practical applications in natural water bodies. Thus, "Reynolds
averaging", a statistical approach was taken by decomposing the flow variables into
a mean and a fluctuating part and averaging the equations over a period of time that
is large compared to the turbulent time scale. The resulting equations thus produced
are called Reynolds averaged equations. The Reynolds averaged mean flow equations
contain terms involving correlations of fluctuating flow variables (i.e., second-order
correlations) that represent fluxes of momentum or scalar quantities caused by tur-
bulent motion. The task of turbulence modeling is to parameterize these unknown
correlations in terms of known quantities.
Numerous turbulence models were developed for the parameterization of turbu-
lence. Some models are empirical, while others are based on more rigorous turbulence
theory. In this section, some available turbulence models are briefly reviewed with an
emphasis on the simplified second-order closure model.
2.4.1 Eddy Viscosity/Diffusivity Concept
By analogy with molecular transport of momentum, the turbulent stresses are
assumed proportional to the mean-velocity gradients. This can be expressed as
P = Vt( + (2.9)
Ox1 Oxi
14
where uu and u. are fluctuating velocity components in the x, and xj directions, while
u; and uj are mean velocity components, and vt is the turbulent eddy viscosity.
Similarly, an eddy diffusivity Kt can be defined:
-uV = K (2.10)
where b' and 0 are the fluctuating and mean temperature or salinity or scalar con-
centration, and Kt is the turbulent eddy diffusivity.
In order to close the system of mean flow equations for ui, it is necessary to obtain
an expression for vt in terms of known mean flow variables. Several options are given
in the following:
2.4.2 Constant Eddy Viscosity/Diffusivity Model
Earlier models used constant eddy viscosity/diffusivity models. Although this
allows easy determination of analytical solutions for the 1-D equation of motion and
easy programming, there are many disadvantages. Turbulence is spatially and tempo-
rally varying, hence constant eddy viscosity is not realistic. It is difficult to calibrate
the constant eddy coefficient model even if extensive field data exists.
2.4.3 Munk-Anderson Type Model
Prandtl (1925) assumed that eddy viscosity is proportional to the product of a
characteristic fluctuating velocity, V, and a mixing length, L. He suggested that
V = Imu and A, = I~2 The only parameter to be specified is length scale, Im,
which is assumed to be a linear function of z.
Following Prandtl, we can define the "neutral" vertical eddy viscosity as follows:
AvO [(=] + (2.11)
H Tz Tz
where Ao is assumed to be a linear function of z increasing with distance above the
bottom or below the free surface and with its peak value at mid-depth, while not
exceeding a certain fraction of the local depth. In the presence of strong waves,
15
turbulence mixing in the upper layers may be significantly enhanced. In such case,
the length scale Ao throughout the upper layers may be assumed to be equal to the
maximum value at mid-depth (Sheng, 1983).
When a lake is stratified, vertical turbulence is affected by buoyancy induced
by the vertical non-homogeneity. In this situation, vertical eddy coefficients should
be modified to account for this effect. This is parameterized by introducing the
Richardson number:
R = p T2 + -a 2 (2.12)
p 5z Qz) 9z
Ri is positive when flow is stable (a < 0) and when Ri is negative when flow is
unstable (2 > 0). Generally, eddy viscosity and eddy diffusivity are expressed as
follows:
A, = A,,0i(Ri); K,, = K,,o2(Ri) (2.13)
where q1 and 02 are stability functions and A,, and Kvo are eddy viscosity and eddy
diffusivity when there is no stratification. Stability functions have the following forms:
01 = (1 + i Ri)m'; 2 = (1 + a2Ri)"2 (2.14)
Based on comparing model results with field data, Munk and Anderson (1948) devel-
oped the following formula:
S= (1 + lORi)-1/2; 2 = (1 + 3.33Ri)-3/2 (2.15)
Many similar equations with different coefficients were suggested based on numerous
site-specific studies. These coefficients, however, are not universal, and care must be
taken when applying these formulae to a new water body where little data exist.
2.4.4 Reynolds Stress Model
One can obtain an equation for the time-averaged second-order correlations by
following the procedure: (i) decompose the dependent variables into mean compo-
nents and fluctuating components, (ii) substitute the decomposition into continuity,
16
momentum equation and heat equation, and (iii) take time-average of all equations.
For example, the resulting time-averaged equation for u'u' (e.g., Donaldson, 1973;
Sheng, 1982) is
au'u a9uu --Ouj --9, u ujp up
+ uk = -iUk -k UjUk +t 9g
,t k k k Po Po
'2ek [ -2- (aukkuiuj) (2.16)
2e~iktfkutuj 2EjkiLZUkUi Xk
ui ap uj 9p a u uuj u u
+v 2v
p zxj p Oxi aXkaxk Oxk 9zk
Similar equations for up' and p'p' can be obtained. Unresolved third-order cor-
relations and pressure correlations are modeled using the simplest possible forms
(Donaldson, 1973).
u, Ou, bS.q3 avuuiu
Uvxu;xk +A (2.17)
ZXk OXk 3A A2
p ( Su u q ar q2
(x+-2(u j6 8 -) (2.18)
(,,, up Pu p aqA uiuA 3
k u + up OqAcuu, (2.19)
-5+ 7 =va-k---.
where q is the total fluctuating velocity and A is the turbulence macroscale. The
model constants (a, b, and v,) are determined from a wide variety of laboratory data
(Lewellen, 1977). Thus, a full Reynolds stress model consists of six equations for
velocity fluctuations u u', three equations for the scalar fluxes, u p', and one equation
for the variance, p'p'. Considering the required computer storage and CPU time for
the turbulence models, it is desirable to use a simplified form of the Reynolds stress
model.
2.4.5 A Simplified Second-Order Closure Model: Equilibrium Closure Model
The complete second-order closure model is too complicated to be used in a three-
dimensional model. A simplified second-order closure model can be developed with
17
the following assumptions: 1) Second-order correlations have no memory effect. That
means correlations at the previous time have no effect on correlations at the next time.
Therefore, D = 0. 2) Correlations at a point do not affect the value at another point.
Therefore, all the diffusive terms are dropped. These conditions are approximately
true if the time scale of turbulence is much less than the mean flow time scale and
the turbulence does not vary significantly over the macroscale, i.e., the turbulence is
in local equilibrium. Then the remaining equations become as follows (Sheng, 1983):
9uj ui ujp u'p
0 = -uiuk5 UUk 9X 9i 9j-
=k -uk k Po Po
2eikefkUtU'j Cjklfet'ku (2.20)
q q2 q3
A(u "U 6s'T
S p -,-, ui giP'P'
0 = -ujiuj 7 .P Zj P
'Oxj u Po
2Eijk jUkp' -0.75q (2.21)
'- P 0.45qp'p'
0 = 2uP + (2.22)
These algebraic equations can be solved with ease, once the mean flow conditions are
known. In order to complete the system of equations, q and A need to be solved
following the procedure described in Sheng et al., 1989b.
The above "equilibrium closure" model was applied to the Atlantic Ocean (Sheng
and Chiu, 1986), Chesapeake Bay (Sheng et al., 1989a) and the James River (Sheng
et al., 1989b). More details of the model will be given in Chapter 8.
2.4.6 A Turbulent Kinetic Energy (TKE) Closure Model
To introduce some dynamics of turbulence into the simplified second-order closure
model, one can add a dynamic equation for q2(q2 = u'u' + v'v' + w'w'), which is twice
the turbulent kinetic energy (Sheng and Villaret, 1989). This TKE closure model has
18
been applied to James River (Choi and Sheng, 1993) and Tampa Bay (Schoellhamer
and Sheng, 1993).
Two other turbulence models, which are based on the so-called k model (Rodi,
1980), are described in the following:
2.4.7 One-Equation Model (k Model)
Using the eddy viscosity/diffusivity concept, the choice of velocity scale can be
v/k, where k = (u2 + v2 + w2)/2 is the kinetic energy of the turbulent motion. When
this scale is used, the eddy viscosity is expressed as
vt = c'VkL (2.23)
where c' is an empirical constant and L is the length scale. To determine k, an
equation is derived from the Navier-Stokes equation as:
Ok Ok a u u' p' u, l u, u,
+ = ~j[( + -)] uiU Aflgiui<- v2 (2.24)
2t axi axi 2 p Oxj Oxyj Oxj
To obtain a closed set of equations, diffusion term and dissipation term must be
modelled. The diffusion flux is often assumed proportional to the gradient k as
uj' p vt Ok
ui( ) = (2.25)
2 p
where Oak is an empirical diffusion constant. The dissipation term e, which is the last
term of Eq. 2.24, is usually modelled by the expression
k3/2
S= CD- (2.26)
L
The length scale, L, needs to be specified to complete the turbulence model. Usually
L is determined from empirical relations.
2.4.8 Two-Equation Model (k e Model)
To avoid the empirical specification of length scale, another equation for the
dissipation e is needed. Then eddy viscosity and eddy diffusivity are expressed as
Vt = c,- (2.27)
f
19
Ft = -- (2.28)
at
where c, is an empirical constant and at is the turbulent Prandtl/Schmidt number.
An equation for e is derived from the Navier-Stokes equation and is
de c ae 0 t Oe e 62
+ u,- = ( ) + c (Pc3,G) c2c (2.29)
(9t xi Oxi 9ao l c+x k k
where P and G are the stress and buoyancy term in Equation 2.24, respectively, and
cie, c2c and c3, are empirical model coefficients.
2.5 Previous Lake Okeechobee Studies
Numerous studies on Lake Okeechobee have been performed, but most of them
focused on water quality.
Whitaker et al. (1975) studied the storm surges in Lake Okeechobee while con-
sidering the vegetation effect in the western marsh area, by using a two-dimensional,
vertically integrated model. They simulated the seiche in the lake during the 1949
and 1950 hurricanes. In their study, the bottom friction coefficient was parameterized
as a function of depth to achieve better agreement of storm surge height.
Schmalz (1986) investigated hurricane-induced water level fluctuations in Lake
Okeechobee. His study consisted of two parts: a hurricane submodel and a hydrody-
namic submodel. The hurricane submodel used hurricane parameters such as central
pressure depression, radius to maximum winds, maximum wind speed, storm track,
storm forward speed, and azimuth of maximum winds, and determined the wind and
pressure field that were used as forcing terms for the hydrodynamic submodel.
The hydrodynamic submodel solved the depth-averaged momentum equations
and continuity equation. Finite-difference method was used for the numerical solution.
For treatment of marsh area, an effective bottom friction which relates the Manning's
n to water depth and canopy height was used. The 2-D hydrodynamic model can
resolve flooding and drying: during strong wind conditions such as a hurricane, a
portion of the lake can become dried because of the excessive setdown by wind, while
20
other portion of the lake can become flooded because of excessive setup by wind.
A three-dimensional Cartesian-grid hydrodynamic and sediment transport model
for Lake Okeechobee was recently developed (Sheng et al., 1991a; Sheng, 1993). In
addition, these models were extended to produce a three-dimensional phosphorus
dynamics model (Sheng, et al., 1991c). These models use the simplified second-order
closure model and the sigma- stretched grid, however, did not consider the effects of
vegetation and thermal stratification
2.6 Present Study
The present work focuses on the study of effects of vegetation and thermal strat-
ification on wind-driven circulation in Lake Okeechobee. As will be shown later, the
three-dimensional curvilinear-grid model (CH3D) will be significantly enhanced to
allow accurate simulation of the observed circulation. Model features are compared
with model features of some previous lake studies in Tables 2.1 and 2.2. It is apparent
that the 3-D model developed in this study is more comprehensive than those used
in previous lake studies.
Table 2.1: Selected features of lake models.
Author Dimensionality Type of Temporal Turbulence Advection
Model Dynamics
Welander 3-D AN T.D. A No
1957
Liggett 3-D F.D. T.D. A No
1969
Lee + Liggett 3-D F.D. S.S. A No
1970
Liggett + Lee 3-D F.D. S.S. A No
1971
Gedney + Lick 3-D F.D. T.D. A No
1972
Goldstein + Gedney 3-D A.N. B No
1973
Sengupta + Lick 3-D F.D. T.D. D Yes
1974
Simons 3-D F.D. T.D. B
1974
Sheng 3-D F.D. S.S. A No
1975
Thomas 3-D F.D. S.S. B No
1975
Whitaker et al. 2-D F.D. T.D. Yes
1975
Witten + Thomas 3-D F.D. S.S. C No
1976
Lien + Hoopes 3-D F.D. S.S. A No
1978
Schmalz 2-D F.D. T.D. Yes
1986
Sheng + Lee 3-D F.D. T.D. E Yes
1991a _
* F.D. : Finite difference
* AN : Analytic
* S.S. : Steady state
* T.D. : Time dependent
* A: Constant
* B : Dependent on wind
* C : Exponential form
* D : Munk-Anderson type
* E : Simplified second-order closure model
Table 2.2: Application features of lake models.
Author Basin Dimen- Mean Forcing Vege- Grid
tion Depth W H R station Hor. ver.
Liggett Idealized Yes No No No U
1969 Basin
Lee + Liggett Idealized Yes No No No U
1970 Basin
Liggett + Lee Idealized Yes No No No U
1971 Basin
Gedney + Lick Lake Erie 400 km 20 m Yes No Yes No U
1972 100 km
Sengupta + Lick Squire 1.89 m Yes Yes No No N
1974 Valley
Simons Lake Yes Yes No No U
1974 Ontario
Sheng Lake Erie 400 km 20 m Yes No Yes No N oa
1975 100 km
Thomas Idealized Yes No No No U
1975 Basin
Whitaker et al. Lake 57 km 2.5 m Yes No No Yes U
1975 Okeechobee 60 km
Witten + Thomas Idealized 300 x Max Yes No No No
1976 Basin 87 km 180 m
Lien + Hoopes Lake Yes No No No U
1978 Superior
Schmalz Lake 57 km 2.5 m Yes No No No U
1986 Okeechobee 60 km
Sheng + Lee Lake 57 km 2.5 m Yes No No Yes C a
1991a Okeechobee 60 km
* W: Wind
* H: Heating
* R: River
* U : Uniform Cartesian grid
* C : Curvilinear grid
* N : Non-uniform Cartesian grid
* : Vertically stretched grid
CHAPTER 3
GOVERNING EQUATIONS
3.1 Introduction
This chapter presents the basic equations which govern the water circulation in
lakes, reservoirs, and estuaries. Because the details can be found in other references
(e.g., Sheng, 1986; Sheng, 1987; Sheng et.al., 1989c), the governing equations are
presented here without detailed derivations.
3.2 Dimensional Equations and Boundary Conditions in a Cartesian Coordinate System
The equations which govern the water motion in the water bodies consist of
the conservation of mass and momentum, the conservation of heat and salinity, and
the equation of state. Because Lake Okeechobee is a fresh water lake, the salinity
equation is not considered. The following assumptions are used in the Curvilinear
Hydrodynamic Three-dimensional Model (CH3D) model.
(1) Reynolds averaging: Three components of velocity, pressure, and temperature
are decomposed into mean and fluctuating components and time-averaged.
(2) Hydrostatic assumption: Vertical length scale in lakes is small compared to
the horizontal length scale, and the vertical acceleration is small compared with the
gravitational acceleration.
(3) Eddy viscosity concept: After time-averaging, the second-order correlation
terms in the momentum equation are turbulence stresses, which are related to the
product of eddy viscosity and the gradient of mean strain.
(4) Boussinesq approximation: Density variation of water is small, and variable
density is considered only in the buoyancy term.
3.2.1 Equation of Motion
With above assumptions, the equations of motion can be written in a right-handed
Cartesian coordinate system as follows:
Ou Ov Ow
+ -+-=0 5(3.1)
9x 9y 9z
Ou Ou2 Ouv Ouw 1 9p 0 (A u\
+ -- + + -- =f- +- AoH
9t y P o x O Ox
9 9 9 ((o 9
+ (AH ) + aZ (3.2)
9v Ouv 9v2 Ovw 1 p O ( av\
+- + -+ --fu--+ AH
at Ozx y o z Po Y Ox Ox
+ AH-) +- A, O (3.3)
ap
z -pg (3.4)
OT OuT OvT OwT (KOT
S+ + + y+ z -
+ KH + K, a z (3.5)
p = p(T, S) (3.6)
where (u, v, w) are velocities in (x, y, z) directions, f is the Coriolis parameter defined
as 2f sine where 0 is the rotational speed of the earth, q is the latitude, p is density,
p is pressure, T is temperature, (AH, KH) are horizontal turbulent eddy coefficients,
and (A,, K,) are vertical turbulent eddy viscosities.
For the equation of state, Eqn. 3.6, there are many different formulae that can be
used. For the present study, the following equation given by Eckart (1958) is used:
p = (1 + P)/(a + 0.698P)
(3.7)
25
P = 5890 + 38T 0.375T2 + 3S (3.8)
a = 1779.5 + 11.25T 0.0745T2 (3.8 + 0.01T)S (3.9)
where T is water temperature in degree C, S is salinity in ppt, and p is density in
g/cm3.
Besides governing equations, boundary conditions should be specified.
3.2.2 Free-Surface Boundary Condition (z = 77)
(1) Kinematic boundary condition:
w = + u- + v (3.10)
Ot lo dy at
(2) Surface heat flux:
9T
q = K ,,z = K(T,-T,) (3.11)
where T, is the lake surface temperature, Te is the equilibrium temperature, K is a
heat transfer coefficient, and q is positive upward. (3) Surface stress:
9u 9v
Az= rx, Az = (3.12)
where the wind stresses rT and r, must be specified.
3.2.3 Bottom Boundary Condition (z = -h)
(1) Heat flux is specified as zero,i.e., T = 0
(2) Quadratic bottom friction law is used, i.e.,
Tx = pCdu,2 + 12u1,r = pCd Ul2 +v12v1 (3.13)
where ul and vl are velocity components at the first grid point above the bottom.
3.2.4 Lateral Boundary Condition
(1) Heat flux is assumed zero,i.e., T = 0
(2) No flow through boundary,i.e., u = 0 or v = 0
26
3.3 Vertical Grid
No natural water bodies have strictly flat bottoms. Therefore to represent the
variable bottom topography, a- stretching is used by defining a new variable o- :
z- (x, y,t)
a = (3.14)
h(x, y) + ((x, y, t) (
The advantage of a-stretching is that the same vertical model resolution can be main-
tained in both shallow and deeper parts of a lake. The disadvantage is that it in-
troduces additional terms in the equations. Details of r-stretching can be found in
Sheng and Lick (1980) and Sheng (1983).
3.4 Non-Dimensionalization of Equations
By introducing reference values, the governing equations can be non-dimensionalized.
The purpose is to make it easier to compare the relative importance of each term.
The following relations were used (Sheng, 1986).
(u*, v',w*) = (u,v, wX,/Z) /U,
(x',y',z*) = (X,y,zX,r/Z,) /X,
(7;,T;) = (r", -r')/lpfZU,
t* = tf
q: = To/(T, T) q/poc,fzTo
(* = g/lfUX, = C/S, (3.15)
P" = (P po)/(Pr Po)
T* = (T- To)/(T, To)
AH = AH/AHr
A* = A,,/A,,
KI = KH/KHr
K, = Kv/Kv
w* = wXIU,
27
where variables with asterisks are non-dimensional variables and variables with r are
the reference values.
3.5 Dimensionless Equations in a-Stretched Cartesian Grid
The transformation relations from a Cartesian coordinate (x, y, z) to a vertically
stretched Cartesian coordinate (x,y, c) may be found in Sheng (1983). Using the
relation presented in the previous section, the following dimensionless equations are
obtained:
ac 8Hu 8Hv 9w
- + 5 + --y + H- =0
at Bz 9y 8a
1 8Hu ( E, (A u
H at xz H- a \ a +)
H \Ox + ay + ao
+ E A, + a + H.O.T.
Ox Eu 9y ay)
Ro r oOrp 0H I]
[H -da + aj pda + ap
Fr aQx ax
ac E + a ,( au
ax H2a ka, a+
1 aHv (_ E, a (Av \
H 9t ~ Ty TH2,a u
Ro (Huv aHvv aHvw
H x + ay a
+ EH AH + y aA + H.O.T.
ax x) ay 5ay
Ro [H -d + aH o pdH+ +p)]
Fr2[ J ay "y a r
yo E a o 9a
a +E, a av + By
ay H2 aa ao)
1 aHT E, a T\
H at PrH, H2+a Iv a
Ro aHuT 9HvT aHowTN
H ax +y + a
(3.16)
(3.17)
(3.18)
(3.19)
28
+ Kg + ay H ) a + H.O.T.
Pr x ax ay a 9
p = p(T, S) (3.20)
where H is total depth, / = gZr/X2f2 and H.O.T. is higher order terms.
3.5.1 Vertically-Integrated Equations
The CH3D model can solve the depth-integrated equations and the three-dimensional
equations. The vertically integrated momentum equations are obtained by integrating
the three-dimensional equations from bottom to top.
-+( -+ -= 0 (3.21)
at a9x ay
-a H +T~ -nb + V
at ax
[ f UU ) a UV\1
o [a + (3.22)
EL9 AOU a aLU
+ EH [ AH + a AmHa
a x T x ay ay)
Ro H2 ap
Fr~ 2 6x
-V a(
O--x
1= H +r,,- 7,by U
at ay
I 9 UV\ 9 rVV\
[- Ro \+ H o(3.23)]
+ EH AH- + AH
[ax ax] Ty ay
Ro H2 ap
Fr2 2 Oy
S- H- + D,
9y
3.5.2 Vertical Velocities
(1 + a) 9( 1 (OHu OHv
S= H Ot H d-i + a (3.24)
w = Hw (1 + ) d( ( h dh
w =dt H + v (3.25)
) dt 8x T y
where U = f01 udo, V = f1 vdo, r,,, r., are wind stresses at the surface and Trb, Tby
are bottom stresses.
3.6 Generation of Numerical Grid
3.6.1 Cartesian Grid
In order to numerically solve the governing equations, finite difference approxima-
tions are introduced to the original governing equations, and solutions are obtained
at discrete points within the domain. Therefore, a physical domain of interest must
be discretized. When a simple physical domain is considered, cartesian grid can be
used and hence grid generation and development of finite-difference equations are
relatively easy.
Unfortunately, most physical domains in lakes or estuaries are complex. Pre-
viously rectangular grid was widely used. This method has such disadvantages as
inaccuracies at boundaries and complications of programming due to unequal grid
spacing near boundaries.
3.6.2 Curvilinear Grid
To better resolve the complex geometries in the physical domain, boundary-fitted
(curvilinear) grid can be used. In general, a curvilinear grid can be obtained by use
of (1) algebraic methods, (2) conformal mapping, and (3) numerical grid generation.
Algebraic grid generation uses an interpolation scheme between the specified
boundary points to generate the interior grid points. This is simple and fast com-
putationally, while the smoothness and skewness are hard to control. Conformal
30
mapping method is based on complex variables, so the determination of mapping
function is a difficult task. Therefore, many practical applications rely on numerical
grid generation techniques.
3.6.3 Numerical Grid Generation
Partial differential equations are solved to obtain the interior grid points with
specified boundary points. Thompson (1983) developed an elliptic grid generation
code (WESCORA) to generate a two-dimensional, boundary-fitted grid in a complex
domain.
To help understand the physical reasoning of this method, consider a rectangular
domain. When the temperature is specified along the horizontal boundary, then
the temperature distribution inside can be obtained by solving the heat equation.
Therefore, isothermal lines can be drawn. Also, other isothermal lines can be obtained
with the specified temperature in the vertical direction. By superimposing these
isothermal lines, intersection points of isothermal lines can be considered as grid
points.
WESCORA solves Poisson equations with same idea in a complex domain. Con-
sider the following set of equations (see Figure 3.1):
G, + ~ = P (3.26)
'rz+ + 7yy = Q (3.27)
with the following boundary conditions:
S= (x,y) on 1 and 3
7 = constant (3.28)
S= constant on 2 and 4
q = r7(x,y) (3.29)
31
where the functions P and Q may be chosen to obtain the desired grid resolution and
alignment. In practice, one actually solves the following equations which are readily
obtained by interchanging the dependent and independent variables in Eqs. 3.26 and
3.27:
ax-2 203x, + xz,, + aPzx + ^Qx, = 0 (3.30)
aye 2/3y, + yy, + aPyW + ^Qy, = 0 (3.31)
where
2 2
/ = x7 +y Y
1 = 2
P = (x7 + y2) (3.32)
1 2
Q = 7(x +yi)Q
J = xzy, Xzy
with the transformed boundary conditions:
x = fi(,rl 7) on i = 1 and 3
Y = gi(, ri) (3.33)
x = fi(,i, 7) on i = 2 and 4
Y = gi(i, 7) (3.34)
3.7 Transformation Rules
Generations of a boundary-fitted grid is an essential step in the development of
a boundary-fitted hydrodynamic model. It is, however, only the first step. A more
important step is the transformation of governing equations into the boundary-fitted
coordinates. A straightforward method is to transform only the independent variables,
C ((.y) y,
or
q (,.y) yS
\ ---3.
Y. Z,
PROTOTYPE
Tf
T1RAS. Y M
TRANSFORMED
Figure 3.1: A computational domain and a transformed coordinate system.
* (I ,Xa)
yjt( It. -I a)
S 3'
33
i.e., the coordinates, while retaining the Cartesian components of velocities. Johnson
(1982) developed such a 2-D vertically-integrated model of estuarine hydrodynamics.
The advantage of the method is its simplicity in generating the transformed equations
via chain rule. The dimensional forms of the continuity equation and the vertically-
integrated momentum equations are shown by Eqs. (20) and (21) in Appendix A of
Sheng (1986). The resulting equations, however, are rather complex. Even when an
orthogonal or a conformal grid is used, the equations do not become any simpler. Ad-
ditional disadvantages are (1) the boundary conditions are quite complicated because
the Cartesian velocity components are generally not aligned with the grid lines, (2)
the staggered grid cannot be readily used, and (3) numerical instability may develop
unless additional variables (e.g., surface elevation or pressure) are solved at additional
grid points, (Bernard, 1984, cited in Sheng (1986)).
To alleviate the problems mentioned in the previous paragraph, Sheng (1986)
chose to transform the dependent variables as well as the independent variables.
Equations in the transformed coordinates (, 7) can be obtained in terms of the con-
travariant, or covariant, or physical velocity components via tensor transformation
(e.g., Sokolnikoff, 1960). As shown in Fig. 8 of Appendix A of Sheng (1986), the
contravariant components (u') and physical components u(i) of the velocity vector
in the non-Cartesian system are locally parallel or orthogonal to the grid lines, while
the covariant components (ui) are generally not parallel or orthogonal to the local
grid lines. The three components are identical in a Cartesian coordinate system. The
following relationships are valid for the three components in a non-Cartesian system
are
ui = (gi)-1/2u(i) (no sum on i) (3.35)
ui = (gii)-'/2lgi(j) (no sum on i) (3.36)
u(i) = iifi
(3.37)
34
where g is the diagonal element of the metric tensor gij:
9 xj 6Se.6 (3.38)
which for the two-dimensional case of interest is
9i + yC XCX17 + yy07 911 912
g7XC + Y1 x2 +2 921 922( (3.39)
(- x+ + y~ye x y Y- g21 g22
The three components follow different rules for transformation between the prototype
and the transformed plane:
U' = -- (3.40)
Oxi
Ui = uj (3.41)
(i) = ( 2 uz) (3.42)
where the unbarred quantities represent the components in the prototype system,
while the barred quantities represent the components in the transformed system.
3.8 Tensor-Invariant Governing Equations
Before transforming the governing equations, it is essential to first write them in
tensor-invariant forms, i.e., equations which are independent of coordinate translation
and rotation. For simplicity, unbarred quantities are used to denote the variables in
the transformed system unless otherwise indicated.
Following the rules described in the previous paragraph, the following equations
are obtained (Sheng, 1986):
P a
+ o xk (V Huk) = 0 (3.43)
1 dHuk
H !k g "jun
H 19t ~ gne _
Ro HhukW
o ( Hueuk), + hw (3.44)
H I \
E, a uk ak
+ AH2 + EAHUk
Ro F 'M.
iH p!da +H!k (f pd + ap)
where O/Oxk is the partial derivative, ge, is the metric tensor while go = J =
zy, x,y is the determinant of the metric tensor, Uk is the contravariant veloc-
ity, ( ),t represents the covariant spatial derivative, !k represents the contravariant
spatial derivative, and ekj is the permutation tensor and
12 1
Vjo
621 = (3.45)
Vo
11 22 = 0
The covariant and contravariant differentiations are defined by
: = + Di ju (3.46)
S!k = gkmS,m (3.47)
where :j represents partial differentiation and D,, represents the Christoffel symbol
of the second kind:
D k = giDnk (3.48)
where g'" represents the inverse metric tensor, hi,, and Dnjk is the Christoffel symbol
of the first kind:
1
Dijk = (gij:k + gik:i gjk:i) (3.49)
3.9 Dimensionless Equations in Boundary-Fitted Grids
Expanding the equation 3.41 and 3.42, the following equations are obtained:
(t +
V9.
1 9Hu
H at
+
(3.50)
-( + 9 120 ~ g U12 922
9 + O ) u + -0v
H [(Huu) + (Huv) + (2D11 + D 2)Huu
+ (OHuw
+ (3D12 + D22)Huv + D2gHvv + 9H
TUZZdcr
H2 Eo (A OU)
+ H--a--a)o, O
Ro r, f H 19p
FrLH-f t, \Y I
(3.51)
+p- da
12i d
all
[/ a lH aH\
+ AH(Horizontal Diffusion)
+ EHAH(Horizontal Diffusion)
(g21
Ro
H [
~IIH
+ 22
077)
(H + (H ) + D
(Huv) + -(Hvv) + D1 Huu
a( d77
+ (D1, + 3D2)Huv + (D12 + 2D2)Hvv]
+ -- 2 A 2
H2 m I o-( v
H2O* O9aj
Ro H\
F,2 I
21 aH
+ (9 -+
( 21 0 1
22H
g9 a
-22 )P da
( 9+
pda +ap]
+ EHAH(Horizontal Diffusion)
where the horizontal diffusion terms are listed in Sheng(1986). The temperature
equation can be obtained according to the same procedure as
1 OHT
H at
E,, 0
Pr, H2 r,
K-}
FaTJ
1 OHv
H Ot
(3.52)
g"
921
a( Hv) =0
0
Ro 1 9 8 Ro8HwT
Ro [ (v [ g(HuT) + (1HvT)\ RoH
H rgL;r VO H 0a
+ H-[I T,11 + 12T,1,2
PrH
+ g2'T2,1 + g2T2,2] (3.53)
3.10 Boundary Conditions and Initial Conditions
3.10.1 Vertical Boundary Conditions
The boundary conditions at the free surface (a = 0) are
( uOu -H
A, = u,('7, ),=) (3.54)
aT HPr,
=q-
5a E,
The boundary conditions at the bottom (a = -1) are
(au Qv H
A, a = -(rb, 7rbn)
= HZCd [gnu2 + 2gt2uiv1 + g22v (l, Vl)
Avr
9T
= 0 (3.55)
ao"
where ul and vl are the contravariant velocity components at the first grid point
above the bottom.
3.10.2 Lateral Boundary Conditions
Due to the use of contravariant velocity components, the lateral boundary condi-
tions in the ({, 77, a) grid are similar to those in the (x, y, a) system. Along the solid
boundary, no-slip condition dictates that the tangential velocity is zero, while the slip
condition requires that the normal velocity is zero. When flow is specified at the open
boundary or river boundary, the normal velocity component is prescribed.
3.10.3 Initial Conditions
Initial conditions on vectors, if given in the Cartesian or prototype system, such
as the velocity and the surface stress, must be first transformed before being used
38
in the transformed equations. Thus, the surface stress in the transformed coordinate
system is given by
S= 1+ ar2 (3.56)
-2 = + + 2 (3.57)
where T1, T2 are the contravariant components of the stress in the transformed system
and rT, r2 are the contravariant components in the Cartesian system. Note that in
the Cartesian system, the contravariant, covariant, and physical components of a
vector are identical. The contravariant components of the initial velocity vectors can
be transformed in the same manner to obtain the proper initial conditions for the
transformed momentum equations.
CHAPTER 4
VEGETATION MODEL
4.1 Introduction
The western portion of Lake Okeechobee is covered with an extensive amount
of vegetation. The vegetation can affect the circulation in several different ways.
First of all, wind stress over the emergent vegetation is reduced below that over the
open water. Furthermore, the submerged vegetation introduces drag force to the
water column. Because most of the vegetation stalks are elongated cylinders without
large leaf areas, the drag force is primarily associated with the profile drag (or form
drag) instead of the skin friction drag. The profile drag can reduce the flow and is
proportional to the "projected area" of vegetation in the direction of the flow.
The presence of vegetation also can affect the turbulence in the water column.
The characteristic sizes of the horizontal and vertical eddies generally are reduced by
the vegetation. This usually leads to a reduction of turbulence, although some wake
turbulence may be generated on the downstream side of vegetation.
In order to simulate the effects of vegetation, several approaches have been under-
taken in previous investigations. For example, Saville (1952) and Sheng et al. (1991b)
used an empirical correction factor to simulate the reduction of wind stress over the
vegetation area. Sheng et al. (1991b) also adjusted the bottom friction coefficient
over the vegetation area. For simplicity, however, Sheng et al. (1991b) did not include
the effect of vegetation on mean flow and turbulence in the water column, because the
primary focus of that study was the internal loading of nutrients from the bottom sed-
iments in the open water zone. Whitaker et al. (1975) developed a two-dimensional,
vertically-integrated model of storm surges in Lake Okeechobee. The profile drag cre-
40
ated by the vegetation was included in the linearized vertically-integrated equations of
motion, which did not contain the nonlinear and diffusion terms. Sheng (1982) devel-
oped a comprehensive model of turbulent flow over vegetation canopy by considering
both the profile drag and the skin friction drag in the momentum equations in addi-
tion to the reduction of turbulent eddies and the creation of turbulent wake energy.
Detailed vertical structures of mean flow and turbulence stresses were computed by
solving the dynamic equations of all the mean flow and turbulent quantities. Model
results compared well with available mean flow and turbulence data in a vegetation
zone.
For the present study, due to the lack of detailed data on vegetation and mean
flow and turbulence in the vegetation zone, a relatively simple vegetation model which
is more robust than Whittaker et al.'s model yet simpler than Sheng's 1982 model is
developed. Due to the shallow depth in the vegetation zone, it is feasible to treat the
water column with no more than two vertical layers. When the height of vegetation is
greater than 80% of the total water depth, the flow is considered to be one-layer flow,
i.e., the entire water column is considered to contain uniformly distributed vegetation.
When the height of vegetation is between 20% and 80% of the total water depth, the
flow is considered to be two-layer flow, i.e., the water column consists of a water layer
on top of a vegetation layer. The vegetation effect is neglected when the height of
vegetation is less than 20% of the total depth. The profile drag introduced by the
vegetation can be formulated in the form of a quadratic stress law:
Canopy = pCdU1 Uil AN (4.1)
where u is the vertically averaged velocity in the vegetation layer (layer 1), p is the
density of water, A is the projected area of vegetation in the direction of the flow, N is
the number of stalks per unit horizontal area, and ed is an empirical drag coefficient.
Tickner (1957) performed a laboratory study. Strips of ordinary window screen 0.1
foot in height were placed across a channel to simulate a vegetative canopy. Using
41
Tickner's experimental results, Whitaker et al. (1975) calculated cd(1.77) which was
used in this study. Roig and King (1992) showed Cd is a function of Froude Number,
Reynolds Number, vegetation height, spacing, and diameter of vegetation. As the
water level changes, the flow regime over a vegetation area may change from one-
layer to two-layer flow, and vice versa.
4.2 Governing Equations
Let us consider an x, y, z coordinate system with the velocity components in the
(x, y, z) directions as (u, v, w). The lower layer (layer I) of the water column is covered
with vegetation, while the upper layer (layer II) is vegetation-free (Figure 4.1).
Flow in the vegetation layer (Ul, vi) and flow in the vegetation-free layer (u2, v2)
both satisfy the equations of motion.
4.2.1 Equations for the Vegetation Layer (Layer I)
Oul Ou2 Ouivvl Oulwl 1 api 0 A i
-ti +x+ fvl + AH
at xz y + z p x ax Ox
9 Buil 9 9u
4+ A- A + ar A- (4.2)
Ovi Ouivl Ov 9wi vw1 -pi iA
ty foy O z z J
a' + O ( + + = fu p + AH
0A 0 l (4.3)
8p
S= -Pg (4.4)
Integrating Equation (4.4) vertically:
P = Pa + pg(( z) (4.5)
Integrating Equations (4.2) and ( 4.3) vertically from z = -h to z = -f:
a-, + a (L1] U y L) a + gLl(C = f V + (,ri 7bx Fcx)
at axL ay L,
+ aZ[-[A,,, aU1 + a [AHr2-F]
+ A A 1 + AH ] (4.6)
+ x azx Ty yJ
Y
S(x,y,t)
Figure 4.1: Schematics of flow in vegetation zone.
at V x+ L (
T Tx LL /
1
gL1(y = -f U + -(ri, Tby Fcy)
p
oa [A oy [A 'y
O+ '9X TY ayj
(4.7)
where L1 = h i and U1 and V1 are the vertically-integrated velocities within the
vegetation layer:
(U1, Vi) f (ul,v I)dz
(4.8)
7bT and Tb, are bottom stresses, Ti, and ri, are interfacial stresses between layer I and
layer II, and F, and F, are form drags due to the vegetation canopy.
4.2.2 Equations for the Vegetation-Free Layer (Layer II)
0u2 Ou2 Ou2v2
at+ ,x+ oy
at 8x By
Ov2 au2vy
V2 + 22
't ax
vy
+ O+
'y
OU2W2
Oz
aV2W2
Ou9Zw
1op2 a A _u21
Sfv2 P+ AHf
p yy -+ zAS
+ AH- +- A-j-
9y L y 9z 9zV
1 0P2 0 Fv2l
--fu2 TA-,
p dy Oz L[ ax
+ -[ AH + A--]
(4.9)
(4.10)
Integrating Equation (4.9) and (4.10) vertically from z = -e to z = ( and defining
L2 = + C:
a
Y
U2)V2
L2
+ gL2( = fV2 + -(5r. r7.)
P
8 AU, [AU2
9az 8[ Dy J i/ y
(4.11)
OV2 + 1 U2V2
at Ox L 2
O (V2\
2 +
+ T(L2r-I
9y L2
1
gL2(y = -fU2 + -(r, Ty)
P
9 axQ T ay
S \A-- +-A -
(4.12)
where U2 and V2 are vertically-integrated velocities within the vegetation-free layer:
(U2,V2)= f(u2,v2)dz (4.13)
aU2
Ot
S( V,_ +
+ Txz ,L2
44
4.2.3 Equations for the Entire Water Column
Instead of solving the above equations for the vegetation layer and the vegetation-
free layer, it is more convenient to solve the vertically-integrated equations for the
entire water column, which can be readily derived by combining the equations for the
two layers. First, the vertically-integrated velocities over the entire water column in
the vegetation zone, U and V, can be defined as
(U, V) =- (U1 + U2, Vi + V2) = (Liii + L2ii2, L1i1 + L202) (4.14)
where uz u2, iU and U2 are vertically-averaged velocities within layer I and layer II,
respectively, while L1 and L2 are the thicknesses of layer I and layer II, respectively.
Adding the U1 equation and the U2 equation leads to
aU 9 (1U U\2 U Ux 2V2 .
+ + + L + gH(.
at Tx Li L2 dy L, L2
= fV + -(r, rb Fx)
P
r au] a a u
+ -AH- AHU] (4.15)
while the summation of the Vi equation and the V2 equation results in
av U V2 Vul
1
-9- + a(-n + --2 + 9y + -2 + gH(y
7t Tx L, L2 Ty L, L2 )
= -fU+-(r,, -Fy )
P
+ aAnV] + AHv] (4.16)
x a x Ty [A y
All the stress terms are computed as the quadratic power of the flow velocity. For
example, 7,, and r,y are computed as quadratic functions of the wind speed, rb, and
rby are computed as quadratic functions of U and V, and r, and riy are quadratic
functions of (U2 U1) and (V2 VI). The form drags associated with the vegetation
are:
S A2 ) (4.17
F, = pcdAN, (4.17)
Lx \LI/
FU, + v12 1V1
F, = pcdAN L L (4.18)
Ti, = PCdai -L + V- )2 (U (4.19)
GL2 L,1 L2 L, L2 TL
U2_ U2 + V2 V2 L 2 VI
7T = PCdi (4.20)
(4.21)
where Cd and Cdi are empirical drag coefficients.
An additional equation that must be satisfied is the continuity equation:
a( au ov
+ + 0 (4.22)
4.2.4 Dimensionless Equations in Curvilinear Grids
The above dimensional equations were presented to illustrate the development of
the vegetation model. In the curvilinear-grid model, however, dimensionless equations
in curvilinear grids are solved. These dimensionless equations are presented in the
following in terms of the contravariant velocity components in two layers:
+ L,[g"~ + gl2]
912 Ulr + U"22 + Ti F*
+ (Horizontal Diffusion), (Ui, U")
SRo [Nonlinear Terms(U, U")] (4.23)
aut
S+ L1[g21 +g22 1
-911 U 921 U+ 7r F*
+ (Horizontal Diffusion); (Uf, U")
+Ro [Nonlinear Terms(Uf, Ul)] (4.24)
+ i L-'"
Ot
+ L2[g1 g'12C,
912 U2 + U9 4+ '*
+ (Horizontal Diffusion), (U2, U2)
+Ro [Nonlinear Terms(Uf, U")] (4.25)
S+ L2 [21( +g22(]
= 9 U2 921 U" + r* '*
+ (Horizontal Diffusion), (Uf, U;7)
+R [Nonlinear Terms (U2, U) (4.26)
where
7,, = :r*. + (4.27)
7,*, = q71*. + T,7,*, (4.28)
7 = ri I[pfU (4.29)
ri4 = pfU I + ^ pfUrZ
7* + ,nri (4.30)
L" pfU~Z, pfU,. Z,
F: = ( + p ,Z (4.31)
P U'. Z, p f Uz, Z
F = ,~ ,z]+ 7Y yz (4.32)
Defining
U4 = U + U$ (4.33)
U'" = Ux + U (4.34)
we can obtain the following equations for the entire water column:
-t-+ + g[g11 +g2 + U T+ F
+ (Diffusion), [U14, U']
+ (Diffusion), [U, U2]
Ro
+ [Nonlinear Terms (Uf, Ur)]
Ro
+ -2[Nonlinear Terms (U2, U2")] (4.35)
aU7+ H [gl2 2 + q(7 2911 UU ..
a t [21 +g22] = + /:r,,7 b,7 F;,7
+ (Diffusion),[Uf, U1]
+ (Diffusion),[U, UV1)
Ro
+ ([Nonlinear Terms (Uf, U7)]
Ro
+ -2[Nonlinear Terms (Ut, U7)] (4.36)
The dimensionless continuity equation in the curvilinear grids becomes:
Ct + (-U) + ( U1) = 0 (4.37)
CHAPTER 5
HEAT FLUX MODEL
5.1 Introduction
All water bodies exchange heat with the atmosphere at air-sea interface. To
estimate the net heat flux at the air-sea interface, it is necessary to consider seven
processes: short-wave solar radiation, long-wave atmosphere radiation, reflection of
solar radiation, reflection of atmospheric radiation, back radiation, evaporation, and
conduction.
One way to estimate the net heat flux is to combine the seven complicated pro-
cesses to an equation in terms of an equilibrium temperature and a heat exchange
coefficient (Edinger and Geyer, 1967) as a boundary condition for the temperature
equation as follows:
0T
PoK, z = q, = K(T Te) at z = 9 (5.1)
where T,, the equilibrium temperature, is defined as the water surface temperature
at which there is no net heat exchange. This method will be called the "Equilib-
rium Temperature Method." Method to estimate the heat flux includes the "Inverse
Method" developed by Gaspar et al, 1990.
In the following, two methods for defining the boundary condition are described.
5.2 The "Equilibrium Temperature" Method
To determine an equilibrium temperature, following the procedure first devel-
oped by Edinger and Geyer (1967), meteorological data and empirical formulas are
required. Data from the South Florida Water Management District (SFWMD) and
"Climatological Data" published by the National Oceanic Atmospheric Administra-
tion (NOAA) were used.
5.2.1 Short-Wave Solar Radiation
The amount of short-wave radiation reaching the earth's atmosphere varies with
latitude on the earth, time of day, and season of the year. However, the amount
of short-wave radiation is reduced as it pass through the atmosphere. Cloud cover,
the sun's altitude, and the atmospheric transmission coefficient affect the amount.
Empirical formulas are used (Huber and Perez, 1970) to compute the amount of
short-wave radiation.
However, the short-wave solar radiation is more easily measured than computed
(Edinger and Geyer, 1967). SFWMD measured the solar radiation at the Station
L006, which is located at the south of Lake Okeechobee (See figure 5.1). The unit is
mLy/min.
5.2.2 Long-Wave Solar Radiation
The magnitude of the long-wave radiation may be estimated by use of empirical
formula. Brunt formula (Brunt, 1932) was used.
H, = a(T,- + 460)4(C + 0.031Ve-) [BTUFt-2Day-1] (5.2)
Ha = Long-wave atmospheric radiation,
Ta = Air-temperature in OF measured about six feet
above the water surface,
where e = Air vapor pressure in mmHg measured about six feet
above the water surface, and
C = A coefficient determined from the air temperature and
ratio of the measured solar radiation to the clear sky
solar radiation.
5.2.3 Reflected Solar and Atmospheric Radiation
The fractions of the solar and atmospheric radiant energy reflected from a water
surface are calculated by using the following reflectivity coefficients :
Hgs
R,, = H7 (5.3)
H,
Har
Rar = Ha (5.4)
Ha
50
Solar Radiation
1500
41000
C
-'J
S500
0 2 4.5
100 l Relative Humidity
90
80
70
60
50
32.5 Air Temperature
30.0
27.5 W
S25.0
22.5
Julian Day
Figure 5.1: Meteorological data at Station L006.
51
The solar reflectivity, R,,, is a function of the sun's altitude and the type and amount
of cloud cover. Because there was no cloud data, 0.1 was assumed for R,,. The
atmospheric reflectivity, Rar, 0.03, was assumed following the study of Lake Hefner.
5.2.4 Back Radiation
Water sends energy back to the atmosphere in the form of long-wave radiation.
This can be calculated by the Stephan-Boltzman fourth-power radiation law (Edinger
and Geyer, 1967 and Harleman, et al., 1973):
Hbr = y7,(T, + 460)4 (5.5)
Hbr = Rate of back radiation in BTUFt-2Day-',
where Yw = Emissivity of water, 0.97,
a = Stephan-Boltzman constant (4.15 x 10-OBTUFt-2Day-1 R-4), and
T, = Water temperature, OF.
5.2.5 Evaporation
Heat is lost from a body of water to the atmosphere through evaporation of
the water. Frequently, evaporation is related to meteorological variables (Brutsaert,
1982). The most general form is
He = (a + bW)(e, eC) (5.6)
a, b = coefficients depending on the evaporation formula employed,
W = Wind speed in miles per hour,
where ea = air vapor pressure in mmHg, and e, and
e, = Saturation vapor pressure of water determined
from the water surface temperature, T,.
e, are related to the temperature of air and relative humidity ( Lowe, 1977). NOAA
data include daily-averaged values of evaporation, wind speed, and air and water tem-
perature. These data are measured at the Belle Glade Station near the southeastern
shore of Lake Okeechobee. Those data are averaged daily. By using the least square
method, a and b were determined to be 5.663 and 296.36, respectively.
5.2.6 Conduction
Water bodies can gain or lose heat through conduction due to temperature differ-
ence between air and water. Heat conduction is related to evaporation by the Bowen
Ratio (Bowen, 1926).
B= H! (5.7)
He
B C(T, Ta) P
B = CT T) P (5.8)
(e, ea) 760
where
P = barometric pressure in mmH,,
C' = a coefficient determined from experiments = 0.26.
Thus, conduction is related to the other parameters as follows:
H, = 0.26(a + bW)(T, Ta) [BtuFt-2day-1] (5.9)
5.2.7 Equilibrium Temperature
Of the seven processes, four processes are independent of surface water temper-
ature: short-wave solar radiation, long-wave atmospheric radiation, reflected solar
radiation, and reflected atmospheric radiation. The sum of these four fluxes is called
absorbed radiation (HR). Thus the net heat flux can be written as follows:
AH = HR Hb, H, H, (5.10)
When the net flux AH is zero, HR becomes
HR = 7y(Te + 460)4 + 0.26(a + bW)(Te Ta) + (a + bW)(e, e,) (5.11)
The net heat flux can be expressed as follows:
AH = -K(T, Te) (5.12)
where K is the heat exchange coefficient.
5.2.8 Linear Assumption
Vapor pressure difference, e, ea, is assumed to have a linear relationship with
temperature increment as
e, e, = P(T, Te) (5.13)
53
Also, the fourth-power radiation term can be approximated by a linear term with less
than 15% error (Edinger and Geyer, 1967). Therefore, AH become as follows:
AH = -15.7 + (0.26 + #)(a + bW)(T, T,) + 0.051(T T,2) (5.14)
= -K(T, T)
Neglecting the quadratic term,
K = 15.7 + (0.26 + /)(a + bW) (5.15)
Using the above relation, an equation for Te can be derived as follows:
0.051T,2 H 1801 K 15.7 ea c(3) 0.26Ta
Te + + [06 + a] (5.16)
K K K 0.26 + 0 0.26 + P
where c(/) is intercept for the temperature and vapor pressure approximation.
5.2.9 Procedure for an Estimation of K and T,
Step 1. Compute HR.
Step 2. Assume T,.
Step 3. Find K for given wind and temperature.
Step 4. Compute the right hand side of Eq. 5.16.
Step 5. Compute the left hand side of Eq. 5.16.
Step 6. Compute the difference of Step 4 and Step 5.
Step 7. If error is not within error limit, go to step 2.
54
Step 8. If error is within error limit (0.5 C), K and Te are correct estimated
values.
An actual equilibrium temperature file was created using SFWMD data, which
were measured at 15-minute intervals. Wind speeds at Station A,B,C,D,E were used
for the computation.
5.2.10 Modification of the Equilibrium Temperature Method
By using the equilibrium temperature method, model-predicted temperature in
Lake Okeechobee was found to be unrealistic. Therefore, b of evaporation formula was
multiplied by a factor of 0.1. Further, K, the heat exchange coefficient, at Station C
was multiplied by 10 to give stronger stratification.
There are many uncertain empirical coefficients in the computation of an equilib-
rium temperature. First, evaporation data are averaged daily, but model time step
is 5 minutes. That means the estimation of evaporation data at a short interval is
difficult. Second, wind speed is also averaged daily, and evaporation is correlated with
this average wind speed. In actual computation, wind speed at 15-minute intervals
was used. Surface water temperature data are uncertain. Considering the sharp gra-
dient of water temperature that usually exists near the water surface, the error can be
large. Third, all the meteorological data used are from L006 station. Considering the
spatial variation of meteorological condition over the lake, error can be large. Most
other thermal models simulate the long-term variation of temperature with a time
step of one day. Therefore, it seems that a short-term variation of meteorological
data did not create serious problems.
5.3 The "Inverse" Method
When there are insufficient meteorological data, the errors in the estimation of
total heat flux at the air-sea interface can be large. To better estimate the total flux,
the so-called inverse method (Gaspar et al., 1990) was used in this study.
Total flux (qt) can be divided into two parts: solar (q,o1ar) and nonsolar (qnosolar).
55
While incoming solar radiation data are usually available, the nonsolar part is esti-
mated by solving the vertical one-dimensional temperature equation coupled with the
momentum equation.
5.3.1 Governing Equations
T 0 T
T a OuT
(K ) (5.17)
at az Oz
av a av
+ fu = a(Av) (5.19)
8t Tz 8z
5.3.2 Boundary Conditions
At the free surface
aT qt
K, = (5.20)
9z p
au 7,
A z (5.21)
az p
A, p = (5.22)
where K, is eddy diffusivity, A, is eddy viscosity, qt is the total heat flux, and 7, and
7y are wind stresses.
At the bottom
aT
S= 0 (5.23)
az
rTb = pcd Iu + V1ui (5.24)
Tby = pcd J + V1u v (5.25)
56
where Cd is the drag coefficient and ul and vl are velocities at the lowest grid point
above bottom.
Total flux qt cannot be specified a priori because of unknown nonsolar flux,
qnonsolar- Therefore, a value of qnonsoiar is first guessed and then corrected until the
calculated and measured water temperatures are within an error limit. In this way,
the total flux can be determined by summing up the solar and nonsolar parts. This
total flux was later used as a boundary condition of temperature equation in the
three-dimensional simulation.
5.3.3 Finite-Difference Equation
Treating the vertical diffusion term implicitly, Eq. 5.18 is written in the finite-
difference form as follows:
At the interior points,
(n+l + ,+1 n1n+l
U"n+ U" A j '+' 1 A )-
at aZ
At f.- =X Az (5.26)
where At is time increasement and Az is vertical spacing.
Applying free surface boundary condition,
U+l un+l
Un+1 -n U Avim-1 .zm
IMt fvi, m Az (5.27)
At IM Az
Where im is the index of surface layer.
Applying bottom boundary condition,
<1- ... ^^-Qyn+l n+1
Un+1 U1 A1 (n + Cd = n+1
At fv = Az (5.28)
Similar form for 5.19 is as follows:
At the interior points,
v1 A (VI+IV (- _+1)
'+1 + fu! = -' A Av (5.29)
At Az
Applying free surface boundary condition,
1,,n+1 _n+l
n+l ~ Avi m-1 IM im--l
-+1 -Vim--
Vim -- 2 + fu!n = p v (5.30)
At t Az
Applying bottom boundary condition,
n+1- +n
V1 V + fUn =
At 1
Av1 -Cd + + V+l
A,, z C
(5.31)
5.17 can be written as follows:
At the interior points,
Tin+1 Tn
At
T+ ( -T') (T"+ -Ti+1 )
SKAz.1 Az
Applying free surface boundary condition,
S1Tn+lTn+1
Tn+l __ Tm q I viml ~'im--l
A- TqjM1 Az
At Az
Applying bottom boundary condition,
T,+I T K (T+ -Tn)
Az
At Az
5.3.4 Procedure for an Estimation of Total Heat Flux
Step 1. Solve for u for given wind stress.
Step 2. Solve for v for given wind stress.
Step 3. Guess nonsolar heat flux and solve for T.
Step 4. compare the computed temperature to measured temperature.
Step 5. If error is not within error limit (0.5 OC), go to step 3.
Step 6. If error is within error limit (0.5 C), go to step 1.
(5.32)
(5.33)
(5.34)
CHAPTER 6
FINITE-DIFFERENCE EQUATIONS
This chapter describes the finite-difference equations for the governing equations
in (, t7, a coordinate system.
6.1 Grid System
Earlier models used a non-staggered grid so that all the variables were calculated
at the same point (center point). This has a disadvantage. When centered difference
scheme is used, one-sided difference scheme near the boundary should be used to
maintain the same order of accuracy. Therefore, it is inconvenient. CH3D uses a
staggered grid as shown Figure 6.1. Surface elevation and temperature are computed
at the center of a cell, while the velocities are computed at the face of a cell. A
vertical grid is shown in Figure 6.1, and all the variables are computed at the middle
of the layer.
6.2 External Mode
The external mode equations consist of the equations for surface displacement
c and the vertically-integrated velocities U and V. Treating the wave propagation
terms in the finite-difference equations implicitly and factorizing the matrix equation,
the following equations for the i-sweep and r7-sweep are obtained:
i-Sweep Equations
+ Tia~~(v~S U*) = C a((1 T1)8~(vS U") (6.1)
-a (T V")
Tli f"6S( .) [I + At(TC, _T g _l
0 U, u
a V, v
A C w. T,
o Ao 04
0 0 0
a a
o o
a a
0 a0
0a a
0 0A 0
u
Aw
0 T, p
A
0 0
*o
00o
A
0 0
A
0 0
nnc}
z-+ h
+h h
Figure 6.1: Horizontal and vertical grid system.
0
0 o
a 0
A 0A
60
= Tv1 12,((" ) + vIU" (6.2)
(1 T1)vfo1y(8"n) (1 TI) V 7126,((")
V U" At(1 T2)C, At(1 T) + D"
77-Sweep Equations
(+f + Tla,, ( V+) = ( + Tia,6,(x/V")) (6.3)
Tv1 ,22S((n+1) + i+ t (T2Cr, + T (6.4)
= TVy712S (*) (1 Tl)V 7126(") (1 Ti)4V;7'22S(nC)
+ V V V [At( T2)C7 + At( )1 + D
Ti 71286(Cn+1) + vU'n+ = v-U* + Tiv~ 12(C") (6.5)
where T, T2, and T3 are all constants between 0 and 1, superscripts indicate the
time level, and C, and C,, are the bottom friction terms. For example, the wave
propagation terms are treated explicitly if Ti = 0, but implicitly if T1 = 1. D" and D'"
are explicit terms in the U and V equations excluding the Coriolis, bottom friction,
and wave propagation terms. Additional parameters appearing in the equations are
aat
gdA(
1Atg"
A( (6.6)
12 HAtg12
7A AC
,22 HAtg2
= A7
21 HAt2
A77
In the Lake Okeechobee application, the external mode is first solved over the
entire lake. For the open water region, the above i-sweep and 7-sweep equations
61
are used. For the vegetation zone, however, the -sweep and rl-sweep equations are
modified by the presence of vegetation and are derived from Equations (4.23) through
(4.26). However, it is only necessary to solve the finite-difference equations for the
integrated velocities in the entire water column, (U, V), and the velocities in one
layer, (UI, Vi) or (U2, V2). The velocities in the other layer can be readily obtained
by subtracting the one-layer velocities from the total velocities.
6.3 Internal Mode
The internal mode is defined by the equations for the deficit velocity fi and v,
(l, i) (u -, v -) The equations for ii and i are obtained by subtracting the
vertically-averaged equations from the three-dimensional equations for u and v:
a g12 22
Hii = Hi + -H + F3
at I- -o
E, a L) (6.7)
+ H2 (A)A, (r rb) F2 (6.7)
a g1 g21
-Hi = g HL --Hv + G3
at -/go .go
+ H2 (Ao a ) (r., b,) G2 (6.8)
where ii and i are the deficit velocities in the (7, 7) directions, F3 and G3 indicate all
the explicit terms in the u and v equations, respectively, while F2 and G2 indicate all
the explicit terms in the U and V equations, respectively.
Applying a two-time-level scheme to the above equations leads to the following
finite-difference equations:
(Hfi)"+1 = (Hii)" + At- Ha+n+l
V9_
22
+ At-/ -H"f" + At(F3- F2)
+ At(H E ) A (H"+1in+l\
(rs _bE)n+I (6.9)
11
(Hv)n+l = (Hfi)" At f- Hn"+
-At H+I'v1 + At(G3 G2)
+4At E- A "(H n+1p+1)
(Hin)2 [Tva v
(r, rb)n+1 (6.10)
For the open water zone, the above internal mode equations are solved after the
external mode solutions are obtained. For the vegetation zone, no internal mode
equations are solved. This is consistent with the assumption that, in the vegetation
zone, the velocities are fairly uniform within the vegetation layer and the vegetation-
free layer.
6.4 Temperature Scheme
This section describes the finite difference equation which is used for solving the
temperature equation. Equation 3.53 is written in finite difference form using the
forward scheme in time and the centered difference in the vertical diffusion term.
At Ro [_
H+Tk1 -- H"T n T (HuTVr) (6.11)
a a
+ -(HvT + (Hw)"
E,(At) 1 Dv+ (n+1 +
+ Hn+lT. Ak Aa+ k+1
v- (Tn+l rpn+l
.Aa Ik t I-lI
HEH 2At 2T 82T 22T nT
+ TH g + 2g2 + 22 (6.12)
Dividing the above equation by H"+1 and collecting all the unknown terms in the
left-hand side and the known terms in the right-hand side, and writing advection
terms and diffusion terms separately,
E, At (D, D,+ ,n-+1 (6.1
(H+ A ,k-1 A J+l (6.13)
(Hn+l)2P,u aak sTyv- 1 sa,y~
63
E, At D,, D,
(Hn+1)2P,,aEOk D+o A o jk
H "n Ti A t n + V g '
6. 1~P, , + A+ (HwT)n
= Hn+1 (HuT + (HvT + io (HwT)"
H" EH'At n11 a2 2T+ 22 2Tn
+ +l+ 7 + g22 (6.14)
6.4.1 Advection Terms
Many different schemes can be used for the advection terms. When there is a
sharp discontinuity, it is difficult to model the convection without numerical diffu-
sion. Leonard (1979) introduces the QUICKEST( Quadratic Upstream Interpolation
for Convective Kinematics with Estimated Stream Terms), which gives good results
without excessive numerical diffusion.
This QUICKEST scheme treats the advection terms in the i-direction as follows:
( (,IogHuT)+ (VJgHuT)_
0(x~HuT) = (6.15)
where the first term in the right hand side is the flux at the right face of the cell
and the second term is at the left face. These two terms are difference differently
depending on the direction of current as follows:
At the right face of a cell :
When u is positive :
(v'lHuT)+ = (vgH);+,jui+j,i,k[ (Ti,j,k +Ti+,j,k)
(1 (i+1, )j)(Ti+lj,k 2Ti,k + Ti-1,j,k)
1 Ui+1,j,kAt T+f
1- (Ti+1,j,k Tij,,k)] (6.16)
When u is negative :
(xv HuT)+ = ( .\H)i+,,Iui+,,k [1(Ti,j,k + Ti+l,j,k)
1 ui+Iklst
(1 ( ui k ) (Ti+2j,k 2Ti+l,,k + Tj,k)
6 g
1 ui+lkAtk
il (Ti+1,j,k T,,k)] (6.17)
At the left face of a cell :
When u is positive :
( HuT)_ = (vH)ii,j,k [(T-Ij,k + Tij,k)
S (1 (i )2)(Ti,j,k 2Ti-,,k + Ti-2,j,k)
1 ui,j,kAt
-2 A U (Tij,k -i)] (6.18)
When u is negative :
1
(. THuT)_ = (v/H);,jui,j,k[ i(Ti-,j,k + Tid,k)
1 Uiij,klt
(1 u' ))(Ti+l,j,k 2Ti,j,k + Ti-1,j,k)
1 Ui,j,kAt
2 A (Ti,j,k i-,j,k)] (6.19)
The QUICKEST method treats the t7 direction advection term as follows:
(vHvT) ( vHvT)_
( HvT) (6.20)
At the top face of a cell :
When v is positive :
1
(v-HvT)+ = (gH)ij+li,j+l,k[-(T ,j,k+ Ti,j+l,j,k)
2
(1 (Vi+kt )2)(i,+,k 2Ti,j,k + Ti,j-1,k)
6 A-
1 vij+,kAt
-2 A (Ti,j+,k Ti,k)] (6.21)
When v is negative :
(J.HvT)+
=- (TH)i,j+lvi,j+l,k[1(Ti,j,k + Ti,j+l,j,k)
S(1- (Vi'j+'^At)')(Ti,j+2,j,k 2Ti,j+l,k + Tij,k)
1 vi,i+1,kat
S(Ti,+1,k ,k)]
-2-- -^Jl A.J,
At the bottom face of a cell :
When v is positive :
(\yHvT)
= ( vH)i,ji,,k[(Ti,j-,k + Tij,k)
(1 ( kAt)2)(Tj,k 2Tij-l,k + T,-2,k)
1 vi,j,kAt
2 A (T,,,k Tj-l1,)]
When v is negative :
(1gHvT)_ =
H),k
(V-/H),,jvj,k[I (Tij-l,k + Tij,k)
1 v i,j,k At
2 ( ( Vi ) 2) l, 2T + T ,)
I vij,kAt(T..i j, j-l,k)]
2 A ,,:k
The QUICKEST method treats the r-direction advection term as follows:
.(9 (HwT)+ (HwT)_
HwT) =
9r Aa
At the top of a cell :
When w is positive :
(HwT)+ = Hi,~w,,k[(Ti,j,k+l + Tj,k)
(6.22)
(6.23)
(6.24)
(6.25)
66
1 ( I,, 2 ) j,k+ 2Ti,j,k + Ti,,k-1)
W, (Tik+ Ti,j,k)] (6.26)
When w is negative :
(HwT)+ H= Hiwi,j,k Ti,jk+ + Ti,,k)
S(1 ( ,kA )2)(Ti,k+2 2Ti,j,k+l + Ti,,k)
2 ik (Ti,k+l Ti,j,k)] (6.27)
At the bottom of a cell :
When w is positive :
(HwT). = Hi,Iw,ij,k- l[(Ti,j,k + Ti,j,k-)
(1 wi'i'k-lAt)2T.
(1-('i- ) )(Ti,j,k 2Ti,j,k-1 + Ti,j,k-2)
6
1 Wi,j,k-1L~t
2Wi t (Ti,jk Ti,j,k-1)] (6.28)
2 a
When w is negative :
(HwT)_ = Hi,jwi,j,k-1 [(Ti,j,k + Ti,j,k-1)
1 wi,i,k-lat
(1 )2)(Ti,k+l 2T,j,k + Ti,,k-l)
1 oi1jk-lAt (Ti,j,k Ti,I,k-1)] (6.29)
2 Au
6.4.2 Horizontal Diffusion Term
The centered difference scheme was used for the horizontal diffusion. For the
mixed derivative term, temperature and depth at the corners of a cell are obtained
by averaging the values using four neighboring points at the center points of a cell.
1182(HT) 11 (HT)i+l,j,k 2(HT)i,j,k + (HT)i-,j,k
9 12 gi,j (-2 (6.30)
67
12 02(HT) 1
912 = g ,j[(HT)"i+l/2,j+l/2,k (H)i+l/2,j-1/2,k
(HT)i-1/2j+1/2,k + (HT)i-1/2,_-1/2,k]/A Aq (6.31)
22 2(HT) _2 (HT)ij+,k 2(HT)i,j,k + (HT),j,k32)
9 2 ( 2 (6.32)
CHAPTER 7
MODEL ANALYTICAL TEST
The purpose of model analytical test is to examine a model's capability to re-
produce well-known physical phenomena for which the model is designed for, by
comparing model results with analytical solution.
7.1 Seiche Test
The CH3D model has been tested for wind-driven circulation in an idealized
enclosed lake which is 11 km long and 11 km wide with a uniform depth of 5 m. A
uniform rectangular grid of 1 km grid spacing was used. To perform the seiche test,
the initial surface elevation was given as ( = (o cos(2rx/) where C( is an amplitude
and I is a wave length. In the test, (o was set to 5 cm and was set to 10 km.
Since the lake is of homogeneous density and without bottom friction and diffu-
sion, seiche period can be calculated as T = 2 where i is the basin length and h
is the mean depth. For the test basin, the seiche period is 0.87 hours. The simulated
surface elevation in the test basin over a 12-hour period is shown in Figure 7.1. The
result shows that the surface elevation was not damped and the seiche period agrees
with analytical seiche period.
7.2 Steady State Test
When a uniform wind blows in the same direction over a rectangular lake with
same magnitude over a long period, the lake circulation eventually reaches steady
state. Neglecting advection, horizontal diffusion, and bottom friction, the setup equa-
tion can be obtained as follows:
g h(7.1)
8 x- =h
69
Sufac El atten at Notth End
6
4)
6!-- 2----
-2
Ti.. l.ew I Tim. I- W .J
-4
S-
6 "2
Sirf.e Eleti.n aIt West bn Surface Elevatlon at Cantlr Surfce Etevlation at Cot End
6. 6. 6
.-2-- -2 -2
44 -4
-6 -6-- -6
b 6 '2 12 I 6 12
Surfagc EIlevaten at Sauth End
6
aI
-2
-4
-8
0 6 12
T(ia IHeural
Figure 7.1: Model results of a seiche test.
70
where p is density of water, 77 is surface elevation, rr, is wind stress, and h is water
depth. Using the same rectangular grid in the seiche test, a uniform wind stress of
1 dyne/cm2 from east to west was imposed. After 48 hours, a steady state is reached.
As shown in Figure 7.2, the surface elevation has a setup in the western part and
setdown in the eastern part. The setup across the lake is 1.12cm, which is exactly
the same as given by the analytical solution Eqn. 7.1.
7.3 Effect of Vegetation
In order to investigate the ability of vegetation model to represent the effect of
vegetation, CH3D was applied to a rectangular lake with a constant depth of 1 m and
horizontal dimensions of 4 km by 9 km. At first, vegetation was not considered, and
a wind stress of 5 dyne/cm2 was imposed. Then, vegetation canopy with width of 1
cm and density of 500 stalks / m2 was added to the western half of the lake. After
that, vegetation density was increased to 5000 stalks / m2. Vegetation height was
assumed to be the same as water depth.
With a time step of 5 minutes, the model was run for 24 hours. Time history
of surface elevation in the northern end of the vegetation area was plotted in Figure
7.3. As expected, surface elevation rises slowly for the second case and reaches steady
state after 5 hours. With high-density vegetation, surface elevation rises at a slower
rate compared to the second case. When wind blows uniformly over long time, vege-
tation effect disappears and reaches steady state. Additional resistance term due to
vegetation becomes smaller because the currents also become smaller at the steady
state and wind stress, and pressure gradient and bottom friction are balanced.
7.4 Thermal Model Test
The purpose of this test is to demonstrate how the velocity can be changed with
the consideration of thermal stratification. Surface heat flux was idealized using the
sine function as follows:
Tp
pC,K, = K(T Te) (7.2)
71
L.
I I I
I I
J I i 1
CONTOUR FRI -1.0000 TO 1.0000 CTOUR INTERYV O 0.10000 PTI3.31= 0.61442
SI I "
Figure 7.2: Surface elevation contour when the lake is steady state with uniform wind
stress of -1 dyne/cm.
Time history of surface elevation
0 10 0
0 5 10 15 20
Hour
Figure 7.3: Effect of vegetation on surface elevation evolution in a wind-driven rectan-
gular lake. Solid line is without vegetation, broken line is with low vegetation density,
and dotted line is with high vegetation density.
73
2drt
Te = 26 + 10 sin( -) (7.3)
where p is density of water, Cp is specific heat of water, K is heat exchange coefficient,
Te is an equilibrium temperature in 0 C, and P is period of 24 hours. Also, wind stress
was idealized as shown in Figure 7.4. Time history of currents at all five layers are
shown in Figures 7.4 and 7.5. It is apparent that when thermal stratification is
considered, currents at the surface layer are much stronger during increasing wind
condition because initial momentum is confined to a thinner surface layer.
74
Wind Stress at center
0. 12. 24. 36.
48. 60. 72. 84. 96.
Velocity at center without temperature
0. 12. 2'4. 36. 48.
Hours
60. 72. 84.
Figure 7.4: Time history of wind stress and currents at the center of lake at all five
levels. Thermal stratification is not considered.
E
- 0.0
C
-1.5
I I
-10.
-20.
i i
' ''''
i\
Wind Stress at center
I i T 1
/I I
0. 12. 24. 36. 48.
72. 84. 96.
Velocity at center with temperature
I I S' I
I I]
K
0. 12. 24. 36. 48. 60. 72. 84. 96.
Hours
Figure 7.5: Time history of wind stress and currents at the center of lake at all five
levels. Thermal stratification is considered.
ClJ
E
C
-1 .5
10. -
L
IT. "
CHAPTER 8
MODEL APPLICATION TO LAKE OKEECHOBEE
8.1 Introduction
Before the description of the application of CH3D to simulate the wind-driven
circulation in Lake Okeechobee, it is worth investigating the characteristics of the
lake.
8.1.1 Geometry
Lake Okeechobee, located between latitudes 27012'N and 26040'N and longitudes
80037'W and 81008'W, is the largest freshwater lake in America, exclusive of the
Great Lakes. With an average depth of approximately 3m, and the deepest part less
than 5m deep, the basin is shaped like a very shallow saucer. The western part of
the lake contains a great deal of emergent and submerged vegetation. According to
satellite photos, marsh constitutes 24% of the lake surface area.
8.1.2 Temperature
Due to the location of the lake in sub-tropical latitude, the annual fluctuations
of water temperatures are relatively small. The mean lake temperature based on
SFWMD monitoring in the 1970s and 1980s ranges from 15C to 34C (Dickinson et
al., 1991).
8.2 Some Recent Hydrodynamic Data
During the fall of 1988 and the spring of 1989, field data were collected by the
Coastal and Oceanographic Engineering Department, University of Florida (Sheng et
al., 1991a). Details of the field experiments and field data are described by Sheng
et al. (1991a). For completeness, some 1989 field data are described briefly in this
section.
Six platforms were set up in Lake Okeechobee. Locations of the six platforms are
shown in Figure 8.1. Station A was located in the northern portion of the lake, east
of north lake shoal, during the 1988 deployment. This Station A was moved south of
the rocky reef area during the 1989 deployment. The other five platforms remained at
the same locations during the 1988 and the 1989 deployments. Station B was located
near the Indian Prairie Canal. Station C was located in the center of the lake. Station
D was located 1.5 km west of Port Mayaka Lock. Station E was located about 1 mile
from the boundary between the vegetation zone and the open water, i.e., the littoral
zone. Station F was located within 30 m from the littoral zone.
Station A was selected to quantify the flow system in the northern zone during
the 1988 field survey. Because there was no measurement in the southern zone, it was
moved south to quantify the flow during the 1989 field survey. Station B was selected
to quantify the flux in the northern littoral zone. Stations C and D were selected to
calibrate the model and quantify the flux in the mud zone. Stations E and F were
selected to help the computation of phosphorus flux between the vegetation and the
open water.
Measured data at these six platforms include wind, current, water temperature,
wave, and turbidity. In this study, wind data were used to compute the wind stress
field, which is an essential boundary condition for the simulation of the wind-driven
circulation. Current data were used to calibrate and validate the 3-D hydrodynamics
model. The installation dates and locations of the platforms during the two deploy-
ment periods are shown in Table 8.1.
8.2.1 Wind Data
Wind speed and direction data averaged over 15-minute intervals were collected
from five stations in Lake Okeechobee during the spring of 1989. Because Station F
was close to Station E, wind data were, not collected at Station F. The data collection
E Emergent Vegetation
] Submerged Vegetation
0 5 10 Miles
i I i
81000'
80045'
Figure 8.1: Map of Lake Okeechobee.
1 r
I I
Table 8.1: Installation dates and locations of platforms during 1988 and 1989.
TIME OF LOCATION DATE LATITUDE LONGITUDE DEPTH
YEAR (cm)
Site A 09-20-88 27 06.31 80 46.21 396.0
Site B 09-17-88 27 02.78 80 54.31 274.0
FALL Site C 09-21-88 26 54.10 80 47.36 518.0
Site D 09-21-88 26 58.47 80 40.34 457.0
Site E 09-18-88 26 52.81 80 55.96 274.0
Site F 09-19-88 26 51.90 80 57.09 183.0
Site A 05-16-89 26 45.67 80 47.83 183.0
Site B 05-18-89 27 02.78 80 54.31 152.0
SPRING Site C 05-20-89 26 54.10 80 47.36 366.0
Site D 05-20-89 26 58.47 80 40.34 335.0
Site E 05-19-89 26 52.81 80 55.96 152.0
Site F 05-18-89 26 52.03 80 56.91 91.0
started on Julian Day 136.708. However, the direction of the anemometer was not
properly oriented until Julian Day 141.5. The location and height of the anemometer
are shown in Table 8.2.
As described in Sheng et al. (1991a), the measured wind over Lake Okeechobee
often exhibited significant diurnal variations associated with the lake breeze. During
relatively calm periods, significant spatial variation is often found in the wind field.
Water motion in the lake is significantly influenced by the wind. Figure 8.2
shows the wind rose diagram at Station C between Julian days 147 and 161. The
number inside the triangle indicates the percentage of wind data in that direction.
For example, wind from east to west is 27%. Wind speed between 4-6 m/sec is
about 45%. The governing wind direction is from east to west due to the location
of Lake Okeechobee. The surface area of Lake Okeechobee is big enough to create
its own lake breeze. During the daytime, wind blows from lake to land because the
air over the land is warmer than that over the lake. Because the Florida peninsula
is located between the Atlantic Ocean and the Gulf of Mexico, sea breeze affects
the wind direction. As Pielke (1974) indicated, the typical summer wind direction
Table 8.2: Instrument mounting, spring deployment.
SITE ARM CURRENT ELEV AZIM TEMP ELEV OBS ELEV
___(cm) (deg) (cm) (cm)
A Pressure Sensor 55695, Elev. 96 cm
Wind Sensor 5200, Elev. above MWL 670 cm minus water depth
S1 80673 71 342 07 86 0076 86
B Pressure Sensor 55696, Elev. 86 cm
Wind Sensor 5202, Elev. above MWL 518 cm minus water depth
1 80674 25 204 04 43 0078 43
2 80675 114 181 02 132 0075 132
C Pressure Sensor 48228, Elev. 297 cm
Wind Sensor 5203, Elev. above MWL 884 cm minus water depth
1 80679 61 330 01 79 0079 79
2 80681 123 333 09 140 0077 140
3 80680 284 342 06 302 0084 302
D Pressure Sensor 55694, Elev. 254 cm
Wind Sensor 5200, Elev. above MWL 883 cm minus water depth
1 80672 79 270 08 97 0083 97
2 80677 241 255 03 259 0082 259
E Pressure Sensor 55699, Elev. 104 cm
Wind Sensor 5199, Elev. above MWL 518 cm minus water depth
1 80671 36 59 11 53 0081 53
2 80678 116 72 05 135 0080 135
F Pressure Sensor 55697, out of water
No Wind Sensor
S1 80676 25 305 10 43 0085 43
81
is from southeast or southwest in south Florida. Those wind systems determine the
dominant wind direction in Lake Okeechobee. For the synoptic study of spring 1989,
the mean wind speed is about 5.1 m/sec, but it can exceed 11 m/sec. The wind field
is characterized by temporal and spatial non-uniformities. However, during strong
wind periods, the wind tends to be more uniform.
8.2.2 Current Data
Current data were collected at 15-minute intervals at the six locations. The data
collection started on Julian Day 136.708. The number of instrument arms at each
station depended upon the water depth and how many vertical levels of data were
desired. At station C, which is located in the center of lake, three current meters were
installed to measure the vertical variation of currents. Two sensors were installed at
Stations B, D and E. One sensor was installed at Stations A and F. The location and
height of each of the current meters are shown in Table 8.2.
As discussed in Sheng et al. (1991a), current data showed significant diurnal
variations in direct response to the wind. During a period of significant change in
wind direction, which usually follows a peak wind period, seiche oscillation generally
leads to significant current speed over several seiche periods (multiples of 5 hours).
Magnitudes of currents are very small at all stations. At Station C, mean magni-
tudes of u and v at arm 3 between Julian day 147 and 156 are 2.09 and 1.52 cm/sec,
respectively. Maximum magnitudes of u and v at arm 3 are 11.7 and 6.8 cm/sec,
respectively. Considering the accuracy of current meter, 2-3 cm/sec, currents are
very small.
8.2.3 Temperature Data
The lake temperature data showed that a significant vertical temperature gradi-
ent can be developed during part of a day when wind is low and atmospheric heating
is high. However, over the relatively shallow littoral zone and transition zone, tem-
perature appeared to be well mixed vertically much of the time.
82
LEGEND
: /
A
WIND ROSE AT STATION C
Figure 8.2: Wind rose at Station C.
83
Lake Okeechobee does not have a strong thermal stratification. But during the
daytime, lake becomes stratified. Temperature field data show the water temperature
difference between the near top and near bottom at Station C can reach 4 C. However
this stratification disappears in a short time due to the strong wind. This stratification
can affect the eddy viscosity significantly.
8.2.4 Vegetation Data
The distribution of aquatic vegetation in Lake Okeechobee was determined by the
use of recent satellite imageries and ground truth data obtained during field surveys.
Satellite imageries from the SPOT satellite were received and processed to create a
vegetation map of Lake Okeechobee (J.R.Richardson, personal communication, 1991).
This map was overlayed onto a curvilinear grid created by us. The number of pixels
with the same color was counted to classify the vegetation type on each grid cell. The
vegetation data include the vegetation class, the height of vegetation, the number of
stalks per unit area, and the diameter of each class. A total of more than 25 kinds of
vegetation were identified. The range of vegetation height is between 0.5 m and 4 m.
The density ranges from 10 stalks per m2 to 2000 stalks per m2 of bottom area, while
the diameter of the stalks ranges from 0.25 cm to 15 cm. The most popular type
of vegetation was cattail. These vegetation data were provided by John Richardson
from the Department of Fisheries and Aquaculture, University of Florida.
8.3 Model Setup
The following describes the procedure for simulation of wind-driven circulation in
Lake Okeechobee during spring 1989.
8.3.1 Grid Generation
The first step is to select the grid and grid size. Rectangular Cartesian grids were
widely used in hydrodynamic models. These grids are easy to generate but have a
disadvantage. Because these grids have to represent the boundaries in a stair-stepped
fashion, they cannot represent the complex geometry accurately unless a large number
|
RSS
HIDE
