THREEDIMENSIONAL CURVILINEARGRID MODELING OF BAROCLINIC
CIRCULATION AND MIXING IN A PARTIALLYMIXED ESTUARY
By
JEIKOOK CHOI
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OFTHE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1992
ACKNOWLEDGEMENTS
This dissertation could not have been completed without the advice, ideas, and
financial support of my committee chairman, Dr. Y. Peter Sheng throughout this
study. I would also like to express my thanks to my dissertation committee member,
Dr. Robert G. Dean, Dr. Max Sheppard, and Dr. Louis H. Motz, for their patience in
reviewing this dissertation. I would also like to thank Dr. Albert Y. Kuo of Virgina
Institute of Marine Science, for providing the field data obtained from the James
River. Initial funding for this research was provided by the state of Virginia through
Virginia Institute of Marine Science.
I would like to dedicate this dissertation to my parents who have been praying
for me to be faithful, kind, and honest all the time.
Last but not least, I would like to thank my wife, Mija for her encouragemnet
and support, and my son, Jungmin for his sincere prayer.
TABLE OF CONTENTS
ACKNOWLEDGEMENTS ............................ ii
LIST OF FIGURES .... ............................ iv
LIST OF TABLES ................................. v
ABSTRACT .................................... vi
CHAPTERS
1 INTRODUCTION .................... ........... 1
1.1 Background ....................... ......... 1
1.2 Review of Previous Works ........................ 4
1.3 Modeling Strategy ................... ......... 7
1.4 Scope of This Study .................... ....... 10
1.5 Objective of This Study ................... ...... 13
1.6 James River ............ ...... .............. 13
2 GOVERNING EQUATIONS AND SOLUTION TECHNIQUE ...... 16
2.1 Differential Equations .................... ....... 16
2.2 Governing Equations in Stretched Coordinates ............. 20
2.3 Dimensionless Equations in Stretched Coordinates ........... 24
2.4 Dimensionless Vertically Integrated Equations ............. 29
2.5 BoundaryFitted Coordinate Generation ................ 30
2.6 Transformation of Governing Equations Into BoundaryFitted Grids 32
2.6.1 Transformation Rules . . . . . 32
2.6.2 Dimensionless Equations in BoundaryFitted Grids . 33
2.7 Boundary Conditions and Initial Conditions . . .... 35
2.7.1 Vertically Stretched System . . . . .
2.7.2 Vertically Stretched and Laterally BoundaryFitted System
2.8 Solution Technique ..........................
2.9 Finite Difference Equation . . . . . .
2.9.1 External Mode in (x, y, a) (TwoTimeLevel Scheme) .
2.9.2 Internal Mode in (x, y, a) (TwoTimeLevel Scheme) ....
2.9.3 External Mode in (7, 7, a) (TwoTimeLevel Scheme) .
2.9.4 Internal Mode in (n, r, a) (TwoTimeLevel Scheme) .
2.9.5 Salinity Scheme (TwoTimeLevel Scheme) . . .
2.9.6 Advection Terms .......................
2.9.7 Horizontal Diffusion Terms . . . . .
2.9.8 Numerical Stability ......................
2.10 Turbulence Model for Mass and Momentum Exchange Coefficients
3 MODEL ANALYTICAL TEST: COMPARISON BETWEEN MODEL
SULTS AND ANALYTIC SOLUTIONS ...................
3.1 Forced Constant Surface Slope ....................
3.2 Seiche Test .... .. .. ... .. .. .... ... .. .. ...
3.3 Tidal Forcing .............................
3.4 W ind Forcing .............................
3.5 DensityDriven Currents .......................
4
103
3.6 Summ ary ................................
ANALYSIS OF OBSERVED HYDRODYNAMIC DATA FROM JAMES R
3
Water Level Data .....
Tidal Currents ......
Circulation Pattern ..
SlackWater Salinity Data
RE
IVER
104
107
107
108
.
.
.
.
.
.
...........
...........
. . . .
. . . .
4.5 Front Dynamics ................... ......... 108
5 APPLICATION OF THREEDIMENSIONAL MODEL TO JAMES RIVER 117
5.1 Procedure of Model Simulation ....................... 117
5.2 Simulation of Tidal Circulation ................. ... 123
5.3 Sensitivity Tests .............................. 133
5.3.1 Effect of Bottom Friction ..................... 134
5.3.2 Effect of Bathymetry ....................... 139
5.3.3 Effect of Horizontal Grid ..................... 139
5.3.4 Effect of Time Step ........................ 140
5.4 NeapSpring Stratification and Destratification Cycle . ... 140
5.5 Salinity Intrusion ................... .......... 151
5.6 Frontal Dynamics ................... .......... 158
5.7 Model Calibration and Validation . . . ..... 229
5.8 Comparison of a and zbaroclinic approaches . . .. 234
6 CONCLUSIONS .... .............. ........... 252
BIBLIOGRAPHY ................................. .256
APPENDICES
A HIGHERORDERTERMS (H.O.T.) IN aGRID SYSTEM ....... ..263
B TRANSFORMATION OF GOVERNING EQUATIONS .......... 264
B.1 Transformation Rules ................... ........ 264
B.2 Dimensionless Tensor Invariant Governing Equations . ... 265
B.3 Transformation of stress components . . . ... 266
C HORIZONTAL DIFFUSION TERMS IN BOUNDARYFITTED GRIDS 267
D MORE COMPACT MOMENTUM EQUATIONS .............. 269
E FINITE DIFFERENCE EQUATIONS IN (, 7r, a) (THREETIMELEVEL) 271
E.1 External Mode .................... ......... 271
E.2 Internal Mode .................... ........... 274
LIST OF FIGURES
1.1 The James River Estuary from Old Point Comfort, at the mouth,
to Jamestown Island (upper left corner). Seed oyster beds are
shown as irregularly shaped darkened areas (Byrne et al., 1987). 15
2.1 The righthanded Cartesian coordinate system. . ... 18
2.2 at density field calculated from Eckart formula. . ... 19
2.3 The vertically stretched grid (From Sheng and Butler, 1982). 22
2.4 Field transformation (From Thompson et al., 1977). . ... 31
2.5 Staggered numerical grids in lateral and vertical directions (From
Sheng, 1983) ... .. .. ... .. .. .. .. ... .. .. .. .. 41
2.6 Comparison of stability functions for (a) some epirical formulae
versus (b) the simplified secondorder closure model (From Sheng,
1983) . . . . . . . . 67
3.1 A nonuniform rectangular grid system for testing the forced sur
face slope test .................. ........ 73
3.2 Simulated steadystate verticallyaveraged velocity field in a rect
angular basin forced by a constant surface slope. Horizontal grid
is shown in Figure 3.1.............. .......... 74
3.3 Same as Figure 3.2, except that the basin is rotated by 300 in the
counterclockwise direction. . . . .. 75
3.4 An irregular grid system for rectangular basin with irregular grids 76
3.5 Same as Figure 3.2, except that the irregular curvilinear grid as
shown in Figure 3.4 is used ................... .. 77
3.6 Comparison between simulated surface elevation (diamond sym
bol) and analytic solution (solid lines) for seiche oscillation in a
closed basin (T is seiche period) . . . .. 79
3.7 Comparison between simulated verticallyaveraged velocities (di
amond symbol) and analytic solution (solid line) for seiche oscil
lation in a closed basin (T is seiche period) . . .... 80
3.8 Maximum velocities for seiche case. . . . ... 81
3.9 Comparison between simulated surface elevations and velocities
(*'s) and analytic solutions for a tidally forced flatbottom basin
Elevation and velocity are in phase) . . . .... 85
3.10 Comparison of simulated maximum surface elevations at the closed
end and maximum velocities at the open end with the analytic
solution for a tidally forced sloping bottom basin .. . .... 86
3.11 Grid system for the annular section . . . ... 87
3.12 Comparison between simulated surface elevations and velocities
(*'s) and analytic solutions for a tidally forced flatbottom annular
section... ............................ 88
3.13 Maximum ebb and flood velocities in a tidally forced flatbottom
annular section at the end of tidal 10 cycles run. . ... 89
3.14 Comparison between simulated maximum surface elevation and
analytic solution for a sloping bottom basin forced by wind (ao is
the setup in the flat bottom).................. .. 93
3.15 Comparison between simulated vertical velocity structure and an
alytic solution for a constant wind in a closed basin. ...... ..94
3.16 Grid system for wind setup test in a bowtieshaped basin. .... 95
3.17 Steady state setup in cm for (a) 1 dyne/cm2 eastward wind, and
(b) 1 dyne/cm2 northward wind. . . . ... 96
3.18 Comparison between simulated results and analytic solution for
the vertical structure of horizontal velocity induced by density
gradient. . . .. . . . . .. 100
4.1 Lower James River and Hampton Roads (After Hepworth and Kuo
,1989) . . . . . . . .. 105
4.2 Observed Eulerian residual velocity over transect 'JB' (After Ham
rick et al. ,1989) ............................ 106
4.3 Surface elevation at the open boundary between June 19 and July
17, 1985. . . . . . . .. 109
4.4 Observed longitudinalvertical salinity at slackwater before flood
on June 19, 1985. .......................... 110
4.5 Observed longitudinalvertical salinity at slackwater before flood
on July 3, 1985. .. ............... ........ 111
4.6 Observed longitudinalvertical salinity at slackwater before flood
on July 9, 1985. .................. ........ 112
4.7 Observed longitudinalvertical salinity at slackwater before flood
on July 17, 1985 ................. ......... 113
4.8 Configuration of fronts at 1118 hours. From aerial photographs
of August 29, 1985. Generalized flow direction shown by arrows.
Dotted lines are depth contours in ft. (After Byrne, 1987) . 114
4.9 Configuration of fronts at 1400 hours. From aerial photographs
of August 29, 1985. Generalized flow direction shown by arrows.
Dotted lines are depth contours in ft. (After Byrne, 1987) . 115
5.1 Curvilinear grid for the James River coarse grid system (physical
dom ain) .. .. .. .. ... .. .. .. .. .. .. .. ... 119
5.2 Curvilinear grid for the James River fine grid system (physical
dom ain) .. .. .. .. ... .. .. .. .. ... .. ... 120
5.3 Curvilinear grid for the James River coarse grid system (compu
tation domain) ......................... 121
5.4 Curvilinear grid for the James River fine grid system (computation
dom ain) ... .......................... 122
5.5 Eulerian residual current at surface in the James River between
June 19 and July 17, 1985. ................... .. 124
5.6 Eulerian residual current at middepth in the James River between
June 19 and July 17, 1985. ................... .. 125
5.7 Eulerian residual current near bottom in the James River between
June 19 and July 17, 1985. ..................... 126
5.8 Vertical velocity profile at spring near the James River Bridge.
The numbers on top of plot indicate the time evolution in hours. 127
5.9 Vertical velocity profile at neap near the James River Bridge. The
numbers on top of plot indicate the time evolution in hours. 128
5.10 Computed nearsurface tidal currents near the James River Bridge
between June 19 and July 17, 1985. . . . ... 129
5.11 Contours of 29day residual longitudinal velocity at a transect near
the James River Bridge (along i=65) . . . .. 130
5.12 Simulated longitudinalvertical salinity at slackwater before flood
on July 3, 1985. ........ .......... ..... ... 142
5.13 Simulated longitudinalvertical salinity at slackwater before flood
on July 9, 1985. ........................... 143
5.14 Computed elevation, and velocity and salinity at surface, mid, and
bottom level for 29day simulation at Station 1 . .... 145
5.15
5.16
5.17
5.18
5.19
5.20
5.21
5.22
5.23
5.24
5.25
5.26
Instantaneous nearsurface horizontal
River at 1
horizontal
River at 2
horizontal
River at 3
horizontal
River at 4
horizontal
River at 5
horizontal
River at 6
horizontal
River at 7
distribution in lower James
5.27 Instantaneous nearsurface
distribution in lower James
5.28 Instantaneous nearsurface
distribution in lower James
Instantaneous nearsurface
distribution in lower James
5.30 Instantaneous nearsurface
distribution in lower James
5.31 Instantaneous nearsurface
distribution in lower James
5.32 Instantaneous nearsurface
distribution in lower James
velocity and salinity (ppt)
a.m. 1 on July 3, 1985. .
velocity and salinity (ppt)
a.m. on July 3, 1985. .
velocity and salinity (ppt)
a.m. on July 3, 1985. .
velocity and salinity (ppt)
a.m. on July 3, 1985. .
velocity and salinity (ppt)
a.m. on July 3, 1985. .
velocity and salinity (ppt)
a.m. on July 3, 1985. .
velocity and salinity (ppt)
a.m. on July 3, 1985. .
Computed elevation, and velocity and salinity at surface, mid, and
bottom level for 29day simulation at Station 2 . . .
Computed elevation, and velocity and salinity at surface, mid, and
bottom level for 29day simulation at Station 3 . . .
Model predicted longitudinal surface salinity distribution with 3000
ft3/sec of river discharge on July 3, 1985. . . . .
Model predicted longitudinal surface salinity distribution with 3000
ft3/sec of river discharge on July 9, 1985. . . . .
Model predicted longitudinal surface salinity distribution with 3000
ft3/sec of river discharge on July 17, 1985. . . .
Model predicted longitudinal surface salinity distribution with 1000
ft3/sec of river discharge on July 3, 1985. . . . .
Model predicted longitudinal surface salinity distribution with 1000
ft3/sec of river discharge on July 9, 1985. . . . .
Model predicted longitudinal surface salinity distribution with 1000
ft /sec of river discharge on July 17, 1985. . . .
Model predicted longitudinal surface salinity distribution with 9000
ft3/sec of river discharge on July 3, 1985. . . . .
Model predicted longitudinal surface salinity distribution with 9000
ft3/sec of river discharge on July 9, 1985. . . . .
Model predicted longitudinal surface salinity distribution with 9000
ft3/sec of river discharge on July 17, 1985. . . .
146
147
148
149
150
152
153
154
155
156
157
163
164
165
166
167
168
169
5.29
1
5.33 Instantaneous nearsurface horizontal velocity and salinity (ppt)
distribution in lower James River at 8 a.m. on July 3, 1985. 170
5.34 Instantaneous nearsurface horizontal velocity and salinity (ppt)
distribution in lower James River at 9 a.m. on July 3, 1985. 171
5.35 Instantaneous nearsurface horizontal velocity and salinity (ppt)
distribution in lower James River at 10 a.m. on July 3, 1985. 172
5.36 Instantaneous nearsurface horizontal velocity and salinity (ppt)
distribution in lower James River at 11 a.m. on July 3, 1985. 173
5.37 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 1 a.m. on July 3, 1985. 174
5.38 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 2 a.m. on July 3, 1985. 175
5.39 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 3 a.m. on July 3, 1985. 176
5.40 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 4 a.m. on July 3, 1985. .. 177
5.41 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 5 a.m. on July 3, 1985. 178
5.42 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 6 a.m. on July 3, 1985. 179
5.43 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 7 a.m. on July 3, 1985. 180
5.44 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 8 a.m. on July 3, 1985. 181
5.45 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 9 a.m. on July 3, 1985. 182
5.46 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 10 a.m. on July 3, 1985. 183
5.47 Instantaneous nearbottom horizontal velocity and salinity (ppt)
distribution in lower James River at 11 a.m. on July 3, 1985. 184
5.48 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 1 a.m. on July 3, 1985. . . .... 185
5.49 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 2 a.m. on July 3, 1985. . . .... 186
5.50 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 3 a.m. on July 3, 1985. . . .... 187
5.51 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 4 a.m. on July 3, 1985 . . .. 188
5.52 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 5 a.m. on July 3, 1985. . . ... 189
5.53 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 6 a.m. on July 3, 1985. . . ... 190
5.54 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 7 a.m. on July 3, 1985. . . ... 191
5.55 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 8 a.m. on July 3, 1985 . . .. 192
5.56 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 9 a.m. on July 3, 1985. . . ... 193
5.57 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 10 a.m. on July 3, 1985. . . ... 194
5.58 Instantaneous velocity in (x, z) domain along the channel in lower
James River at 11 a.m. on July 3, 1985. . . ... 195
5.59 Instantaneous velocity in (x, z) domain along the onegrid away
from the channel in lower James River at 1 a.m. on July 3, 1985. 196
5.60 Instantaneous velocity in (x, z) domain along the onegrid away
from the channel in lower James River at 2 a.m. on July 3, 1985. 197
5.61 Instantaneous velocity in (x, z) domain along the onegrid away
from the channel in lower James River at 3 a.m. on July 3, 1985. 198
5.62 Instantaneous velocity in (x, z) domain along the onegrid away
from the channel in lower James River at 4 a.m. on July 3, 1985. 199
5.63 Instantaneous velocity in (x, z) domain along the onegrid away
from the channel in lower James River at 5 a.m. on July 3, 1985. 200
5.64 Instantaneous velocity in (x, z) domain along the onegrid away
from the channel in lower James River at 6 a.m. on July 3, 1985. 201
5.65 Instantaneous velocity in (x,z) domain along the onegrid away
from the channel in lower James River at 7 a.m. on July 3, 1985. 202
5.66 Instantaneous velocity in (x, z) domain along the onegrid away
from the channel in lower James River at 8 a.m. on July 3, 1985. 203
5.67 Instantaneous velocity in (x, z) domain along the onegrid away
from the channel in lower James River at 9 a.m. on July 3, 1985. 204
5.68 Instantaneous velocity in (x, z) domain along the onegrid away
from the channel in lower James River at 10 a.m. on July 3, 1985. 205
5.69 Instantaneous velocity in (, z) domain along the onegrid away
from the channel in lower James River at 11 a.m. on July 3, 1985. 206
5.70 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 1 a.m. on July 3, 1985. . . ... 207
5.71 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 2 a.m. on July 3, 1985. . . ... 208
5.72 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 3 a.m. on July 3, 1985. . . ... 209
5.73 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 4 a.m. on July 3, 1985. . . ... 210
5.74 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 5 a.m. on July 3, 1985. . . ... 211
5.75 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 6 a.m. on July 3, 1985. . . ... 212
5.76 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 7 a.m. on July 3, 1985. . . .... 213
5.77 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 8 a.m. on July 3, 1985. . . ... 214
5.78 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 9 a.m. on July 3, 1985. . . ... 215
5.79 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 10 a.m. on July 3, 1985 . . ... 216
5.80 u at surface, w at middepth, salinity at surface, mid, and bottom
depth, and elevation along the onegrid away from the channel in
lower James River at 11 a.m. on July 3, 1985 . . .... 217
5.81 Salinity contours (ppt) in (x, z) domain along the onegrid away
from the channel in lower James River at 1 a.m. on July 3, 1985. 218
5.82 Salinity contours (ppt) in (x, z) domain along the onegrid away
from the channel in lower James River at 2 a.m. on July 3, 1985. 219
5.83 Salinity contours
from the channel
5.84 Salinity contours
from the channel
5.85 Salinity contours
from the channel
5.86 Salinity contours
from the channel
5.87 Salinity contours
from the channel
5.88 Salinity contours
from the channel
5.89 Salinity contours
from the channel
5.90 Salinity contours
from the channel
5.91 Salinity contours
from the channel
(ppt) in (x, z) domain along the onegrid away
in lower James River at 3 a.m. on July 3, 1985.
220
(ppt) in (x, z) domain along the onegrid away
in lower James River at 4 a.m. on July 3, 1985. 221
(ppt) in (x,z) domain along the onegrid away
in lower James River at 5 a.m. on July 3, 1985. 222
(ppt) in (x, z) domain along the onegrid away
in lower James River at 6 a.m. on July 3, 1985. 223
(ppt) in (x,z) domain along the onegrid away
in lower James River at 7 a.m. on July 3, 1985. 224
(ppt) in (x,z) domain along the onegrid away
in lower James River at 8 a.m. on July 3, 1985. 225
(ppt) in (x,z) domain along the onegrid away
in lower James River at 9 a.m. on July 3, 1985. 226
(ppt) in (x,z) domain along the onegrid away
in lower James River at 10 a.m. on July 3, 1985. 227
(ppt) in (x, z) domain along the onegrid away
in lower James River at 11 a.m. on July 3, 1985. 228
5.92 Time series of uvelocity and salinity in zbaroclinic with diffusion
at grid point (65,16) ........................
5.93 Time series of uvelocity and salinity in cbaroclinic with diffusion
at grid point (65,16) ........................
5.94 Time series of uvelocity and salinity in abaroclinic without dif
fusion at grid point (65,16)......................
5.95 Time series of uvelocity and salinity in zbaroclinic with diffusion
at grid point (65,17) ........................
5.96 Time series of uvelocity and salinity in rbaroclinic with diffusion
at grid point (65,17) ........................
5.97 Time series of uvelocity and salinity in abaroclinic without dif
fusion at grid point (65,17) . . . . . .
5.98 Time series of uvelocity and salinity in zbaroclinic with diffusion
at grid point (65,18) ........................
5.99 Time series of uvelocity and salinity in cbaroclinic with diffusion
at grid point (65,18) ........................
5.100 Time series of uvelocity and salinity in cbaroclinic without dif
fusion at grid point (65,18)......................
236
237
238
239
240
241
242
243
244
5.101 Time series of uvelocity and salinity in zbaroclinic with diffusion
at grid point (65,19). ........................ 245
5.102 Time series of uvelocity and salinity in abaroclinic with diffusion
at grid point (65,19). ..... .............. ..... .. 246
5.103 Time series of uvelocity and salinity in cbaroclinic without dif
fusion at grid point (65,19). . . . ..... 247
5.104 Time series of salinity in zbaroclinic with diffusion, abaroclinic
with diffusion, and cbaroclinic without diffusion at grid point
(65,16) . . . . . . . . 248
5.105 Time series of salinity in zbaroclinic with diffusion, rbaroclinic
with diffusion, and rbaroclinic without diffusion at grid point
(65,17).... ............... ............ 249
5.106 Time series of salinity in zbaroclinic with diffusion, ubaroclinic
with diffusion, and rbaroclinic without diffusion at grid point
(65,18)... ................ ............ 250
5.107 Time series of salinity in zbaroclinic with diffusion, rbaroclinic
with diffusion, and abaroclinic without diffusion at grid point
(65,19)..................... ........... .. 251
LIST OF TABLES
2.1 Characteristic Time Scales of Various Physical Processes in Estu
aries and Lakes. .......................... .. 27
3.1 Characteristics of model analytical tests. . . ... 102
4.1 Tidal amplitude harmonic constants obtained from June 19 to July
17, 1985. (period in hour, amplitude in cm, phase lag in degree) 104
5.1 Comparison of the model predicted harmonic constants of water
level with the measured data. . . . ..... 131
5.2 Comparison of the model predicted harmonic constants of velocity
with the measured data (St. 4, Depth=8.5 m) . . ... 132
5.3 Test cases performed. ........................ 135
5.4 Comparison between observed and simulated tidal elevation con
stituents at Station 1 (amplitude in cm, phase lag in degree). 136
5.5 Comparison between observed and simulated tidal elevation con
stituents at Station 2 (amplitude in cm, phase lag in degree). 137
5.6 Comparison between observed and simulated tidal elevation con
stituents at Station 3 (amplitude in cm, phase lag in degree). 138
5.7 Comparison of data and model tidal constituents of surface u
velocity at Station 4 (amplitude in cm, phase lag in degree, 6/19
7/17/85) . . . . . . .. 140
5.8 Comparison of data and model tidal constituents of surface v
velocity at Station 4 (amplitude in cm, phase lag in degree, 6/19
7/17/85) . .. . . . . .. 141
5.9 Model skill parameters for the sensitivity test to surface elevation
at Station 1 ................... ......... 232
5.10 Model skill parameters for the sensitivity test to surface elevation
at Station 2. ...... ....... ...... ......... 232
5.11 Model skill parameters for the sensitivity test to uvelocity at Sta
tion A . . . . . . . .. 233
5.12 Model skill parameters for the sensitivity test to vvelocity at Sta
tion A . . . . . . . .. 233
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
THREEDIMENSIONAL CURVILINEARGRID MODELING OF BAROCLINIC
CIRCULATION AND MIXING IN A PARTIALLYMIXED ESTUARY
By
JeiKook Choi
December 1992
Chairman: Dr. Y. Peter Sheng
Major Department: Coastal and Oceanographic Engineering
Threedimensional curvilineargrid hydrodynamic models for partiallymixed es
tuaries developed by Sheng (1989) and Sheng et al. (1989) were significantly en
hanced, applied to the James River, and tested to determine the sensitivity of model
results to numerous model parameters. The model computes the intertidal and
longterm dynamics of surface elevation and threedimensional field of velocity, and
salinity in partiallymixed estuaries with complex geometry and bathymetry.
The momentum and continuity equations in the curvilinear coordinates are solved
on a spacestaggered grid system using an implicit finitedifference algorithm for the
"external" verticallyintegrated flow and vertically implicit procedure for the "inter
nal" flow. Details of the numerical procedures are presented. The numerical accuracy
of the model is verified by comparing model results with analytic solutions for a vari
ety of idealized problems describing circulations in simple basins forced by constant
surface slope, seiche, tide, wind, and density gradient.
Model application to the James River reproduces the observed densitydriven
current. Effect of springtoneap tidal dynamics on vertical stratification and density
driven currents is simulated by the model. Longitudinal position of the head of
xvii
salinity intrusion was found to move up and down the estuary due to change in river
discharge and tidal condition. Increasing river discharge caused an increase in relative
stratification. Spring tide produced more mixing and enhanced salt intrusion while
neap tide produced more stratification. The observed frontal structure near Newport
News Point was successfully simulated as a narrow zone of surface flow convergence
and strong vertical transport. Simulation of the formation and evaluation of the front
show reasonable agreement with available field data. Sensitivity of the model results
to such parameters as turbulence model, advection scheme, horizontal grid, vertical
grid, and bottom boundary condition was investigated and the results were presented.
xviii
CHAPTER 1
INTRODUCTION
1.1 Background
Estuaries are the transition regions where the fresh water from the land mixes
with the saline water from the sea. Cameron and Pritchard (1963) have defined
an estuary as "a semienclosed coastal body of water which has a free connection
with the open sea and within which sea water is measurably diluted with fresh water
derived from land drainage" (p. 306). This definition implies two important dynamics
in an estuary: a horizontal density gradient between the fresh water and the sea
water, and a generation of kinetic energy caused by tidal currents. These dynamics
control the circulation pattern, the mixing processes, and the salinity and temperature
distributions that are the major physical problems that need to be understood in an
estuary. Ecologically, most estuaries are controlled by physical factors (Kinne, 1967)
such as tide, river discharge, wind, salinity, and morphology of the estuarine area.
Vertical mixing caused by the bottom and/or surface friction, tidal currents, and
the river discharge produces vertical and longitudinal gradients of salinity in the
mixing zone. Large river discharge and small mixing energy form sharp vertical
gradients in an estuary, which is referred to as "stratified". Saltwedge is the limiting
case of the stratified estuary. If the river discharge is small and the tidal mixing
energy large, slight vertical gradients of salinity will be developed. In the limiting
case, a "wellmixed" estuary is formed.
Partiallymixed estuaries are intermediate between saltwedge and wellmixed es
tuaries. James River is a partiallymixed, coastal plain type estuary like the majority
of estuaries along the east coast of the United States. These estuaries are charac
2
terized dynamically by vertical salinity variations of 110 /oo (Bowden, 1967) and
a net circulation pattern which moves seaward in the upper layers and landward in
the lower layers due to the pressure gradient caused by the density difference of the
salt and fresh waters. Rattray and Hansen (1962) showed that residual currents were
directly proportional to the longitudinal density gradient and inversely proportional
to the vertical eddy viscosity. Stratification exists in partiallymixed estuaries, and
this changes in various ways according to the relative amount of river discharge and
available tidal mixing energy. If the mixing energy is cyclical, like a springneap tidal
variation, alternate stratification and destratification may be evident (Haas, 1977;
Cerco, 1982). Currentinduced turbulent mixing plays an important role in affect
ing the stratification and destratification in a partiallymixed estuary. Vertical salt
flux occurs through vertical advection and entrainment of salinity across the interface
between fresh water and salt water. Tidal motion enhances the vertical turbulent ex
change of momentum and salt. A twolayered flow exists in a partiallymixed estuary
with the surface of no net motion near the maximum salinity gradient (halocline). In
addition to the longitudinal and vertical salinity gradients, a lateral salinity gradient
also occurs because of the Coriolis force and complex shorelines. Thus, motions in
partiallymixed estuaries are basically threedimensional and timedependent because
of irregular shape of the basin, densityinduced baroclinicity and the Coriolis force.
Salinity and temperature influence the density of estuarine waters, but salinity is
more effective than temperature because the range of salinity variation is much larger
than that of temperature variation. Salinity ranges from freshwater to saltwater in
the estuary. Wind also influences the vertical salinity and velocity structures while
interacting with internal friction and bottom friction. Surface currents are generally
in the same direction as the local wind stress even though they are affected by the
Coriolis force. Thus downestuary winds increase stratification while upestuary winds
decrease it. Nonlocal wind effects can affect stratification through change of sea level
3
at the mouth of the estuary (Wang and Elliott, 1978).
River or open sea boundary conditions vary with time, and the location of salinity
intrusion changes accordingly. If the head of salt intrusion is defined as the point
where freshwater meets saltwater, the head of salt intrusion is dependent on the
amount of river discharge, the salinity at the open sea boundary, tidal currents, wind
stresses, and the geometry of the basin. They are all timedependent except the
geometry of the basin. Accurate specification of the boundary conditions is essential
to accurate prediction of circulation and salinity distribution.
Turbulence models have been developed to simulate the mean and turbulent mo
tions in most stratified flow situations. A turbulence model is defined as a system
of equations which calculate Reynolds stress terms in the meanflow equations and
similar terms in temperature and salinity equations. They parameterize the effect of
turbulence on the mean flow. Through the Reynolds decomposition of the governing
equations and Reynolds averaging, Reynolds stress equations are developed. However,
the system of equations is not closed since every equation contains terms consisting
of higher order correlations which are unknown. The system can only be closed by
using a semiempirical (e.g. various levels of secondorder closure model) or adhoc
empirical models (e.g. simple eddyviscosity models). An equilibrium closure model
(Sheng et al., 1989) is included in the threedimensional model used for this study.
To reduce computational effort, the system of equations can be simplified with such
assumptions as local equilibrium of turbulence. Although Reynolds stress models can
provide rather accurate results (e.g. Sheng, 1984; Sheng and Villaret, 1989), Reynolds
stress models are quite complicated, and require substantial computer resources.
At present, hydrodynamic models of circulation are important tools used by both
engineers and scientists attempting to quantify and understand complex physical,
chemical and biological processes in estuaries and coastal areas. Historically, engineers
and scientists preferred physical modeling to mathematical modeling for studying cir
culation and transport in most of the complicated environmental systems because
of difficulties in finding mathematical solutions. However, physical modeling has in
herent scaling problems and is very costly. As computers become faster and more
accurate, mathematical models have replaced physical modeling as the major tool for
studying circulation and transport. Some differences exist between numerical results
and analytical solutions of mathematical problems because of the digital nature of
the computations. To reduce these differences and enhance efficiency of the compu
tations, numerous numerical schemes and solution techniques have been developed.
The numerical models developed must be checked against analytical solutions before
applying them. Agreement between numerical results and analytical solutions sug
gests accurate and consistent numerics of the numerical model. To test the model's
ability to simulate such physical processes as turbulent mixing and bottom friction, it
is necessary to calibrate and verify the numerical model with laboratory and/or field
data on circulation and transport.
1.2 Review of Previous Works
Modeling of circulation and salinity transport dynamics in the James River has
been active over the last decade. Circulation modeling has included the flow induced
by tide, wind, and density gradients, with the primary emphasis on the tidally driven
circulation. Cerco (1982) applied a twodimensional laterally averaged model to the
James River to predict intertidal variation of surface level, velocity, and salinity.
This model could simulate the destratificationstratification cycle coincident with the
springneap tidal cycle in the James River Estuary. Heltzel and Granat (1988) applied
a twodimensional verticallyaveraged model to the lower James River to study the
impact of artificial islands for dredge spoil disposal, channel deepening, and port
expansion. This model could simulate the convergence phenomena in the tidal flow,
but could not show the complex vertical structure of the flow field.
Most of the early work on the numerical modeling of estuaries used vertically
integrated equations (Hansen, 1956; Leendertse, 1967; Blumberg, 1975; Wang, 1975;
Johnson, 1980). This approach was found to be useful for studying tidal dynamics
and storm surges although it cannot simulate the internal features of the estuar
ine circulation. In other cases, laterally averaged equations were used to study the
vertical structure of currents and species (Spaulding, 1974; Blumberg, 1976; Cerco,
1982), although no lateral flow structure could be provided. The development of
threedimensional numerical models started in the late 1960's, but have progressed
significantly during the past decade because of the availability of high speed comput
ers. Details of a few threedimensional numerical estuarine models can be found in
Caponi (1976), Sheng (1983,1987), Lynch (1983), Blumberg and Mellor (1983), and
Oey et al. (1985).
The early threedimensional models were developed by dividing the vertical water
column into layers and solving each layer with a vertically averaged twodimensional
model (e.g., Laevastu, 1974). Laevastu (1974) introduced a twolayer model consist
ing of a stack of two vertically averaged models and a staggered finite difference grid
system. Wang and Connor (1975) also used the layer concept using a finite element
method. This layered approach has the disadvantage that a layer thickness can disap
pear in simulations with strong vertical velocities such as downwelling and upwelling.
To alleviate this weakness, many researchers have used the leveled approach. In a
leveled model, the level interfaces are fixed, and free exchanges of water mass, mo
mentum, salt and heat are allowed across the interfaces. The leveled approach is used
by Leendertse and Liu (1975), Caponi (1976), and Sheng (1983, 1987).
Sheng (1986c, 1986d, 1987) developed a threedimensional hydrodynamic model
which can be used in nonorthogonal boundaryfitted grid as well as orthogonal or
Cartesian grid. As in the earlier Cartesian grid model (Sheng, 1983), the curvilinear
grid model uses similar modelsplitting technique, vertical turbulence closure scheme,
modified alternating direction implicit (ADI) scheme for the external mode, and
verticallyimplicit numerical schemes. The governing equations in the curvilinear grid
system were derived by tensor transformation in terms of the contravariant velocity
components.
Swanson (1987) developed a threedimensional boundaryfitted hydrodynamic
model in terms of Cartesian velocity components. He developed a Helmholtz equation
for surface elevations, which is derived by eliminating the velocities from the continu
ity equation and the momentum equations, and solved it with the Gaussian elimina
tion or successive over relaxation (SOR) method. This method has an advantage in
handling boundary conditions since velocity components are known simultaneously
at the new time step compared with an ADI method. However, the advantage of an
ADI scheme is far greater in terms of its efficiency. The tridiagonal matrix associated
with the ADI method is very efficient to solve on any computer, while SOR needs
more computing time because of the repeated iterations required. Swanson used con
stant vertical eddy coefficients, which is incorrect (turbulence is timedependent and
spatiallyvarying) and is difficult to adjust in each new model.
The development of turbulence closure models has been of major significance be
cause of the need to accurately simulate the vertical fluxes of salt and momentum
and their dependence on stratification. A number of turbulence models have been de
veloped. For example, Donaldson (1973) and his colleagues at Aeronautical Research
Associates of Princeton developed a hierachy of turbulence models. Sheng (1982,
1984) modified the Reynolds stress model and Sheng and Chiu (1986) and Sheng
and Villaret (1989) modified the equilibrium closure and turbulent kinetic energy
(TKE) closure models for marine and freshwater applications. Mellor and Yamada
(1982) described a similar turbulence model which was applied by Blumberg and
Mellor (1983) and Oey et al. (1985). They classified the turbulence model as level
4 (Reynolds stress), level 3 (TKE closure), level 2 (equilibrium closure), and level
2.5 models according to the complexity of the turbulence equations. The turbulence
closure models have been applied to simulate complex mixing processes (Sheng, 1982,
1985; Blumberg and Mellor, 1983; Oey et al., 1985). In addition to the Reynolds stress
equations of turbulent motion, equations for turbulent kinetic energy and turbulent
macroscale had to be added to the system of equations to complete the turbulence
model. Although various levels of turbulence closure models have been used by vari
ous modelers, little has been done to compare the results obtained with various levels
of closure models versus realistic data from estuaries.
Conventional finitedifference models use only Cartesiangrid and have shortcom
ings for applications where physical boundaries are very complex and thus may require
excessively fine grid in order to resolve the boundary in the model. To allow a more
flexible arrangement and spatial variation of the grids, one can use curvilineargrid
finitedifference models which were discussed above or finiteelement models (Wang,
1975; Wang and White, 1976). The finiteelement models allow flexibility in the dis
cretization of the spatial domain. However, because the matrix of the discretized
system is irregular, direct inversion method which requires very extensive compu
tational effort, is usually used. In the finitedifference models, however, tridiagonal
matrix system generally exists and can be efficiently solved. The curvilineargrid
finitedifference models (e.g. Spaulding, 1984; Sheng, 1986c, 1986d, 1987). This al
low extensive user flexibility in spacing of the horizontal grid structure while still
maintaining a rather simple matrix system and allowing efficient solutions. However,
there are some disadvantages of this approach such as more computational efforts
due to the extra terms in the transformed equations and potential numerical dissipa
tion/dispersion problems caused by skewness and rapid change in grids.
1.3 Modeling Strategy
The examination and analysis of available data, such as river flow, water eleva
tions, wind, current, and salinity, should be performed to classify the dynamics of
the water body along with its range of variability. Models that do not adequately
resolve the depths of navigation channels of an estuarine system may have difficulty
in predicting the correct degree of stratification and the strength of tidally averaged,
nearbottom upstream flow. The Hansen and Rattray theory (1965) indicates the
importance of tidal dispersion in maintaining the salt balance in an estuary.
The tide in an estuary may have the characteristics of a progressive, standing, or
mixed wave. Flood and ebb tide are not only unequal in velocity at a fixed point in
most estuaries, but also do not always flow opposite to each other. Ludwick (1975)
attributed this inequality in tidal flow to three factors: (1) mixing between river
water and sea water produced by turbulent exchange and vertical entrainment across
a density interface, (2) the Coriolis force which deflects the ebb and flood dominant
flows along different perimeters of the estuary, and (3) evasion of ebb and floodtidal
pathway. Although all three factors may occur concurrently, the effect of the Coriolis
force is usually overshadowed by the other two factors.
Exchange of various forms of energy produces the mixing in a vertically stratified
system. For example, turbulent kinetic energy extracted from the mean flow in a
frictional boundary layer produces mixing when exchanged with potential energy.
Various parameters have been introduced to classify the degree of stratification. Ippen
(1966) introduced the stratification parameter defined as the ratio of rate of turbulent
energy dissipation divided by rate of gain of potential energy. Hansen and Rattray
(1965) developed the dimensionless parameter by which the degree of stratification
can be related to buoyancy and mixing energy.
Festa and Hansen (1976) discussed the effects of depth and river discharge on
estuarine circulation and salinity in their twodimensional model based on vorticity
and saltbalance equations. They found that decreasing the discharge allowed the
head of salt intrusion to move upstream. As the discharge is decreased, estuarine
circulation is also weakened, although it occurs over the larger area. Greater depths
enhanced circulation and moved the head of salt intrusion upstream.
9
Geyer and Farmer (1989) measured the time variation of salt intrusion in an estu
ary during a variety of tidal and runoff conditions using highresolution instruments
and sampling methods. The salt wedge was found to vary considerably in position
and vertical structure throughout the tidal cycle due to interaction of tidal flow and
densityinduced motion. During the flood cycle, salt wedge progresses up the estuary
as a gravity current, while during the ebb the salinity structure was eroded by shear
instability, in which vertical shears become so large that the pycnocline becomes un
stable and causes the salt wedge to break down. Here the most important dynamic
variable is the gradient Richardson number (Ri), which reflects the stability of the
shear layer and thus determines whether momentum and salt can be exchanged across
the pycnocline.
Tides in estuaries are modulated by riverforced estuarine plumes. A numerical
study of this phenomena was done by Chao (1990). In a typical midlatitude estuary
shelf environment, the riverforced circulation can be divided into three regions: a
seaward surface flow and a landward undercurrent inside the estuary, a strong anti
cyclonic surface gyre overlying a weak cyclon off the estuary mouth, and a coastal jet
flowing in the direction of Kelvin wave propagation farther downstream. The tidal
residual eddies enhance the plume growth but tend to trap the plume.
Frontal systems in estuaries are also very important in estuarine circulation and
transport. Limited observations by Byrne et al. (1987) suggested that oyster larvae
remain in the water column for a period of two to three weeks and are carried by water
movement associated with tidal front and tidal residual circulation. A mathematical
model that describes the formation and dilution of a frontally bounded river plume is
presented by O'Donnell (1990). The behavior of plumes in timedependent flows like
tidal current were examined in his numerical test. The interaction of the plume front
and the ambient flow was found to have major influence on the spreading and mixing
of the buoyant water. River plumes may only mix significantly with the underlying
10
ambient fluid at particular phase of the tide, but the mixing generally occurs over
a large area. Thus the boundary conditions applied in estuarine circulation models
may need to be carefully reconsidered.
In any case, the density distribution determines the stratification. Therefore, a
synoptic densityfield survey should be made. To remove tidal effects as much as
possible, it is desirable to sample at slack. During slack water, stratification may be
stronger than during the other stages, making water column sampling easier in low
currents. The usual way to sample at slack is that the boat is following slack water.
1.4 Scope of This Study
Modeling of a partially mixed estuary is one of the most challenging areas in com
putational coastal hydrodynamics. Partial stratification takes place in many coastal
plain estuaries like the James River (Figure 1.1), where freshwater inflows and tidal
mixing interact significantly. Partial stratification requires that vertical variation of
flow variables be adequately resolved such that the model can simulate the density
driven currents. Thus, in contrast to traditional barotropic tidal hydrodynamics, the
simulation of baroclinic circulation requires the accurate representation of vertical ve
locity and density structures, and the use of threedimensional models (Sheng, 1983;
Oey et al.). Sheng et al. (1989a) simulated the tidal, winddriven, and baroclinic cir
culation in Chesapeake Bay, while Sheng et al. (1989b) simulated the tidal circulation
and baroclinic circulation in the James River using a very coarse grid. This study
provide more detailed simulations of tidal and baroclinic circulation in the James
River.
Traditionally, Cartesian grids have been used for the estuarine hydrodynamic
model. For most estuaries, the geometry and bathymetry are quite complex and ir
regular so that the generation of an orthogonal grid is very difficult. Models developed
by Blumberg and Mellor (1987) require an orthogonal grid and cannot accept a non
orthogonal grid, and hence cannot accurately resolve the complex shorelines. On the
11
other hand, boundaryfitted grid hydrodynamic models appear to be an excellent tool
to allow the user to accurately resolve complex geometries for estuarine and coastal
applications. Few studies focused on model development and validation against an
alytic solutions. For example, Sheng and Hirsh (1984) developed a twodimensional
curvilineargrid model for coastal ocean circulation. Spaulding (1984) applied a two
dimensional model to simulate tides in the North Sea. Practical applications generally
require a fully threedimensional model analysis because of stratification, complex ge
ometry, and topography. Although the threedimensional boundaryfitted grid hydro
dynamic model developed by Sheng has been tested against some analytic solutions
(Sheng and Hirsh, 1984; Sheng, 1987), this study provides more detailed analytical
tests of the threedimensional model.
Despite the complexity of the Reynolds stress models, such models do contain
more physics and hence the model constants remain unchanged for all model ap
plications (Sheng, 1982). A Reynolds stress model was used to produce successful
simulations of the wave boundary layer (Sheng, 1982; Sheng, 1984; Sheng, 1986a;
Sheng and Villaret, 1989), wavecurrent boundary layer (Sheng, 1984), and vegeta
tion boundary layer (Sheng, 1982). However, threedimensional models generally do
not use Reynolds stress models to represent the turbulent mixing, because of the com
plexity. Instead, simplified secondorder closure models have been developed and used
with threedimensional models. These include the turbulent kinetic energy closure
model (Sheng and Villaret, 1989) and the equilibrium closure model (Sheng and Chiu,
1986; Sheng et al., 1989a). Both of these models contain the assumption that the
turbulence time scale is much less than the mean flow time scale. The equilibrium
closure model contains the additional assumption that turbulence does not change
significantly over the turbulent macroscale. The TKE closure model resolves the
turbulent dynamics by including the differential transport equation for q2, while the
other secondorder correlation equations are simplified to algebraic equations. The
1
12
TKE closure model is similar to the Level3 model of Mellor and Yamada (1982) and
the k e model of Rodi (1980), but uses an equation for the macroscale (A) instead of
the equation for turbulence dissipation (e) used by the other two models. The equilib
rium closure model is slightly more simplified than the TKE closure model, since only
algebraic equations are used to resolve the secondorder correlations. The equilibrium
closure model is similar to Mellor and Yamada's Level2 model. Instead of solving
the differential transport equations for A which consists completely of modeled terms,
simplified secondorder closure models often solve a set of integral constraints for A
(Sheng and Chiu, 1986).
Results of the simplified secondorder closure models have been quite successful.
Sheng and Villaret (1989) used the TKE closure model to simulate a wave boundary
layer, a thermally stratified boundary layer and a sedimentladen boundary layer.
On the other hand, the equilibrium closure model has also been found to be very
successful for simulating the vertical mixing in estuaries (Sheng et al., 1989b, 1989c;
Johnson et al., 1991), despite the further simplification. There exists no unequivocal
evidence to suggest that the TKE closure model is superior to the equilibrium closure
model. This subject will be investigated in this study by comparing results of TKE
closure model and equilibrium closure model versus field data from the James River.
Recent field studies in the James River have obtained quantitative or qualitative
information on the temporal and spatial variations of salinity stratification, frontal
structure, and gyre formation (Bryne, 1987; Hepworth and Kuo, 1989). These field
data provide a good opportunity for detailed calibration and validation of the three
dimensional curvilinear grid hydrodynamic model developed by Sheng (1989). Model
calibration and validation will be conducted in terms of detailed intertidal dynamics,
water level vertical flow and salinity structures, springtoneap and residual longterm
circulation and salinity transport, and formation and evolution of tidal fronts in the
Hampton Roads area. Sheng et al. (1989c) reported that the use of ogrid may lead
13
to numericallyinduced salinity transport from the navigation channel to the shallow
shoal, unless special care is taken to improve the evaluation of baroclinic terms and
the advection terms. This subject will also be further investigated in this study.
1.5 Objective of This Study
The objectives of this study are (1) to further enhance the threedimensional
hydrodynamic curvilineargrid model developed by Sheng (1989) and Sheng et al.
(1989) for simulating circulation and transport patterns in coastal areas and estu
arine systems, (2) to apply the model to the James River to simulate the observed
shortterm and longterm circulation features, (3) to determine the model sensitiv
ity to various model features including horizontal grid resolution, bottom friction,
turbulence closure scheme, and vertical grid.
1.6 James River
The Chesapeake Bay estuarine system was formed by the drowning of a cascading
series of river valleys, the most seaward of which was the James. Estuaries formed
in this way are called coastal plain estuaries. In this type of estuary, the primary
forcing functions for salinity and current distribution are tidal mixing which depends
on the tidal range and basin configuration, gravitational circulation due to freshwater
discharge, and Coriolis force due to the earth rotation. An estuary where gravita
tional circulation is dominant is called a salt wedge estuary. But the tidal oscillation
increases mixing between the salt and freshwater and reduces the effectiveness of the
gravitational circulation. James River is called a partially mixed estuary because
both the gravitational circulation and the tidal mixing play an important role in the
estuarine dynamics.
The James River (Figure 1.1) originates in the Appalachian mountains of Virginia
at the confluence of the Jackson and Cowpasture Rivers. The river runs approximately
530 km into Chesapeake Bay southeastward. About one third of it, the 160 km (from
the mouth at Sewell's Point to the fallline at Richmond) is a tidal zone, and about
14
one third of this tidal zone is the estuary proper. The central half of the estuarine
portion of the James River contains some of the most productive seed oyster beds in
the world. These beds cover an area of 25,000 acres, of which 16,150 acres have been
determined to be higher or moderately productive (Haven and Whitcomb, 1983).
The average tidal range varies from 76 cm at Sewell's Point to 58 cm near the
Chickahominy River mouth to 98 cm at Richmond. Mean annual freshwater dis
charge at Richmond is 215 m3/sec and in the summer months flow is more typically
in the range of 35 to 75 m3/sec. (Cerco, 1982). Salt generally intrudes upstream to
the Jamestown Island, but the intrusion may approach Hopewell or be forced down
stream as far as Newport News due to low flow and heavy rainfall. The lower James
River near Hampton Roads generally shows an extremely complex threedimensional
and timevarying density field due to complicated geometry, highly irregular bottom
topography, tidal flow, wind effect, and variable river discharge.
As shown in Figure 1.1, the James River is channelized as far upstream as Rich
mond with a minimum depth at midchannel of approximately 8 m. Maximum depth
in the estuary is 10 15 m and average depth is about 6 m. Crosssectional area
ranges from approximately 27,000 m2 near the James River Bridge to 1,000 m2 at
Richmond. In general, the basin is on the order of 7 km wide from the mouth at
Old Point Comfort to a constricted region at Mulberry Point, just upstream of the
broader (10 km wide) region of Burwell Bay. The remaining portion of the estuary,
from Mulberry Point to Jamestown Island is from 4 to 6 km wide.
15
ac
I. *. :
C,
or ca`,o
0 r
~\ ~L,
I i c,'L

Cd
/'I
QCd
.0.
.Y~ Qi:a
.~~0. I
ri ~ ~ *a
C I)
EEE
NCSJCJB(
CHAPTER 2
GOVERNING EQUATIONS AND SOLUTION TECHNIQUE
In this chapter, the basic mathematical model for estuarine circulation is de
scribed. First of all, assumptions and differential equations for estuarine circula
tion models are described in terms of a Cartesian (x, y, z) coordinate system. The
equations are then transformed into a verticallystretched (x, y, c) coordinate system.
Following that, equations in dimensionless form are presented. Anticipating model
applications to estuaries with complex geometry, a method for generating boundary
fitted grid is briefly discussed. Equations in the boundaryfitted (, 7, r) coordinate
system are then presented. Following a presentation of boundary conditions, initial
conditions, solution technique, and finitedifference equations are presented.
2.1 Differential Equations
The four basic assumptions used here for deriving the governing equations are:
(1) The flow is incompressible so that the flow is nondivergent. (2) The vertical
accelerations are negligible compared to gravity, thus, the vertical pressure distribu
tion satisfies hydrostatic approximation. (3) Density is approximated by its mean
value except when multiplied by gravity, i.e., the Boussinesq approximation is valid.
(4) The concept of turbulent (or eddy) viscosity and diffusivity is introduced. The
smallest scales of motion in real fluids, i.e., the dissipative scales, are always several
orders of magnitude larger than the molecular scales (Tennekes and Lumley, 1972).
In contrast to their molecular parts, turbulent viscosity and diffusivity are not fluid
characteristics but depend strongly on the state of turbulence. They are properties of
the flow instead of the fluid and hence are temporally and spatially dependent. Us
ing a righthanded Cartesian coordinate system (x, y, z) with the origin at the water
surface and z measured upward as shown in Figure 2.1, the equations may be written
as follows (Sheng, 1989)
Ou 9v Ow
T+ + 0=
ox dy dz
Ou 9u2 Ouv Ouw
S+ + y + oz
at 8Z By 8za
1 p
Sfv +
Po8OX
+ 9, T TY
by \ 9y
o (A o'u)
\ Oxz
+ 9z Tz
(2.1)
(2.2)
OT OuT OvT dwT
o + +N + z
OS 9uS OvS owS
9t Ox 9y 9z
f 1 p (Av\
 f + Ar I
Po Oy z ox]
0 (A, av\ 9 (A v\
+ ATy + A
9y \ 9y 8z \ 9z)
0p
9z
Oa (K T T
+ K + z
+ OD a
0 B
P = P(T, S) (2.7)
where (u, v, w) are mean fluid velocities in the (x, y, z) directions, f is the Coriolis
parameter defined as 20 sin 0 where Q is the rotational speed of the earth and 0 is
the latitude, g = 9.81ms2 is the Earth's gravitational acceleration, p is density, p is
pressure, T is temperature, S is salinity, (AH, KH, DH) are horizontal turbulent eddy
coefficients, and (A,, K,, D,) are vertical turbulent eddy coefficients. Variable eddy
Ov Ouv Ov2 Ovw
ot + x + y + O
DT tz ay az
(2.3)
(2.4)
(2.5)
(2.6)
Displaced Water
Surface (x= )
y,v
z,w
Bottom
(z=h)
Nominal Water
Surface (z=0)
Figure 2.1: The righthanded Cartesian coordinate system.
30.
0D 20.
CCO
10.
0.
20. 25. 30. 35. 40.
Salinity (0/,)
Figure 2.2: at density field calculated from Eckart formula.
20
coefficients are allowed in the present model. In addition to the above equations, the
program code contains an equation for dissolved species concentration similar to Eq.
(2.6).
Equation (2.1) represents the continuity equation as conservation of water mass.
Equations (2.2) and (2.3) are conservation of momentum in the x and y directions.
Eq. (2.4) expresses the hydrostatic balance. Equations (2.5) and (2.6) represent a
conservation of energy and conservation of salt, respectively. An equation of state for
seawater, Eq. (2.7), relating salinity and temperature to density closes the system.
In the present model, the equations given by Eckart (1958) are used
p = (P+ 1)/(a + 0.698P)
P = 5890 + 38T 0.375T2 + 3S (2.8)
a = 1779.5 + 11.25T 0.0745T2 (3.8 + 0.01T)S
where T is in C, S is in /oo, and p is in g cm3. Figure 2.2 shows at density field
calcualted from Eckart formula.
2.2 Governing Equations in Stretched Coordinates
To simulate circulation and transport in water bodies with gradual bathymetric
variation, it is common to use utransformation such that both the free surface and the
bottom become the coordinate surfaces with an equal number of coordinate surfaces
in between. This transformation, i.e., the so called astretching, leads to a smooth
representation of the topography with the same order of vertical resolution for the
shallow and deeper parts of the water body (Figure 2.3). The astretching introduces
extra terms to the original Cartesian equations of motion. However, most of these
extra terms appear in the horizontal diffusion terms, which are generally not very im
portant compared to the other terms. The agrid model is suitable for simulating flow
and salinity transport in regions of gradual bathymetric variation and gives a more
accurate estimation of bottom stress than the zgrid model, which resolves the depth
L
21
with "stairstep" grids. Recent studies by Sheng et al. (1989c) and Haney (1990)
found that, when there are sufficient grid points across regions of sharp bathymetric
change, the agrid model can yield good results. In case of insufficient grid points,
special care must be exercised. Sheng et al. (1989c) suggested direct evaluation of
the horizontal density gradient terms along the constant zplane and avoiding the use
of higherorder advective schemes along the sharp bathymetric variation to reduce
numerical error. This issue will be further discussed later.
The transformation from the Cartesian grid (x, y, z, t) to the verticallystretched
grid (x', y', a, t') is defined as
z (x,y,t)
= x, y' = y, ax, y, z, t) = (,+ t' = t (2.9)
h(x, y) + ((x, y, t)
where h is the water depth relative to mean sea level, and ( is the free surface ele
vation relative to mean sea level, and a is the transformed vertical coordinate such
that a = 0 at the free surface and a = 1 at the bottom. Using the chain rule, the
following partial differential expression can be derived
OF F Ox' OF Oy' OF Oa OF Ot'
Ox x + + + (2.10)
where F(x, y, z, t) = F(x', y', a, t'). The firstorder differential operators in the (x, y, z, t)
coordinate system are related to those in the (x', y', a, t) coordinate system as
a a a
9y 9y' 'Y9a
a 1 a
(2.11)
az H Oo
a a 1 aa 8
H[(1 + ) ]
at at' H Ot' a
where
1 8 aOh
,x = H[(1 + a) +a ]
OX, O]
j
22
zC
+h
Figure 2.3: The vertically stretched grid (From Sheng and Butler, 1982).
1
23
S= [(1 + 0) +a ] (2.12)
H = h+(
The hydrostatic equation in the transformed coordinate system is
ap = pg9z = pgHao" (2.13)
Integrating from the surface (a = 0) to some depth a,
p(a) = pa g pH d (2.14)
where pa is atmospheric pressure. Taking the xderivative while applying the coordi
nate transformation results in
Op apa ap 1 a9C H 9p H(
Ox 9x'g J [ H( + j ]H da + gpoa (2.15)
where po is the mean density. After some algebraic manipulation, Eq. (2.15) can be
written as
9p 9p, a OpH dH a(
Ox' 9 x' da + gap + gz (2.16)
x a9x' Jo 9x' _x1 ox'
Now we will drop all the primes(') for convenience. Using the flux conservative
form for the nonlinear advective terms and neglecting the atmospheric pressure gra
dients, complete equations of motion can be written as follows (higher orderterms
(H.O.T.) are shown in Appendix A)
Conservation of momentum in x direction
9Hu d9Hu2 aHuv a8Huw _H ppH o H) (.
S+ + + = fvH + g( 'da pa (2.17)
at ax dy a po J x ax
( u ( u\ 1 9 (Au\
gH + H[+ AR + A + a A,, + H.O.T.
9x Tx ( Tx) Ty \ 9y H2 9a 50a
24
Conservation of momentum in y direction
OHuv OHv2
+ +
ox ay
 gH + H[
By Tz
OHvw
+ 
Oa
(AH av
AOj
\z 9
Conservation of water mass
OC OHu
t +
fuH + gH(
po Jo
a (A,. v\ )
+ ATHT +
9y yj H
Conservation of heat
OHT
at
OHuT
Ox
OHvT
Oy
OHwT 0 f T\
+ a =Hw H= aKH 
0e O T\ ox)
+ H KH + aK,,T
Ty Ty H TO ( TO,
Conservation of salt
OHuS
+ +
8O
OHvS
Oy
OHwS
+ O 
00
H DH
OF Tx
where the vertical velocity, w, is
doH w
dt H
1 aO
H(1+ )  T7ru 7,v
which is a result of the a coordinate transformation.
2.3 Dimensionless Equations in Stretched Coordinates
The dimensionless equations make it easy to compare the relative importance of
various terms in the equations. In a properly nondimensionalized equation to be
OHv
at
pH OH
da pa)
Qy By
(2.18)
 \ A + H.O.T.
5a \ Da
OHv
+ y
By
OHw
+ 0
Oa
(2.19)
(2.20)
OHS
at
+ H IDH
Oy \ yJ
(2.21)
1 0
H aT
OSa )
( D.}
Oc9a
(2.22)
25
shown later, the part of each term contained within a parenthesis or bracket should
be of order unity, while the part of that term containing the dimensionless numbers)
should indicate the order of the term. This should allow one to compare the rela
tive importance of all the terms and to simplify the model equations under certain
circumstances. The following nondimensionalization procedure follows that used by
Sheng (1986b) and Sheng et al. (1991). The added insight justifies the slightly more
complicated form of the dimensionless equations.
Dimensionless Variables
(u*, vI*, W*)
(*, y*, z*)
t*
C*
p*
T*
A*
KH
K
D*
D:
W*
= (u, v, wX,/Z,) /U
= (x,y, zX,/Z) /X,
= (7,, r)/(p fZ, U,) = (r, )/r,
= tf
= To/(Tr To) q/(pocfZrTo)
= gc/(fUrXr)= (/Sr
= (P o)/(r Po)
= (T To)/(T, To)
= AH/AHr
= AV/AV,
= KH/KHr
= K,/K,,
= DHIDHr
= D/D,,
= wX,/U,
(2.23)
Dimensionless Parameters
Vertical Ekman
Lateral Ekman
Vertical Prandtl
Lateral Prandtl
Vertical Schmidt
Lateral Schmidt
Froude
Rossby
Number:
Number:
Number:
Number:
Number:
Number:
Number:
Number:
Densimetric Froude Number:
Ev = A,,/(fZ,) = tilt dm
EH = AH,/(fX,) = tilttdm
Pr, = A,,/Kvr = tvdh/tvdm
PrH = AHr/KHr = ttdh/ttdm
SC. = Avr/Dvr = tvdsltvdm
SCH AHr/DHr = ttds/ttdm
Fr = U,/(gZr)1/2 = tgelte
Ro = U,/(fX,) = ti/t
S= gZr/(f2X) (Ro/Fr)2 = (t/tge)2
S, = Sr/Z
= (pr po)/Po
FrD = Fr/^/e = tgi/tc
where the various t's appearing in the above expressions represent the characteristic
time scales associated with various physical processes in lakes and estuaries (Ta
ble 2.1).
Utilizing the dimensionless variables and parameters defined above and suppress
ing the asterisks for clarity, one can derive the following dimensionless equations
o( dHu MHv Ow
+0 + + PH = 0o
ot ax oy 00
1 oHu o8 E, a du
H t x H2 T, O +v
Ro (Huu o Huv +Huw
H a + By + oa
(2.25)
(2.26)
(2.24)
Table 2.1: Characteristic Time Scales of Various Physical
Processes in Estuaries and Lakes.
Physical Time Order of
Process Scale Magnitude
Periodic Forcing tf 1/w
Convection tc X /Ur
Inertia Oscillation t 1/ f
Vertical Turbulent
Diffusion tvdm, tvdh, tvds Z~ /Av, Z /Kr, Z, /Dv
Lateral Turbulent
Diffusion ttdm,ttidh, tds X2/AHr, X2/KHr, X2/DHr
Gravity Wave tge X, / r
Internal
Gravity Wave tgi X/ /gZxAp/po
+ E ) + ( aH + H.O.T.y
Ro aO Dp DH f O \
H d + O pda + ap
Fri Ja Dx Ja +B
 AE a+B
S x H2 DTo a)Do
(
dy
E, a
H2 da
(A oV) 
\ Qa
Ro (Huv DHvv
UH \ + Dy +
dHvw
To,)
+ EAH + ) AH + H.O.T.
 a( + A,, + B
y H2 d D(a
E. (
Pr, H2Da
Ro (aHuT
PH [X (
PrH I TX
aT
aa)
DHvT
ay
E, a aS
= D ( DS
Sc, H2da \ a
Ro (HuS DHvS
H x + Dy
(2.28)
DHwT
+D )
(2.29)
DHwS
a~o~)
EH
+ SC
SCH
p = p(T, S)
where H = h(x, y) + ((x, y, t) is the total depth. The H.O.T.s are omitted here for
simplicity.
1 dHv
H t
(2.27)
1 DHT
H at
1 OHS
H 9t
(2.30)
)H + a KH + H.O.T.
 O a y 8 OyS
D a a D + H.O.T.
TXDz Ox +y Ty Ds
29
2.4 Dimensionless Vertically Integrated Equations
The external mode consists of the surface elevation and the verticallyintegrated
velocities. The vertically integrated equations represent the "external mode" of the
3D model. The dimensionless vertically integrated equations are
0+ P + =0 (2.31)
at ax ay,
U H + r. 7 + V (2.32)
Ot Ox
Ro [ UU) + Uy
+ EH k AH) + (AI a
TX TX Ty ay
Ro H2 p
FrD2 2 ax
S H + D,
Ox
O H + ,y 7n U (2.33)
at ay
+ [a (AH ) + (AHV)I
T T9 T 9V\9j 9V
Ro H2Op
FrD2 2 y
H + D,
Oy
where (U, V) = f0 H(u, v)d are the vertically integrated velocities. The nonlinear
inertia, lateral diffusion, and baroclinic pressure gradient terms in Equations (2.32)
and (2.33) are obtained by vertically integrating the corresponding terms in Equations
(2.26) and (2.27), respectively. However, these terms contain idealized assumptions
30
about the vertical profiles of horizontal velocities and density, i.e., the horizontal
velocities and density are assumed to be uniform in the vertical direction, and hence
are not exact.
2.5 BoundaryFitted Coordinate Generation
The basic idea of the boundaryfitted coordinate system approach is to make
computations on curvilinear grids so that one of the curvilinear coordinates always
follows a boundary. According to the work of Thompson et al. (1974), the natural
coordinates ( and qr are taken as solutions of an elliptic boundary value problem with
one of the coordinates being constant on the boundaries. The same methodology was
used to generate the curvilinear grids for the present study. Detail discussions are
given Thompson et al. (1974, 1976, 1977a, b), Thames (1975), and Thames et al.
(1977).
The curvilinear coordinates are determined by solving an elliptic system
S + P = P(, 7) (2.34)
X + ryy = Q(r, 1) (2.35)
with Dirichlet boundary conditions,i.e., the values of and 7 are specified as constants
on the boundaries. (x, y) is coordinate system in the physical domain, and (, q) is
coordinate system in the transformed domain. Subscripts indicate differentiation. To
obtain the grid locations in the (x, y) plane corresponding to known uniform rectan
gular grid in the (, 7) plane, the dependent and independent variables in Eq. (2.34)
and Eq. (2.35) are interchanged. This results in a coupled system of quasilinear
elliptic equations for determining x(,, n7) and y(, 77) in the transformed plane. Fig
ure 2.4 shows the relation between physical domain and transformed domain.
axf 2/xa, + yx,7, + ax~P + yx,Q = 0
(2.36)
Physical Plane
n=2
1 ,
n=1
3
1
TVJiiili ILIZLLEEIIA21
'It 13 __
Transformed Plane
Figure 2.4: Field transformation (From Thompson et al., 1977).
Region D*
I + i
, I H
ayt 2/y,7 + yy,, + aytP + yy,Q = 0 (2.37)
where the coefficients are
a = x2 + y2
3 = xsxJ+y1y7
7 = x +y2 (2.38)
1 
P(M, ) = .aP
1 
Q(V, 77) yQ
J = x y, z y
Here P and Q are the control functions that serve to concentrate the grid lines.
With a weight function w, P = wI/w, Q = w,/w.
2.6 Transformation of Governing Equations Into BoundaryFitted Grids
2.6.1 Transformation Rules
Generation of a boundaryfitted grid is an essential step in the development of
a boundaryfitted hydrodynamic model. It is, however, only the first step. A more
important step is the transformation of governing equations into the boundaryfitted
coordinates. A straightforward method is to only transform the dependent variables,
i.e., the coordinates, while retaining the Cartesian components of velocities. John
son (1982) developed such a twodimensional verticallyintegrated model of estuarine
hydrodynamics. The advantage of such method is its simplicity in generating the
transformed equations via the chain rule. The dimensional forms of the continuity
equation and the vertically integrated momentum equations are shown by Equations
(2.17), (2.18) and (2.19). The resulting equations, however, are rather complex. Even
when an orthogonal or a conformal grid is used, the equations do not become any
simpler. Additional disadvantages are: (1) the boundary conditions are quite compli
cated because the Cartesian velocity components are generally not aligned with the
33
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 (Sheng, 1989).
To alleviate the problems mentioned in the previous paragraph, we have chosen
to transform the dependent variables as well as the independent variables. Equations
in the transformed coordinates ( 77) can be obtained in terms of the contravariant,
or covariant, or physical velocity components via tensor transformation (Sokolnikoff,
1960). Rules of transformation and dimensionless tensor invariant governing equations
are described in Appendix B.
2.6.2 Dimensionless Equations in BoundaryFitted Grids
It is clear that if a proper transformation is given by Equations (2.36), (2.37),
(2.38) and (2.39), one can compute the metric tensor and Christoffel symbols accord
ingly, and then expand Equations (2.48) and (2.49) in terms of the two contravariant
velocity components and sum up all the terms over the range of indices. This process
is rather tedious and highly prone to human error. To alleviate these problems, we
used a symbolic manipulator to expand the tensorinvariant equations. After rear
ranging the resulting terms of expansion, we obtain
C + L (Vrg/Hu) + 1(,?Hv) + w= 0 (2.39)
1 aHu 11 1 2 912 2 2
H 9t 9+ _L?_U + 2V
H (Huu) + (Huv) + (2D11 + D2)Huu
dOHuw
+ (3D12 + D22)Huv + D22Hvv +
+ a(A, a (2.40)
Ro H 0 g 119P 12'P d
Fr2(H g + g (dj )
+ g + g1O) 0 pdo + ap
34
+ EHAH(Horizontal Diffusion)
= g +g g ) u + v
Ro a a2
S (Huv) + (Hvv) + D,1Huu
+ (Di, + 3D2)Huv + (D12 + 2D2)Hvv]
+ L A a A, a (2.41)
H2, ( a9 f
Ro H 218P +22 p d
Fr2 a)r
Sg21H g22 Hd + ap
+ EHAH(Horizontal Diffusion)
where the horizontal diffusion terms are listed in Appendix C. The temperature and
salinity equations can be obtained according to the same procedure as
1 aHT E, a ( T\
H 9t Pr, H2a ~" a)
Ro 1 a H ) + Ro HwT
Y 0(9 voHuT) + ( HvT)i H 19a
+ EH [giT1,1 + g12T1,2
PrH
+ g21T2,1 g$22T,2] (2.42)
1 HS E, a D, S\
H at Sc, H2aa \ 9au
[~o ( HuS) + (VHvS)
RoaHwS EH
Roll awS +EH [g 11S,1,1 + g12S,1,2
H Be SCH
+ g21S,2,1 + g22S,2] (2.43)
Recently, the nonlinear terms in the momentum equations have been replaced by a
more compact formulation (Sheng, 1989). The uequation and vequation are shown
35
in Appendix D. The more compact equations do yield better numerical stability than
the earlier equations shown in Equations (2.40) and (2.41).
2.7 Boundary Conditions and Initial Conditions
Boundary conditions are required for computing the six dependent variables on
or near the boundary. If the closed boundary is fixed in space, the component of
fluid velocity normal to it must be zero (V n = 0, where n is the unit normal vector
to the boundary). If the fluid is viscous, an additional condition, i.e., the noslip
condition (V. g = 0, where sis the unit tangential vector to the boundary) is applied.
Application of the above closed and noslip conditions is made here which result in
the fluid at the boundary having zero velocity with respect to the boundary itself. At
the closed boundary, no normal flux of momentum, salt, or energy is allowed since
the velocity is zero. All the variables shown below are dimensionless based on the
nondimensionalization in section 2.3 except for S.
2.7.1 Vertically Stretched System
The boundary conditions at the free surface (a = 0) are
S(u 9v\ H
9T HPr,
q =  (2.44)
as E,
dS
=0
The boundary conditions at the bottom (a = 1) are
(u 9v\ H
A, 0'a TE(rE.,'by)
= Hr ZrCd 1/211) (2.45)
A,,
OT
= 0 (2.46)
8cr
36
dS
= 0 (2.47)
where ul and v1 represent the velocity components at the first grid point above the
bottom.
Along the shoreline where river inflow or outflow may occur, the conditions are
generally
u = u(x,y,a,t)
v = v(x,y,a,t)
w = 0 (2.48)
T = T(x,y,a,t)
S = S(x,y,a,t)
The kinematic boundary condition is also imposed at the surface.
aC aC aC
w = + u + v (2.49)
at 9x ay
The boundary condition required for the pressure is P = 0 at a = 0 which
expresses the assumption that either the contribution of the atmospheric pressure
to the pressure of the fluid column is negligible, or that the atmospheric pressure
gradients are unimportant.
The open boundary condition is much more complicated, since momentum and
mass transport occur through the boundary. Salinity and temperature need to be
specified with C or velocity at the open boundary. For the elevation C, there are
currently several options
'ma' /72rt
C=C(x,y,t)= n A, cos +s ,) (2.50)
n=1 T.
=0 or = 0 (2.51)
Dx ay
L
37
where An, T, and ,n are the amplitude, period and phase angle of the tabular tidal
constituents. When open boundary conditions are given in terms of C, the normal
velocity component is assumed to be of zero slope while the tangential velocity com
ponent may be either (1) zero, or (2) zero slope, or (3) computed from the momentum
equations. But in reality the velocity gradients are impossible to obtain from field
observation. In the present study, it is assumed that the velocity gradients are zero
at the open boundary.
The salinity along an open boundary or river entrance is generally prescribed
throughout the tidal cycle. The present model allows the salinity at the open bound
ary to be computed from a 1D advection equation during the outflow. For example,
along an open boundary on the left or right (Sheng, 1986)
9HS OHuS
+ Ro = 0 (2.52)
where the spatial salinity flux is evaluated from the salinity values at the boundary
and the interior grid point via a onesided differencing scheme. During the inflow,
however, the salinity value at the open boundary can either take on a prescribed value
or be determined from the 1D advection equation while using the boundary salinity
value and the prescribed salinity value to evaluate the spatial flux term.
To initiate a run, the initial spatial distributions of C, u, v, w, T and S need to
be specified. When initial data are unknown, one could start with zero initial fields.
When initial data are known at a limited number of locations, an initial field can
be generated by some interpolation scheme. It is desirable if the interpolated field
satisfies the conservation law governing that field variable.
For practical simulations, the present model usually assumes zero initial flow if
little initial data are known. In case of salinity simulation, the spinup time is longer
and an interpolation routine is provided to produce a reasonable initial field from
limited data points.
38
2.7.2 Vertically Stretched and Laterally BoundaryFitted System
The boundary conditions at the free surface (a = 0) are
(Ou Ou H
A,, = (7,, 7,,) (2.53)
OT HPr,
= q
0o E,,
dS
a 0
The boundary conditions at the bottom (a = 1) are
(Ou Ov\ H
a,, = =,
A,,
 HrZrCd [gllu2 + 2g12U1v1 + 2212] (Ui, V1)
Avr
9T
T 0
Oa
dS
o 0 (2.54)
where ul and vl are the contravariant velocity components at the first grid point
above the bottom.
Due to the use of contravariant velocity components, the lateral boundary condi
tions in the (,, 7r, a) grid are similar to those in the (x, y, a) system. Along the solid
boundary, the noslip 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.
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
in the transformed equations. Thus, the surface stress in the transformed coordinate
system is given by
I = +1 i 2 (2.55)
2 = T 1 0 71 2 (2.56)
ax 1y
39
where 71, 72 are the contravariant components of the stress in the transformed system
and Tr, 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.
2.8 Solution Technique
The solution of the difference equations may differ from the solution of the differ
ential equations because of a variety of errors (truncation error, dispersion error, and
roundoff error) even if the scheme is stable. For explicit methods, the stability con
dition is that the CFL (CourantFriedrichsLewy) criteria based on the shallow water
wave speed (Roach, 1976) should be satisfied. The basic numerical technique to en
hance the model efficiency is the modesplitting technique (Sheng et al., 1978; Sheng,
1985; Madala and Piascek, 1977), which basically solves the equations of motion by
splitting them into an "external" or verticallyaveraged mode and an "internal" or
vertical structure mode. The external mode is solved implicitly, and the internal
mode is solved explicitly except for the vertical diffusion terms. Semiimplicit meth
ods (Sheng, 1983; Sheng and Choi, 1989) can be introduced to take advantage of the
merits of both the explicit method and the implicit method. This approach solves
implicitly those terms in the governing equations which possess rather stringent sta
bility criteria if explicit schemes were used. These terms include those that govern the
propagation of surface gravity waves: local acceleration, and surface slope terms in
the momentum equations and all terms in the continuity equation. Furthermore, ver
tical diffusion terms in the momentum and transport equations are solved implicitly.
All other terms in the equations are solved explicitly.
The governing equations are solved by a finite difference scheme along with proper
boundary conditions. A spacestaggered grid system shown in Figure 2.5, is used to
40
discretize the differential equation (Sheng, 1983). In the staggered grid, elevation,
vertical velocity, temperature, salinity, and density are defined at the center of a cell,
and uvelocity at the left and right sides of a cell, and vvelocity at the top and bottom
sides of a cell. A staggered grid system has an accuracy of O(Ax2) with the same
mesh system compared with an accuracy of O(Ax) in a nonstaggered grid system.
For the time marching, forward time difference (twotimelevel) or leapfrog technique
(threetime level) can be used. The present study uses a twotimelevel scheme.
Both the leapfrog technique and spatial difference are centered differences, and the
accuracy in space and time is secondorder. The separation of solutions caused by
the leapfrog method is corrected by using a periodic time smoothing. The diffusion
terms are evaluated at the previous time level to avoid numerical instabilities.
In the computational procedure, first the spatial derivatives x, x,, yY, y,, and Ja
cobian J at each grid point are computed using the transformation equations (Eq.
(2.34) Eq. (2.38)). These derivatives are used as input to the circulation and trans
port model. The "external mode" equations are solved using an alternatingdirection
implicit (ADI) method. This solution technique consists of a twostep calculation
which advances the solution first from time level n to (an intermediate step) and
then from to n+1. In the first step, C and U are advanced from n to by an implicit
calculation along each row in the ( direction, while treating V and yderivative terms
explicitly. In the second step, C and V are advanced from to n + 1 by an implicit
calculation along every column in the ] direction, while treating U and iderivatives
explicitly. Then deviations of the velocities from the vertically averaged velocities are
calculated, while treating implicitly the vertical diffusion terms in the equations for
the deficit velocities. The deficit velocities are then added to the verticallyaveraged
velocities to produce the threedimensional velocity field. Conservation equations for
salinity and heat are then solved explicitly, except for the vertical diffusion term.
The implicit schemes used in the external mode and the internal mode require the
o U, u u
0 Vv w
Sw, T, p T, p
Oo A o AO o o
o n D A
0 0 0 0
Q0A 0n 0 0 0
0 0 0 0o
oA o A o A
oo
A o& Ao
a  7777'77
Figure 2.5: Staggered numerical grids in lateral and vertical directions (From Sheng,
1983).
42
inversion of only tridiagonal matrix systems.
2.9 Finite Difference Equation
In the following, a set of twotimelevel difference equations in (x, y, a) system
are presented. A set of threetimelevel difference equations in (, ri, a) system are
presented in Appendix E. The results presented in this study are obtained with the
twotimelevel equations. The finite difference representation of a dependent variable
can be written
'jk = O Y, a, t = (iAx,jAy, kAcr, nAt)
Sx. = (0+1/2 1/2
Finitedifference equations for the external mode and the internal mode are listed
in the following. The basic twotimelevel schemes are similar to those presented in
Sheng (1983) and Sheng (1987). For threetimelevel schemes, refer to Appendix E.
2.9.1 External Mode in (x, y, a) (TwoTimeLevel Scheme)
The external mode consists of the surface displacement C and the vertically
integrated velocities U and V and is described by Equations (2.31) through (2.33).
The finitedifference schemes treat all of the terms in the continuity equation implic
itly. The time derivatives and surface slopes in the momentum equations are generally
treated implicitly, while the bottom stresses are computed explicitly from the latest
vertical profiles of horizontal velocities assuming the presence of a turbulent bottom
boundary layer. In that case, one can either specify a drag coefficient or a rough
ness height, instead of the conventional ChezyManning coefficients which are used
in many two dimensional models. When used in twodimensional mode, however,
the model computes implicitly the bottom stresses in the u and v equations .'
with ChezyManning formulation. The dimensionless finitedifference equations in
the following are
AF"+1 = F" + AtD"
(2.57)
43
where superscripts n + 1 and n denote the present and previous timestep of integra
tion, A t is the timestep, and
1 cS, ay,
A = 7A 1 0
7A8 0 1
F= U)
D= D)
Dj
with D, and Dy as defined by Equations (2.32) and (2.33), and
PAt
x = A (2.58)
At (2.59)
Ht= (2.60)
HAt
S A (2.61)
where Ax and Ay are grid spacings in the x and y directions, respectively.
For increased efficiency of computation as previously discussed, the ADI method
factors the matrix A into the sum of three matrices: the identity matrix I, a matrix
AX, and a matrix A,:
0 a 0 (.
A. = 7,5= 0 0 (2.62)
0 0 0
Ay = 0 0 0 (2.63)
0 0
To solve the factored version of Equation (2.57), the ADI method proceeds by first
computing in the xdirections (the "xsweep"), and then computing in the ydirection
(the "ysweep"). Carrying out the factorization of Equation (2.57), and neglecting
terms of order At" yields
xsweep:
[I + A,]F* = [I A,]F" + AtD" (2.64)
ysweep:
[I + Ay]F"+I = F* + AyF" (2.65)
where F* represents the intermediate solution
F* = U* (2.66)
V*
The more detailed equations for the xsweep and ysweep are
xsweep:
(* + acSU* = (" Cae6V" (2.67)
U* + 7,C(* = U" + AtD' (2.68)
V* = V" 7SyC + AtDx (2.69)
The (* and U* are two unknowns in the xdirection. The above equation can be
readily solved by simple inversion of tridiagonal matrices along all the grid rows.
ysweep:
cn+l + aySyV+l = (* + a6yyVn (2.70)
V"+1 + y6,C+l = V* + 7,5,C" (2.71)
Un+1 = U* (2.72)
The ("+1 and V"+1 are two unknowns in the y direction. Again, inversion of
tridiagonal matrices along all the grid columns are performed to solve the above
equation by substituting the equation for V* from the xsweep into the equation in
the ysweep.
2.9.2 Internal Mode in (x, y, a) (TwoTimeLevel Scheme)
Defining deficit velocities as fi = u u and f v u, we obtain the differential
equations for the internal mode by subtracting the verticallyaveraged momentum
equations from the threedimensional momentum equations multiplied by H (Sheng,
1986)
1 t B x + I A, (1+ fi) (2.73)
H Ot H H 2 aa 00,
1 H D E, (U f)
= By +H ( + 6,) (2.74)
H at H H2a 1\"da
The finite difference equations are
(Hfi)"+1 = (Hii)" + At(HB, Dx)"
E a.
+ At a ((V + iL)n+) (2.75)
HH" an +
(H6)"+1 = (H.)" + At(HB, D,)"
+ At A (( + )n+) (2.76)
The vertical diffusion terms in the momentum equations are treated implicitly to
ensure numerical stability in shallow waters. It is also important to ensure that the
46
verticallyintegrated perturbation velocities are always equal to zero
kmax
Sfii,j, = 0 (2.77)
k=l
kmax
E i,,k = 0 (2.78)
k=1
In order to ensure that Equations (2.77) and (2.78) are satisfied, the nonlinear
inertia and turbulent diffusion terms in the verticallyintegrated Equations (2.32) and
(33) must be evaluated by summing the corresponding terms in the threedimensional
equations at all vertical levels.
Once i"n+1 and un+1 are obtained, u and v can be obtained from
un+l = jin+l + Un+l/Hn+l (2.79)
n+l = n+1 + Vn+l/Hn+1 (2.80)
Following these, the vertical velocities wn+1 and Wn+1 can be computed from
=A 8 n+l A [ MHu "+1 + Hv n1
w/+1 W 31l Ot Ub x k ] (2.81)
1 PH at H k ay k
n+11 + o d(n+1 + O A hA dh +
wn+l = Hn+lwn+l + + + v (2.82)
Sp dt x ay k(2
The vertical velocity w should be nearly zero at the free surface.
2.9.3 External Mode in (,, r,, a) (TwoTimeLevel Scheme)
The present study generally uses the following twotimelevel scheme instead of the
twotimelevel scheme described earlier. Verticallyintegrated differential equations
are written as follows (Sheng, 1986)
( + (v0U) + ( )= 0 (2.83)
Ut + Hgll" + Hg12 = U C,U + D, (2.84)
Vt + Hg'12 + Hg227 = 2V + D (2.85)
V>9o v9_o
1
47
where Ct is drag coefficient in idirection, and D, and D' are all the terms in the
xdirection and ydirection which are not shown in the momentum equations. The
above equations are written in finite difference forms as follows
n+ n + PAT [TX16(vAU)n+' + (1 T1)68(vs U")]
+ A [T1,(vV+') + (1 T1)6,(VrV")] = 0
(2.86)
Un+1 
Vn"+ 1
Instead of solve
ADI scheme as fol
.sweep:
U"
(2.87)
+ AT [T1(Cn+') + (1 T,)6s(")]
Hg12 TJ8,((n+l) + (1 _T) ,(
+ [Tc(+ (n)]}
nT _T 912 rn+1
 AT [T2CUn+ + (1 T2)CTU" T 3
(1 912 +
 (1 T3) U]I +Dx
g12o ] D
Vn + AT {H' [T6(Cn+) + (1 TI)6(n)] (2.88)
+ g22 [T6,(n+) + (1 Ti)6,)]
g12 LVn+l
= AT [T2C.V"+ + (1 T2)C.V" 1 + T3V1
+ (1 T3) 912 V] +Dgn'
ing the above equations directly, however, they can be written in
lows
(* + Tia6A(vU*) = C" a(1 Ti)6S(vU) )
78,,(,\yv")
(2.89)
48
Tiv ll6~((*) + [1 + At (C T3 ] (/fU*)
= Ti '26,(("C) + ( VgUn) (2.90)
(1 Tr)v 71r((") (1 Ts)g x26,g(")
(v nU) AT( T)C T At(1 T3) ] + D
77sweep:
C"+1 + Tras,((yv +) = (* + Tia,1S,(VyV") (2.91)
Ti 1((n1) + 1 + At f(T2C + T, 1 [avn+1] (2.92)
= Tlv/i126y *) (1 Ti)Vy7128,(c) (1 T) 7226 ")
+ (VVo) (VV') [At(1 T2)C + At(1 T3) ] + D
Ti 712(sn+( ) + (VUV"+~) = (sU*) + TiV, 126,((n1) (2.93)
where T1, T2 and T3 are all constants between 0 and 1. For example, the wave
propagation terms are treated explicitly if T1 = 0, but implicitly if T1 = 1. m can be
n 1 or n. Additional parameters appearing in the equations are
jfat
i goA
pat
11 HA (2.94)
12 HAtg12
S22 HAtg22
A77
21 HAtg12
A77
The system of equations can be rewritten in the following matrix form
[A]F"+1 = AFn1 + D
[I + AJ]F* = [B A,7]Fn1 + D
[I + A,]F"+1 = F* + A,F"1
ere
1 aTib6 a yTi
A] = V "Tis T 1 + AT (C,T, 2T3) 0
0. / .22T, 1, 0 1 + AT (C,.T.2 + 912T3
L
\ "' /go
0 a.Ti15
11"T AT (C,~T 12 T)
V90 v9_0
[Ad =
[A,] =
0 0
0 0
ayT6~
(2.100)
0 AT (C,, T2 + 91T3)
1 0 0
[B]= 0 1 0
0 0 1
F= vgVo
vsfg
wh
(2.95)
(2.96)
(2.97)
(2.98)
(2.99)
(2.101)
(2.102)
I
}
50
aY(1 Tib1(VgU")) ay(1 Ti)S(VV")
JU"AT ((1 T2)C7 + (1 T3) ) /"(1 Ti)6C
D = v,Y21C, + ATD' No (2.103)
gV"AT ((1 T2)C,, (1 T3) _) 712
NO722(1 Ti),"Cn + ATD'/g
2.9.4 Internal Mode in (, 77, a) (TwoTimeLevel Scheme)
S12 g22
Hi Hf + H + F3
E, 1 OHii
+ A, o (7e (r ) F2 (2.104)
S gll ~ g21
Hi = Hii HO + G3
E, 8 8Hi3
+ 02 AV a ) (r, rb,) G2 (2.105)
and F2 and G2 indicate all the explicit terms in U and V equations, respectively, and
F3 and G3 indicate all the explicit terms in the u and v equations, respectively.
(Hi)n+1 = (Hii)" + 2At Hn+lin+l
22
+ 2At H"" + 2At(F3 F2)"
+ 2At () A, (H+n+1)
(Hn)2d pToI ]
(Ts t b)n+l (2.106)
11
(Hu)" = (Hu)" 2At H"n"n+
21
2At H+1iin+1 + 2At(G3 G2)"
+ 2At "a [AV (Hn+n+)
(H")2 1 (2.10
(rs s)"n+l (2.107)
51
2.9.5 Salinity Scheme (TwoTimeLevel Scheme)
The finite difference equation for salinity are similar to the internal mode solution.
The vertical diffusion terms are treated implicitly, while the horizontal diffusion and
advection terms are treated explicitly. The salinity equation in the ( 0, 7, a) coordinate
system is shown in the following (Sheng, 1986)
HS E, ( as\
at HS a da) + G (2.108)
where G is advection plus horizontal diffusion terms, and
A,
D,?
S= AHR
DHR
Ro a a a8 1
G = (HuS/) + (HvSV ) + (HS)
HEH ( 18S d 8S 228S
+ SH 1 + 2' 2+
Using i, j, and k as grid indices in the , i, and adirections, respectively, the
finite difference equation for salinity can be written as (Sheng and Choi, 1989)
Hn+l1n~l n n At Ro,
H Sk HSk = A R (HuS )" (2.109)
+ (HvSg)' + V, (HwS)]
E,(At) 1 [ D+ ,n+l S1
+ Hn+S Aak A+ (S+1
DV /on+l C_+l t
a \ (sk l ]
HEH. At Qg S 2g S a Ss
+ g + 2g12 + a 2s I7
Sc (9 aQ,977 Oq)
Dividing by H"+1 gives
E, At (D ,_ D,+ (+1
(H+12S,. Aak Aa_ 1+ Siak+l (2.110)
+ 1 + ( E, At D,+ + D7. S"i
Si +(Hn+1)2S 1, k (aa + o )] S'Vk
n Sn
fHn+l 'jk
At Ro
H"+1 Vga
Hn EH At ( 1192S 2g2 2S
+ n S 2g9 [ 2
Hn+l SH a52 8489
(HuSn + S(HvS ) + H vS(HwS)
a~ d0,
g2 S)
+ g2
asz
Rearranging Eq. (2.110) gives
7k k1 S~Sk1 +
[1 + 7k(ak_ + ak+ )] n+lk 7_kk in+1
Hn
= SHn+S,k + Gn
This is applied for 1 < k < km: where
EAt
k (H+,) 2 ScAk
D,+
ak+l/2 =
Dv
ak1/2 
Gk = advection + diffusion
For k = 1, i.e., the lowest alevel:
[1 + 713/2] Sa,3 7l3/2S
Hn
=H.+I S1i + G
(2.111)
(2.112)
(2.113)
(2.114)
(2.115)
(2.116)
53
For k = km, i.e., the uppermost alevel:
n+l [l l (2.117)
i7kmakm_ Ss i +1 + 7kmIakm~ Sijk (2.117)
H"
Hn+1 Sj,km + Gkm
Equations (2.111), (2.116), (2.117), can be written as the following tridiagonal system
of algebraic equations
A1(K) S(K 1) + A2(K) S(K) + A3(K) S(K + 1) = A4(K) (2.118)
where Al, A2, A3, and A4 are defined as
A0K) = 0 K=1
A() 7k k1/2 1 < K < KM
S 1 + 71a3/2 K = 1
A2(K) = 1 + 7k(akml/2 + k+1/2) 1 < K < KM
1 + 7fkm1/2OCkm1/2 K = KM
(2.119)
A3(K) = 'ykak+1/2 1 < K < KM
0 K=KM
A4(K) = S,k + G 1 < K
The above tridiagonal system of equations given by Equation (2.118) is efficiently
solved by a tridiagonal solver first developed by Thomas (1954), which is often called
"The Thomas Algorithm".
2.9.6 Advection Terms
1. Central Difference Scheme
The central difference option has the advantage of secondorder accuracy, but a
disadvantage of numerical dispersion. When a large concentration gradient exists in
54
the model domain with insufficient grid points across the gradient, negative concen
trations may sometimes develop in the simulation. The finitedifference form is
S(HuS ) (HS/^u)+ (HS.u)_
(HuSqy") =(2.120)
1 r(HS~Jf)i,j,k + (HSV )i+l,j,k
= A 2 ui+l,i,k
(HSvj)ij,k + (HSv),)i1,,, Ui, k
2 ijk
8 HS (HSv v)+ (HSv vv)
(HvS ) )(2.121)
1 "(HSV )i,j,k + (HSv)i,j+l,k
= 2H vijv+lk
An 2
(HSNo/)i,j,k + (HSvo)ij1,k
2 i
(HwS) = wk 1 (Sij,k + Si,j,k+1) Wki (j, + i,k1)] (2.122)
2. Upwind Difference Scheme
The upwind difference option has the advantage of increased stability versus the
central difference option, but at a cost of increased numerical diffusion. The upwind
difference option guarantees positive salinity concentrations. The finitedifference
form is
(HuSV~ ) = (S+Ui+lI,k SujI,k)/A (2.123)
where uT = Huv/g.
if Uf+l,,k > 0
S+ = Si,j,k (2.124)
if +,j,k < 0
S+ = si+lj,,k
(2.125)
Thus
S+i,j,k = (+ijk1 + juTi + 1,j, k ) Si,j,k (2.126)
+ 1(U T JUT
+ 2 +1,j,k I+1j,k S i+l,j,k
Likewise,
Su, = (T,k + IUJ,k) Si,l,k (2.127)
+ (U.,k Iui,kl Sij,k
3. Combined CentralUpwind Difference Option
The combined central and upwind option gives the advantages of central and
upwind options together, but minimizes the disadvantages. The combined option is,
however, more computationally demanding. The finitedifference form is
a(HuSv ) = (S+uT+l,k Sui,,k)/Al (2.128)
where uT = HuVo.
if ,,k > 0 and Sij,k > Si+1,,k S+ = ci,k + Ci+,j,k (2.129)
Uj3 2
if uj,k > 0 and Sij,k < Si+lj,k S+ = Sij,k (2.130)
if uTjk < 0 and Si,j,k > Si+,j,k S+ = Si+lj,k (2.131)
if u,,k < 0 and Sij,k < Si+lj,k S+ = Sij,k + Si+,i,k (2.132)
4. Higher Order Advective Schemes
To improve the accuracy of advective terms in the numerical model, it is possible
to use higher order advective schemes which have formal truncation errors that are
third order or fourth order in Ax. For example, Sheng (1983) applied the Flux
CorrectedTransport (FCT) scheme (Zalesak, 1979) to simulate the advection and
diffusion of salinity and suspended sediments by winddriven currents, and compared
56
the performance of the FCT scheme versus the three advective schemes discussed
above. The performance of the combined central upwind option was found to be
quite comparable to that of the FCT scheme. More recently, Johnson et al. (1989)
applied the QUICKEST scheme developed by Leonard (1979) to the curvilinear grid,
and obtained improved results when the zgrid version of CH3D was used. Sheng et
al. (1989) found that, when a agrid model is used, higherorder advective schemes
should be avoided in regions of sharp bathymetric gradient when insufficient grid
points exist to resolve the sharp gradient.
The higherorder schemes have been primarily applied to solve the advection
and diffusion equations but not to solve the NavierStokes equations. Several other
higherorder schemes can be used to improve the accuracy of advective terms. The
application of higherorder advective schemes requires additional programming and
computational effort, due to the fact that a differencing formula involving four or five
points is needed to represent the firstorder derivative term.
The QUICKEST scheme used for this study has third order accuracy. This
method uses cubic upstreamweighted interpolation with five points of i2, i1, i, i+
1,and i + 2. Because this is a fivepoint interpolation the fluxes near boundaries must
be calculated using lower order schemes.
The finite difference equations for advection terms using the QUICKEST scheme
are in idirection,
8 1
S(HuSV/a) = (Ui+l,j,k Iui+1,j,k) (Si,j,k + Si+,j,k)
(1 (ui+,,kAt)2) (Si+2,j,k 2Si+lj,k + Si,j,k)
11
~n+lj~kat (si+l,,k Si,j,k)
+ (ui+1,i,k + IUi+1,j,k) (Si,j,k + Si+1,j,k)
(l (uitlj.,kAt)) (Si+l,,k 2Si,k + Si1,j,k)
1 1
ui+1,j,kAt (Si+1,j,k Si,j,k) 2 vgwinuHi+1,j,u
57
1
(uj,k 2 (SSijk) (2.133)
1( (ui,j,kAt)2) (Si+1,Ik 2Siik + Si1,j,k)
ui,j,kAt (Si,j,k Si,j,k)
11
+ (.,+j,k + Ij,skl) 2 (Si1,k + Sij,k)
1 (1 (ui,jkAt)2) 2Sil,i,k + Si2,j,k)
1 1 1
i,j,kAt (Si,j,k Sil,j,k) jg jHi ,y 
2 2 Hij,< "'jc
and in the qdirection,
(HvS =g) (vi+x,j,k Ili+lj,kl) (Si,j,k + Si+lj,k)
(1 (i+l,,k At)2) (Si+2,j,k 2Si+1,j,k + Sij,k)
1 1
i+l,j,kAt (Si+l,j,k Si,j,k) / H i+,jV
1
(vi+l,j,k Ivj,kI) (Si,j, + Si,,j ,k) (2.134)
(1 (vil,kA t)2) (Si+l,j,k 2Si,j,k + Si1,,k)
Vi,,kAt (Si,j,k, Sil,,k)
+ (ijk Ii,j,kl) 2 (Sil,j,k + Sidj,k)
(1 (vi,j,kAt)2) (Si,k 2Sii,,k + Si2,j,k)
1 1 1
Vi,,kAt (Si,j,k Siij,k) ~ v Hij,v H
2 (VijjHk +* Vj+J,<
2.9.7 Horizontal Diffusion Terms
The finitedifference form of the horizontal gradient in the horizontal diffusion
terms is the standard secondorder divided difference. There is no need to treat the
horizontal diffusion terms implicitly because the time step constraint associated with
58
them are usually much larger than that allowed by the advection terms.
11 2S 11 Si+l,j,k 2Sij,k + Si,j,k (2.135)
g g ( (2.135)
(2.136)
82S
g12 g12 [(Si+2,i+1/2,k Si+1 1/ 2) (2.137)
(Si.1/2,i+1/2,k) Si1/2j1/2,k Si1/2,I1/2,k)] /A, A1
22 2S 22. Si,j+l,k 2Sij,k + Si,jl,k (18
gJk = gf, (2.138)
a772 (A77)2
2.9.8 Numerical Stability
Since the numerical scheme used here is not fully implicit, the time step of numer
ical computation cannot exceed the limits associated with the advection terms, the
Coriolis terms, the baroclinic pressure gradient terms, and the horizontal diffusion
terms
(At)advection <
(At)coriolis < 1/f
(At)b.oinic 1 (2.139)
2 gHmaX( +Ly)max
(At)diff.sion < I 
If the allowable internal time step according to the above equations is much (10
to 100 times) larger than the at allowed by the surface gravity wave as shown in
Equation (2.139), then it is feasible to use a small time step for the internal mode
(At,). If the difference between the allowable internal time step is only a few times
larger than the external step, then it is better to simply use the same Ate and At,.
The selection of time step is still very much an "art" because all the time step limits
59
are based on stability analysis which is strictly valid only in a linear system without
considering the effects of nonlinearity and boundary conditions.
2.10 Turbulence Model for Mass and Momentum Exchange Coefficients
Various models are available to parameterize the vertical turbulent mixing of
momentum, heat, and salinity in the water column. For reviews of turbulence models,
the readers are referred to Lewellen (1977), Rodi (1980), Mellor and Yamada(1982),
and Sheng (1986). The starting point of a turbulence model is to derive the mean flow
equations from the conservation equations for mass, momentum, energy, and species.
For incompressible flows, conservation equations can be written in tensor notation as
aui
S= 0 (2.140)
+ U + gi Po 2eijkfjUk + (2.141)
xat Po dOz Po x, axj
+ = + 9 k (2.142)
where q is either density or temperature or salinity, p is density, Po is reference density,
eijk is the permutation tensor, fj is the rotational speed of earth, and v and k are
the molecular viscosity and diffusivity.
When the variables are separated into mean and fluctuating quantities, the mean
flow equations can be written as
= 0 (2.143)
dx1
_8i Nuug 1 p (p Po)
D" + ,77 1 (P+ po (2.144)
at Oxj oxj Po Di po
26ijkj'Uk + ( (2j
a0 + 7 x D + a x j (2.145)
9t 9xj 9xj x j
60
where ui, uj, and u are the mean velocity components, u u and u' are the fluctuating
velocity components, 4 and 0' are the mean and fluctuating temperature or salinity,
and p is the mean pressure. This system of equations is not closed because of the
Reynolds stress terms.
Basically, turbulence models can be classified as eddy viscosity models and second
order closure models. The eddy viscosity models assume the turbulent Reynolds
stresses to be the products of mean velocity gradients and "eddy viscosities," and
the turbulent heat and mass fluxes to be the products of mean temperature and
concentration gradients and "eddy diffusitivities"
uu = h x + 6*S* (2.146)
u = Kt (2.147)
_xi
where q2 = Uui,, and At and Kt are turbulent eddy viscosity and diffusivity, respec
tively.
The eddy coefficients are assumed to be certain known functions of, for example,
mean flow parameters (Munk and Anderson, 1948), depths and wind speeds. These
functions, however, are difficult to determine and are often adhoc in nature, and
hence may have to be adjusted for application to a new site or a new flow situation.
Obviously, this limits the "predictability" of the hydrodynamics model.
The secondorder closure models, on the other hand, resolve the dynamics of
turbulence by including the differential transport equations for the turbulence vari
ables; i.e., the secondorder correlations (uu,, u', and V'70'). In the most compli
cated case, i.e., the Reynolds stress model, differential transport equations are solved
(Lewellen, 1977; Sheng,1982; Sheng and Villaret, 1989). The time averaged mo
mentum equations are subtracted from the instantaneous NavierStokes equations to
produce the primedvariable equations. The resulting equation for the icomponent
is multiplied by u', and j component equation by ut. Adding two equations and
61
averaging by time, gives the u uequation.
at
au;.u _ uuu u; u
u'' F BY U U '1s
+ auk U k k +9i +9j
Szk Z k k k __
2eiklku'u 2ejlc2kl kU (k U)
u' ap' ui p' a2Uuj au aOu'
Po xj X+ pv 2vxk x
Po 82j Po 8k; 4 dz~k azk
 + =a u u + g 2icjkjU'I
at u aj# a71 ip' 2 ,9
OuU0 T 19j7 a277
v V kua
axj Po 9xi ax 9x
a7fol a~o a0; a
S= 2u' (u'')
at a7  Xk
+ ka 2k
9x x. 9xj
(2.148)
(2.149)
(2.150)
where xi are coordinate axes, t is time, (ui, ij,, uk) are the mean velocity components,
(u', u', u') are the fluctuating velocity components, j and 0' are the mean and fluc
tuating values of density or temperature or salinity, Cijk is the alternating tensor, and
Q is earth rotation.
The above system of equations is not closed since each equation contains unknown
terms. Following the secondorder closure model described in Lewellen (1977) and
62
Sheng (1982), the above equations are replaced by the following equations.
ot + 17 = (Evolution)
 Uuk aO j Uk aXk
+ i,7 l
 2eikIfkIuU 2ejlknlUkU(
+ vC (qAS
xk ( k
A U;; 35
_ 2b 3 2avuu
"3 12A A2
at j a9xi
 u ui ~u rj
+ gi
 2ejfk(juh')
+ vC a(qA2E i
A Au 
(E
(P
(Production)
(Buoyancy)
(Rotation)
(Diffusion)
(Tendency)
(Dissipation)
evolution)
reduction)
(Buoyancy)
(Rotation)
(Diffusion)
(Dissipation)
(Evolution)
2u W'9 (Production)
+ vcL (qA E) (Diffusion) (2.53)
2bg'' (Dissipation)
where Sij is the Kronecker delta, q = (uu1)1/2 is the total rms fluctuating velocity,
and A is the turbulence macroscale which is a measure of the average turbulent eddy
size. A total of five model coefficients: a, A, b, ve, and s are contained in the above
equations. By matching model results to data from critical laboratory experiments
(2.151)
(2.152)
t , +79a
at a a xi
63
where only one or two of the turbulent transport processes are important, these model
constants were found to be (Lewellen, 1977) a = 2.5, A = 0.75, b = 0.125, vc = 0.3
and s = 1.8. An additional equation for A is needed to close the system.
The present study tested several methods for calculating the vertical turbulent
eddy coefficients: (1) constant eddy coefficients, (2) MunkAnderson type eddy coef
ficients, (3) a simplified secondorder closure model of turbulent transport, i.e., the
equilibrium closure model, and (4) another simplified secondorder closure model, i.e.,
the turbulent kinetic energy (TKE) closure model. These methods are discussed in
the following.
Constant Eddy Coefficients
Constant eddy coefficients are generally used in preliminary model runs. The
rule of thumb is that the vertical eddy coefficients should increase with increasing
wind stress and/or tidal amplitude. Coefficients can be adjusted until a reasonable fit
exists between model results and data. In the absence of data, however, it is difficult
to adjust the constant eddy coefficients to give good results.
MunkAnderson Type Eddy Coefficients
In nonstratified flow situations, a variable eddy coefficient of the following form
is used
AW r = + (2.154)
v H [9a) ) (
where u and v are dimensionless horizontal velocities, the turbulence macroscale Ao
is assumed to be a linear function of o, increasing with distance above the bottom or
below the free surface, and with its peak value at middepth not exceeding a certain
fraction of the local depth. In the presence of strong waves, turbulent mixing in
64
the upper layers may be significantly enhanced. In such a case, the length scale Ao
throughout the upper layers may be assumed to be equal to the maximum value at
middepth.
In stratified flow situations, the effect of buoyancy on eddy coefficients can be
parameterized in terms of several Richardson numberdependent stability functions
A, = A,,o1(Ri); K, = K,,o2(Ri); D, = D,,02(Ri) (2.155)
where
gHZ,e Op [/uu (Ov 1
Ri + e=)( a [lap +a +ep) (2.156)
where A,,, Kvo, and D,, are eddy coefficients in the absence of any density stratifica
tion, and 1 and 02 are stability functions traditionally assumed to be of the following
forms
q1 = (1 + al Ri)~; q2 = (1 + a2 Ri)O2 (2.157)
As discussed by Sheng (1983), there exist great discrepancies among the many
existing empirical forms of the stability functions. To determine the coefficients in
these stability functions, one needs sufficient data to achieve the "best fit" between
modeled and measured results. As such, the specific formulae and coefficients will
vary with the problem, the environment, and the available data. Munk and Anderson
(1948) used the following formula in a study of the marine thermocline
01 = (1 + lORi)1/2; q2 = (1 + 3.33Ri)3/2 (2.158)
Others (Kent and Pritchard, 1959; Blumberg, 1975; Bowden and Hamilton, 1975)
used 01 and 02 in the form of Equation (2.157) but with coefficients different from
those in Equation (2.158). Not only the form of the stability function may vary
with the problem, but different Richardson numbers may also be used for different
types of problems. For example, the formation and deepening of the thermocline in
65
a relatively shallow basin depends strongly on the relative importance of wind stress
and heat flux at the free surface. In such a case, the following Richardson number
could be used
R= 2gHZa2e 8p
Ri = U2( + pa (2.159)
U,2n2(1 + ep) aa
where K is the vonKarman constant, and u. is the dimensionless friction velocity at
the free surface.
A Simplified SecondOrder Closure Model: The Equilibrium Closure Model
An equilibrium closure model (Sheng, 1985; Sheng and Chiu, 1986; Sheng et al.,
1989) has been used to compute the vertical turbulence. As discussed before, the local
equilibrium condition is valid when the time scale of mean flow is much larger than
that for turbulence and when turbulence varies little over the turbulence macroscale.
In this case, the evolution terms and the diffusion terms in Equations (2.151)(2.153)
become negligible with respect to other terms. The equilibrium closure model is
significantly simpler than the complete secondorder closure model (Reynolds stress
model) and has been found to give very good results in mean flow, salinity and
temperature with little or no tuning of model coefficients.
In addition to the mean flow equations, a set of algebraic equations are solved for
the secondorder correlation quantities to obtain the stability functions q1 and q2 in
terms of the mean flow variables. These equations, when written in dimensional and
tensor forms, are
0= u + S + (2.160)
zxk dexk po po
2eiktIkU't 2EjtkStU kU
q 4q q
A 3 12A
0 = U , ug' (2.161)
axj 9Xj PO
2eijk~Iukp 0.75q
0 2uW 0.45qp (2.162)
0 = 2up' + A
where the subscripts i,j, k can take on the values of 1, 2, or 3. Hence Equation
(2.160) represents six equations and Equation (2.161) represents three equations in
general. A total of five model coefficients are contained in the above equations.
These coefficients were determined from comparing model results with data from
critical laboratory experiments where only one or two of the turbulent transport
processes are dominant. These coefficients remained "invariant" in application of
the equilibrium closure model, the kinetic energy closure model (Sheng and Villaret,
1989) and the Reynolds stress model (Sheng, 1982; Sheng, 1984; Sheng and Villaret,
1989). A graphical comparison of this formulation versus some empirical formulae is
shown in Figure 2.10.
As shown in Sheng et al. (1989), q2 can be determined from the following dimen
sional equations when mean flow variables are known
3A2b2sQ4 + A[(bs + 3b + 7b2s) Ri Abs(1 2b)]Q2 (2.163)
+ b(s + 3 + 4bs) Ri2 + (bs A)(1 2b) Ri = 0
where A, b, and s are model constants, Ri is the Richardson number and
q=QA + \ (2.164)
Once q2 are computed, the vertical eddy coefficients can be computed from
:b)
4.0
3.2
,Blumberg
Kent & Pritchard
r Bowden & Hamilton
S2(Ri)
~1 (F
I I
0 2.5 5.0
Ri
_________i
10 .5 .4 2
.2 .4 .b
Figure 2.6: Comparison of stability functions for (a) some epirical
formulae versus (b) the simplified secondorder closure model (From
Sheng, 1983).
c.50
_C_
2.4
1.6
*K.
A+w w'w'
A, = + wW Aq (2.165)
A wq2
K, = b W Aq (2.166)
(bs w)A q2
where
w = Ri/(AQ2) (2.167)
iw = w/(1 w/bs) (2.168)
1 2b
w'w = )2 (2.169)
S 3(1 2)
In the complete secondorder closure model, a differential transport equation for
the turbulent macroscale A is derived. The equation is
9A 9A A Bui 8 9 A
+ = 0.35 uu + 0.75q + 0.3 ( )
et Ox, q2 q Ox1j axi Ox
0.375 OqA 0.8A u'.
q( ) + gi (2.170)
q 9x q O
However, the A equation contains four model coefficients that must be calibrated
with experimental data. For ease of application but without loss of generality, the
turbulent macroscale A is assumed to satisfy a number of integral constraints. First
of all, A is assumed to be a linear function of vertical distance immediately above
the bottom or below the free surface. In addition, the turbulent macroscale A must
satisfy the following relationships (Sheng and Chiu, 1986; Sheng et al., 1989)
A < 0.65 (2.171)
dz
A < C H (2.172)
A < C.Hp (2.173)
A < C2 q2 (2.174)
A < (2.175)
where C1 is a number usually within the range of 0.1 to 0.25; H is the total depth; Hp is
the depth of the pycnocline; C2, which is between 0.1 and 0.25, is the fractional cutoff
limitation of turbulent macroscale based on 8q2, the spread of turbulence determined
from the turbulent kinetic energy (q2) profile; and N is the BruntVaisila frequency
defined as
N = p 1/2 (2.176)
Although the simplified secondorder closure model presented above is strictly
valid when the turbulence time scale (A/q) is much less than the mean flow time
scale and when turbulence does not change rapidly over A, it has been found to be
quite successful in simulating vertical flow structures in estuarine and coastal waters.
Recent simulations of windinduced mixing of salinity in Chesapeake Bay (Sheng et
al., 1989a; Johnson et al., 1989) and winddriven circulation and sediment transport
in Lake Okeechobee (Sheng et al., 1989b) gave good results with little or no tuning
of model coefficients.
Turbulent Kinetic Energy (TKE) Closure Model
A simplified TKE closure model (Sheng and Villaret, 1989) has been used to compute
the vertical turbulence. When the time scale of the mean flow evolution is much
greater than the time scale of turbulence (A/q), the evolution terms and the diffusion
70
terms in the secondorder closure model become negligible with respect to other terms.
Instead of calculating q2 by algebraic equations, q2 is obtained in a dynamic equation
9q2 8_q2 u'q
+ uk = 2uIu + 2g,
9t OXk k k
+ 0.3 a(qA (2.177)
aXk ak 4A(2.177)
In the onedimensional vertical case of the ogrid, the equation is
q2 1 9Huq2 wq2 1 8 Aq
Ot H Ox H =ao H2 ao q a)
Oi+ + HOa H 2 (0.3qA )
A, Ou 2gK, dp q3
+ 2 (a)2 + g O 4A (2.178)
H2 90 (H aa 4A
Applying the rules in section 2.3, nondimensionalization of the above equation gives
(all the variables are now nondimensional)
9q2 Ro Huq2 Ro,8wq2 1 0 Oq2
+ + = . 3qA )
at 'H ax Haa H2 da da
E,A, Ou 2E,K, 8p q3
+ 2:E (0 )2 + F H a Ro, (2.179)
H2 Tb Fr2H Ob 4A
where A, and K, are nondimensionalized with A,,, Ro, = U,/fZ, is a vertical type
Rossby number, and E, = A,,,/fZ2.
Again, integral constraints for A are used as the equilibrium closure model.
CHAPTER 3
MODEL ANALYTICAL TEST: COMPARISON BETWEEN MODEL RESULTS
AND ANALYTIC SOLUTIONS
When numerical models are used to simulate circulation and transport in realistic
estuarine environments, it is common to find discrepancies between numerical results
and data. The discrepancies often result from a variety of errors: errors associated
with measurement, errors associated with the formulation of the physical problem
(e.g., dimensionality, turbulence closure scheme, moving boundary scheme, etc.), er
rors associated with the formulation of the numerical problem (e.g., finitedifference
schemes or finiteelement schemes for various terms of the governing equations), and
errors associated with the computer and programming (e.g., computer roundoff error
and programming error). Before the model is applied to a realistic estuary, it is nec
essary to test and improve the numerical accuracy of the various numerical schemes
used in the model. This chapter presents a study on the numerical accuracy of the
threedimensional hydrodynamic model by comparing model results with a number of
analytic solutions describing the circulation in idealized basins forced by a constant
surface slope, seiche, tide, wind, and density gradients.
3.1 Forced Constant Surface Slope
The first test considers the steadystate response of a water body in a rectangular
basin forced by constant surface slope with uniform density. At steady state, the
surface slope balances the bottom stress. Assuming onedimensional flow with the
quadratic bottom friction, the governing equation is given as
gh = Cd u I U (3.1)
with boundary conditions ( = a at x = 0, and ( = a at x = I (Figure 3.1),
where u is the verticallyaveraged velocity, C is surface elevation above mean sea level,
h is water depth, g is gravitational acceleration (= 981 cm/s2), 1 is the basin length,
and Cd is the bottom friction factor. A nonuniform rectangular grid as shown in
Figure 3.1 was used. The basin is open at the left and right ends. Model simulation
was started with zero elevation and velocity initially. Model tests were performed
for the cases of linear bottom friction and quadratic bottom friction. The following
parameter values were used
a = 25 cm
I = basin length = 62 km
basin width = 14 km (3.2)
h = 10 m
Cd 0.003 for quadratic bottom friction
0.3 for linear bottom friction
According to Equation (3.1), the steady state velocity was found to be u =
52.17 cm/sec for the quadratic bottom friction case, and u = 27.2 cm/sec for the
linear bottom friction case.
Using a time step of 360 seconds, model results reached steady state in about 300
time steps. The steadystate verticallyaveraged velocities in the xy plane as shown
in Figure 3.2 agree well with the analytical results. Figure 3.3 shows the velocity
distribution in a basin which is rotated in the xy plane by 300 in the counterclockwise
direction. The results appear to be the same. The model was then tested using an
irregular curvilinear grid shown in Figure 3.4. The same steady velocity is reproduced
(Figure 3.5).
Figure 3.1: A nonuniform rectangular grid system for testing the forced surface slope
test.
50 cm/sec
Figure 3.2: Simulated steadystate verticallyaveraged velocity field in a rectangular
basin forced by a constant surface slope. Horizontal grid is shown in Figure 3.1.
50 cm/sec
Figure 3.3: Same as Figure 3.2, except that the basin is rotated by 300 in the coun
terclockwise direction.
Figure 3.4: An irregular grid system for rectangular basin with irregular grids
50 cm/sec
  
 . 
Figure 3.5: Same as Figure 3.2, except that the irregular curvilinear grid as shown in
Figure 3.4 is used.
78
3.2 Seiche Test
The second test was designed to test the model's ability to simulate a seiche
(standing wave) in a closed basin. The local acceleration term balances the pressure
term. Neglecting vertical and horizontal diffusion, bottom and sidewall friction,
convective acceleration, and Coriolis acceleration, the equations of motion reduce to
the following onedimensional equations
u 0( (3.3)
at ax
h =a (3.4)
ax at
with boundary conditions u = 0 at x = 0, 1, where u is velocity, ( is surface elevation,
h is water depth, g is gravitational acceleration (= 981 cm/s2), and I is the basin
length. The analytical solutions of the lowest mode for the above two equations are
( = a cos wt cos kx (3.5)
gak
u = sin wt sin kx (3.6)
where a is the wave amplitude, w is the circular frequency, and k is the wave number.
The grid system used here is the same as that shown Figure 3.1. The following values
were used for this test case
a = 15 cm
I = basin length = 62 km = A/2 (3.7)
basin width = 14 km
h = 5m
t= diamond: model
Z).50 0.50
O
S t= T/4
0.00 Y < 0.00
0.50 0.50
t=3T/8
t=T/
1.00 1.00
0.00 0.25 0.50 0.75 1.00
BASIN LENGTH FROM ONE END (xil)
Figure 3.6: Comparison between simulated surface elevation (diamond symbol) and
analytic solution (solid lines) for seiche oscillation in a closed basin (T is seiche period).
SEICHE TEST
S1.00
0.50
0.00
0.50
1.00
1.00
1.00
S
d
0.50
0
._j
w
0.00
0.50
1.00 
0.00
Figure 3.7: Comparison between simulated verticallyaveraged velocities (diamond
symbol) and analytic solution (solid line) for seiche oscillation in a closed basin (T is
seiche period).
0.25 0.50 0.75
BASIN LENGTH FROM ONE END (xll)
20 cm/sec
20 cm/sec
Figure 3.8: Maximum velocities for seiche case.
 e L 
 * ) ) 
 
  
 C *  
t=0
t=T/2
* 4  
C * 4 C 
 
 c4 
c 4 
The wave characteristics are calculated as
k = 27r/A = 5.0671 107/cm
T = A/Vgh= 4.92 hours (3.8)
w = 27r/T = 3.5488 104/s
The above values and a time step of 300 seconds were used in the numerical simulation,
using the verticallyintegrated version of the threedimensional model. Figures 3.6
and 3.7 show the comparison between model results and analytic solution after the
model has been run for ten tidal cycles. Figure 3.6 shows the surface elevations
distribution along the channel at increments of one eighth tidal period. Figure 3.7
shows the velocity distribution at corresponding times. The model predictions match
the analytic solution very well with a maximum error of 3 % in amplitude of both
the elevation and the velocity. Figure 3.8 shows the instantaneous velocity fields at
two instants of time when the velocity is maximum. The maximum velocity occurs
at the center and zero velocity occurs at both closed ends.
3.3 Tidal Forcing
One of the most common applications of an estuarine hydrodynamic model is
tidal simulation. Therefore, before applying the model to a real field, the model result
should be compared with the analytic solution describing the response of a rectangu
lar estuary to forced oscillation at the open boundary. Lynch and Gray (1978) derived
the analytic solutions for a tidally forced estuary with a flat bottom and a sloping
bottom. Neglecting nonlinear, Coriolis, horizontal diffusion, and ydirection terms
and assuming linear bottom friction, one obtains the following verticallyaveraged
equations
9( 9hu
+ = 0 (3.9)
9t + =0
