UFL/COEL 2000/007
USE OF THREEDIMENSIONAL 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 THREEDIMENSIONAL 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 ULTIMATEQUICKEST 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 ULTIMATEQUICKEST Algorithm ....................... 34
4.4 CH3D with ULTIMATEQUICKEST 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 Twodimensional 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 Threedimensional Test ................ ....................... 68
5.3.1 Test 5: Advection of a Cubic Block in 3D Uniform Flow ...... 68
5.4 ULTIMATEQUICKEST 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 9Segment 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 BoundaryFitted 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 Onedimensional 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 eightsegment scheme by Sheng et al.............. 16
3.2 The Indian River Lagoon 18segment 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
ninesegment scheme. ............................................ 23
4.1 Control volume stencil for QUICKEST scheme......................... 29
4.2 Local monotonicity criteria and universal limiters ....................... 31
5.1 The 100x100 nonorthogonal 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
ULTIMATEQUICKEST 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. 2D projection along the x axis through the center of the analytical
solution. The ULTIMATEQUICKEST scheme was used.................. 51
5.7 Comparison of numerical and analytical solutions for the advection of a cubic hill in
uniform flow. 2D projection along the y axis through the center of the analytical
solution. The ULTIMATEQUICKEST 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
ULTIMATEQUICKEST 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. 2D projection along the x axis through the center of the analytical
solution. The ULTIMATEQUICKEST scheme was used ................. 57
5.13 Comparison of numerical and analytical solutions for the advection of a cone hill in
uniform flow. 2D projection along the y axis through the center of the analytical
solution. The ULTIMATEQUICKEST 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
ULTIMATEQUICKEST 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. 2D projection along the x axis through the center of the analytical
solution. The ULTIMATEQUICKEST scheme was used ................. 63
5.19 Comparison of numerical and analytical solutions for the advection of a Gaussian hill
in uniform flow. 2D projection along the y axis through the center of the analytical
solution. The ULTIMATEQUICKEST 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
ULTIMATEQUICKEST 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. 2D projection along the x axis through the center of the analytical
solution. The ULTIMATEQUICKEST scheme was used ................. 68
5.25 Comparison of numerical and analytical solutions for the advection of a column in
circular flow. 2D projection along the y axis through the center of the analytical
solution. The ULTIMATEQUICKEST 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 ULTIMATEQUICKEST 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 nonsplitting ULTIMATEQUICKEST 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 ULTIMATEQUICKEST 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 8721147 (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 8721456 (Titusville Brewer Causeway) from September
1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. ........... 86
6.7 Wind data at FDEP station 8721789 (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 8722213 (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 ULTIMATEQUICKEST 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 ULTIMATEQUICKEST 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 ULTIMATEQUICKEST 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: ULTIMATEQUICKEST advection scheme with no diffusion, A=0,
Test 2: ULTIMATEQUICKEST advection scheme with moderate diffusion, Ai=1m2/s,
Test 3: Upwind advection scheme with moderate diffusion, A=lm2/s,
Test 4: ULTIMATEQUICKEST advection scheme with higher diffusion, A=5m2/s.
..............................................................99
7.1 Maps of Scheme A, an 8segment scheme and Scheme B, an 18segment scheme.
................................................ ............102
7.2 Indian River Lagoon 9segment 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 oneyear simulation using splitting ULTIMATEQUICKEST
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 oneyear simulation using
ULTIMATEQUICKEST 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 (19431992) within each of the
18 segments according to SJRWMD ............................... 117
7.10 Contour map of R50 value within each of the nine segments (ninetracer 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 THREEDIMENSIONAL 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 seagrass loss. For
example, in the Indian River Lagoon (IRL), the average phosphorus and nitrogen loadings
have been doubled, and the total seagrass coverage has decreased about 50% during the past
50 years. In order to improve the IRL water quality, restore the seagrass 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 IRLPLR (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
(ULTIMATEQUICKEST) for the modeling of tracer advection process. The second
component is the application of the threedimensional hydrodynamic numerical model CH3D
with the ULTIMATEQUICKEST scheme to the Indian River Lagoon flushing and
segmentation study. A boundaryfitted nonorthogonal 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 nonorthogonal
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 ninesegment scheme has
been proposed. Flushing simulations based on the entire Indian River Lagoon and each of
the nine segments were carried out separately using 1tracer method and 9tracer 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
Seagrass 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 photosynthetically active radiation (PAR) through the water column is the
primary determinant of the depth limitation of seagrass 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 seagrass 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 nontidal 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 (IRLPLR) 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 3D
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 thirdorder accurate and oscillationfree 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 verticallyintegrated 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 seagrass
perspective (Virstein et al. 1994). A ninesegment scheme has been provided and more
accurate flushing calculations for the entire lagoon and each individual segment will be
conducted by using the fully threedimensional 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 ULTIMATEQUICKEST advection scheme (Leonard 1991) to the three
dimensional nonorthogonal curvilinear coordinate system.
* Verify the scheme in nonorthogonal curvilinear coordinate system by conducting
various numerical tests.
* Develop CH3D flushing model with the ULTIMATEQUICKEST 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 ULTIMATEQUICKEST 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 semienclosed coastal
body of water which has a free connection with the open sea and within which sea water is
measurably diluted with fresh water derived from land drainage." 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 shoreparallel."
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:
(SVS)
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 onedimensionally 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+,+(1a)VnCL; (n>2) (2.4)
At low water,
V,,CL=(a V,+R)CI+ [(1a)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 EquationBased 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 backflushing from the mouth based on the following
equation,
Vb=P(1), (2.8)
35
where, Vb is the fresh water volume reentered to the estuary by backflushing, 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 nontidal 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 crosssection 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 backflushing; (5) These
models are onedimensional 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 NavierStokes
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 advectiondiffusion 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 onedimensional, two
dimensional or threedimensional model can be used. The flushing results are very sensitive
to the numerical schemes used for solving the advectiondiffusion 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 oscillationfree advection scheme is needed for
flushing calculation, especially in those regions where the flows are advectiondominated.
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 0D instant and not not
Estuarine cycle ("box" complete considered considered
volumes, Tidal model) in estuary
prism
Fraction of Harmonic tide, 1 tidal 0D 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 1D 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 0D 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 3D 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 24 km wide) and shallow (typically 13
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 semidiurnal 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 intercoastal 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 eightsegment 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 56 free 67 free 78 free 
Boundary shore shore shore highway highway connection connection connection
South 405 3 516 56 free 67 free 78 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 (IRLPLR) 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 onedimensional 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 semidiurnal 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 onedimensional flow between the segments and then
solved a onedimensional advectiondiffusion 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 subbasins
of Indian River Lagoon 50% renewal times were about 5 and 12 days respectively while the
northern subbasin 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 onedimensional lagoonwide model and a three
dimensional lagoon wide model to calculate the flushing in the Indian River Lagoon. They
found that the onedimensional 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 threedimensional curvilineargrid hydrodynamic model CH3D
was used to calculate the IRL flushing. Because of lack of data, they assumed semidiurnal
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 eightsegment 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 subbasin
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 eightsegment (Sheng et al. 1995) and eighteen
segment schemes (SJRWMD 1994). The ninesegment 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 18segment scheme by SJRWMD.
600000
iiE, 1.
0)
C
0
CM
I
2 "
zE
010
C0) ) *
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.
CJ
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 threedimensional 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):
au+ + 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 verticallystretched and horizontally boundaryfitted 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, thirdorder 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 socalled "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 onedimensional and twodimensional 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 twodimensional 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 ULTIMATEQUICKEST 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 costeffective thirdorder upwind (QUICKEST) ULTIMATE scheme gives
excellent practical results (Leonard 1991). The following sections will only present the
formulation of the ULTIMATEQUICKEST scheme for the threedimensional non
orthogonal curvilinearcoordinate numerical model CH3D.
4.3.1 Third Order Transient Interpolation Model
Considering unsteady onedimensional pure advection of a scalar 4 at constant
velocity u, the equation is
a. +u = 0. (4.9)
at ax
Applying centerspace 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 wellknown schemes is the LaxWendroff
scheme. The right face value for the LW 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 LaxWendroff value and an upwind biased
curvature (Leonard 1991),
(3)( =(2)(0 (1c2) 1+SGN(c)( n 2 ) 1 SGN(c) (n +1n)]. (4.15)
r(0=< (rl+22j+ 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 4i2 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 rc 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.
i1 i i+l i1 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 lefthand side inequality in equation (4.18) is less restrictive than
4b<'n+l. (4.19)
For the righthand 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 nonuniform 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 Onedimensional 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 ULTIMATEQUICKEST Algorithm
Based on the knowledge of the third order transient interpolation model and the
universal limiters, the ULTIMATEQUICKEST 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= Id24 + nI for each
face. If ACURV>ADEL(nonmonotonic), 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 ULTIMATEQUICKEST 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
nondimensionalized 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+ cnl c b
i:+1fRn+l r[
Rbz 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)
r2 )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 onedimensional limiter method
described above can be directly applied to threedimensional control volume by taking three
direction sweeps ((, rl, a) separately.
First considering t sweep only, we take the first and second terms on the righthand
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 onedimensional 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:
(1c) 1+SGN(c) 1SGN(c,)
S2 2n n
6 2 2 (4.35)
in 5 direction,
(j) n n n
YW>=<~(4j+1 ; Cy (;+i
(1c2) 1+SGN( n 1SGN(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 onedimensional
ULTIMATEQUICKEST scheme to the threedimensional CH3D. Simply implement three
sweeps in E, rl, and a directions consecutively while updating the central value within each
sweep. This method solves only onedimensional advection equation at every intermediate
time step, and can make better use of the onedimensional limiters. The advantage of this
method is its accuracy and its applicability on relative large Courant number conditions
(C < 1).
Nonsplitting method. The basic procedure is similar but with only one sweep. 1)
Apply onedimensional 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 threedimensional advection
equation at a time using the onedimensional 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 onedimensional pure advection tests which showed
the advantages of the ULTIMATE algorithm against other schemes. This chapter will focus
on two and threedimensional tests in a nonorthogonal 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 twodimensional uniform flow, advection of a column in twodimensional
circular flow, and advection of a cubic block in threedimensional uniform flow) were
conducted. Results of the splitting and nonsplitting ULTIMATEQUICKEST 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 gridlines are established in order to test the new scheme in the nonorthogonal
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+(51i)/20], 32i<51;
xi=x51+&x(i51), 52 ik 101, (5.1)
(2) Along the northern boundary (Y=Y ):
x,=Ax(i1), 15 i 31;
x,=x,_1 +ax[0.5+(i32)/20], 32 i 51;
xi=x51+Ax(i51), 52< i101, (5.2)
(3) Along the western boundary (X=0):
yj=Ay( 1), 1j 3 1;
yj=yj_ +Ay[0.5+(5lj)/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+(f032)/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 twodimensional and one threedimensional bench mark tests are selected to
verify the new scheme.
25 5000 75 00 10000 X (m)
Figure 5.1 The 100x100 nonorthogonal curvilinear test grid system.
5.2 Twodimensional 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, Ixxol\dx and yy0
( =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 nonsplitting
ULTIMATEQUICKEST) 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 ULTIMATEQUICKEST 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
ULTIMATEQUICKEST 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. 2D 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. 2D 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, (xxf)2 +(yyFo) 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
nonsplitting ULTIMATEQUICKEST) 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
ULTIMATEQUICKEST 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
ULTIMATEQUICKEST 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. 2D 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. 2D 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
nonsplitting ULTIMATEQUICKEST) 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 ULTIMATEQUICKEST 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
202
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 ULTIMATEQUICKEST 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. 2D
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. 2D
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+(Yyo0)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(yy,),
v=((xxc). Three numerical schemes (the upwind, QUICKEST and nonsplitting
ULTIMATEQUICKEST) 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 nonuniform flow. Again the test
results show that the ULTIMATEQUICKEST 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
ULTIMATEQUICKEST 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. 2D 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. 2D projection
along the y axis through the center of the
analytical solution. The ULTIMATE
QUICKEST scheme was used.
5.3 Threedimensional Test
5.3.1 Test 5: Advection of a Cubic Block in 3D Uniform Flow
The initial concentration distribution for this test is:
f(=1.0, Ixxol dx, lyyo01 dy, zzo 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 nonsplitting ULTIMATEQUICKEST) 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 ULTIMATEQUICKEST 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
ULTIMATEQUICKEST schemes. A detailed analysis, however, does show significant
advantage of the ULTIMATEQUICKEST 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 ULTIMATEQUICKEST 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 ULTIMATEQUICKEST 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 MassTotal 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 ULTIMATEQUICKEST Scheme With Splitting
The bench mark tests using the nonsplitting ULTIMATEQUICKEST 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 nonsplitting 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 nonsplitting 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 nonsplitting ULTIMATEQUICKEST 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 nonsplitting ULTIMATEQUICKEST 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 ULTIMATEQUICKEST method. (Cx)max=(Cy)max=0.45.
Based on these tests, It seems that there is little difference between the results of
splitting and nonsplitting 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 nonsplitting method.
5.5 Summary of Tests
This chapter has described several numerical test problems: (1) Advection of a cubic
hill in 2D uniform flow; (2) Advection of a cone hill in 2D uniform flow; (3) Advection
of a Gaussian hill in 2D uniform flow; (4) Advection of a column in 2D circular flow; and
(5) Advection of a cubic block in 3D 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) ULTIMATEQUICKEST schemes were compared. The test grid is
specially designed for verification of the schemes based on the nonorthogonal 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 ULTIMATEQUICKEST scheme seems to have the best performance. For the
ULTIMATEQUICKEST scheme, both the splitting and nonsplitting method were
implemented. Further tests with different Courant numbers show that both the Splitting and
nonsplitting 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 nonsplitting 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 ULTIMATEQUICKEST 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 semidiurnal 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 ULTIMATEQUICKEST 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 boundaryfitted grid is necessary. By applying
Thompson's (1983) grid generation method, Sheng et al (1995) generated a nonorthogonal
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 northsouth 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 onedimensional 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 oneyear 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 hourlyaveraged
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 River1 269 3
6 Eau Gallie River2 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 169173 21
2 Northern Indian River 176 199 20
3 203 212 20
4 215 234 20
5 Banana.River 177 199 28
6 204213 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
8721456'.
3.15E+06
SFDEP 8721789
S3.1 E+06
3.05E+06
FDEP 8722213
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 8721 47 (NorthSouth)
10
5
5o
10
265 0 100 200 300
Julian Day 1997 Julian Day 1998
Figure 6.4 Wind data at FDEP station 8721147 (Ponce Inlet) from September 1 (Julian Day
244), 1997 December 31 (Julian Day 365), 1998.
station 8721147x (EastWest) I
I / I FFDEPE
