Use of three-dimensional hydrodynamic model and accurate advection scheme for modeling flushing dynamics and developing segmentation schemes in Indian River Lagoon

UFL/COEL -2000/007

USE OF THREE-DIMENSIONAL HYDRODYNAMIC
MODEL AND ACCURATE ADVECTION SCHEME FOR
MODELING FLUSHING DYNAMICS AND DEVELOPING
SEGMENTATION SCHEMES IN INDIAN RIVER LAGOON

by

Haifeng Du

Thesis

2000

USE OF THREE-DIMENSIONAL HYDRODYNAMIC MODEL AND
ACCURATE ADVECTION SCHEME FOR MODELING FLUSHING DYNAMICS
AND
DEVELOPING SEGMENTATION SCHEMES IN INDIAN RIVER LAGOON

By

HAIFENG DU

A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE

UNIVERSITY OF FLORIDA

2000

Copyright 2000

by

Haifeng Du

ACKNOWLEDGMENTS

First, I would like to express my gratitude to my advisor, Dr. Sheng, for his direction

and financial support in making this study a success. I would also want to thank the other two

members of my committee, Dr. Dean and Dr. Thieke, for carefully reviewing this thesis.

My appreciation is extended to the St. Johns River Water Management District. As

part of the Indian River Lagoon Pollutant Load Reduction Model Development Project, this

study could not be fulfilled without the funding from them.

I would like to thank all the colleagues in room 429: Justin, Detong and Chenxia, for

their help and assistance on the model application, developing and debugging; Dave, for

being a good buddy on field work; Vadim and Jack, for much advice on making the slides;

Kijin and Jun, for being good studying and working partners.

Special thanks go to Sidney, Vik, Vernon, J. J., and Chuck, all the coastal lab

specialists, for their guidance and assistance with field work. Grateful thanks are also

extended to Becky, Lucy, Sandra, Helen and Doretha, for the administrative work and help.

Finally I would like to thank my parents, my wife and my daughter, for their love and

support throughout my graduate research.

TABLE OF CONTENTS

ACKNOW LEDGMENTS ................................................ 111

LIST OF TABLES ................................................... vii

LIST OF FIGURES .................................................... viii

ABSTRACT ..................................................... ..... xiv

CHAPTERS

1 INTRODUCTION ................................................... 1
1.1 Background ................................................1
1.2 Scope of The Study ............................................ 2

2 FLUSHING IN ESTUARIES ......................................... 5
2.1 Definitions and Classifications .................................. 5
2.2 Review of Estuarine Flushing Methods .............................. 6
2.2.1 Tidal Prism M ethod .................................. 6
2.2.2 Fraction of Fresh W ater Method ............................ 7
2.2.3 Modified Tidal Prism Method ............................ 8
2.2.4 Difference Equation Method ........................ ... 9
2.2.5 Hydrodynamic M ethod ................................. 11

3 INDIAN RIVER LAGOON SEGMENTATION .......................... 14
3.1 Describing Indian River Lagoon .................................. 14
3.2 PLRG and PLR Model ......................................... 18
3.3 Previous Indian River Lagoon Segmentation Studies .................. 19

4 NUMERICAL MODEL .............................................. 24
4.1 Governing Equations .......................................... 24
4.2 Review of Advection Schemes ................................. 25
4.3 ULTIMATE-QUICKEST Scheme ............................... 27
4.3.1 Third Order Transient Interpolation Model .................. 27
4.3.2 Universal Limiter ................................... 30

Monotonic Criteria ................ .................. 30
Universal Limiters ............... .................. 31
4.3.3 ULTIMATE-QUICKEST Algorithm ....................... 34
4.4 CH3D with ULTIMATE-QUICKEST Scheme ....................... 34
4.4.1 Dimensionless Advection Equation in CH3D ................ 34
4.4.2 Universal Limiters for CH3D ............................ 36
4.4.3 Third Order Transient Interpolation Method (QUICKEST) for CH3D
.......................................... .........40
4.4.4 Applying ULTIMATE Algorithm to CH3D .................. 41

5 NUMERICAL TESTS OF ADVECTION SCHEMES ....................... 43
5.1 Test Grid System ............................................. 43
5.2 Two-dimensional Tests ........................................ 45
5.2.1 Test 1: Advection of A Cubic Hill in Uniform Flow ........... 45
5.2.2 Test 2: Advection of A Cone Hill in Uniform Flow ............ 51
5.2.3 Test 3: Advection of A Gaussian Hill in Uniform Flow ......... 57
5.2.4 Test 4: Advection of A Column in Circular Flow ............. 63
5.3 Three-dimensional Test ................ ....................... 68
5.3.1 Test 5: Advection of a Cubic Block in 3-D Uniform Flow ...... 68
5.4 ULTIMATE-QUICKEST Scheme With Splitting ..................... 72
5.5 Summary of Tests ................ ........................... 75

6 INDIAN RIVER LAGOON FLUSHING CALCULATION ................... 76
6.1 Flushing Time Calculation Using Simple Methods .................... 76
6.2 IRL Flushing Calculation Using CH3D ........................... 77
6.2.1 IRL Grid System ..................................... 77
6.2.2 Forcing Terms / Boundary Conditions ...................... 79
6.2.3 Method for Calculating Flushing Rate and Time .............. 92
6.3 Sensitivity Study of Flushing Rate and Time Using CH3D .... ........ 93

7 CALCULATION OF INDIAN RIVER LAGOON FLUSHING TIME .......... 101
7.1 Indian River Lagoon Segmentation Schemes ....................... 101
7.2 IRL Flushing Results with 9-Segment Scheme ...................... 105
7.2.1 Hydrodynamic Calibration .............................. 105
7.2.2 Flushing Rate Calculation ............................... 107
7.2.3 50% Renewal Time ................ .................. 112
7.2.4 Flushing Maps ................ ...................... 113
7.2.5 M ass Conservation Check ............................... 115
7.2.6 Discussion .......................................... 116

8 CONCLUSION AND RECOMMENDATION ........................... 120

APPENDICES

A TRANSFORMED EQUATIONS IN CURVILINEAR GRID SYSTEM ........ 124
Coordinate Transformation .................................. 124
Dimensionless Equations in Boundary-Fitted Coordinates .......... 126
Numerical Algorithms of CH3D .............................. 130

B OPEN BOUNDARY INFLOW/OUTFLOW CONDITIONS ................. 131

LIST OF REFERENCES ............................................... 132

BIOGRAPHICAL SKETCH ............................................ 136

LIST OF TABLES

Table page

2.1 Comparison between different estuarine flushing calculation methods .......... 13

3.1 Outline of the general information about Indian River Lagoon ................ 17

4.1 One-dimensional universal limiters ................................... 33

4.2 Universal limiters on faces normal to direction ........................... 37

4.3 Universal limiters on faces normal to r1 direction .......................... 38

4.4 Universal limiters on faces normal to a direction .......................... 39

5.1 Comparisons of numerical schemes for Test 5: The advection of a cubic block in
uniform flow ..................................................71

6.1 (I, J) locations of fresh water input .................................... 81

6.2 50% renewal time at each of the eight segments of IRL ...................... 99

7.1 The boundaries of each of the nine segments of Indian River Lagoon .......... 105

7.2 50% renewal time at each of the nine segments and the entire Indian River Lagoon
...................................................... 113

LIST OF FIGURES

Figure page

3.1 The Indian River Lagoon eight-segment scheme by Sheng et al.............. 16

3.2 The Indian River Lagoon 18-segment scheme by SJRWMD ............... 22

3.3 Changes of seagrass coverage and nutrient loadings in each of the 18 segments by
SJRWMD over the past 50 years (1943 1992). The # are segment numbers for the
nine-segment scheme. ............................................ 23

4.1 Control volume stencil for QUICKEST scheme......................... 29

4.2 Local monotonicity criteria and universal limiters ....................... 31

5.1 The 100x100 non-orthogonal curvilinear test grid system ................. 45

5.2 Analytical solution for the advection of a cubic hill in uniform flow after 100 hours.
The center of the cubic hill is initially located at xo=1.26x104m and yo=0.95x104m,
which is close to the origin. ........................................ 47

5.3 Numerical solution for the advection of a cubic hill in uniform flow using the upwind
scheme after 100 hours. The center of the cubic hill is initially located at
x0=1.26x104m and y0=0.95x104m, which is close to the origin. ............. 48

5.4 Numerical solution for the advection of a cubic hill in uniform flow using the
QUICKEST scheme after 100 hours. The center of the cubic hill is initially located
at xo=1.26x104m and yo=0.95x104m, which is close to the origin. ........... 49

5.5 Numerical solution for the advection of a cubic hill in uniform flow using the
ULTIMATE-QUICKEST scheme after 100 hours. The center of the cubic hill is
initially located at xo=1.26x104m and yo=0.95x104m, which is close to the origin.
..................................................... ........ ..50

5.6 Comparison of numerical and analytical solutions for the advection of a cubic hill in
uniform flow. 2-D projection along the x axis through the center of the analytical

solution. The ULTIMATE-QUICKEST scheme was used.................. 51

5.7 Comparison of numerical and analytical solutions for the advection of a cubic hill in
uniform flow. 2-D projection along the y axis through the center of the analytical
solution. The ULTIMATE-QUICKEST scheme was used ................. 51

5.8 Analytical solution for the advection of a cone hill in uniform flow after 106.7 hours.
The center of the cone hill is initially located at xo=1.26x104m and yo=0.95x104m,
which is close to the origin. ........................................ 53

5.9 Numerical solution for the advection of a cone hill in uniform flow using the upwind
scheme after 106.7 hours. The center of the cone hill is initially located at
xo=1.26x104m and yo=0.95x104m, which is close to the origin. ............ 54

5.10 Numerical solution for the advection of a cone hill in uniform flow using the
QUICKEST scheme after 106.7 hours. The center of the cone hill is initially located
at xo=1.26x104m and yo=0.95x104m, which is close to the origin. ........... 55

5.11 Numerical solution for the advection of a cone hill in uniform flow using the
ULTIMATE-QUICKEST scheme after 106.7 hours. The center of the cone hill is
initially located at x0=1.26x104m and y0=0.95x104m, which is close to the origin.
...... ......................................... . ..56

5.12 Comparison of numerical and analytical solutions for the advection of a cone hill in
uniform flow. 2-D projection along the x axis through the center of the analytical
solution. The ULTIMATE-QUICKEST scheme was used ................. 57

5.13 Comparison of numerical and analytical solutions for the advection of a cone hill in
uniform flow. 2-D projection along the y axis through the center of the analytical
solution. The ULTIMATE-QUICKEST scheme was used ................. 57

5.14 Analytical solution for the advection of a Gaussian hill in uniform flow after 100
hours. The center of the Gaussian hill is initially located at xo=1.62x104m and
yo=1.35x104m, which is close to the origin ............................ 59

5.15 Numerical solution for the advection of a Gaussian hill in uniform flow using the
upwind scheme after 100 hours. The center of the Gaussian hill is initially located at
x0=1.62x104m and y0=1.35x104m, which is close to the origin. ............. 60

5.16 Numerical solution for the advection of a Gaussian hill in uniform flow using the
QUICKEST scheme after 100 hours. The center of the Gaussian hill is initially
located at xo=1.62x104m and yo=1.35x104m, which is close to the origin.
............................................ ................. .. 61

5.17 Numerical solution for the advection of a Gaussian hill in uniform flow using the
ULTIMATE-QUICKEST scheme after 100 hours. The center of the Gaussian hill is
initially located at x0=1.62x104m and yo=1.35x104m, which is close to the origin62

5.18 Comparison of numerical and analytical solutions for the advection of a Gaussian hill
in uniform flow. 2-D projection along the x axis through the center of the analytical
solution. The ULTIMATE-QUICKEST scheme was used ................. 63

5.19 Comparison of numerical and analytical solutions for the advection of a Gaussian hill
in uniform flow. 2-D projection along the y axis through the center of the analytical
solution. The ULTIMATE-QUICKEST scheme was used ................. 63

5.20 Analytical solution for the advection of a column in circular flow after 100 hours.
The center of the column is initially located at x0=7.08x104m and y0=6.9x104m. 64

5.21 Numerical solution for the advection of a column in circular flow using the upwind
scheme after 100 hours. The center of the column is initially located at x=7.08x 104m
and yo=6.9x104m. ................. ................................ 65

5.22 Numerical solution for the advection of a column in circular flow using the
QUICKEST scheme after 100 hours. The center of the column is initially located at
x0=7.08x104m and yo=6.9x104m. .................................. 66

5.23 Numerical solution for the advection of a column in circular flow using the
ULTIMATE-QUICKEST scheme after 100 hours. The center of the column is
initially located at xo=7.08x104m and y=6.9x104m. ...................... 67

5.24 Comparison of numerical and analytical solutions for the advection of a column in
circular flow. 2-D projection along the x axis through the center of the analytical
solution. The ULTIMATE-QUICKEST scheme was used ................. 68

5.25 Comparison of numerical and analytical solutions for the advection of a column in
circular flow. 2-D projection along the y axis through the center of the analytical
solution. The ULTIMATE-QUICKEST scheme was used ................. 68

5.26 Numerical solution after 86.7 hours for the advection of a cubic block in uniform
flow using the upwind scheme. The cubic block on the left front corer is the initial
condition. ...................................................... 69

5.27 Numerical solution after 86.7 hours for the advection of a cubic block in uniform
flow using the QUICKEST scheme. The cubic block on the left front corer is the
initial condition. ................................................. 70

5.28 Numerical solution after 86.7 hours for the advection of a cubic block in uniform
flow using the ULTIMATE-QUICKEST scheme. The cubic block on the left front
corer is the initial condition. ..................................... 70

5.29 Numerical solution for the advection of a cubic hill in uniform flow (Test 1 case)
using the non-splitting ULTIMATE-QUICKEST method. (Cx)max=(Cy)max=O.45.
................................................................. 73

5.30 Numerical solution for the advection of a cubic hill in uniform flow (Test 1 case)
using the splitting ULTIMATE-QUICKEST method. (Cx)max=(Cy)max=O.45.
................................................................. 74

6.1 The Indian River Lagoon grid system ....... ......................... 78

6.2 Tide data (relative to NAVD88) at Ponce, Sebastian, Ft. Pierce and St. Lucie inlets
from September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998.
...................................................82

6.3 Location map of five wind stations.................................... 83

6.4 Wind data at FDEP station 872-1147 (Ponce Inlet) from September 1 (Julian Day
244), 1997 December 31 (Julian Day 365), 1998 ....................... 84

6.5 Wind data at USGS station 02248380 (Haulover Canal) from September 1 (Julian
Day 244), 1997 December 31 (Julian Day 365), 1998 ................... 85

6.6 Wind data at FDEP station 872-1456 (Titusville Brewer Causeway) from September
1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. ........... 86

6.7 Wind data at FDEP station 872-1789 (Banana River, south end) from September 1
(Julian Day 244), 1997 December 31 (Julian Day 365), 1998 ............. 87

6.8 Wind data at FDEP station 872-2213 (Ft. Pierce Inlet) from September 1 (Julian Day
244), 1997 December 31 (Julian Day 365), 1998 ....................... 88

6.9 River discharge data at Eau Gallie River, Crane Creek and Turkey Creek from
September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 89

6.10 River discharge data at Goat Creek, Trout Creek and Fellsmere Canal from
September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 90

6.11 River discharge data at South Prong river, North Canal and Main Canal from

September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 91

6.12 River discharge data at South Canal, North Fork -1 -2 and South Fork -1 -2 from
September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 92

6.13 Relative remaining tracer mass (in percentage of the initial mass) in each of the eight
segments of IRL, obtained using ULTIMATE-QUICKEST advection scheme with
AH= AV=O. ................... ............................ ... 95

6.14 Relative remaining tracer mass (in percentage of the initial mass) in each of the eight
segments of IRL, obtained using ULTIMATE-QUICKEST advection scheme with
moderate diffusion of A,=lm2/. ............................... . 96

6.15 Relative remaining tracer mass (in percentage of the initial mass) in each of the eight
segments of IRL, obtained using the upwind advection scheme with moderate
diffusion of AH=lm2/s ................................................ 97

6.16 Relative remaining tracer mass (in percentage of the initial mass) in each of the eight
segments of IRL, obtained using ULTIMATE-QUICKEST advection scheme with
higher diffusion of AH=5m2/s. ..................................... 98

6.17 Comparison of 50% renewal time at each of the eight segments of IRL:
Test 1: ULTIMATE-QUICKEST advection scheme with no diffusion, A-=0,
Test 2: ULTIMATE-QUICKEST advection scheme with moderate diffusion, Ai=1m2/s,
Test 3: Upwind advection scheme with moderate diffusion, A-=lm2/s,
Test 4: ULTIMATE-QUICKEST advection scheme with higher diffusion, A=5m2/s.
..............................................................99

7.1 Maps of Scheme A, an 8-segment scheme and Scheme B, an 18-segment scheme.
................................................ ............102

7.2 Indian River Lagoon 9-segment scheme. ................. ............ 104

7.3 One year measured and simulated water surface elevation at FDEP Station 872-
21251 (Vero Bridge) from 1997.9.1 ................................106

7.4 One year measured and simulated water surface elevation at USGS station 02248380
(Haulover Canal) from 1997.9.1. ................................... 107

7.5 Flushing results of nine tracer method: Relative remaining mass in each of the nine
segments based on one-year simulation using splitting ULTIMATE-QUICKEST
advection scheme and normal diffusion coefficient (AH =lm2/s). The lower figure is
a closer look during the first 100 days. .............................. 110

7.6 Flushing results of one tracer method: Relative remaining mass in each of the nine
segments and the entire Indian River Lagoon based on one-year simulation using
ULTIMATE-QUICKEST advection scheme and normal diffusion coefficient (AH
=lm2/s). The upper figure includes results in segment 1, 2, 3, 4 and 5; the lower
figure includes results in segment 6, 7, 8, 9 and the entire Indian River Lagoon
(TOTL). .......................................................... 111

7.7 IRL flushing maps: (a) P=0.9, a contour map of 10% renewal time; (b) P=0.75, a
contour map of 25% renewal time. .................................. 114

7.8 IRL flushing maps: (c) P=0.5, a contour map of 50% renewal time; (d) P=0.37, a
contour map of 63% renewal time. ................................ 115

7.9 Change of seagrass coverage during the past 50 years (1943-1992) within each of the
18 segments according to SJRWMD ............................... 117

7.10 Contour map of R50 value within each of the nine segments (nine-tracer method).
...................................................... 118

7.11 Interpolated map of the Indian River Lagoon bottom sediment size (D50).
...................................................... 119

A.1 Transformation from z coordinate to a coordinate.
........................................................125

A.2 Transformation from (x, y) coordinates to ( ,ir) coordinates.
........................................................126

Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science

USE OF THREE-DIMENSIONAL HYDRODYNAMIC MODEL AND
ACCURATE ADVECTION SCHEME FOR MODELING FLUSHING DYNAMICS
AND
DEVELOPING SEGMENTATION SCHEMES IN INDIAN RIVER LAGOON

By

Haifeng Du

August 2000

Chairperson: Dr. Y. Peter Sheng
Major Department: Civil and Coastal Engineering

Recent studies show that many Florida estuaries are facing the problems of excessive

nutrient loading and consequent effects of degenerated water quality and sea-grass loss. For

example, in the Indian River Lagoon (IRL), the average phosphorus and nitrogen loadings

have been doubled, and the total sea-grass coverage has decreased about 50% during the past

50 years. In order to improve the IRL water quality, restore the sea-grass beds and protect the

estuarine ecosystem, a goal of maximum loadings that IRL can sustain (Pollutant Load

Reduction Goal) needs to be established.

As part of the IRL-PLR (Pollutant Load Reduction) Model Development Project

under the supervision of Dr. Y. P. Sheng, this study focuses on two major components. The

first component is the description and the application of an accurate numerical scheme

(ULTIMATE-QUICKEST) for the modeling of tracer advection process. The second

component is the application of the three-dimensional hydrodynamic numerical model CH3D

with the ULTIMATE-QUICKEST scheme to the Indian River Lagoon flushing and

segmentation study. A boundary-fitted non-orthogonal curvilinear grid has been deployed

in order to best represent the complex geometry of Indian River Lagoon. Vertically 0-

stretching technique is employed. Various numerical tests and comparisons have been,

conducted to verify the superior accuracy of the new advection scheme in non-orthogonal

curvilinear grid. Numerical sensitivity studies were also performed to see how different

advection schemes and model parameters could influence the flushing results. Previous IRL

segmentation schemes have been reviewed and compared and a nine-segment scheme has

been proposed. Flushing simulations based on the entire Indian River Lagoon and each of

the nine segments were carried out separately using 1-tracer method and 9-tracer methods.

The simulation results show that the Banana River, which is almost an enclosed water body,

has the lowest flushing rate with R50 value (50% renewal time) equal to about 200 days in

the northern section of the Banana River, and R50 value equal to about 100 days in the

southern section of the Banana River. However, some segments, which have inlets directly

connected to the ocean, have very high flushing rates, e.g., Segments 7 (including Sebastian

Inlet) and 9 (including Ft. Pierce Inlet and St. Lucie Inlet). The R50 values of these two

segments are about 2 days and 11 days, respectively. The flushing rates of all other segments,

e.g., Segments 4 and 5 (the northern Indian River), Segment 8 (Vero Beach area) and

Segment 1 (Mosquito Lagoon), are more or less a month. The total tracer mass was also

checked to ensure the mass conservation of the model.

CHAPTER 1
INTRODUCTION

1.1 Background

Sea-grass beds provide habitats and nursery areas for estuarine animals, enhance

water quality by removing nutrients and stabilizing sediments, and are one of the most

productive ecosystems on the planet (Zieman 1982). Recent studies showed that the

penetration of photo-synthetically active radiation (PAR) through the water column is the

primary determinant of the depth limitation of sea-grass system (Kenworthy et al. 1991,

Morris and Vimstein 1993). With increasing population and development, many Florida

estuaries such as Indian River Lagoon, are facing the problems of declining sea-grass beds

due to excessive nutrient/pollutant loadings, degrading water quality and increasing light

attenuation. While some of the nutrients/pollutants can be quickly flushed out of the estuary

through inlets, some of them will remain in the estuary for a long time. If the

nutrient/pollutant loading rate is faster than the estuarine natural flushing rate, the estuarine

water quality will degenerate gradually. Hence, in order to improve estuarine water quality,

we need to first understand estuarine hydrodynamics and flushing characteristics.

Generally speaking, flushing is the capability of an estuary or a part of an estuary to

replace its old water mass with new ocean water mass. It is conceptually very simple but

2

complex to calculate in reality. There are many methods that have been used to quantify

estuarine flushing characteristics. Traditionally volumetric driven methods were used to

calculate flushing characteristics (Ketchum 1951, Stommel and Arons 1951, Dyer and Taylor

1973). Despite their simplicity, these methods have several drawbacks and are not reliable

to represent the real world situation. First, most of these methods assume complete mixing.

Second, they consider a single principal tidal constituent but exclude water movements

caused by other tidal constituents and non-tidal forces. Third, estuarine hydrodynamics are

not incorporated. More recently the hydrodynamic driven method (Edinger et al. 1998, Sheng

et al. 1995), which considers the estuarine hydrodynamics, is being used in favor of the

traditional volumetric driven methods.

1.2 Scope of The Study

As part of the Indian River Lagoon Pollutant Load Reduction (IRL-PLR) model

project sponsored by the St. John's River Water Management District (SJRWMD), this study

focuses on the quantification of the Indian River Lagoon flushing characteristics, using a 3-D

hydrodynamic model CH3D (Sheng 1986a, 1986b, 1987, 1989). Here the Indian River

Lagoon (IRL) is a linked system of three water bodies: Mosquito Lagoon, Banana River

Lagoon and Indian River Lagoon, not only the Indian River itself. A relatively fine boundary-

fitted curvilinear grid (477x43x4) for the Indian River Lagoon was used in order to fit the

complex geometry and bathymetry. A third-order accurate and oscillation-free advection

scheme has also been developed for CH3D to calculate the flushing time of a passive tracer.

3

Because of IRL's large size and complex geometry, flushing characteristics are very

different from place to place. Some regions in the lagoon may have good flushing and some

may have very poor flushing. A reasonable segmentation scheme, which reflects the spatial

variation in flushing characteristics in IRL, is necessary. Sheng et al. (1995) calculated the

IRL flushing characteristics by using the vertically-integrated version ofCH3D and provided

a preliminary segmentation scheme with eight segments. However, Sheng et al. (1995) only

considered the flushing due to tidal circulation forced by simple harmonic tides at the tidal

inlets. This segmentation scheme will be compared with that produced from the sea-grass

perspective (Virstein et al. 1994). A nine-segment scheme has been provided and more

accurate flushing calculations for the entire lagoon and each individual segment will be

conducted by using the fully three-dimensional version of CH3D with a highly accurate

advection scheme and considering realistic circulation driven by tide, wind and density over

a one year period.

In summary, the scope of the study includes:

* Apply the ULTIMATE-QUICKEST advection scheme (Leonard 1991) to the three-

dimensional non-orthogonal curvilinear coordinate system.

* Verify the scheme in non-orthogonal curvilinear coordinate system by conducting

various numerical tests.

* Develop CH3D flushing model with the ULTIMATE-QUICKEST advection scheme.

* Apply CH3D flushing model to the study of Indian River Lagoon flushing

characteristics.

* Analyze flushing results and compare different Indian River Lagoon segmentation

4

schemes and propose an optimal segmentation scheme.

Chapters two and three present the background knowledge for this study. The

definition of flushing, different flushing methods, and previous flushing and segmentation

studies of Indian River Lagoon will be reviewed in these two chapters. Chapters four and five

cover the numerical model and test results, including CH3D hydrodynamic model, flushing

model, the ULTIMATE-QUICKEST scheme with CH3D, and scheme verification. Chapters

six and seven present various IRL segmentation schemes and flushing results. Conclusions

and recommendations are presented in chapter eight.

CHAPTER 2
FLUSHING IN ESTUARIES

2.1 Definitions and Classifications

An estuary is where the river meets the sea (Fischer et al. 1979). Pritchard (1967),

based on his study on James River in Virginia, 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." This definition, however,

does not work very well with all kinds of estuaries in the world. In Florida, for example,

most of the estuaries are in the form of coastal lagoons, e.g., Indian River Lagoon, Sarasota

Bay, and Biscayne Bay, etc. Kjerfve (1994) defined the coastal lagoon as, "a shallow coastal

water body separated from the ocean by a barrier, connected at least intermittently to the

ocean by one or more restricted inlets, and usually oriented shore-parallel."

Since almost half of the world's population are now living around estuaries, growing

human activities, (e.g., effluent disposal, dredging of navigational channels, reclamation of

wetlands or open water, agriculture runoff, etc.) are causing adverse impacts on the estuarine

environments.

One of the major concerns nowadays is the impact of pollutant loadings on estuarine

water quality, habitat, and living resources. Understanding estuarine flushing is essential to

6

understanding the impact of pollutant loadings on the estuarine ecosystem. Flushing

characteristics are generally described in terms of "flushing time," which is related to

"residence time" or "detention time." Fischer et al. (1979, p. 274) defined flushing time as

"the mean time that a particle of tracer remains inside an estuary." Dyer (1973, p. 109) and

Officer (1976, p. 155) defined the flushing time as "the time required to replace the existing

fresh water in the estuary at a rate equal to the river discharge." Pritchard (1960, p. 56)

defined the flushing time as the "50% renewal time of new ocean water." Although there is

not a standard definition of "flushing time," it is useful to examine the relative flushing time

in different estuaries or various parts of an estuary.

2.2 Review of Estuarine Flushing Methods

There are many methods to calculate estuarine flushing time. The traditional methods

are all relatively simple volumetricc" methods which do not consider the estuarine

hydrodynamics, while most recent methods are hydrodynamicc" methods, which explicitly

consider estuarine hydrodynamics.

2.2.1 Tidal Prism Method

This method assumes that the sea water entering an estuary during the flood tide is

fully mixed with the estuarine water, and the volume of the sea water and river water

introduced during a tidal cycle is equal to the tidal prism, the volume difference between the

high tide volume and the low tide volume (Dyer 1997). Hence the flushing time t can be

calculated in the following formula:

t= T, (2.1)
P

where: V is the estuarine water volume at high tide, P is the tidal prism and T is the tidal

period. This method has been proved to be not reliable since the flushing time calculated

using this method is always considerably lower than that calculated using other methods due

to the incorrect assumption of complete mixing.

2.2.2 Fraction of Fresh Water Method

This method uses fresh water rather than salt, as a tracer to calculate the flushing

time. According to Fischer et al. (1979), "freshness" is defined as:

(SV-S)
f= (- (2.2)
so

where, So is the salinity of the ocean water andS is the estuarine average salinity. Pure fresh

water has a freshness of one and pure ocean water has a freshness of zero. Assuming the total

volume of water in an estuary is V, the amount of fresh water in an estuary will be V,=fV,

then the flushing time can be described as:

t Qf (2.3)

where, Qf is the river inflow (Fischer et al. 1979). However this method seems to be too

strictly defined to be used in systems that lack river or river discharge is much less significant

8

than the tidal prism. For example, in Indian River Lagoon, most rivers are very small creeks

or canals and some segments even do not have any river inflow, i.e., Qf is very small or even

equals to zero in some segments. According to this method, the flushing time could reach

infinity, which is not true in the real world.

2.2.3 Modified Tidal Prism Method

This method was originally introduced by Ketchum (1951) and has been further

developed by Dyer and Taylor (1973) and Wood (1979). Like the fraction of fresh water

method, this method uses fresh water as the tracer. It divides an estuary into several

segments, 0, 1, 2, .... considering one-dimensionally from the head of the estuary down

to the mouth of the estuary, and assumes that the high tide and low tide are in phase for all

segments. The inner segment (segment 0) is defined based on the assumption that the

intertidal volume (tidal prism) in segment 0, Po, is supplied by the river flow during the flood

tide, i.e., Po = R/2, where R is the river inflow during one tidal cycle. Second, segment 1 is

placed so that Vi = R, where V, is the low tide volume of segment 1. Then the tidal prism in

segment 1, P,, is defined as a portion a of the low tide volume of segment 2, i.e., aV2 = P1,

and subsequently, in general, aV,, = aV, +P, for the rest of the segments. High and low

water concentrations of fresh water CH, CL are defined for each segment. At high water,

(Vn+P,)CH=aV,,+,C+,+(1-a)VnCL; (n>2) (2.4)

At low water,

V,,CL=(a V,+R)CI+ [(1-a)V,-R]CH, (n>2) (2.5)

where, V, is the low tide volume of segment n, P, is the intertidal volume between high tide

and low tide in segment n. These two equations can be solved and CH, CnL can be obtained.

Then the estuarine flushing time t can be calculated in the following formula:

t= H (2.6)
(Vn +P)Cn(

Although this method considers incomplete mixing in each segment by adjusting

parameter a, for the same reason of "Fraction of Fresh Water Method," it seems only

applicable to the flushing calculations of river estuaries but not coastal lagoons, especially

those with relatively small rivers or even no river inflow, e.g., Indian River Lagoon.

2.2.4 Difference Equation Method

More recently, Lowery (1998) developed a "Difference Equation-Based Estuarine

Flushing Model" and applied it to a few Gulf of Mexico estuaries. Focusing on "tidal

flushing" rather than "river flushing," this method incorporates the functionality of the "tidal

prism method," though it still uses fresh water as the tracer. By tracking a single input of

fresh water over several time steps (tidal periods) until it is flushed out of the estuary, the

flushed fresh water volume is defined as:

V_ ,(n.p
V+1 (2.7)
; V +P

10

where, yin"1 is the flushed fresh water volume at time step n+l; Vr" is the remaining fresh

water volume at time step n, V"=VV-^-Vyn; P is the tidal prism, which is the volume

between high tide and low tide; V, is the low tide estuarine volume. Lowery calculated the

Gulf of Mexico's estuaries' low tide volumes and intertidal volumes from NOAA's Physical

Environmental Characterizations Branch's (PECB) unpublished data bases. The long term

fresh water inflows were obtained from United States Geologic Survey (USGS) and Texas

Department of Water Resources' (TDWR) data bases. After this, he then calculated the

estuarine flushing time. Lowery's model provides some improvements over previous

volumetric methods in: 1) considering the estuarine stratification by adjusting the low tide

estuarine volume to 75%, 50% or 25% depending on the extent of the vertical stratification;

and, 2) considering the fresh water back-flushing from the mouth based on the following

equation,

Vb=P(1--), (2.8)
35

where, Vb is the fresh water volume re-entered to the estuary by back-flushing, P is the tidal

prism, So is the salinity at the mouth, the ocean water salinity is assumed as 35 ppt. This

method is simple and can apply to the flushing calculations of coastal lagoons. However, it

does not include the estuarine hydrodynamic process. The narrow focus on tidal flushing also

excludes many other non-tidal flushing features.

2.2.5 Hydrodynamic Method

The volumetric methods are simple and direct for calculating the estuarine flushing

characteristics. One doesn't need to understand the estuarine hydrodynamics before applying

these methods. However, these methods have severe restrictions due to their narrow focus

and unrealistic assumptions: (1) Tidal prism method only focuses on the tidal flushing

without considering fresh water input. It uses one simple harmonic tidal constituent and

considers only one tidal period; Or in terms of fresh water input, the fresh water flushing

method calculates the estuarine flushing without explicit consideration of exchange with

ocean water; (2) Most of these methods use only the tidal prism, estuarine volumes at high

tide and low tide, and cross-section averaged velocities to drive the model; (3) Besides tide

and/or river discharge (fresh water input), none of other major driving forces, such as wind

and density gradients, is considered in these methods; (4) Most of these methods assume

steady and uniform estuarine conditions, and instant and complete mixing in estuarine

sections; (5) Some of the older methods ignore the effect of tidal back-flushing; (5) These

models are one-dimensional models along the head of the estuary to the mouth of the estuary;

(6) Estuarine hydrodynamics are not incorporated in either of these methods. Because of the

narrow focus and unrealistic assumptions, the results based on some of theses methods

become very unreliable. For example, according to the "freshwater flushing rate method,"

it takes 715 days to replace Tampa Bay's fresh water and in contrast, only 5 days are required

to flush Tampa Bay according to the "tidal prism flushing method" (Biggs 1986).

More recently, with the development of more powerful computers, the hydrodynamic

12

method is becoming more popular (Sheng et al. 1995, Edinger et al. 1998). This method

requires a quantitative understanding of the estuarine hydrodynamics under the influence of

multiple forcing. Usually the continuity equation and the Reynold's averaged Navier-Stokes

equations of motion need to be solved in order to obtain the surface elevation and the flow

field. Then the estuarine flushing time (e.g., 50% renewal time) can be determined by

tracking the time history of relative remaining tracer mass contained in any segment of the

estuary (or the entire estuary) which is initially filled with a uniform tracer concentration

(e.g., 100). The conservative advection-diffusion equation needs to be solved along with the

equations of motion until the total tracer mass in the segment (or the entire estuary) drops to

below 50% of its initial total mass. Rather than considering one tidal cycle as the time step

for estimating the flushing time, this method solves the transport equation, the equations of

motion, and the continuity equation simultaneously at a much smaller time step (on the order

of seconds or minutes). At each time step the tracer concentration will be calculated

according to the updated new surface elevation and flow field. The one-dimensional, two-

dimensional or three-dimensional model can be used. The flushing results are very sensitive

to the numerical schemes used for solving the advection-diffusion equation. Traditional first

order upwind scheme usually overestimates the flushing rate due to excessive numerical

diffusion. Certain higher order schemes may give negative concentration because of

numerical dispersion. An accurate and oscillation-free advection scheme is needed for

flushing calculation, especially in those regions where the flows are advection-dominated.

Table 2.1 summarizes all of the methods for calculating flushing time in estuaries.

13

Table 2.1 Comparison between different estuarine flushing calculation methods

Methods Driving Forces Time Dimension Mixing Tidal Back- Estuarine
step Type flushing hydro-
dynamics
Tidal Prism Harmonic tide, 1 tidal 0-D instant and not not
Estuarine cycle ("box" complete considered considered
volumes, Tidal model) in estuary
prism

Fraction of Harmonic tide, 1 tidal 0-D instant and not not
Fresh Water Estuarine cycle ("box" complete in considered considered
volumes, Tidal model) estuary
prism, Fresh
water input

Modified Harmonic tide, 1 tidal 1-D considering not not
Tidal Prism Estuarine cycle (from head incomplete considered considered
volumes, Tidal to mouth) mixing by
prism, Fresh adjusting
water input empirical
parameters

Difference Harmonic tide, 1 tidal 0-D considering considered not
Equation Estuarine cycle ("box" vertical considered
volumes, Tidal model) stratification
prism, Fresh by simple
water input adjustment
of
coefficient

Hydro- Real tide, seconds 1, 2 or 3-D Eddy considered considered
dynamic wind, river or viscosity
discharge, minutes and
density diffusivity
gradient ____ concept ___

CHAPTER 3
INDIAN RIVER LAGOON SEGMENTATION

3.1 Describing Indian River Lagoon

The Indian River Lagoon complex system, which lies on the east coast of Florida,

including Mosquito Lagoon, Banana River Lagoon and Indian River Lagoon (Figure 3.1),

is a long (more than 200 km long), narrow (typically 2-4 km wide) and shallow (typically 1-3

m deep) estuary with four inlets connecting it to the Atlantic Ocean and numerous rivers

and canals discharging fresh water into it. The northern end of Mosquito Lagoon is

connected to the ocean through Ponce de Leon Inlet. The southern part of Indian River

Lagoon is connected to the ocean by Ft. Pierce Inlet and St. Lucie Inlet. Sebastian Inlet

connects the middle part of Indian River Lagoon to the ocean. The southern end of Mosquito

Lagoon is connected to Indian River Lagoon through Haulover Canal. Banana River is an

almost enclosed water body with only two narrow connections with the Indian River Lagoon.

Indian River Lagoon is a micro tidal system with the semi-diurnal M2 tide as the principal

constituent (Smith 1987). Tide is the major force causing flushing for most of the Indian

River Lagoon, except in the northern IRL and Banana River, where tidal action is weaker.

The geometry of the Indian River Lagoon is very complicated because of highly irregular

boundaries, barrier islands, and numerous causeways. Indian River Lagoon has a mean depth

15

of about two meters, with a maximum depth of nearly ten meters and a minimum depth less

than half meter. Dredging of the inter-coastal waterway has added to the complexity of the

bathymetry. Due to the complex geometry and bathymetry, the flushing characteristics are

very different within different segments of the Indian River Lagoon. Segments with strong

tidal action and less constriction usually have higher flushing rates, while segments with

weak tidal action and more constriction have slower flushing rates. Based on the natural

boundaries and the preliminary flushing study (Sheng et al. 1995), Indian River Lagoon was

divided into eight segments (Figure 3.1). Table 3.1 gives an outline of some general

information about Indian River Lagoon and its eight segments.

3.2E+06 -

(2)

3.15E+06 -

(3)

-
S (4)\

3.1 E+06 -(5)

3.05E+06 -

3E+06 '---l--- l --l -
500000 525000 550000 575000 600000
X UTM (m)

Figure 3.1 The Indian River Lagoon eight-segment scheme by Sheng et al..

17

Table 3.1 Outline of the general information about Indian River Lagoon

Segment 1 2 3 4 5 6 7 8 Total

North 405 516 5-6 free 6-7 free 7-8 free -
Boundary shore shore shore highway highway connection connection connection

South 405 3 516 5-6 free 6-7 free 7-8 free -
Boundary shore highway highway highway connection connection connection shore

West
Boundary shore shore shore shore shore shore shore shore

East
Boundary shore shore shore shore shore shore shore shore

Ponce de Sebastian Ft Pierce
Tidal Inlet Leon Inlet --- --- --- --- Inlet --- Inlet
St Lucie
Inlet

Annually Tumbull Sykes Cr.: Horse Cr.: Crane Cr.: Fellsmere North South
Averaged - Cr.: 3.61. 0.60. 0.15; 1.41; Canal: Canal: Canal:
River Big Eau Gallie TurkeyCr.: 3.28; 1.17; 1.28; 53.17
Discharge Founder: River_l + 9.52; South Main North Fork:
(m'/s) 0.04. Eau Gallie Goat Cr.: Prong: Canal: 12.7;
River_2: 0.78; 2.97. 1.86. South Fork:
0.46. Trout Cr.: 12.7.
0.64.

Surface
Area 134.9 163.2 183.2 137.3 48.36 47.72 25.43 153.1 893.21
(km2)

Minimum
Depth 0.20 0.44 0.60 0.67 0.80 0.66 0.71 0.47 0.20
(m) __

Maximum
Depth 9.02 3.79 7.17 6.14 3.70 4.39 4.01 11.14 11.14
(m)

Average
Depth 1.47 1.53 1.69 2.11 2.06 1.62 1.42 1.82 1.71
(m)

Minimum
Water 0.1813 0.2097 0.2859 0.2816 0.0897 0.0663 0.0299 0.2433 1.3877
Volume
(km3)

Maximum
Water 0.2075 0.2765 0.3214 0.3119 0.1067 0.0938 0.0441 0.3245 1.6864
Volume
(km3)

Average
Water 0.1978 0.2498 0.3093 0.2900 0.0995 0.0771 0.0361 0.2786 1.5382
Volume
(km')

Tidal Prism -- --- --- -- --- - 0.1493
(km3) -- 019

3.2 PLRG and PLR Model

In order to improve the Indian River Lagoon water quality and restore the sea grass

bed, Florida Department of Environmental Protection (FDEP), with St. John's River Water

Management District's (SJRWMD's) agreement, requires that numerical Pollution Load

Reduction Goals (PLRG's) be established for the major sources of pollution in Indian River

Lagoon. A PLRG specifies the maximum amount of a pollutant (contaminated sediments,

pesticides, nutrients, and toxic chemicals) that will be permitted to enter into Indian River

Lagoon or a segment of IRL through various tributaries. In order to develop scientifically

defensible PLRG's, an Indian River Lagoon Pollutant Load Reduction (IRL-PLR) model is

being developed by the University of Florida (Dr. Y. P. Sheng is the Principal Investigator),

with funding from the St. Johns River Water Management District (SJRWMD). The project

aims to develop a Pollution Load Reduction (PLR) model which incorporates a circulation

model, a sediment transport model, a water quality model, a light model and a sea grass

model. An important first step is to develop a segmentation scheme and an understanding of

the flushing characteristics. As part of the PLR model development effort, this study will

provide an accurate calculation of flushing characteristics within the entire lagoon and each

of the nine segments determined in this study. The flushing model sensitivity study was also

conducted.

3.3 Previous Indian River Lagoon Segmentation Studies

Smith (1993, 1996) used a one-dimensional lagoon wide model to calculate the

flushing time in Indian River Lagoon. He divided Indian River Lagoon into 16 segments

(excluding Mosquito Lagoon and Banana River) and Banana River into 3 segments. He then

quantified the intertidal volume by using harmonic constituents (principal semi-diurnal and

diurnal) recorded at 31 NOS stations. He multiplied the tidal amplitudes by the surface areas,

and incorporated the phases to represent the tidal volumes and the tidal movements. The

model first estimated the averaged one-dimensional flow between the segments and then

solved a one-dimensional advection-diffusion equation. In his 1996 version model, he added

"Flexible Grid" technique which allows control of the extent to which segment boundaries

move along the longitudinal axis of the lagoon with the flow. If the boundaries can move

freely with the flow, the transport equation is solved with diffusion only, and if the

boundaries are fixed, then the transport equation is solved with advection only. Based on one

year simulation he found that 50% renewal time is about 140 days if transport by advection

(fixed grid), and 50% renewal time is not reached until the end of the simulation if transport

by diffusion (fully flexible grid). The results suggest that it tends to underestimate flushing

if using fully flexible grid and overestimate flushing if using fixed grid. Another series of

simulations using completely flexible grid show that in the southern and central sub-basins

of Indian River Lagoon 50% renewal times were about 5 and 12 days respectively while the

northern sub-basin of Indian River Lagoon and the Banana River reached only about 30%

20

renewal by the end of one year simulation. Although the method he used is comparable to

the traditional modified tidal prism method, his results still can be used as a rough indication

of the IRL flushing characteristics.

Sheng et al. (1993a, 1993b) used a one-dimensional lagoon-wide model and a three-

dimensional lagoon wide model to calculate the flushing in the Indian River Lagoon. They

found that the one-dimensional model tended to overestimate the flushing rate. Further, as

part of the preliminary study for IRLPLR model, Sheng et al. (1995) conducted a preliminary

study of the Indian River Lagoon segmentation scheme. A vertically integrated two-

dimensional version of the three-dimensional curvilinear-grid hydrodynamic model CH3D

was used to calculate the IRL flushing. Because of lack of data, they assumed semi-diurnal

M2 tides with amplitudes of 5 cm, 15 cm and 30 cm at Ponce, Sebastian and Ft. Pierce inlets

respectively. Without considering any other forcing functions, they simulated the flushing

characteristics corresponding to four segmentation schemes. Based on the flushing results

and the hydrodynamic and water quality data from 24 FDEP or USGS monitoring stations,

they found that the eight-segment scheme (Figure 3.1) appeared to be the best scheme.

In 1994, St. Johns River Water Management District (SJRWMD) conducted a spatial

analysis of seagrass changes versus changes in nutrient loadings over the past 50 years

(Virnstein et al. 1994). Based on changes of seagrass coverage and drainage sub-basin

boundaries, they divided Indian River Lagoon into 18 segments (Figure 3.2). Then they

calculated change of sea grass coverage and changes of nutrient loadings for each segment

from 1943 to 1992, and regressed percentage change of seagrass coverage on percentage

changes of nutrient loadings. The results as shown in Figure 3.3 indicate that most of the

21

segments have significant increases in nutrient loadings and corresponding decreases in

seagrass coverage during the past 50 years except segment 15 where the fast flushing rate

appears to have flushed increased nutrient loadings into the ocean.

In this study, simulation of flushing characteristics will be first conducted with eight

segments developed by Sheng et al. (1995), and finally with the nine segments developed

based on the result of comparison between eight-segment (Sheng et al. 1995) and eighteen-

segment schemes (SJRWMD 1994). The nine-segment scheme will be described in

following chapters of this thesis.

3.2E+06-

3.15E+06

E

I--
3.1 E+06

3.05E+06

3E+06
500000

550000
X UTM (m)

Figure 3.2 The Indian River Lagoon 18-segment scheme by SJRWMD.

600000

iiE, 1.

0)
C
0
CM
I
2 "

zE
010
C-0) ) *
LO C
* *
CD I- *_1
*c) 0~ '- Eli*
0 s *E*II

