• TABLE OF CONTENTS
HIDE
 Front Cover
 Title Page
 Acknowledgement
 Table of Contents
 List of Figures
 List of Tables
 Abstract
 Introduction
 Literature review
 Governing equations
 Model analytical test
 Model Application to Lake...
 Appendix A: Simulated currents...
 Bibliography






Group Title: Technical report – University of Florida. Coastal and Oceanographic Engineering Program ; 103
Title: Wind-driven circulation in Lake Okeechobee, Florida
CITATION PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00075325/00001
 Material Information
Title: Wind-driven circulation in Lake Okeechobee, Florida the effects of thermal stratification and aquatic vegetation
Physical Description: xvi, 202 leaves : ill. ; 29 cm.
Language: English
Creator: Lee, Hye Keun, 1955-
Publication Date: 1993
 Subjects
Subject: Okeechobee Lake (Fla.)   ( lcsh )
Coastal and Oceanographic Engineering thesis Ph. D
Dissertations, Academic -- Coastal and Oceanographic Engineering -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
 Notes
Thesis: Thesis (Ph. D.)--University of Florida, 1993.
Bibliography: Includes bibliographical references (leaves 196-201).
Statement of Responsibility: by Hye Keun Lee.
General Note: Typescript.
General Note: Vita.
Funding: Technical report (University of Florida. Coastal and Oceanographic Engineering Dept.) ;
 Record Information
Bibliographic ID: UF00075325
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved, Board of Trustees of the University of Florida
Resource Identifier: aleph - 001934047
oclc - 30811771
notis - AKB0129

Table of Contents
    Front Cover
        Front Cover
    Title Page
        Title Page
    Acknowledgement
        Acknowledgement
    Table of Contents
        Table of Contents 1
        Table of Contents 2
        Table of Contents 3
    List of Figures
        List of Figures 1
        List of Figures 2
        List of Figures 3
        List of Figures 4
        List of Figures 5
        List of Figures 6
        List of Figures 7
    List of Tables
        List of Tables 1
        List of Tables 2
    Abstract
        Abstract 1
        Abstract 2
    Introduction
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
    Literature review
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
    Governing equations
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
    Model analytical test
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
    Model Application to Lake Okeechobee
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
        Page 173
        Page 174
        Page 175
    Appendix A: Simulated currents by inverse method
        Page 176
        Page 177
        Page 178
        Page 179
        Page 180
        Page 181
        Page 182
        Page 183
        Page 184
        Page 185
        Page 186
        Page 187
        Page 188
        Page 189
        Page 190
        Page 191
        Page 192
        Page 193
        Page 194
        Page 195
        Page 196
    Bibliography
        Page 197
        Page 198
        Page 199
        Page 200
        Page 201
        Page 202
Full Text



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




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs