Three-dimensional curvilinear-grid modeling of baroclinic circulation and mixing in a partially-mixed estuary

THREE-DIMENSIONAL CURVILINEAR-GRID MODELING OF BAROCLINIC
CIRCULATION AND MIXING IN A PARTIALLY-MIXED ESTUARY

By

JEI-KOOK 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 Boundary-Fitted Coordinate Generation ................ 30

2.6 Transformation of Governing Equations Into Boundary-Fitted Grids 32

2.6.1 Transformation Rules . . . . . 32

2.6.2 Dimensionless Equations in Boundary-Fitted Grids . 33

2.7 Boundary Conditions and Initial Conditions . . .... 35

2.7.1 Vertically Stretched System . . . . .

2.7.2 Vertically Stretched and Laterally Boundary-Fitted System

2.8 Solution Technique ..........................

2.9 Finite Difference Equation . . . . . .

2.9.1 External Mode in (x, y, a) (Two-Time-Level Scheme) .

2.9.2 Internal Mode in (x, y, a) (Two-Time-Level Scheme) ....

2.9.3 External Mode in (7, 7, a) (Two-Time-Level Scheme) .

2.9.4 Internal Mode in (n, r, a) (Two-Time-Level Scheme) .

2.9.5 Salinity Scheme (Two-Time-Level 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 Density-Driven Currents .......................

4

103

3.6 Summ ary ................................

ANALYSIS OF OBSERVED HYDRODYNAMIC DATA FROM JAMES R

3

Water Level Data .....

Tidal Currents ......

Circulation Pattern ..

Slack-Water Salinity Data

RE-

IVER

104

107

107

108

.

.

.

.

.

.

...........

...........

. . . .

. . . .

4.5 Front Dynamics ................... ......... 108

5 APPLICATION OF THREE-DIMENSIONAL 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 Neap-Spring 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 z-baroclinic approaches . . .. 234

6 CONCLUSIONS .... .............. ........... 252

BIBLIOGRAPHY ................................. .256

APPENDICES

A HIGHER-ORDER-TERMS (H.O.T.) IN a-GRID 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 BOUNDARY-FITTED GRIDS 267

D MORE COMPACT MOMENTUM EQUATIONS .............. 269

E FINITE DIFFERENCE EQUATIONS IN (, 7r, a) (THREE-TIME-LEVEL) 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 right-handed 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 second-order closure model (From Sheng,
1983) . . . . . . . . 67

3.1 A non-uniform rectangular grid system for testing the forced sur-
face slope test .................. ........ 73

3.2 Simulated steady-state vertically-averaged 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 vertically-averaged 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 flat-bottom 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 flat-bottom annular
section... ............................ 88

3.13 Maximum ebb and flood velocities in a tidally forced flat-bottom
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 bowtie-shaped 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 longitudinal-vertical salinity at slackwater before flood
on June 19, 1985. .......................... 110

4.5 Observed longitudinal-vertical salinity at slackwater before flood
on July 3, 1985. .. ............... ........ 111

4.6 Observed longitudinal-vertical salinity at slackwater before flood
on July 9, 1985. .................. ........ 112

4.7 Observed longitudinal-vertical 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 mid-depth 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 near-surface tidal currents near the James River Bridge
between June 19 and July 17, 1985. . . . ... 129

5.11 Contours of 29-day residual longitudinal velocity at a transect near
the James River Bridge (along i=65) . . . .. 130

5.12 Simulated longitudinal-vertical salinity at slackwater before flood
on July 3, 1985. ........ .......... ..... ... 142

5.13 Simulated longitudinal-vertical 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 29-day 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 near-surface 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 near-surface
distribution in lower James

5.28 Instantaneous near-surface
distribution in lower James

Instantaneous near-surface
distribution in lower James

5.30 Instantaneous near-surface
distribution in lower James

5.31 Instantaneous near-surface
distribution in lower James

5.32 Instantaneous near-surface
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 29-day simulation at Station 2 . . .

Computed elevation, and velocity and salinity at surface, mid, and
bottom level for 29-day 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 near-surface horizontal velocity and salinity (ppt)
distribution in lower James River at 8 a.m. on July 3, 1985. 170

5.34 Instantaneous near-surface horizontal velocity and salinity (ppt)
distribution in lower James River at 9 a.m. on July 3, 1985. 171

5.35 Instantaneous near-surface horizontal velocity and salinity (ppt)
distribution in lower James River at 10 a.m. on July 3, 1985. 172

5.36 Instantaneous near-surface horizontal velocity and salinity (ppt)
distribution in lower James River at 11 a.m. on July 3, 1985. 173

5.37 Instantaneous near-bottom horizontal velocity and salinity (ppt)
distribution in lower James River at 1 a.m. on July 3, 1985. 174

5.38 Instantaneous near-bottom horizontal velocity and salinity (ppt)
distribution in lower James River at 2 a.m. on July 3, 1985. 175

5.39 Instantaneous near-bottom horizontal velocity and salinity (ppt)
distribution in lower James River at 3 a.m. on July 3, 1985. 176

5.40 Instantaneous near-bottom horizontal velocity and salinity (ppt)
distribution in lower James River at 4 a.m. on July 3, 1985. .. 177

5.41 Instantaneous near-bottom horizontal velocity and salinity (ppt)
distribution in lower James River at 5 a.m. on July 3, 1985. 178

5.42 Instantaneous near-bottom horizontal velocity and salinity (ppt)
distribution in lower James River at 6 a.m. on July 3, 1985. 179

5.43 Instantaneous near-bottom horizontal velocity and salinity (ppt)
distribution in lower James River at 7 a.m. on July 3, 1985. 180

5.44 Instantaneous near-bottom horizontal velocity and salinity (ppt)
distribution in lower James River at 8 a.m. on July 3, 1985. 181

5.45 Instantaneous near-bottom horizontal velocity and salinity (ppt)
distribution in lower James River at 9 a.m. on July 3, 1985. 182

5.46 Instantaneous near-bottom horizontal velocity and salinity (ppt)
distribution in lower James River at 10 a.m. on July 3, 1985. 183

5.47 Instantaneous near-bottom 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 one-grid 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 one-grid 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 one-grid 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 one-grid 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 one-grid 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 one-grid 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 one-grid 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 one-grid 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 one-grid 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 one-grid 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 one-grid away
from the channel in lower James River at 11 a.m. on July 3, 1985. 206

5.70 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid away from the channel in
lower James River at 1 a.m. on July 3, 1985. . . ... 207

5.71 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid away from the channel in
lower James River at 2 a.m. on July 3, 1985. . . ... 208

5.72 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid away from the channel in
lower James River at 3 a.m. on July 3, 1985. . . ... 209

5.73 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid away from the channel in
lower James River at 4 a.m. on July 3, 1985. . . ... 210

5.74 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid away from the channel in
lower James River at 5 a.m. on July 3, 1985. . . ... 211

5.75 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid away from the channel in
lower James River at 6 a.m. on July 3, 1985. . . ... 212

5.76 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid away from the channel in
lower James River at 7 a.m. on July 3, 1985. . . .... 213

5.77 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid away from the channel in
lower James River at 8 a.m. on July 3, 1985. . . ... 214

5.78 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid away from the channel in
lower James River at 9 a.m. on July 3, 1985. . . ... 215

5.79 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid away from the channel in
lower James River at 10 a.m. on July 3, 1985 . . ... 216

5.80 u at surface, w at mid-depth, salinity at surface-, mid-, and bottom-
depth, and elevation along the one-grid 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 one-grid 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 one-grid 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 one-grid away
in lower James River at 3 a.m. on July 3, 1985.

220

(ppt) in (x, z) domain along the one-grid away
in lower James River at 4 a.m. on July 3, 1985. 221

(ppt) in (x,z) domain along the one-grid away
in lower James River at 5 a.m. on July 3, 1985. 222

(ppt) in (x, z) domain along the one-grid away
in lower James River at 6 a.m. on July 3, 1985. 223

(ppt) in (x,z) domain along the one-grid away
in lower James River at 7 a.m. on July 3, 1985. 224

(ppt) in (x,z) domain along the one-grid away
in lower James River at 8 a.m. on July 3, 1985. 225

(ppt) in (x,z) domain along the one-grid away
in lower James River at 9 a.m. on July 3, 1985. 226

(ppt) in (x,z) domain along the one-grid away
in lower James River at 10 a.m. on July 3, 1985. 227

(ppt) in (x, z) domain along the one-grid away
in lower James River at 11 a.m. on July 3, 1985. 228

5.92 Time series of u-velocity and salinity in z-baroclinic with diffusion
at grid point (65,16) ........................

5.93 Time series of u-velocity and salinity in c-baroclinic with diffusion
at grid point (65,16) ........................

5.94 Time series of u-velocity and salinity in a-baroclinic without dif-
fusion at grid point (65,16)......................

5.95 Time series of u-velocity and salinity in z-baroclinic with diffusion
at grid point (65,17) ........................

5.96 Time series of u-velocity and salinity in r-baroclinic with diffusion
at grid point (65,17) ........................

5.97 Time series of u-velocity and salinity in a-baroclinic without dif-
fusion at grid point (65,17) . . . . . .

5.98 Time series of u-velocity and salinity in z-baroclinic with diffusion
at grid point (65,18) ........................

5.99 Time series of u-velocity and salinity in c-baroclinic with diffusion
at grid point (65,18) ........................

5.100 Time series of u-velocity and salinity in c-baroclinic without dif-
fusion at grid point (65,18)......................

236

237

238

239

240

241

242

243

244

5.101 Time series of u-velocity and salinity in z-baroclinic with diffusion
at grid point (65,19). ........................ 245

5.102 Time series of u-velocity and salinity in a-baroclinic with diffusion
at grid point (65,19). ..... .............. ..... .. 246

5.103 Time series of u-velocity and salinity in c-baroclinic without dif-
fusion at grid point (65,19). . . . ..... 247

5.104 Time series of salinity in z-baroclinic with diffusion, a-baroclinic
with diffusion, and c-baroclinic without diffusion at grid point
(65,16) . . . . . . . . 248

5.105 Time series of salinity in z-baroclinic with diffusion, r-baroclinic
with diffusion, and r-baroclinic without diffusion at grid point
(65,17).... ............... ............ 249

5.106 Time series of salinity in z-baroclinic with diffusion, u-baroclinic
with diffusion, and r-baroclinic without diffusion at grid point
(65,18)... ................ ............ 250

5.107 Time series of salinity in z-baroclinic with diffusion, r-baroclinic
with diffusion, and a-baroclinic 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 u-velocity at Sta-
tion A . . . . . . . .. 233

5.12 Model skill parameters for the sensitivity test to v-velocity 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

THREE-DIMENSIONAL CURVILINEAR-GRID MODELING OF BAROCLINIC
CIRCULATION AND MIXING IN A PARTIALLY-MIXED ESTUARY

By

Jei-Kook Choi

December 1992

Chairman: Dr. Y. Peter Sheng
Major Department: Coastal and Oceanographic Engineering

Three-dimensional curvilinear-grid hydrodynamic models for partially-mixed 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 inter-tidal and

longterm dynamics of surface elevation and three-dimensional field of velocity, and

salinity in partially-mixed estuaries with complex geometry and bathymetry.

The momentum and continuity equations in the curvilinear coordinates are solved

on a space-staggered grid system using an implicit finite-difference algorithm for the

"external" vertically-integrated 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 density-driven

current. Effect of spring-to-neap 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 semi-enclosed 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". Salt-wedge 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 "well-mixed" estuary is formed.

Partially-mixed estuaries are intermediate between salt-wedge and well-mixed es-

tuaries. James River is a partially-mixed, 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 1-10 /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 partially-mixed 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 spring-neap tidal

variation, alternate stratification and destratification may be evident (Haas, 1977;

Cerco, 1982). Current-induced turbulent mixing plays an important role in affect-

ing the stratification and destratification in a partially-mixed 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 two-layered flow exists in a partially-mixed 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

partially-mixed estuaries are basically three-dimensional and time-dependent because

of irregular shape of the basin, density-induced 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 time-dependent 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 mean-flow 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 semi-empirical (e.g. various levels of second-order closure model) or ad-hoc

empirical models (e.g. simple eddy-viscosity models). An equilibrium closure model

(Sheng et al., 1989) is included in the three-dimensional 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 two-dimensional laterally averaged model to the

James River to predict intertidal variation of surface level, velocity, and salinity.

This model could simulate the destratification-stratification cycle coincident with the

spring-neap tidal cycle in the James River Estuary. Heltzel and Granat (1988) applied

a two-dimensional vertically-averaged 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

three-dimensional 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 three-dimensional 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 three-dimensional models were developed by dividing the vertical water-

column into layers and solving each layer with a vertically averaged two-dimensional

model (e.g., Laevastu, 1974). Laevastu (1974) introduced a two-layer 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 three-dimensional hydrodynamic model

which can be used in non-orthogonal boundary-fitted grid as well as orthogonal or

Cartesian grid. As in the earlier Cartesian grid model (Sheng, 1983), the curvilinear-

grid model uses similar model-splitting technique, vertical turbulence closure scheme,

modified alternating direction implicit (ADI) scheme for the external mode, and

vertically-implicit 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 three-dimensional boundary-fitted 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 time-dependent and

spatially-varying) 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 finite-difference models use only Cartesian-grid 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 curvilinear-grid

finite-difference models which were discussed above or finite-element models (Wang,

1975; Wang and White, 1976). The finite-element 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 finite-difference models, however, tridiagonal

matrix system generally exists and can be efficiently solved. The curvilinear-grid

finite-difference 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,

near-bottom 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 flood-tidal

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 two-dimensional model based on vorticity

and salt-balance 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 run-off conditions using high-resolution 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

density-induced 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 river-forced estuarine plumes. A numerical

study of this phenomena was done by Chao (1990). In a typical midlatitude estuary-

shelf environment, the river-forced 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 down-stream. 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 time-dependent 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 density-field 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 three-dimensional models (Sheng, 1983;

Oey et al.). Sheng et al. (1989a) simulated the tidal, wind-driven, 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, boundary-fitted 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 two-dimensional

curvilinear-grid 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 three-dimensional model analysis because of stratification, complex ge-

ometry, and topography. Although the three-dimensional boundary-fitted 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 three-dimensional 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), wave-current boundary layer (Sheng, 1984), and vegeta-

tion boundary layer (Sheng, 1982). However, three-dimensional models generally do

not use Reynolds stress models to represent the turbulent mixing, because of the com-

plexity. Instead, simplified second-order closure models have been developed and used

with three-dimensional 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 second-order correlation equations are simplified to algebraic equations. The

1

12
TKE closure model is similar to the Level-3 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 second-order correlations. The equilibrium

closure model is similar to Mellor and Yamada's Level-2 model. Instead of solving

the differential transport equations for A which consists completely of modeled terms,

simplified second-order closure models often solve a set of integral constraints for A

(Sheng and Chiu, 1986).

Results of the simplified second-order 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 sediment-laden 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, spring-to-neap and residual long-term

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 o-grid may lead

13
to numerically-induced 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 three-dimensional

hydrodynamic curvilinear-grid 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

short-term and long-term 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 fall-line 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 three-dimensional

and time-varying 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 mid-channel of approximately 8 m. Maximum depth

in the estuary is 10 15 m and average depth is about 6 m. Cross-sectional 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 vertically-stretched (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 boundary-fitted (, 7, r) coordinate

system are then presented. Following a presentation of boundary conditions, initial

conditions, solution technique, and finite-difference 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 right-handed 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.81ms-2 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 right-handed 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 cm-3. 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 u-transformation 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 a-stretching, 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 a-stretching 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 a-grid 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 z-grid model, which resolves the depth

L

21
with "stair-step" 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 a-grid 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 z-plane and avoiding the use
of higher-order 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 vertically-stretched
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 first-order 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

z-C

+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 x-derivative while applying the coordi-
nate transformation results in

Op apa ap 1 a9C H 9p H(
-Ox -9x'g J [- H(- + j- ]H da + gpo-a (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 + g--a-p + g-z (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- order-terms
(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 non-dimensionalized 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 non-dimensionalization 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 H-2 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,-)

+ E-AH + ) 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 vertically-integrated
velocities. The vertically integrated equations represent the "external mode" of the
3-D 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 Boundary-Fitted Coordinate Generation

The basic idea of the boundary-fitted 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 quasi-linear
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 Boundary-Fitted Grids

2.6.1 Transformation Rules

Generation of a boundary-fitted grid is an essential step in the development of
a boundary-fitted hydrodynamic model. It is, however, only the first step. A more
important step is the transformation of governing equations into the boundary-fitted

coordinates. A straightforward method is to only transform the dependent variables,

i.e., the coordinates, while retaining the Cartesian components of velocities. John-

son (1982) developed such a two-dimensional vertically-integrated 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 Boundary-Fitted 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 tensor-invariant 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 u-equation and v-equation 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 no-slip

condition (V. g = 0, where sis the unit tangential vector to the boundary) is applied.
Application of the above closed and no-slip 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 1-D 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 one-sided differencing scheme. During the inflow,

however, the salinity value at the open boundary can either take on a prescribed value

or be determined from the 1-D advection equation while using the boundary salinity

value and the prescribed salinity value to evaluate the spatial flux term.
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 spin-up 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 Boundary-Fitted 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 no-slip condition dictates that the tangential velocity is zero while the
slip condition requires that the normal velocity is zero. When flow is specified at the
open boundary or river boundary, the normal velocity component is prescribed.
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

round-off error) even if the scheme is stable. For explicit methods, the stability con-