IC

gCo

6^ C
C,

- c
c0
0 .

* C

C
C'

00.11

#2
W, #3

Figure 3.3 Changes of seagrass coverage and nutrient loadings in each of the 18 segments
by SJRWMD over the past 50 years (1943 1992). The # are segment numbers for the nine-
segment scheme.

C-J
0 to
._ l i. I
rI-

i i li i
.UEEIIIIIIEE .

0)

i0
E

I)
a

iii-

CHAPTER 4
NUMERICAL MODEL

4.1 Governing Equations

The governing three-dimensional momentum equations can be derived from

Newton's second law of motion which is based on the concept of conservation of

momentum. The continuity and scalar transport equations can be derived from the mass

conservation criteria. With the assumptions of incompressible flow, vertically hydrostatic

water column, Boussinesq approximation, and turbulent eddy viscosity concept, these

equations were expressed in the following forms (Sheng 1989):

a-u+ + aw = 0(4.1)
ax ay az

au auu +uv+auw=fv_ 1 p a au a au a au
-+ --+--+-(AH) + --(AH-+-(A-) (4.2)
at ax ay az pO ax ax ax ay +y -z vz )

av avu avvw vw 1 9p 8 v a av .v
--+---+--+- = -fu- -- + -(A )+--(A )+--(A ) (4.3)
at ax ay az pa y ax Hax +y H( y az A (z3)

(p4
-= -Pg (4.4)
8z

as us avs aws a as a as a as
+ + = f-(D )+ (D )+-(Dz) (4.5)
at ax ay az ax ax ay ay az az

aT auT avTawTaw aT a aT a aT(
-+--+- +- =-(KH -)+-(KH )+-(KV-) (4.6)
at ax ay az ax ax ay ay az az

S + v+ aw=(KH )+ (K _)+-(K a) (4.7)
at ax ay az ax x ay ay az az

p = p(S,T) (4.8)

where: t is the time, (u, v, w) are the mean velocities in the (x, y, z) directions, f=2 sin( is

the Coriolis parameter where 0 is the rotational speed of the earth and ( is the latitude, p

is the density of water, p is the pressure, S is the salinity, T is the temperature,

A and D, KH are the horizontal turbulent eddy viscosity and diffusivities,

A, and D, Kr are the vertical turbulent eddy viscosity and diffusivities, g=9.8 lm/s', is the

gravitational acceleration, ( is any kind of dissolved conservative species concentration.

Instead of dealing with these equations directly, CH3D solves the dimensionless

equations in a vertically-stretched and horizontally boundary-fitted system. These equations

were derived by Sheng (1986) and are shown in Appendix A.

4.2 Review of Advection Schemes

Numerical simulation of circulation and transport processes is very important for

many engineering and scientific studies in coastal waters, lakes and estuaries. In particular,

when quantifying the transport of pollutants and water quality parameters, the accuracy of

the advection scheme is very crucial.

In the past decades, many numerical schemes have been developed and applied to the

study of tracer advection. Probably the simplest scheme is the first order upwind scheme.

However, fewer people are using it now because of its high numerical diffusion. Another

popular scheme is the Leap Frog scheme which is used in the princeton Ocean Model

(POM). This scheme does not have as much dissipation as the up Wind scheme but suffers

a strong dispersive error. Numerical filter has to be used to suppress the strong oscillation,

thus increasing the diffusion in the combined scheme. More recently QUICKEST scheme

(Leonard 1979) has been applied to CH3D (Sheng 1989), because of its simplicity,

efficiency, third-order accuracy with little dissipation and dispersion.

Many numerical schemes have been developed recently to eliminate spurious

oscillations while maintain high order accuracy. One of the schemes is the TVD (Total

Variation Diminution) scheme which used so-called "flux limiting" technique, e.g., Roe

(1985). Although TVD scheme is totally free from numerical overshoots and undershoots

near discontinuities, some studies showed that the scheme also has many disadvantages such

as "clipping" of local extrema and steepeningg" of (what should be) gentle gradients

(Leonard 1991). Leonard (1991) proposed the concept of "universal limiter" and developed

an algorithm, ULTIMATE, which can be used to arbitrary high order explicit conservative

finite difference schemes. By comparing the most popular 36 advection schemes (including

TVD and ULTIMATE), Cahyono (1993) concluded that the ULTIMATE scheme seemed to

be the most "attractive." He applied these advection schemes to six numerical test cases for

both one-dimensional and two-dimensional flows and found that the ULTIMATE scheme

27

was more general than the other schemes and was easier to apply. Lin et al. (1997) also

applied this scheme in his two-dimensional transport model and simulated the water quality

and sediment transport in Humber Estuary, U. K..

The next section will describe how the ULTIMATE algorithm is incorporated in

CH3D.

4.3 ULTIMATE-QUICKEST Scheme

The ULTIMATE algorithm can be used to explicit conservative finite difference

schemes of any order accuracy. The higher order of accuracy, the more interpolation points

will be involved and the less efficiency of computation. For most large scale and relatively

smooth flows the cost-effective third-order upwind (QUICKEST) ULTIMATE scheme gives

excellent practical results (Leonard 1991). The following sections will only present the

formulation of the ULTIMATE-QUICKEST scheme for the three-dimensional non-

orthogonal curvilinear-coordinate numerical model CH3D.

4.3.1 Third Order Transient Interpolation Model

Considering unsteady one-dimensional pure advection of a scalar 4 at constant

velocity u, the equation is

a. +u = 0. (4.9)
at ax

Applying center-space explicit method, equation (4.9) can be interpreted into the following

finite difference form

Sl= #-c[_ ,- itL (4.10)

where (, is the value at the center of the control volume, the superscript n+1 means the

current time step and n means the previous time step. 4t and 4r are scalar values at the left

and right faces of the control volume, c is the Courant number where

at
c=u-. (4.11)
AX

As long as we can determine the control volume face values, the center value can be updated

by equation (4.10).

Many interpolation algorithms can be used to determine the control volume face

values. If we consider first order accuracy, a popular scheme is the upwind scheme. For both

positive and negative Courant number c, the right face value can be written as

,)(-1[(4p* +(pn)-SGN(c)(( -4)n], (4.12)
2(4.12)

where SGN(c) is the sign of Courant number. The corresponding left face value can be

obtained from equation (4.12) by reducing the index number by 1 on all involved center

values,

t(i)= 0(i- 1). (4.13)

If we consider second order accuracy, one of the well-known schemes is the Lax-Wendroff

scheme. The right face value for the L-W scheme can be written as

2)= [( n.+ )-c(n 1 )1]. (4.14)

For QUICKEST, the third order upwind interpolation scheme (Leonard 1979), the right face

value can be written in terms of the second order Lax-Wendroff value and an upwind biased

curvature (Leonard 1991),

(3)( =(2)(0 (1-c2) 1+SGN(c)( n -2 ) 1 SGN(c) (n +1n)]. (4.15)
r(0=< (rl+2-2j+ 1)..15)
6 2 2

Similarly the left face value can be determined by equation (4.13). Four grid points are

involved for each face value, as seen from equation (4.15). Substitute equation (4.13) and

(4.15) into equation (4.10), it is clear that for each control volume the QUICKEST scheme

involves a total of five grid points. The universal stencil can be seen from Figure 4.1.

I 4i-2 i i C C i+1 4t;2
-C.

Figure 4.1 Control volume stencil for QUICKEST scheme.

In a practical simulation, the advection velocity u may not be constant. Hence the

Courant numbers at the left and right faces of the control volume have to be considered

separately. Thus equation (4.10) can be rewritten in a more general form:

I=r?- [cr r-c i ,' (4.16)

where cr=ur At and c,= t are the Courant numbers at the right and left faces,
AX AX

respectively. Note that ci)p=c(i- 1).

4.3.2 Universal Limiter

The third order QUICKEST scheme gives better results than the first and second

order schemes. However, numerical undershoots and overshoots still occur near sharp

gradients. In order to remove the spurious oscillations at sharp gradients while still

maintaining high order accuracy at mild gradients, Leonard (1991) introduced the concept

of "Universal Limiter" which can be applied to arbitrarily high order transient interpolation

models.

Monotonic Criteria

The universal limiters are constructed based on the locally monotonic conditions.

Here the limiters are designed to suppress the spurious overshoots and undershoots of the

control volume face values rather than the center value in order to maintain the conservation

property. Three nodal values (, (upstream), d(downstream) and ( 4(centered between the

two) are defined for the convenience of considering both positive and negative velocities, as

shown in figure 4.2.

i-1 i i+l i-1 i i+1
!C.V. Uf C.V.

(a) u>0 (b) u<0
(a) u1>O (b) u1

Figure 4.2 Local monotonicity criteria and universal limiters.

If the right face velocity is positive, i.e., uO0, then define = "= n n n= ; if

uf<0, then define ( j=I c+2 = +, + ,i7. Considering positive right face velocity, i.e.,

upO (Figure 4.2 a), for 4n in the locally monotonic range (p n:,, d,, we get the

monotonic condition for control volume face value:

q ,. (4.17)

In order to maintain monotonicity, the updated center value must also be constrained by

n+l Bn+l1 n+l
nu ++ < d (4.18)

Universal Limiters

The universal limiters are established based on the monotonicity criteria. For pure

advection at constant velocity, if u>0 and (~4"<(,~ (Figure 4.2 a), the concentration

profile will move to the right by the advection of the positive velocity field. The upstream

value ,u at n+1 time step will not be greater than the value at n time step, i.e., ,,+1g 1,.

Hence, the left-hand side inequality in equation (4.18) is less restrictive than

4b<'n+l. (4.19)

For the right-hand side of equation (4.18), the monotonic criteria require the right face

value 4f to be limited by

ln+l .
C :g %"

(4.20)

For the same reason, inequality (4.20) is less restricted than

(4.21)

If advection by non-uniform flow is considered, i.e.,

+ = *C -(cAf cA 1),

then inequality (4.19) and equation (4.22) can be combined into the following form:

which can be written in a more restrictive form as:

which can be written in a more restrictive form as:

d)'j

Combine inequity (4.17) and (4.24), the universal limiter on the right face of control volume

(4.22)

(4.23)

(4.24)

can be described as:

A min(cl: nc c)n+ -
( 4 min[ d, --].

min(c ) ,c )+ )-4
The upper limiter for(f is ,per =min[" ], the lower limiter is
Cf
lower= = n. In similar procedure, it is easy to derive the right face limiters for the rest of the

cases. Table 4.1 lists all the four cases and the right face limiters, in which c.is the control

volume right face Courant number, c, is the next left face Courant number, Cr is the next

right face Courant number.

Table 4.1 One-dimensional universal limiters
Right face Local Right face limiters for C
velocity, monotonicity
uf upper limiter, lower limiter,

C
positive, M n i -min(c c c' + n n

negative, min(cr,,c]n n
<0 4 m 4 max[ &
Cf
negative, n max(c_ ,c+)-4c +
"<0 M min[r, 4 ) r
cf

(4.25)

4.3.3 ULTIMATE-QUICKEST Algorithm

Based on the knowledge of the third order transient interpolation model and the

universal limiters, the ULTIMATE-QUICKEST algorithm can be established as following:

* Designate upstream(u), downstream(d), and central(c) nodes on the basis of SGN(uf)

for each control volume face.

Compute DEL= (-_ ADEL=IDELI and ACURV= Id-24 + nI for each

face. If ACURV>ADEL(non-monotonic), set n.= and proceed to the next face.

* Calculate the limiters( (upper, tower) for each face based on Table 4.1 and set up the

QUICKEST face value 4f based on equation (4.15).

S Limit (f by limiters upper and (jower

* Update the control volume center value c according to equation (4.22).

4.4 CH3D with ULTIMATE-QUICKEST Scheme

4.4.1 Dimensionless Advection Equation in CH3D

Considering the advection terms only in equation (A. 16), the dimensionless three-

dimensional scalar advection equation in CH3D (without diffusion terms) is:

1 a(Rt) + rRrr 11 --Y R (H )
+) R [b [ (H~ )+. ( Vo +)] b (0, (4.26)
H at ai K a a- H 9 aa

H tCOR, -XREF XREF XREF o
where: =tCOR, u=u v=v-- ,a
ZREF UREF UREF UREFv (XREF

35
UREF
non-dimensionalized variables, andRb= (XRE- (COR) is Rossby number.
(XREF)-(COR)
The finite difference equation for (4.26) can be interpreted based on the control

volume method as following:

,)--w _
(4 rH4 -+(
At

R, )x +
Jc +

-RbZ
AU

With some mathematical manipulation, we get:

n+l Hc n A- R R xv/
C+ cn-l c b-
i:+1f-Rn+l- r-[

-Rb-z f V '
b z( y.0

Defining:

A~n
H"
RAHC-
n+'
C+
H
TX F=Rbu
H (Jgo))
TYF=Rbv Y-

TZF=Rb6A ,A
zb AG

Hxf =gl)x,
At

TYL=Rb -2- ,

TZL=R At
b Zo

(4.29)

(4.27)

(4.28)

r--2 )y -{
f]

equation (4.28) can be simplified as:

C'--=(RAHC)-.(
[(T~ )x,- (rxL).)_ [(F)-.4, -(z)).4J [(rzF).4, -(7zo).). (4.30)

4.4.2 Universal Limiters for CH3D

According to Leonard & Niknafs (1990), the one-dimensional limiter method

described above can be directly applied to three-dimensional control volume by taking three

direction sweeps ((, rl, a) separately.

First considering t sweep only, we take the first and second terms on the right-hand

side of equation (4.30):

{' =(RAHC)Q -[(TXF) (TXL).4x) (4.31)

It can be seen that equation (4.31) is in a very similar form to equation (4.22). Using the

same approach, we can get the limiters on faces normal tog direction (Table 4.2):

Table 4.2 Universal limiters on faces normal to direction
Right
face Local Right face limiters for 4X
velocity monotonicity
in ( in g upper lower
direction, direction limiter, limiter,
Uf C f Cx,_

positive, {. n
U? 0 ( I) () () m(n:,,' (j),n
U0 xn d Xdc Xd
min[(MXL)-4,(M)TX) ] +(RAHC)4 "-I

positive, max n
u? 0 n C n maxI Xd{
max[(ra).x ,(TXL). ]+(RAHC). n- -n
Xc M
(TXF)

negative, n x n
u, min[(TXR)., ,(TXR), ] -(RAHC)4, +

negative, n
Uf max[(TXR)_ X_ ]-(RAHC)_ (T

Then, considering Tr and a sweeps respectively,

4C= ^ (RAHCQ-) n- Y -(TYL)- ,

(4.32)

and

S1= (RAHQ-* n -[(TZF).,-(TZL), ,

(4.33)

the limiters on faces normal to r and a direction can be obtained in the same procedure

(Table 4.3 and Table 4.4).

Table 4.3 Universal limiters on faces normal to Ti direction
Right Local Right face limiters for y
face monotonicity
velocity in ir
in TI direction upper lower
direction, limiter, limiter,
Vf

positive, n Yd {4dn'
v + Y n non d
min[(TYL).- ,(TYL)y] ]+(R,,AHC).;- ;

max[(TYL)-. ,(TYL).,Y ] + (RAHC). n.

negative, n{
f<0 n nn g >a nYd
min[(TYR) (YR) ]- (RAHC) +

negative, n n

max[(TrR)- ,(TM R)- ()(]- (RAHC). (l + -Y
I(TMI I

Table 4.4 Universal limiters on faces normal to a direction
Right Local Right face limiters for A
face monotonicity
velocity in a
in a direction upper lower
direction, limiter, limiter,
Cf CIZO,.

positive, n
0 zn< zn,, \n zdm n
min[(TTZL) ,(TZL4 ]+(RAHC). "-

positive, maxn,
)>O0 .n n max I Zd
max[(TZL)-4z,(TZL)-4 ]+(RA C). --

r n n n n max, Zd
0 znegae, Z d C C{z
min[(TZR)-z.,(TZR) ]- (RAHC). 4 +
(TZF)

negative, .
f:<0 o Zc4Zd ltd C
max[(TZR)- Z,,(TZR) zJ-(RAHC)." +- .
(IZF) J

Where: in Tables 4.2, 4.3 and 4.4,

TXR=Rbi Axt xk (V)xr TYR =Rb y Y( r
XtA 'AT)r H,(V^

TZR=Rb. zt.
tAO

(4.34)

40
4.4.3 Third Order Transient Interpolation Method (QUICKEST) for CH3D

The one-dimensional QUICKEST interpolation formula (equation 4.15) can be

directly used in CH3D to obtain the control volume face values if considering three direction

sweeps separately. The only difference is that the Courant numbers have to involve the

dimensionless coefficients because the equations in CH3D are dimensionless.

The QUICKEST right face value of control volume in three directions for CH3D can

be expressed as:

(1-c-) 1+SGN(c) 1-SGN(c,)
S2 2n n
6 2 2 (4.35)

in 5 direction,

(j) n n n
YW>=-<~(4j+1 ;-- Cy (;+i
(1-c2) 1+SGN( n 1-SGN(c)
6 2 +-) 2 (4.36)

in ri direction,

and

2 2

) n (k+12n n Z2+ n n)
6 2 2 (437)-SGN(c)
6 2 2 k--)J (4.37)

in a direction, where:

atc (UREF) At 1 At
c =u- u Rbu--
f AT (XREF) (COR) Al A

At (UREF) At 1 At
c =j = v -Rb '
ao (IREF) (COR) aA Gq

Ao (AREF) (COR) Ao Aa (4.38)

are control volume face Courant numbers in "l, and a directions respectively.

4.4.4 Applying ULTIMATE Algorithm to CH3D

Splitting method. This is the most direct method to extend one-dimensional

ULTIMATE-QUICKEST scheme to the three-dimensional CH3D. Simply implement three

sweeps in E, rl, and a directions consecutively while updating the central value within each

sweep. This method solves only one-dimensional advection equation at every intermediate

time step, and can make better use of the one-dimensional limiters. The advantage of this

method is its accuracy and its applicability on relative large Courant number conditions

(C < 1).

Non-splitting method. The basic procedure is similar but with only one sweep. 1)

Apply one-dimensional ULTIMATE strategy described above in E, Tl, ando directions

respectively to find the limited control volume face values; 2) Update the central value

according to equation (4.30). This method requires solving the three-dimensional advection

equation at a time using the one-dimensional limiters. The Courant number condition for this

method is more restrictive than that for splitting method:

42

t 1
u + Ivl+ L (4.39)
Ax Ay AZ

Both of the methods were tested and the results are included in the next chapter.

CHAPTER 5
NUMERICAL TESTS OF ADVECTION SCHEMES

Leonard (1991) performed many one-dimensional pure advection tests which showed

the advantages of the ULTIMATE algorithm against other schemes. This chapter will focus

on two- and three-dimensional tests in a non-orthogonal curvilinear coordinate system. A test

grid was set up and five bench mark tests advectionn of a cubic hill, a cone hill, and a

Gaussian hill in two-dimensional uniform flow, advection of a column in two-dimensional

circular flow, and advection of a cubic block in three-dimensional uniform flow) were

conducted. Results of the splitting and non-splitting ULTIMATE-QUICKEST methods were

also compared with different Courant numbers.

5.1 Test Grid System

The test domain is a rectangular basin with a uniform depth of 6000 m, and is divided

into 100x100 grid cells with i=1 to 101 in X direction andj=l to 101 in Y direction (Figure

5.1). Instead of a uniform grid spacing for the entire area, a variable grid spacing and non-

orthogonal grid-lines are established in order to test the new scheme in the non-orthogonal

curvilinear grid system. The grid points on the four boundaries are defined as follows:

(1) Along the southern boundary (Y=0):
xi =,x(i- 1), 1 i< 31;
xi=x,_l +x[0.5+(51-i)/20], 32i<51;
xi=x51+&x(i-51), 52 ik 101, (5.1)

(2) Along the northern boundary (Y=Y ):
x,=Ax(i-1), 15 i 31;
x,=x,_1 +ax[0.5+(i-32)/20], 32 i 51;
xi=x51+Ax(i-51), 52< i101, (5.2)

(3) Along the western boundary (X=0):
yj=Ay( 1), 1j 3 1;
yj=yj_ +Ay[0.5+(5l-j)/20], 32j'51;
Y,=y51 +tay(-51), 52j 101, (5.3)

(4) Along the eastern boundary (X=X ):
yJ=Ay(f- 1), 1 yj=yjl+Ay[0.5+(f0-32)/20], 32 j=Y51 +Ay(f -51), 52 where Ax=-1200m and Ay=1000m.

Thus, the grid is orthogonal and uniform with x=-1200m and Ay=1000m in the

following four regions: (1) 1
51
orthogonal with a variable Ax between 600m and 1740m and a variable Ay between 500m

and 1450m. The total length L (in the X direction) and width W (in the Y direction) of the

domain are: L=119400m and W=99500m.

Four two-dimensional and one three-dimensional bench mark tests are selected to

verify the new scheme.

25 5000 75 00 10000 X (m)

Figure 5.1 The 100x100 non-orthogonal curvilinear test grid system.

5.2 Two-dimensional Tests

5.2.1 Test 1: Advection of A Cubic Hill in Uniform Flow

The initial concentration distribution of this test is defined as:

(=1.0, Ix-xol\dx and |y-y0 ( =0.0, elsewhere,

(5.5)

46

is located on the lower left corer of the grid domain shown in Figure 5.1, which is close to

the origin (x=0 and y=0).

The flow field is uniform: u=0.18m/s along the x direction, and v=0.15m/s along the

y direction. Three numerical schemes (the upwind, QUICKEST and non-splitting

ULTIMATE-QUICKEST) were tested. Figures 5.2 5.7 give the results after 1800 time

steps with At=200s (100 hours).

Theoretically, the pure advection of any type of profile will give the exact same shape

of the profile as the initial one. As seen in Figure 5.2, the cubic hill should keep its original

shape after 100 hour advection in uniform flow. Numerically it is hard to obtain the exact

same profile. Different numerical schemes give different accuracy. In this test the upwind

scheme gives the lowest accuracy due to its severe numerical diffusion. Although the

QUICKEST scheme keeps the higher order accuracy, it cannot avoid the spurious oscillations

at sharp gradient regions (shown in figure 5.4). Only the ULTIMATE-QUICKEST scheme

achieves the higher order accuracy without any undershooting and overshooting (Figure 5.5).

100 .

x 10

Y(m)

X (m)

Figure 5.2 Analytical solution for the advection of a cubic hill in uniform flow after 100
hours. The center of the cubic hill is initially located at xo=1.26x104m and yo=0.95x104m,
which is close to the origin.

48

100

80

60

S40-

01

6 10
x10 4 8
2 4 x 10
Y (m) X(m)

Figure 5.3 Numerical solution for the advection of a cubic hill in uniform flow using the
upwind scheme after 100 hours. The center of the cubic hill is initially located at
xo=1.26x104m and yo=0.95x104m, which is close to the origin.

49

1004

80 .

S40-

QUICKEST scheme after 100 hours. The center of the cubic hill is initially located at
S20 n.

x 0 4 6 8
2 4 x .

Figure 5.4 Numerical solution for the advection of a cubic hill in uniform flow using the
QUICKEST scheme after 100 hours. The center of the cubic hill is initially located at
xo=1.26x104m and yo=0.95x104m, which is close to the origin.

x104

x 10

Y(m)

X(m)

Figure 5.5 Numerical solution for the advection of a cubic hill in uniform flow using the
ULTIMATE-QUICKEST scheme after 100 hours. The center of the cubic hill is initially
located at x0=1.26x104m and yo=0.95x104m, which is close to the origin.

Figure 5.6 Comparison of numerical and
analytical solutions for the advection of a
cubic hill in uniform flow. 2-D projection
along the x axis through the center of the
analytical solution. The ULTIMATE-
QUICKEST scheme was used.

Figure 5.7 Comparison of numerical and
analytical solutions for the advection of a
cubic hill in uniform flow. 2-D projection
along the y axis through the center of the
analytical solution. The ULTIMATE-
QUICKEST scheme was used.

5.2.2 Test 2: Advection of A Cone Hill in Uniform Flow

The initial concentration distribution of a cone hill is defined as:

= 1.O- (x -xo)2 + (-yo)/Ro, (x-xf)2 +(y-yFo) R;
S= 0.0, elsewhere,

(5.6)

where xo=12600m, yo=9500m and Ro=8000m.

The flow field is the same as in Test 1: u=0.18m/s along the x direction, and

v=0.15m/s along the y direction. Three numerical schemes (the upwind, QUICKEST and

non-splitting ULTIMATE-QUICKEST) were tested. Figures 5.8 5.13 give the results after

1920 time steps with at=200s (about 106.7 hours).

52

This test was designed to show the performance of the schemes on local extremum.

Due to the high diffusion, the profile is severely damped and spreads around if the upwind

scheme is used (Figure 5.9). The QUICKEST scheme maintains the third order accuracy, but

unfortunately it cannot avoid the undershooting at the base of the profile (Figure 5.10). The

ULTIMATE-QUICKEST scheme seems to have the best performance. It removes the

spurious oscillation while maintaining the same order of accuracy as QUICKEST scheme

(Figure 5.11).

80

S60

S40-

20

0

8 8
6 10
04 6
2 4 x 10
2
Y (m) X (m)

Figure 5.8 Analytical solution for the advection of a cone hill in uniform flow after 106.7
hours. The center of the cone hill is initially located at xo=1.26x104m and yo=0.95x104m,
which is close to the origin.

100

80 .

c60
0

20-.

0O

8
6 10

2 2 < 4 x104

Y (m) X (m)

Figure 5.9 Numerical solution for the advection of a cone hill in uniform flow using the
upwind scheme after 106.7 hours. The center of the cone hill is initially located at
x0=1.26x104m and y0=0.95x104m, which is close to the origin.

100 .

80 :

c 60,

0
20
0-

x10 4W. 6 8
2 2 4 x 10

Y(m) X (m)

Figure 5.10 Numerical solution for the advection of a cone hill in uniform flow using the
QUICKEST scheme after 106.7 hours. The center of the cone hill is initially located at
xo=1.26x104m and yo=0.95x104m, which is close to the origin.

56

1 .. . . .* .
100

80

a 60 ,

20,

0O

2 4 x0

Y(m) X (m)

Figure 5.11 Numerical solution for the advection of a cone hill in uniform flow using the
ULTIMATE-QUICKEST scheme after 106.7 hours. The center of the cone hill is initially
located at xo=1.26x104m and yo=0.95x104m, which is close to the origin.

Figure 5.12 Comparison of numerical and
analytical solutions for the advection of a
cone hill in uniform flow. 2-D projection
along the x axis through the center of the
analytical solution. The ULTIMATE-
QUICKEST scheme was used.

Figure 5.13 Comparison of numerical and
analytical solutions for the advection of a
cone hill in uniform flow. 2-D projection
along the y axis through the center of the
analytical solution. The ULTIMATE-
QUICKEST scheme was used.

5.2.3 Test 3: Advection of A Gaussian Hill in Uniform Flow

The initial concentration distribution of a Gaussian hill is defined as:

(X- Xo)2 +-yo)2
4P =exp[- ], (5.7)
2R2

where, xo=16200m, yo=13500m and Ro=4000m.

The flow field is the same as in Tests 1 and 2: u=0.18m/s along the x direction, and

v=0.15m/s along the y direction. Three numerical schemes (the upwind, QUICKEST and

non-splitting ULTIMATE-QUICKEST) were tested. Figures 5.14 5.19 give the results after

1800 time steps with At=200s (100 hours).

58

Still, in this test the upwind scheme is proved to be with high damping. It seems that

for the smoothly changed profile, both the ULTIMATE-QUICKEST and the QUICKEST

schemes have similarly good performance. Sometimes for a very steep peak value like in this

test case, the "clipping of the peak value" is difficult to avoid when using the ULTIMATE-

QUICKEST scheme.

100

80 .

c 60- .

E40

20 .

0O

6 10
x o' 4" 6
2 2 4 x 10

Y(m) X (m)

Figure 5.14 Analytical solution for the advection of a Gaussian hill in uniform flow after 100
hours. The center of the Gaussian hill is initially located at xo=1.62x 104m and yo=1.35x 104m,
which is close to the origin.

60

100. .

80 8

S60

ED 40-

20

8

x 10 4 6 8
2 4 x 10

Y (m) X (m)

Figure 5.15 Numerical solution for the advection of a Gaussian hill in uniform flow using
the upwind scheme after 100 hours. The center of the Gaussian hill is initially located at
xo=1.62x104m and yo=1.35x104m, which is close to the origin.

61

100

80

S60-

20

0O

6 104

Y (m) X(m)

Figure 5.16 Numerical solution for the advection of a Gaussian hill in uniform flow using
the QUICKEST scheme after 100 hours. The center of the Gaussian hill is initially located
at xo=1.62x104m and yo=1.35x104, which is close to the origin.

62

100

80 -

c 60

40
O
20-2

0O

6 10
x10 4 6
2 2 4 x 10
Y (m) X (m)

Figure 5.17 Numerical solution for the advection of a Gaussian hill in uniform flow using
the ULTIMATE-QUICKEST scheme after 100 hours. The center of the Gaussian hill is
initially located at xo=1.62x104m and yo=1.35x104m, which is close to the origin.

Figure 5.18 Comparison of numerical and
analytical solutions for the advection of a
Gaussian hill in uniform flow. 2-D
projection along the x axis through the center
of the analytical solution. The ULTIMATE-
QUICKEST scheme was used.

Figure 5.19 Comparison of numerical and
analytical solutions for the advection of a
Gaussian hill in uniform flow. 2-D
projection along the y axis through the center
of the analytical solution. The ULTIMATE-
QUICKEST scheme was used.

5.2.4 Test 4: Advection of A Column in Circular Flow

The initial concentration distribution in this test is set as:

= 1.0, ( -xo)2+(Y-yo0)2RO;
S=0.0, elsewhere,

(5.8)

where, xo=70800m, Yo=69000m and Ro=10000m.

The flow field is circular around the center at xc=59700m, yc=49750m. The angular

velocity is set to (0=1/59400 1/s, and the Cartesian velocity can be described as: u= -(o(y-y,),

v=((x-xc). Three numerical schemes (the upwind, QUICKEST and non-splitting

ULTIMATE-QUICKEST) were tested. Figures 5.20 5.25 give the results after 1800 time

steps with At=200s (100 hours).

Since all the tests described above are advections in uniform flow. This test was

designed to see the performance of different schemes in non-uniform flow. Again the test

results show that the ULTIMATE-QUICKEST scheme appears to be the best scheme.

Figure 5.20 Analytical solution for the advection of a column in circular flow after 100
hours. The center of the column is initially located at x0=7.08x104m and y0=6.9x104m.

o40J ,

20

0O

6 10
x W 4 6
2 4 xo10

Y (m) X (m)

Figure 5.21 Numerical solution for the advection of a column in circular flow using the
upwind scheme after 100 hours. The center of the column is initially located at xo=7.08x 10'm
and y0=6.9x 104m.

100,

80.

60-

I 40-

0
20-

O-
0.

x 104

8
6 SEiiBB 10
8

4 v 1

2

Y(m) X (m)

Figure 5.22 Numerical solution for the advection of a column in circular flow using the
QUICKEST scheme after 100 hours. The center of the column is initially located at
xo=7.08x104m and yo=6.9x104m.

04

100-

80-

S60-

E 40-
o
20-

0-

x 10

6- Ertr i85 10
4 6
2 4 x 1

Y (m) X (m)

Figure 5.23 Numerical solution for the advection of a column in circular flow using the
ULTIMATE-QUICKEST scheme after 100 hours. The center of the column is initially
located at xo=7.08x104m and y0=6.9x104m.

04

100

80

1 2 3 4 5 6 7 8 9 10 11
X(m) xa

Figure 5.24 Comparison of numerical and
analytical solutions for the advection of a
column in circular flow. 2-D projection
along the x axis through the center of the
analytical solution. The ULTIMATE-
QUICKEST scheme was used.

Figure 5.25 Comparison of numerical and
analytical solutions for the advection of a
column in circular flow. 2-D projection
along the y axis through the center of the
analytical solution. The ULTIMATE-
QUICKEST scheme was used.

5.3 Three-dimensional Test

5.3.1 Test 5: Advection of a Cubic Block in 3-D Uniform Flow

The initial concentration distribution for this test is:

f(=1.0, Ix-xol dx, ly-yo01 dy, z-zo jdz;
(=0.0, elsewhere,

(5.9)

where, x0=12600m, Yo=9500m, z0= -4700m, dx=10800m, dy=8000m and dz=1000m.

The flow field is uniform: u=0.18m/s along the x direction, v=0.15m/s along the y

direction and w=0.006m/s along the z direction. Three numerical schemes (the upwind,

69

QUICKEST and non-splitting ULTIMATE-QUICKEST) were tested. Figures 5.26 5.28

give the results after 1560 time steps with at=200s (about 86.7 hours).

Figure 5.26 Numerical solution after 86.7 hours for the advection of a cubic block in
uniform flow using the upwind scheme. The cubic block on the left front corer is the initial
condition.

-4000

-6000

100000

75000 100000

X(m)

Figure 5.27 Numerical solution after 86.7 hours for the advection of a cubic block in
uniform flow using the QUICKEST scheme. The cubic block on the left front corer is the
initial condition.

Figure 5.28 Numerical solution after 86.7 hours for the advection of a cubic block in
uniform flow using the ULTIMATE-QUICKEST scheme. The cubic block on the left front
corer is the initial condition.

0 -

-2000
N

100000

-6000 25000
0 25000 50000 75000 100000 0
x (m)

0 -

-2000

-4000
100000
50000 4l
-6000 )
0 25000 50000 75000 100000
x (m)

71

The figures show that there is little difference between the QUICKEST and the

ULTIMATE-QUICKEST schemes. A detailed analysis, however, does show significant

advantage of the ULTIMATE-QUICKEST scheme. Table 5.1 shows the analysis of the

maximum percentage of over and under prediction and the error ratio of mass conservation.

There is no overshooting or undershooting for the ULTIMATE-QUICKEST scheme.

Table 5.1 Comparisons of numerical schemes for Test 5: The advection of a cubic block
in uniform flow
Numerical Scheme Erro Errunder Err
1 upwind -44.10% 0.0% 2.85 x106 %
2 QUICKEST 24.8% -7.7% 2.25 x10- %
3 ULTIMATE-QUICKEST 0.0% 0.0% 2.41 x105 %

Where, in the table, the Percentage of Overshooting ErroVe, is defined as:

(max- (max) )
Err = -xl100%. (5.10)

The Percentage of Undershooting Err der is defined as:

Errnder (mn) (i' a100%. (5.11)

The Percentage of Mass Conservation Error Errm,. is defined as:

Err., = Total Final Mass-Total Initial Mass00% (5.12)
Err = x100%. (5.12)
Total Initial Mass

Where, ((m )a and ((mi.)a are the maximum and minimum values of analytical solution, ( )

72

and (4m)n are the maximum and minimum values of numerical solution.

5.4 ULTIMATE-QUICKEST Scheme With Splitting

The bench mark tests using the non-splitting ULTIMATE-QUICKEST method seem

to give us very encouraging results. Further study shows that they are only applicable to very

small Courant numbers (C<0.1). When the Courant number increases, the test profile could

become distorted when using non-splitting method, although for certain cases, e.g. advection

in circular flow (Test 4), the results still maintain good shape. However, for most cases, e.g.,

advection in diagonal uniform flow across the domain (Test 1), the non-splitting method

gives severely distorted shape at a relatively larger Courant number (C>0.1), while the

splitting method still gives very good results at the same condition. Figures 5.29 and 5.30

show the comparison for the advection of a cubic hill in diagonal uniform flow using non-

splitting and splitting methods at Cmax=0.45.

100.

80

case) using the non-splitting ULTIMATE-QUICKEST method. (60=(C .45.

40,

20,
01

6 10

x l1 4 6

2
Y (m) X (m)

Figure 5.29 Numerical solution for the advection of a cubic hill in uniform flow (Test 1
case) using the non-splitting ULTIMATE-QUICKEST method. (Cx)max=(Cy)max=0.45.

60,

e 40

20,

0O

6 10
x4 6
2 4 x 1

Y (m) X(m)

Figure 5.30 Numerical solution for the advection of a cubic hill in uniform flow (Test 1
case) using the splitting ULTIMATE-QUICKEST method. (Cx)max=(Cy)max=0.45.

Based on these tests, It seems that there is little difference between the results of

splitting and non-splitting schemes when the Courant number is very small (C<0.1). As the

I

75

Courant number becomes larger and approaches unity, the results of splitting method are

much better than that of the non-splitting method.

5.5 Summary of Tests

This chapter has described several numerical test problems: (1) Advection of a cubic

hill in 2-D uniform flow; (2) Advection of a cone hill in 2-D uniform flow; (3) Advection

of a Gaussian hill in 2-D uniform flow; (4) Advection of a column in 2-D circular flow; and

(5) Advection of a cubic block in 3-D uniform flow. These test problems represent different

conditions of the real world. For each test problem, three advection schemes: (1) upwind, (2)

QUICKEST, and (3) ULTIMATE-QUICKEST schemes were compared. The test grid is

specially designed for verification of the schemes based on the non-orthogonal curvilinear

coordinate system. For all the five test problems, the upwind scheme has too much diffusion.

The QUICKEST scheme cannot avoid spurious oscillation at sharp gradient regions. Only

the ULTIMATE-QUICKEST scheme seems to have the best performance. For the

ULTIMATE-QUICKEST scheme, both the splitting and non-splitting method were

implemented. Further tests with different Courant numbers show that both the Splitting and

non-splitting methods work well if the Courant number is less than 0.1. However, when the

Courant number is greater than 0.1 and less than unity, the non-splitting method is proved

to be unacceptable due to the severe deformation of the original test profile, while the

splitting method can still keep the same good result no matter what Courant number it is.

Thus the splitting ULTIMATE-QUICKEST scheme will be used in this study.

CHAPTER 6
INDIAN RIVER LAGOON FLUSHING CALCULATION

6.1 Flushing Time Calculation Using Simple Methods

Several simple methods reviewed in Chapter 2 were used to estimate the Indian River

Lagoon flushing time, and their results are presented in the section.

By using the Tidal Prism Method, the entire Indian River Lagoon flushing time can

be determined by Equation (2.1). Based on the CH3D simulation results, the estuarine

maximum volume is about 1.69x 10 m3. The tidal prism is about 0.15x109 m3 (The estuarine

maximum volume minus the estuarine minimum volume). Assuming semi-diurnal tides with

T= 12.4 hours, the flushing time can be obtained as, t = 5.65 days, thus, the R50 value is less

than 2.83 days.

The second method used here is the Fraction of Fresh Water Method. The total river

discharges into the IRL (the runoff discharge is neglected) is estimated as 53.17 m3/s based

on the summation of averaged one year historical data at the major 17 rivers (will be

discussed later). Assuming the average estuarine salinity is 20 ppt and the ocean salinity is

35 ppt, the flushing time can be determined by Equation (2.3): t = 157.3 days. The R50 value

is 78.7 days.

The large difference between results of the Tidal Prism Method and the Fraction of

77

Fresh Water Method indicates that these simple methods are not very useful when studying

such a complex system like Indian River Lagoon. This points to the need of the

hydrodynamic method which uses a hydrodynamic model to calculate the flushing time.

6.2 IRL Flushing Calculation Using CH3D

6.2.1 IRL Grid System

Before we can apply CH3D with the third order ULTIMATE-QUICKEST advection

scheme to the simulation of the flushing characteristics of Indian River Lagoon, an accurate

grid system has to be set up. In order to represent the Indian River Lagoon's complex

geometry and bathymetry, a high resolution boundary-fitted grid is necessary. By applying

Thompson's (1983) grid generation method, Sheng et al (1995) generated a non-orthogonal

curvilinear grid for Indian River Lagoon (Figure 6.1). The grid is of very high resolution with

horizontally 477x43 cells (I=1 to 477 in the north-south direction and J=1 to 43 in the west-

east direction) and vertically 4 layers (K=1 to 4).

78

Ponce de Leon Inlet

3.2E+06

Tumbull Mosquito
Creek Lagoon
Big Flounder

Indian
River

3.15E+06 Banana
River
Sykes Creek

SHorse Creek
Eau Gallle River- 1
Eau Galle River- 2
Crane Creek
3.1 E+06 -Turkey Creek
> Goat Creek
Trout Creek
Sebastian Inlet
S Fellsmere Canal
South Prong

North Canal
Main Canal
South Canal
3.05E+06

Ft. Pierce Inlet

North Fork 1 2
South Fork -1 -2 St. Luice Inlet
3E+06 I I I
500000 550000 600000
X UTM (m)

Figure 6.1 The Indian River Lagoon grid system.

6.2.2 Forcing Terms / Boundary Conditions

The major forcing functions for IRL circulation need to be determined once the grid

is set up. Four major forcing functions are considered in this study: surface elevations, winds,

river discharges and salinity distributions. The first three forcing functions (surface elevation,

wind and river discharge) are given by time series of measured data at selected locations.

More than one year continuous data were prepared (September 1997 to December 1998) in

order to represent the seasonal variation. The fourth forcing function, salinity distribution,

is taken into account by calculating the baroclinic circulation term in the numerical model.

Other less dominant forcing functions, e.g., precipitation, evaporation, surface runoff (except

for Merritt Island) and ground water seepage, etc., are not considered.

At the open boundary, a simplified one-dimensional advection equation is used. At

tidal boundaries, tidal back flushing was considered by calculating inflow and outflow

advection fluxes based on the direction of velocities. At river boundaries, only inflow

condition is considered, because river discharge is always positive. Details of inflow and out

flow conditions can be found in Appendix B.

Water surface elevation data. The Indian River Lagoon estuarine system has four tidal

inlets connecting the system to the Atlantic Ocean. As one boundary condition, the surface

elevation at an inlet is given by a time series of measured data. These data were collected at

four FDEP tide stations at the inlets and leveled relative to NAVD88. The data used in this

study contains water surface elevations every 30 minutes over a period of more than one year

(September 1997 to December 1998). Figure 6.2 shows the time series of surface elevation

80

data at the four inlets over a one-year period (September 1997 December 1998). It is

interesting to see that there is a declination of mean sea level from north to south along the

shore. According to the data, the one year averaged mean sea level at the Ponce de Leon Inlet

is about 11 cm higher than that at the Ft. Pierce Inlet.

Wind data. The free surface boundary condition requires the wind stress at every grid

cell. In this study, the measured wind speed and direction at selected locations are

interpolated over the entire grid. Five FDEP and USGS wind stations are selected and their

locations can be seen in Figure 6.3. Figures 6.4 6.8 show time series of hourly-averaged

wind speed at X (eastern is positive) and Y (northern is positive) directions respectively at

the five stations over a period of more than one year (September 1997 December 1998).

River discharge data. A total of 17 rivers are considered for the system. The locations

of these rivers can be seen in Figure 6.1 and Table 7.1. The river flow data collected every

30 minutes at FDEP and USGS stations are used for the study. The fresh water input at the

downstream side of the flow meter is not included. Besides these 17 rivers, the runoff from

Merritt Island area into the northern Indian River and the Banana River are also considered,

though they are very small compared to the river discharge. Table 6.1 gives a list of (I, J)

locations of all the river and runoff sites. Figures 6.9 6.12 show the time series river

discharge plots of some of the major rivers. Smaller rivers and runoff from Merritt Island are

not shown in these figures.

Table 6.1 (I, J) locations of fresh water input

(I, J) location
I J

1 Tumbull Creek 101 13

2 Big Flounder 109 13

3 Sykes Creek 216 24

4 Horse Creek 258 7

5 Eau Gallie River-1 269 3

6 Eau Gallie River-2 273 9

7 Crane Creek 280 7

8 Turkey Creek 289 5

h 9 Goat Creek 303 7

10 Trout Creek 309 7

11 Fellsmere Canal 317 9

12 South Prong River 325 1

13 North Canal 365 12

14 Main Canal 377 12

15 South Canal 388 12

16 North Fork 439 11

South Fork

450

1 169-173 21

2 Northern Indian River 176 199 20
3 203 212 20

4 215 234 20

5 Banana.River 177 199 28

6 204-213 28

82

50 26Poncede Leon Inlet

> -50
z

2E -1o

265 0 100 200 300
50 |Sebastian Inlet]

1 0
.50

-100

265 0 100 200 300
265 160 200 300

Julian Day 1997 Julian Day 1998

Figure 6.2 Tide data (relative to NAVD88) at Ponce, Sebastian, Ft. Pierce and St. Lucie
inlets from September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998.

3.2E+06

USGS 02248380

FDEP
872-1456'.

3.15E+06

SFDEP 872-1789

S3.1 E+06

3.05E+06
FDEP 872-2213

3E+06 I III
500000 550000 600000
X UTM (m)

Figure 6.3 Location map of five wind stations.

10

-10

265 0 100 200 300

10
FDEP Station 872-1 47 (North-South)

10

5

-5o

-10
265 0 100 200 300
Julian Day 1997 Julian Day 1998

Figure 6.4 Wind data at FDEP station 872-1147 (Ponce Inlet) from September 1 (Julian Day
244), 1997 December 31 (Julian Day 365), 1998.

station 872-1147x (East-West) I

I / I FFDEPE

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
76	html_echo_mainwriter.add_text_to_page	Finished reading and writing the file