dition is that the CFL (Courant-Friedrichs-Lewy) 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 mode-splitting 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 vertically-averaged 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. Semi-implicit 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 space-staggered 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 u-velocity at the left and right sides of a cell, and v-velocity 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 non-staggered grid system.

For the time marching, forward time difference (two-time-level) or leapfrog technique
(three-time- level) can be used. The present study uses a two-time-level scheme.

Both the leapfrog technique and spatial difference are centered differences, and the

accuracy in space and time is second-order. 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 alternating-direction
implicit (ADI) method. This solution technique consists of a two-step 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 y-derivative 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 i-derivatives
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 vertically-averaged
velocities to produce the three-dimensional 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 two-time-level difference equations in (x, y, a) system

are presented. A set of three-time-level difference equations in (, ri, a) system are

presented in Appendix E. The results presented in this study are obtained with the
two-time-level 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

Finite-difference equations for the external mode and the internal mode are listed

in the following. The basic two-time-level schemes are similar to those presented in

Sheng (1983) and Sheng (1987). For three-time-level schemes, refer to Appendix E.

2.9.1 External Mode in (x, y, a) (Two-Time-Level 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 finite-difference 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 Chezy-Manning coefficients which are used

in many two- dimensional models. When used in two-dimensional mode, however,
the model computes implicitly the bottom stresses in the u and v equations .'

with Chezy-Manning formulation. The dimensionless finite-difference equations in

the following are

AF"+1 = F" + AtD"

(2.57)

43
where superscripts n + 1 and n denote the present and previous time-step of integra-

tion, A t is the time-step, 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 x-directions (the "x-sweep"), and then computing in the y-direction

(the "y-sweep"). Carrying out the factorization of Equation (2.57), and neglecting

terms of order At" yields

x-sweep:

[I + A,]F* = [I A,]F" + AtD" (2.64)

y-sweep:

[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 x-sweep and y-sweep are

x-sweep:

(* + 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 x-direction. The above equation can be

readily solved by simple inversion of tridiagonal matrices along all the grid rows.

y-sweep:

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 x-sweep into the equation in
the y-sweep.
2.9.2 Internal Mode in (x, y, a) (Two-Time-Level Scheme)

Defining deficit velocities as fi = u u and f v u, we obtain the differential
equations for the internal mode by subtracting the vertically-averaged momentum
equations from the three-dimensional 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
vertically-integrated 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 vertically-integrated Equations (2.32) and

(33) must be evaluated by summing the corresponding terms in the three-dimensional
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) (Two-Time-Level Scheme)

The present study generally uses the following two-time-level scheme instead of the
two-time-level scheme described earlier. Vertically-integrated 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 i-direction, and D, and D' are all the terms in the
x-direction and y-direction 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,(v-V+') + (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 + T3-V1

+ (1 T3) 912 V] +Dgn'

ing the above equations directly, however, they can be written in
lows

(* + Tia6A(v-U*) = 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

77-sweep:

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,((n-1) (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 = AFn-1 + D

[I + AJ]F* = [B A,7]Fn-1 + 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)
V9-0-- 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) (Two-Time-Level 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
2A-t 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 (Two-Time-Level 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 a-directions, 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 k-1 S~Sk-1 +

[1 + 7k(a-k_ + 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-
ak-1/2 -

Gk = advection + diffusion

For k = 1, i.e., the lowest a-level:

[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 upper-most a-level:

-n+l [l l (2.117)
i-7kmakm_ Ss i +1 + 7km-Iakm~- 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 k-1/2 1 < K < KM

S 1 + 71a3/2 K = 1
A2(K) = 1 + 7k(akm-l/2 + k+1/2) 1 < K < KM
1 + 7fkm-1/2OCkm-1/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 second-order 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 finite-difference form is

S(HuS ) (HS/^u)+ (HS-.u)_
(HuSqy") =(2.120)
1 r(HS~Jf)i,j,k + (HS-V )i+l,j,k
= A 2 ui+l,i,k
(HSvj)ij,k + (HSv),)i-1,,, 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)ij-1,k
2 i

-(HwS) = wk 1 (Sij,k + Si,j,k+1) Wk-i (j, + i,k-1)] (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 finite-difference
form is

(HuSV~ ) = (S+Ui+lI,k S-ujI,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 Central-Upwind 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 finite-difference form is

a(HuSv ) = (S+uT+l,k S-ui,,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-
Corrected-Transport (FCT) scheme (Zalesak, 1979) to simulate the advection and
diffusion of salinity and suspended sediments by wind-driven 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 z-grid version of CH3D was used. Sheng et

al. (1989) found that, when a a-grid model is used, higher-order advective schemes

should be avoided in regions of sharp bathymetric gradient when insufficient grid
points exist to resolve the sharp gradient.

The higher-order schemes have been primarily applied to solve the advection

and diffusion equations but not to solve the Navier-Stokes equations. Several other

higher-order schemes can be used to improve the accuracy of advective terms. The
application of higher-order 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 first-order derivative term.

The QUICKEST scheme used for this study has third order accuracy. This

method uses cubic upstream-weighted interpolation with five points of i-2, i-1, i, i+
1,and i + 2. Because this is a five-point 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 i-direction,

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 + Si-1,j,k)

ui,j,kAt (Si,j,k Si-,j,k)

11
+ (.,+j,k + -Ij,skl) 2 (Si-1,k + Sij,k)

1- (1 (ui,jkAt)2) 2Si-l,i,k + Si-2,j,k)
1 1 1
i,j,kAt (Si,j,k Si-l,j,k) -jg jHi ,y-- --
2 2 Hij,< -"'jc

and in the q-direction,

(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 + Si-1,,k)

Vi,,kAt (Si,j,k, Si-l,,k)
+ (ijk Ii,j,kl) 2 (Si-l,j,k + Sidj,k)

(1 (vi,j,kAt)2) (Si,k 2Si-i,,k + Si-2,j,k)
1 1 1
Vi,,kAt (Si,j,k Si-ij,k) ~ v Hij,v H
2 (VijjHk +* Vj+J,<

2.9.7 Horizontal Diffusion Terms

The finite-difference form of the horizontal gradient in the horizontal diffusion
terms is the standard second-order 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) Si-1/2j-1/2,k Si-1/2,I-1/2,k)] /A, A1

22 2S 22. Si,j+l,k 2Sij,k + Si,j-l,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

a-0 + 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 ad-hoc 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 second-order 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 second-order 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 Navier-Stokes equations to
produce the primed-variable equations. The resulting equation for the i-component
is multiplied by u', and j- component equation by ut. Adding two equations and

61
averaging by time, gives the u u-equation.

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 a-7 - 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 second-order closure model described in Lewellen (1977) and

62
Sheng (1982), the above equations are replaced by the following equations.

ot + 1-7 = (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) Munk-Anderson type eddy coef-

ficients, (3) a simplified second-order closure model of turbulent transport, i.e., the

equilibrium closure model, and (4) another simplified second-order 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.

Munk-Anderson Type Eddy Coefficients

In non-stratified 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 mid-depth 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

mid-depth.

In stratified flow situations, the effect of buoyancy on eddy coefficients can be

parameterized in terms of several Richardson number-dependent 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 von-Karman constant, and u. is the dimensionless friction velocity at
the free surface.

A Simplified Second-Order 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 second-order 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 second-order 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 second-order 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 second-order 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 cut-off

limitation of turbulent macroscale based on 8q2, the spread of turbulence determined

from the turbulent kinetic energy (q2) profile; and N is the Brunt-Vaisila frequency
defined as

N = p 1/2 (2.176)

Although the simplified second-order closure model presented above is strictly
valid when the turbulence time scale (A/q) is much less than the mean flow time
scale and when turbulence does not change rapidly over A, it has been found to be

quite successful in simulating vertical flow structures in estuarine and coastal waters.

Recent simulations of wind-induced mixing of salinity in Chesapeake Bay (Sheng et
al., 1989a; Johnson et al., 1989) and wind-driven 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 second-order 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 one-dimensional vertical case of the o-grid, the equation is

q2 1 9Huq2 wq2 1 8 Aq
Ot H Ox H =ao H2 ao q a)
O---i-+ + 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., finite-difference

schemes or finite-element schemes for various terms of the governing equations), and

errors associated with the computer and programming (e.g., computer round-off 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

three-dimensional 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 steady-state 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 one-dimensional 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 vertically-averaged 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 non-uniform 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 steady-state vertically-averaged velocities in the x-y 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 x-y 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 non-uniform rectangular grid system for testing the forced surface slope
test.

50 cm/sec

Figure 3.2: Simulated steady-state vertically-averaged 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 side-wall friction,

convective acceleration, and Coriolis acceleration, the equations of motion reduce to

the following one-dimensional 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 vertically-averaged 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 10-7/cm

T = A/Vgh= 4.92 hours (3.8)

w = 27r/T = 3.5488 10-4/s

The above values and a time step of 300 seconds were used in the numerical simulation,

using the vertically-integrated version of the three-dimensional 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 y-direction terms

and assuming linear bottom friction, one obtains the following vertically-averaged
equations

9( 9hu
+ = 0 (3.9)
9t + =0

MILLISECOND	CLASS.METHOD	MESSAGE
0	sobekcm_page_globals.constructor
0	sobekcm_page_globals.constructor	Application State validated or built
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.constructor	Navigation Object created from URI query string
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.display_item	Retrieving item or group information
0	sobekcm_page_globals.get_entire_collection_hierarchy	Retrieving hierarchy information
0	sobekcm_assistant.get_entire_collection_hierarchy
0	cached_data_manager.retrieve_item_aggregation
0	cached_data_manager.retrieve_item_aggregation	Found item aggregation on local cache
0	item_aggregation_builder.get_item_aggregation	Found 'all' item aggregation in cache
0	system.web.ui.page.page_load (ufdc.page_load)
0	sobekcm_page_globals.constructor.on_page_load
0	html_echo_mainwriter.add_style_references	Adding style references to HTML
0	html_echo_mainwriter.add_text_to_page	Reading the text from the file and echoing back to the output stream
62	html_echo_mainwriter.add_text_to_page	Finished reading and writing the file