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

Material Information

Title:
Use of three-dimensional hydrodynamic model and accurate advection scheme for modeling flushing dynamics and developing segmentation schemes in Indian River Lagoon
Series Title:
UFLCOEL-2000007
Creator:
Du, Haifeng, 1968-
Place of Publication:
Gainesville Fla
Publisher:
Coastal & Oceanographic Engineering Program, University of Florida
Publication Date:
Language:
English
Physical Description:
xv, 136 leaves : ill. ; 28 cm.

Subjects

Subjects / Keywords:
Estuaries -- Research -- Florida ( lcsh )
Hydrodynamics -- Mathematical models ( lcsh )
Ocean circulation -- Florida -- Indian River (Lagoon) ( lcsh )
Estuarine pollution -- Florida -- Indian River (Lagoon) ( lcsh )
Sediment transport -- Florida -- Indian River (Lagoon) ( lcsh )
Indian River (Fla. : Lagoon) ( lcsh )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (M.S.)--University of Florida, 2000.
Bibliography:
Includes bibliographical references (leaves 132-135).
Statement of Responsibility:
by Haifeng Du.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. §107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact Digital Services (UFDC@uflib.ufl.edu) with any additional information they can provide.
Resource Identifier:
49575825 ( OCLC )

Full Text
UFL/COEL -2000/007

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

2000




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

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

2000




Copyright 2000
by
Haifeng Du




ACKNOWLEDGMENTS

First, I would like to express my gratitude to my advisor, Dr. Sheng, for his direction and financial support in making this study a success. I would also want to thank the other two members of my committee, Dr. Dean and Dr. Thieke, for carefully reviewing this thesis.
My appreciation is extended to the St. Johns River Water Management District. As part of the Indian River Lagoon Pollutant Load Reduction Model Development Project, this study could not be fulfilled without the funding from them.
I would like to thank all the colleagues in room 429: Justin, Detong and Chenxia, for their help and assistance on the model application, developing and debugging; Dave, for being a good buddy on field work; Vadim and Jack, for much advice on making the slides; Kijin and Jun, for being good studying and working partners.
Special thanks go to Sidney, Vik, Vernon, J. J., and Chuck, all the coastal lab specialists, for their guidance and assistance with field work. Grateful thanks are also extended to Becky, Lucy, Sandra, Helen and Doretha, for the administrative work and help.
Finally I would like to thank my parents, my wife and my daughter, for their love and support throughout my graduate research.




TABLE OF CONTENTS
ACKNOWLEDGMENTS ......................... iii
LIST OF TABLES.................................................. vii
LIST OF FIGURES ................................................. viii
ABSTRACT ...................................................... xiV
CHAPTERS
1 INTRODUCTION .................................................1I
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 Method .................................. 6
2.2.2 Fraction of Fresh Water Method.......................... 7
2.2.3 Modified Tidal Prism Method........................... 8
2.2.4 Difference Equation Method ............................ 9
2.2.5 Hydrodynamic Method................................ 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 NUJMERICAL MODEL............................................ 24
4.1 Governing Equations........................................ 24
4.2 Review of Advection Schemes................................. 25
4.3 ULTIM4ATE-QUICKEST Scheme .............................. 27
4.3.1 Third Order Transient Interpolation Model ..................27
4.3.2 Universal Limiter ................................... 30




M onotonic Criteria .................................... 30
Universal Limiters .................................... 31
4.3.3 ULTIMATE-QUICKEST Algorithm ....................... 34
4.4 CH3D with ULTIMATE-QUICKEST Scheme ....................... 34
4.4.1 Dimensionless Advection Equation in CH3D ................ 34
4.4.2 Universal Limiters for CH3D ............................. 36
4.4.3 Third Order Transient Interpolation Method (QUICKEST) for CH3D ................................................. 40
4.4.4 Applying ULTIMATE Algorithm to CH3D .................. 41
5 NUMERICAL TESTS OF ADVECTION SCHEMES ....................... 43
5.1 Test Grid System .............................................. 43
5.2 Two-dimensional Tests ......................................... 45
5.2.1 Test 1: Advection of A Cubic Hill in Uniform Flow ........... 45
5.2.2 Test 2: Advection of A Cone Hill in Uniform Flow ............ 51
5.2.3 Test 3: Advection of A Gaussian Hill in Uniform Flow ......... 57
5.2.4 Test 4: Advection of A Column in Circular Flow ............. 63
5.3 Three-dimensional Test ......................................... 68
5.3.1 Test 5: Advection of a Cubic Block in 3-D Uniform Flow ...... 68
5.4 ULTIMATE-QUICKEST Scheme With Splitting ..................... 72
5.5 Summ ary of Tests ............................................. 75
6 INDIAN RIVER LAGOON FLUSHING CALCULATION ................... 76
6.1 Flushing Time Calculation Using Simple Methods .................... 76
6.2 IRL Flushing Calculation Using CH3D ............................. 77
6.2.1 IRL Grid System ....................................... 77
6.2.2 Forcing Terms / Boundary Conditions ...................... 79
6.2.3 Method for Calculating Flushing Rate and Time .............. 92
6.3 Sensitivity Study of Flushing Rate and Time Using CH3D .... ......... 93
7 CALCULATION OF INDIAN RIVER LAGOON FLUSHING TIME .......... 101
7.1 Indian River Lagoon Segmentation Schemes ....................... 101
7.2 IRL Flushing Results with 9-Segment Scheme ...................... 105
7.2.1 Hydrodynamic Calibration .............................. 105
7.2.2 Flushing Rate Calculation ............................... 107
7.2.3 50% Renewal Time .................................... 112
7.2.4 Flushing M aps ........................................ 113
7.2.5 M ass Conservation Check ............................... 115
7.2.6 D iscussion ........................................... 116
8 CONCLUSION AND RECOMMENDATION ............................ 120




APPENDICES
A TRANSFORMED EQUATIONS IN CUR VILINEAR GRID SYSTEM ........ 124
Coordinate Transformation................................ 124
Dimensionless Equations in Boundary-Fitted Coordinates .........126 Numerical Algorithms of CH3D)............................ 130
B OPEN BOUNDARY INFLOW/OUTFLOW CONDITIONS............... 131
LIST OF REFERENCES............................................. 132
BIOGRAPHICAL SKETCH.......................................... 136




LIST OF TABLES

Table page
2.1 Comparison between different estuarine flushing calculation methods .........13 3.1 Outline of the general information about Indian River Lagoon ................17
4.1 One-dimensional universal limiters................................... 33
4.2 Universal limiters on faces normal to direction......................... 37
4.3 Universal limiters on faces normal to 71 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 (1, J) locations of fresh water input ................................... 81
6.2 50% renewal time at each of the eight segments of IRL.................... 99
7.1 The boundaries of each of the nine segments of Indian River Lagoon ..........105
7.2 50% renewal time at each of the nine segments and the entire Indian River Lagoon
....................................................... 113




LIST OF FIGURES

Figure page
3.1 The Indian River Lagoon eight-segment scheme by Sheng et al .............. 16
3.2 The Indian River Lagoon 18-segment scheme by SJRWMD ................ 22
3.3 Changes of seagrass coverage and nutrient loadings in each of the 18 segments by
SJRWMD over the past 50 years (1943 1992). The # are segment numbers for the nine-segm ent scheme ............................................. 23
4.1 Control volume stencil for QUICKEST scheme .......................... 29
4.2 Local monotonicity criteria and universal limiters ........................ 31
5.1 The 100x 100 non-orthogonal curvilinear test grid system .................. 45
5.2 Analytical solution for the advection of a cubic hill in uniform flow after 100 hours.
The center of the cubic hill is initially located at x0=1.26x104m and y0=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 x0=1.26x104m and y0=0.95x104m, which is close to the origin ........... 49
5.5 Numerical solution for the advection of a cubic hill in uniform flow using the
ULTIMATE-QUICKEST scheme after 100 hours. The center of the cubic hill is
initially located at x0=1.26x104m and y0=0.95x104m, which is close to the origin.
............................................................. 50
5.6 Comparison of numerical and analytical solutions for the advection of a cubic hill in
uniform flow. 2-D projection along the x axis through the center of the analytical




solution. The ULTIMATE-QUICKEST scheme was used.................. 51
5.7 Comparison of numerical and analytical solutions for the advection of a cubic hill in
uniform flow. 2-D projection along the y axis through the center of the analytical solution. The ULTIMATE-QUICKEST scheme was used.................. 51
5.8 Analytical solution for the advection of a cone hill in uniform flow after 106.7 hours.
The center of the cone hill is initially located at xo=1.26x104m and y0=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 y0=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 y0=0.95x104m, which is close to the origin........... 55
5.11 Numerical solution for the advection of a cone hill in uniform flow using the
ULTIMATE-QUICKEST scheme after 106.7 hours. The center of the cone hill is
initially located at x0=1.26x104m and yo=0.95x104m, which is close to the origin.
. . . . . . . .. . . . . . .. . . .. . . . . . .. . . . . .. .. ..56
5.12 Comparison of numerical and analytical solutions for the advection of a cone hill in
uniform flow. 2-D projection along the x axis through the center of the analytical solution. The ULTIMATE-QUICKEST scheme was used.................. 57
5.13 Comparison of numerical and analytical solutions for the advection of a cone hill in
uniform flow. 2-D projection along the y axis through the center of the analytical solution. The ULTIMATE-QUICKEST scheme was used.................. 57
5.14 Analytical solution for the advection of a Gaussian hill in uniform flow after 100
hours. The center of the Gaussian hill is initially located at xo=1.62x104m and yo=1.35x104m, which is close to the origin............................. 59
5.15 Numerical solution for the advection of a Gaussian hill in uniform flow using the
upwind scheme after 100 hours. The center of the Gaussian hill is initially located at x0=1.62x104m and y0=1.35x104m, which is close to the origin............. 60
5.16 Numerical solution for the advection of a Gaussian hill in uniform flow using the
QUICKEST scheme after 100 hours. The center of the Gaussian hill is initially
located at xo=1.62x104m and y0=1.35x104m, which is close to the origin.
............ .................................................. 6 1




5.17 Numerical solution for the advection of a Gaussian hill in uniform flow using the
ULTIMATE-QUICKEST scheme after 100 hours. The center of the Gaussian hill is initially located at x0=1.62x104m and y0=1.35x104m, which is close to the origin62
5.18 Comparison of numerical and analytical solutions for the advection of a Gaussian hill
in uniform flow. 2-D projection along the x axis through the center of the analytical solution. The ULTIMATE-QUICKEST scheme was used .................. 63
5.19 Comparison of numerical and analytical solutions for the advection of a Gaussian hill
in uniform flow. 2-D projection along the y axis through the center of the analytical solution. The ULTIMATE-QUICKEST scheme was used .................. 63
5.20 Analytical solution for the advection of a column in circular flow after 100 hours.
The center of the column is initially located at x0=7.08x104m and y0=6.9x104m. 64
5.21 Numerical solution for the advection of a column in circular flow using the upwind
scheme after 100 hours. The center of the column is initially located at x0=7.08x 104m and y0=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 y0=6.9x104m ...................................... 66
5.23 Numerical solution for the advection of a column in circular flow using the
ULTIMATE-QUICKEST scheme after 100 hours. The center of the column is initially located at x0=7.08x104m and y0=6.9x104m ....................... 67
5.24 Comparison of numerical and analytical solutions for the advection of a column in
circular flow. 2-D projection along the x axis through the center of the analytical solution. The ULTIMATE-QUICKEST scheme was used .................. 68
5.25 Comparison of numerical and analytical solutions for the advection of a column in
circular flow. 2-D projection along the y axis through the center of the analytical solution. The ULTIMATE-QUICKEST scheme was used .................. 68
5.26 Numerical solution after 86.7 hours for the advection of a cubic block in uniform
flow using the upwind scheme. The cubic block on the left front comer 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 comer is the initial condition ................................................... 70




5.28 Numerical solution after 86.7 hours for the advection of a cubic block in uniform
flow using the ULTIMATE-QUICKEST scheme. The cubic block on the left front com er is the initial condition ....................................... 70
5.29 Numerical solution for the advection of a cubic hill in uniform flow (Test 1 case)
using the non-splitting ULTIMATE-QUICKEST method. (Cx)max=(Cy)max=0.45.
.. .. .. .. .... .... ............. .............. .. ... ...............7 3
5.30 Numerical solution for the advection of a cubic hill in uniform flow (Test 1 case)
using the splitting ULTIMATE-QUICKEST method. (Cx)max=(Cy)max=0.45.
..... .. ... ...... ......................... .... ....... ..........74
6.1 The Indian River Lagoon grid system ................................... 78
6.2 Tide data (relative to NAVD88) at Ponce, Sebastian, Ft. Pierce and St. Lucie inlets
from September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998.
........................................................ 82
6.3 Location map of five wind stations .................................... 83
6.4 Wind data at FDEP station 872-1147 (Ponce Inlet) from September 1 (Julian Day
244), 1997 December 31 (Julian Day 365), 1998 ........................ 84
6.5 Wind data at USGS station 02248380 (Haulover Canal) from September 1 (Julian
Day 244), 1997 December 31 (Julian Day 365), 1998 .................... 85
6.6 Wind data at FDEP station 872-1456 (Titusville Brewer Causeway) from September
1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998 ........... 86
6.7 Wind data at FDEP station 872-1789 (Banana River, south end) from September I
(Julian Day 244), 1997 December 31 (Julian Day 365), 1998 .............. 87
6.8 Wind data at FDEP station 872-2213 (Ft. Pierce Inlet) from September 1 (Julian Day
244), 1997 December 31 (Julian Day 365), 1998 ........................ 88
6.9 River discharge data at Eau Gallie River, Crane Creek and Turkey Creek from
September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 89
6.10 River discharge data at Goat Creek, Trout Creek and Fellsmere Canal from
September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 90
6.11 River discharge data at South Prong river, North Canal and Main Canal from




September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 91
6.12 River discharge data at South Canal, North Fork -1 -2 and South Fork -1 -2 from
September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 92
6.13 Relative remaining tracer mass (in percentage of the initial mass) in each of the eight
segments of IRL, obtained using ULTIMATE-QUICKEST advection scheme with A H = A v= 0. ..................................................... 95
6.14 Relative remaining tracer mass (in percentage of the initial mass) in each of the eight
segments of IRL, obtained using ULTIMATE-QUICKEST advection scheme with moderate diffusion of A,=lm2/s ...................................... 96
6.15 Relative remaining tracer mass (in percentage of the initial mass) in each of the eight
segments of IRL, obtained using the upwind advection scheme with moderate diffusion of AH=lm2/s ............................................. 97
6.16 Relative remaining tracer mass (in percentage of the initial mass) in each of the eight
segments of IRL, obtained using ULTIMATE-QUICKEST advection scheme with higher diffusion of AH=5m2/s .......................................98
6.17 Comparison of 50% renewal time at each of the eight segments of IRL:
Test 1: ULTIMATE-QUICKEST advection scheme with no diffusion, A,=0,
Test 2: ULTIMATE-QUICKEST advection scheme with moderate diffusion, A=1m2/s,
Test 3: Upwind advection scheme with moderate diffusion, A=lm2/s,
Test 4: ULTIMATE-QUICKEST advection scheme with higher diffusion, A=5m2/s.
............ .................................................. 99
7.1 Maps of Scheme A, an 8-segment scheme and Scheme B, an 18-segment scheme.
........... .................................................. 102
7.2 Indian River Lagoon 9-segment scheme............................... 104
7.3 One year measured and simulated water surface elevation at FDEP Station 87221251 (Vero Bridge) from 1997.9.1 .................................106
7.4 One year measured and simulated water surface elevation at USGS station 02248380
(Haulover Canal) from 1997.9.1 ....................................107
7.5 Flushing results of nine tracer method: Relative remaining mass in each of the nine
segments based on one-year simulation using splitting ULTIMATE-QUICKEST advection scheme and normal diffusion coefficient (AH =lm2/s). The lower figure is a closer look during the first 100 days ...............................110




7.6 Flushing results of one tracer method: Relative remaining mass in each of the nine
segments and the entire Indian River Lagoon based on one-year simulation using ULTIMATE-QUICKEST advection scheme and normal diffusion coefficient (AH =lm2/s). The upper figure includes results in segment 1, 2, 3, 4 and 5; the lower figure includes results in segment 6, 7, 8, 9 and the entire Indian River Lagoon (T O TL) ....................................................... 111
7.7 IRL flushing maps: (a) P=0.9, a contour map of 10% renewal time; (b) P=0.75, a
contour map of 25% renewal time ................................... 114
7.8 IRL flushing maps: (c) P=0.5, a contour map of 50% renewal time; (d) P=0.37, a
contour map of 63% renewal time ................................... 115
7.9 Change of seagrass coverage during the past 50 years (1943-1992) within each of the
18 segments according to SJRWMD ................................. 117
7.10 Contour map of R50 value within each of the nine segments (nine-tracer method).
..... .. ... .. ...... ............. ........ .. ... ... ... ........... 118
7.11 Interpolated map of the Indian River Lagoon bottom sediment size (D50).
.. ... . .. ........ ...... ........ .......... ....... .............. 119
A. 1 Transformation from z coordinate to o coordinate.
.. ... .. . .. .. .......... ........ .......... ....... ..... ......... 125
A.2 Transformation from (x, y) coordinates to ( ,Ti) coordinates.
... ... .. ... ....... ............. .......... ... ... ... ... ......... 126




Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science
USE OF THREE-DIMENSIONAL HYDRODYNAMIC MODEL AND
ACCURATE ADVECTION SCHEME FOR MODELING FLUSHING DYNAMICS AND
DEVELOPING SEGMENTATION SCHEMES IN INDIAN RIVER LAGOON By
Haifeng Du
August 2000
Chairperson: Dr. Y. Peter Sheng
Major Department: Civil and Coastal Engineering
Recent studies show that many Florida estuaries are facing the problems of excessive nutrient loading and consequent effects of degenerated water quality and sea-grass loss. For example, in the Indian River Lagoon (IRL), the average phosphorus and nitrogen loadings have been doubled, and the total sea-grass coverage has decreased about 50% during the past 50 years. In order to improve the IRL water quality, restore the sea-grass beds and protect the estuarine ecosystem, a goal of maximum loadings that IRL can sustain (Pollutant Load Reduction Goal) needs to be established.
As part of the IRL-PLR (Pollutant Load Reduction) Model Development Project under the supervision of Dr. Y. P. Sheng, this study focuses on two major components. The first component is the description and the application of an accurate numerical scheme




(ULTIMATE-QUICKEST) for the modeling of tracer advection process. The second component is the application of the three-dimensional hydrodynamic numerical model CH3D with the ULTIMATE-QUICKEST scheme to the Indian River Lagoon flushing and segmentation study. A boundary-fitted non-orthogonal curvilinear grid has been deployed in order to best represent the complex geometry of Indian River Lagoon. Vertically arstretching technique is employed. Various numerical tests and comparisons have been. conducted to verify the superior accuracy of the new advection. scheme in non-orthogonal curvilinear grid. Numerical sensitivity studies were also performed to see how different advection schemes and model parameters could influence the flushing results. Previous IL segmentation schemes have been reviewed and compared and a nine-segment scheme has been proposed. Flushing simulations based on the entire Indian River Lagoon and each of the nine segments were carried out separately using 1-tracer method and 9-tracer methods. The simulation results show that the Banana River, which is almost an enclosed water body, has the lowest flushing rate with R50 value (50% renewal time) equal to about 200 days in the northern section of the Banana River, and R50 value equal to about 100 days in the southern section of the Banana River. However, some segments, which have inlets directly connected to the ocean, have very high flushing rates, e.g., Segments 7 (including Sebastian Inlet) and 9 (including Ft. Pierce Inlet and St. Lucie Inlet). The R50 values of these two segments are about 2 days and 11I days, respectively. The flushing rates of all other segments, e.g., Segments 4 and 5 (the northern Indian River), Segment 8 (Vero Beach area) and Segment 1 (Mosquito Lagoon), are more or less a month. The total tracer mass was also checked to ensure the mass conservation of the model.




CHAPTER 1
INTRODUCTION
1.1 Background
Sea-grass beds provide habitats 'and nursery areas for estuarine animals, enhance water quality by removing nutrients and stabilizing sediments, and are one of the most productive ecosystems on the planet (Zieman 1982). Recent studies showed that the penetration of photo-synthetically active radiation (PAR) through the water column is the primary determinant of the depth limitation of sea-grass system (Kenworthy et al. 1991, Morris and Vimstein 1993). With increasing population and development, many Florida estuaries such as Indian River Lagoon, are facing the problems of declining sea-grass beds due to excessive nutrient/pollutant loadings, degrading water quality and increasing light attenuation. While some of the nutrients/pollutants can be quickly flushed out of the estuary through inlets, some of them will remain in the estuary for a long time. If the nutrient/pollutant loading rate is faster than the estuarine natural flushing rate, the estuarine water quality will degenerate gradually. Hence, in order to improve estuarine water quality, we need to first understand estuarine hydrodynamics and flushing characteristics.
Generally speaking, flushing is the capability of an estuary or a part of an estuary to replace its old water mass with new ocean water mass. It is conceptually very simple but




2
complex to calculate in reality. There are many methods that have been used to quantify estuarine flushing characteristics. Traditionally volumetric driven methods were used to calculate flushing characteristics (Ketchum 1951, Stommel and Arons 1951, Dyer and Taylor 1973). Despite their simplicity, these methods have several drawbacks and are not reliable to represent the real world situation. First, most of these methods assume complete mixing. Second, they consider a single principal tidal constituent but exclude wateI: movements caused by other tidal constituents and non-tidal forces. Third, estuarine hydrodynamics are not incorporated. More recently the hydrodynamic driven method (Edinger et al. 1998, Sheng et al. 1995), which considers the estuarine hydrodynamics, is being used in favor of the traditional volumetric driven methods.
1.2 Scope of The Study
As part of the Indian River Lagoon Pollutant Load Reduction (IRL-PLR) model project sponsoredbythe St. John's River Water Management District (SJRWMD), this study focuses on the quantification of the Indian River Lagoon flushing characteristics, using a 3-D hydrodynamic model CH3D (Sheng 1986a, 1986b, 1987, 1989). Here the Indian River Lagoon (IRL) is a linked system of three water bodies: Mosquito Lagoon, Banana River Lagoon and Indian River Lagoon, not only the Indian River itself. A relatively fine boundaryfitted curvilinear grid (477x43x4) for the Indian River Lagoon was used in order to fit the complex geometry and bathymetry. A third-order accurate and oscillation-free advection scheme has also been developed for CH3D to calculate the flushing time of a passive tracer.




3
Because of IRL's large size and complex geometry, flushing characteristics are very different from place to place. Some regions in the lagoon may have good flushing and some may have very poor flushing. A reasonable segmentation scheme, which reflects the spatial variation in flushing characteristics in IRL, is necessary. Sheng et al. (1995) calculated the IRL flushing characteristics by using the vertically-integrated version ofCH3D and provided a preliminary segmentation scheme with eight segments. However, Sheng et al. (1995) only considered the flushing due to tidal circulation forced by simple harmonic tides at the tidal inlets. This segmentation scheme will be compared with that produced from the sea-grass perspective (Vimstein et al. 1994). A nine-segment scheme has been provided and more accurate flushing calculations for the entire lagoon and each individual segment will be conducted by using the fully three-dimensional version of CH3D with a highly accurate advection scheme and considering realistic circulation driven by tide, wind and density over a one year period.
In summary, the scope of the study includes:
* Apply the ULTIvLATE-QUICKEST advection scheme (Leonard 1991) to the threedimensional non-orthogonal curvilinear coordinate system.
0 Verify the scheme in non-orthogonal curvilinear coordinate system by conducting
various numerical tests.
0 Develop CH3D flushing model with the ULTIMATE-QUICKEST advection scheme.
* Apply CH3D flushing model to the study of Indian River Lagoon flushing
characteristics.
* Analyze flushing results and compare different Indian River Lagoon segmentation




4
schemes and propose an optimal segmentation scheme.
Chapters two and three present the background knowledge for this study. The definition of flushing, different flushing methods, and previous flushing and segmentation studies of Indian River Lagoon will be reviewed in these two chapters. Chapters four and five cover the numerical model and test results, including CH3D hydrodynamic model, flushing model, the ULTIMATE-QUICKEST scheme with CH3D, and scheme verification. Chapters six and seven present various IRL segmentation schemes and flushing results. Conclusions and recommendations are presented in chapter eight.




CHAPTER 2
FLUSHING IN ESTUARIES
2.1 Definitions and Classifications
An estuary is where the river meets the sea (Fischer et al. 1979). Pritchard (1967), based on his study on James River in Virginia, defined an estuary as "a semi-enclosed coastal body of water which has a free connection with the open sea and within which sea water is measurably diluted with fresh water derived from land drainage." This definition, however, does not work very well with all kinds of estuaries in the world. In Florida, for example, most of the estuaries are in the form of coastal lagoons, e.g., Indian River Lagoon, Sarasota Bay, and Biscayne Bay, etc. Kjerfve (1994) defined the coastal lagoon as, "a shallow coastal water body separated from the ocean by a barrier, connected at least intermittently to the ocean by one or more restricted inlets, and usually oriented shore-parallel."
Since almost half of the world's population are now living around estuaries, growing human activities, (e.g., effluent disposal, dredging of navigational channels, reclamation of wetlands or open water, agriculture runoff, etc.) are causing adverse impacts on the estuarine environments.
One of the major concerns nowadays is the impact of pollutant loadings on estuarine water quality, habitat, and living resources. Understanding estuarine flushing is essential to




6
understanding the impact of pollutant loadings on the estuarine ecosystem. Flushing
I
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 riier 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 "volumetric" methods which do not consider the estuarine hydrodynamics, while most recent methods are "hydrodynamic" 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= V 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 h-ixing.
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:
f=(0-S (2.2)
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'jIJ then the flushing time can be described as:
t=f J" (2.3)
Qf Qf
where, Qf is the river inflow (Fischer et al. 1979). However this method seems to be too strictly defined to be used in systems that lack river or river discharge is much less significant




8
than the tidal prism. For example, in Indian River Lagoon, most rivers are very small creeks or canals and some segments even do not have any river inflow, i.e., Qf is very small or even equals to zero in some segments. According to this method, the flushing time could reach infinity, which is not true in the real world.
2.2.3 Modified Tidal Prism Method
This method was originally introduced by Ketchum (1951) and has been further developed by Dyer and Taylor (1973) and Wood (1979). Like the fraction of fresh water method, this method uses fresh water as the tracer. It divides an estuary into several segments, 0, 1, 2, ... considering one-dimensionally from the head of the estuary down to the mouth of the estuary, and assumes that the high tide and low tide are in phase for all segments. The inner segment (segment 0) is defined based on the assumption that the intertidal volume (tidal prism) in segment 0, Pa, is supplied by the river flow during the flood tide, i.e., P0 = R/2, where R is the river inflow during one tidal cycle. Second, segment 1 is placed so that V = 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, V,,+1 = UV, +P,, for the rest of the segments. High and low water concentrations of fresh water C", CL are defined for each segment. At high water, H L L
(V+P,), =V,+C,,+l(n,,C,'2) (2.4)

At low water,




V,,C[, =(a V,,+R)C,_I + [(1-a)V,,-R]C, (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 C H, CL can be obtained. Then the estuarine flushing time t can be calculated in the following formula: t= (2.6)
(V +P)C(2
Although this method considers incomplete mixing in each segment by adjusting parameter c, for the same reason of "Fraction of Fresh Water Method," it seems only applicable to the flushing calculations of river estuaries but not coastal lagoons, especially those with relatively small rivers or even no river inflow, e.g., Indian River Lagoon.
2.2.4 Difference Equation Method
More recently, Lowery (1998) developed a "Difference Equation-Based Estuarine Flushing Model" and applied it to a few Gulf of Mexico estuaries. Focusing on "tidal flushing" rather than "river flushing," this method incorporates the functionality of the "tidal prism method," though it still uses fresh water as the tracer. By tracking a single input of fresh water over several time steps (tidal periods) until it is flushed out of the estuary, the flushed fresh water volume is defined as:
+1 p (2.7)
r;V,+




10
where, vrn+l is the flushed fresh water volume at time step n+l; V" is the remaining fresh water volume at time step n, V"=Vn-I- v; P is the tidal prism, which is the volume between high tide and low tide; V is the low tide estuarine volume. Lowery calculated the Gulf of Mexico's estuaries' low tide volumes and intertidal volumes from NOAA's Physical Environmental Characterizations Branch's (PECB) unpublished data bases. The long term fresh water inflows were obtained from United States Geologic Survey (USGS) and Texas Department of Water Resources' (TDWR) data bases. After this, he then calculated the estuarine flushing time. Lowery's model provides some improvements over previous volumetric methods in: 1) considering the estuarine stratification by adjusting the low tide estuarine volume to 75%, 50% or 25% depending on the extent of the vertical stratification; and, 2) considering the fresh water back-flushing from the mouth based on the following equation,
Vt,=P(1 -)SO (2.8)
35
where, Vb is the fresh water volume re-entered to the estuary by back-flushing, P is the tidal prism, So is the salinity at the mouth, the ocean water salinity is assumed as 35 ppt. This method is simple and can apply to the flushing calculations of coastal lagoons. However, it does not include the estuarine hydrodynamic process. The narrow focus on tidal flushing also excludes many other non-tidal flushing features.




2.2.5 HydrodyLiamic Method
The volumetric methods are simple and direct for calculating the estuarine flushing characteristics. One doesn't need to understand the estuarine hydrodynamics before applying these methods. However, these methods have severe restrictions due to their narrow focus and unrealistic assumptions: (1) Tidal prism method only focuses on the tidal flushing without considering fresh water input. It uses one simple harmonic tidal constituent and considers only one tidal period; Or in terms of fresh water input, the fresh water flushing method calculates the estuarine flushing without explicit consideration of exchange with ocean water; (2) Most of these methods use only the tidal prism, estuarine volumes at high tide and low tide, and cross-section averaged velocities to drive the model; (3) Besides tide and/or river discharge (fresh water input), none of other major driving forces, such as wind and density gradients, is considered in these methods; (4) Most of these methods assume steady and uniform estuarine conditions, and instant and complete mixing in estuarine sections; (5) Some of the older methods ignore the effect of tidal back-flushing; (5) These models are one-dimensional models along the head of the estuary to the mouth of the estuary;
(6) Estuarine hydrodynamics are not incorporated in either of these methods. Because of the narrow focus and unrealistic assumptions, the results based on some of theses methods become very unreliable. For example, according to the "freshwater flushing rate method," it takes 715 days to replace Tampa Bay's fresh water and in contrast, only 5 days are required to flush Tampa Bay according to the "tidal prism flushing method" (Biggs 1986).
More recently, with the development of more powerful computers, the hydrodynamic




12
method is becoming more popular (Sheng et al. 1995, Edinger et al. 1998). This method requires a quantitative understanding of the estuarine hydrodynamics under the influence of multiple forcings. Usually the continuity equation and the Reynold's averaged Navier-Stokes equations of motion need to be solved in order to obtain the surface elevation and the flow field. Then the estuarine flushing time (e.g., 50% renewal time) can be determined by tracking the time history of relative remaining tracer mass contained in any segment of the estuary (or the entire estuary) which is initially filled with a uniform tracer concentration (e.g., 100). The conservative advection-diffusion equation needs to be solved along with the equations of motion until the total tracer mass in the segment (or the entire estuary) drops to below 50% of its initial total mass. Rather than considering one tidal cycle as the time step for estimating the flushing time, this method solves the transport equation, the equations of motion, and the continuity equation simultaneously at a much smaller time step (on the order of seconds or minutes). At each time step the tracer concentration will be calculated according to the updated new surface elevation and flow field. The one-dimensional, twodimensional or three-dimensional model can be used. The flushing results are very sensitive to the numerical schemes used for solving the advection-diffusion equation. Traditional first order upwind scheme usually overestimates the flushing rate due to excessive numerical diffusion. Certain higher order schemes may give negative concentration because of numerical dispersion. An accurate and oscillation-free advection scheme is needed for flushing calculation, especially in those regions where the flows are advection-dominated.
Table 2.1 summarizes all of the methods for calculating flushing time in estuaries.




13
Table 2.1 Comparison between different estuarine flushing calculation methods
Methods Driving Forces Time Dimension Mixing Tidal Back- Estuarine
step Type flushing hydrodynamics
Tidal Prism Harmonic tide, I tidal O-D instant and not not
Estuarine cycle ("box" complete considered considered
volumes, Tidal model) in estuary
Prism
Fraction of Harmonic tide, I tidal O-D instant and not not
Fresh Water Estuarine cycle ("box" complete in considered considered
volumes, Tidal model) estuary
prism, Fresh water input
Modified Harmonic tide, I tidal I-D considering not not
Tidal Prism Estuarine cycle (from head incomplete considered considered
volumes, Tidal to mouth) mixing by
prism, Fresh adjusting
water input empirical
parameters
Difference Harmonic tide, I tidal O-D considering considered not
Equation Estuarine cycle ("box" vertical considered
volumes, Tidal model) stratification
prism, Fresh by simple
water input adjustment
of
coefficient
Hydro- Real tide, seconds 1, 2 or 3-D Eddy considered considered
dynamic wind, river or viscosity
discharge, minutes and
density diffusivity
I gradient I I I concept I I




CHAPTER 3
INDIAN RIVER LAGOON SEGMENTATION
3.1 Describing Indian River Lagoon
The Indian River Lagoon complex system, which lies on the east coast of Florida, including Mosquito Lagoon, Banana River Lagoon and Indian River Lagoon (Figure 3.1), is a long (more than 200 km long), narrow (typically 2-4 km wide) and shallow (typically 1-3 m deep) estuary with four inlets connecting it to the Atlantic Ocean and numerous rivers and canals discharging fresh water into it. The northern end of Mosquito Lagoon is connected to the ocean through Ponce de Leon Inlet. The southern part of Indian River Lagoon is connected to the ocean by Ft. Pierce Inlet and St. Lucie Inlet. Sebastian Inlet connects the middle part of Indian River Lagoon to the ocean. The southern end of Mosquito Lagoon is connected to Indian River Lagoon through Haulover Canal. Banana River is an almost enclosed water body with only two narrow connections with the Indian River Lagoon. Indian River Lagoon is a micro tidal system with the semi-diurnal M2 tide as the principal constituent (Smith 1987). Tide is the major force causing flushing for most of the Indian River Lagoon, except in the northern IRL and Banana River, where tidal action is weaker. The geometry of the Indian River Lagoon is very complicated because of highly irregular boundaries, barrier islands, and numerous causeways. Indian River Lagoon has a mean depth




15
of about two meters, with a maximum depth of nearly ten meters and a minimum depth less than half meter. Dredging of the inter-coastal waterway has added to the complexity of the bathymetry. Due to the complex geometry and bathymetry, the flushing characteristics are very different within different segments of the Indian River Lagoon. Segments with strong tidal action and less constriction usually have higher flushing rates, while segments with weak tidal action and more constriction have slower flushing rates. Based on the natural boundaries and the preliminary flushing study (Sheng et al. 1995), Indian River Lagoon was divided into eight segments (Figure 3.1). Table 3.1 gives an outline of some general information about Indian River Lagoon and its eight segments.




3.2E+06
3.15E+06
- \ (3)
I-I
3.1E+63.05E+06
>.' :(8)
3E+06'
500000 525000 550000 575000 600000
X UTM (m) Figure 3.1 The Indian River Lagoon eight-segment scheme by Sheng et al..




17
Table 3.1 Outline of the general information about Indian River Lagoon
Segment 1 2 3 4 5 6 7 8 Total
North 405 516 5-6 free 6-7 free 7-8 free --Boundary shore shore shore highway highway connection connection connection
South 405 3 516 5-6 free 6-7 free 7-8 free --Boundary shore highway highway highway connection connection connection shore
West
Boundary shore shore shore shore shore shore shore shore
East
Boundary shore shore shore shore shore shore shore shore
Ponce de Sebastian Ft Pierce
Tidal Inlet Leon Inlet --- --- --- --- Inlet --- Inlet
St Lucie
Inlet
Annually Turnbull 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 Gallic TurkeyCr.: 3.28; 1.17; 1.28; 53.17
Discharge Flounder: Riverl + 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
(i) II_ I
Maximum
Depth 9.02 3.79 7.17 6.14 3.70 4.39 4.01 11.14 11.14
(in)
Average
Depth 1.47 1.53 1.69 2.11 2.06 1.62 1.42 1.82 1.71
(in)
Minimum
Water 0.1813 0.2097 0.2859 0.2816 0.0897 0.0663 0.0299 0.2433 1.3877
Volume (kin3)
Maximum
Water 0.2075 0.2765 0.3214 0.3119 0.1067 0.0938 0.0441 0.3245 1.6864
Volume (ki3)
Average
Water 0.1978 0.2498 0.3093 0.2900 0.0995 0.0771 0.0361 0.2786 1.5382
Volume
(km3)
Tidal Prism - 0.1493
(m 3)" -- 019




3.2 PLRG and PLR Model
In order to improve the Indian River Lagoon water quality and restore the sea grass bed, Florida Department of Environmental Protection (FDEP), with St. John's River Water Management District's (SJRWMD's) agreement, requires that numerical Pollution Load Reduction Goals (PLRG's) be established for the major sources of pollution in Indian River Lagoon. A PLRG specifies the maximum amount of a pollutant (contaminated sediments, pesticides, nutrients, and toxic chemicals) that will be permitted to enter into Indian River Lagoon or a segment of IRL through various tributaries. In order to develop scientifically defensible PLRG's, an Indian River Lagoon Pollutant Load Reduction (IRL-PLR) model is being developed by the University of Florida (Dr. Y. P. Sheng is the Principal Investigator), with funding from the St. Johns River Water Management District (SJRWMD). The project aims to develop a Pollution Load Reduction (PLR) model which incorporates a circulation model, a sediment transport model, a water quality model, a light model and a sea grass model. An important first step is to develop a segmentation scheme and an understanding of the flushing characteristics. As part of the PLR model development effort, this study will provide an accurate calculation of flushing characteristics within the entire lagoon and each of the nine segments determined in this study. The flushing model sensitivity study was also conducted.




3.3 Previous Indian River Lagoon Segmentation Studies
Smith (1993, 1996) used a one-dimensional lagoon wide model to calculate the flushing time in Indian River Lagoon. He divided Indian River Lagoon into 16 segments (excluding Mosquito Lagoon and Banana River) and Banana River into 3 segments. He then quantified the intertidal volume by using harmonic constituents (principal semi-diurnal and diurnal) recorded at 31 NOS stations. He multiplied the tidal amplitudes by the surface areas, and incorporated the phases to represent the tidal volumes and the tidal movements. The model first estimated the averaged one-dimensional flow between the segments and then solved a one-dimensional advection-diffusion equation. In his 1996 version model, he added "Flexible Grid" technique which allows control of the extent to which segment boundaries move along the longitudinal axis of the lagoon with the flow. If the boundaries can move freely with the flow, the transport equation is solved with diffusion only, and if the boundaries are fixed, then the transport equation is solved with advection only. Based on one year simulation he found that 50% renewal time is about 140 days if transport by advection (fixed grid), and 50% renewal time is not reached until the end of the simulation if transport by diffusion (fully flexible grid). The results suggest that it tends to underestimate flushing if using fully flexible grid and overestimate flushing if using fixed grid. Another series of simulations using completely flexible grid show that in the southern and central sub-basins of Indian River Lagoon 50% renewal times were about 5 and 12 days respectively while the northern sub-basin of Indian River Lagoon and the Banana River reached only about 30%




20
renewal by the end of one year simulation. Although the method he used is comparable to the traditional modified tidal prism method, his results still can be used as a rough indication of the IRL flushing characteristics.
Sheng et al. (1993a, 1993b) used a one-dimensional lagoon-wide model and a threedimensional lagoon wide model to calculate the flushing in the Indian River Lagoon. They found that the one-dimensional model tended to overestimate the flushing rate. Further, as part of the preliminary study for IRLPLR model, Sheng et al. (1995) conducted a preliminary study of the Indian River Lagoon segmentation scheme. A vertically integrated twodimensional version of the three-dimensional curvilinear-grid hydrodynamic model CH3D was used to calculate the IRL flushing. Because of lack of data, they assumed semi-diurnal M2 tides with amplitudes of 5 cm, 15 cm and 30 cm at Ponce, Sebastian and Ft. Pierce inlets respectively. Without considering any other forcing functions, they simulated the flushing characteristics corresponding to four segmentation schemes. Based on the flushing results and the hydrodynamic and water quality data from 24 FDEP or USGS monitoring stations, they found that the eight-segment scheme (Figure 3.1) appeared to be the best scheme.
In 1994, St. Johns River Water Management District (SJRWMD) conducted a spatial analysis of seagrass changes versus changes in nutrient loadings over the past 50 years (Virnstein et al. 1994). Based on changes of seagrass coverage and drainage sub-basin boundaries, they divided Indian River Lagoon into 18 segments (Figure 3.2). Then they calculated change of sea grass coverage and changes of nutrient loadings for each segment from 1943 to 1992, and regressed percentage change of seagrass coverage on percentage changes of nutrient loadings. The results as shown in Figure 3.3 indicate that most of the




21
segments have significant increases in nutrient loadings and corresponding decreases in seagrass coverage during the past 50 years except segment 15 where the fast flushing rate appears to have flushed increased nutrient loadings into the ocean.
In this study, simulation of flushing characteristics will be first conducted with eight segments developed by Sheng et al. (1995), and finally with the nine segments developed based on the result of comparison between eight-segment (Sheng et al. 1995) and eighteensegment schemes (SJRWMVD 1994). The nine-segment scheme will be described in following chapters of this thesis.




3.2E+06

550000
X UTM (m)

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

-2>
<8 <>k <9
<1 0>
- <4>
<11>
- Ay<5>
<12>
<7>
- <13>
-<14>
<1 6 >
<17>
<1 8>

3.15E+06

E
I
3.1 E+06

3.05E+06

3E+06

500000

I I

600000




% Seagrass Change,
by segment

-80 -40 0

40 80

'4 -TI
-< C/)
rd
0w (D
S 4'q
D0
(D ~
-D
CD

% Change In N Loading
by segment
0 50 100150200250300 12
3
Mm4
5 Me6
-a7
8
9
10
lS
12
13
5 165 116
17
18TOTAL
TOTAL

% Change in P-loading,
by segrhent
0 100 200 300 400 500 12
3

% Change in TSS-loadln
by segment
0 200 400 600
Hi
12

14
- 5
7 7

gIo /11
13
14
15
116
17 = i
TOTAL

18
19
=10 11
12
13
- 14
15
*16 M 17 =18
M TOTAL

TOTAL




CHAPTER 4
NUMERICAL MODEL
4.1 Governing Equations
The governing three-dimensional momentum equations can be derived from Newton's second law of motion which is based on the concept of conservation of momentum. The continuity and scalar transport equations can be derived from the mass conservation criteria. With the assumptions of incompressible flow, vertically hydrostatic water column, Boussinesq approximation, and turbulent eddy viscosity concept, these equations were expressed in the following forms (Sheng 1989):
a-u+ v+ aw =0 (4.1)
ax ay az
au auu +auv+auw=fv 1 ap a au a au a au
- +-+ -+-++-(AH-) + -(AH-) + (A ) 4.2)
at ax ay az p0 axx a x ay ay +az A vz)
av avu av 8vw 1 8p 8 ( v 8 8v 8 8v
-+.+-+ -fu -- + -(A (A + (A-)(43
at ax ay az po ay ax Hax + (Ay )+a z(A ) (4.3)z z
0 = -P9 (4.4)
8z




as aus+avs aws a as + a as a as
- + + -= -(DH )+ ( -)+-(D-,) (4.5)
at ax ay az ax ax ay ay az az aT auT avT awT a aT a aT a aT
-+--+-+ -(KH-)+-(KH )+-(Kr-) (4.6)
at ax ay az ax ax ay ay az az
a +p + au + av + aw (KH2 )+ (K _)+a(Kr ) (4.7)
at ax ay az ax ax ay ay az az
p = p(S,2') (4.8)
where: t is the time, (u, v, w) are the mean velocities in the (x, y, z) directions, f=20 sin D 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 H 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, 4 is any kind of dissolved conservative species concentration.
Instead of dealing with these equations directly, CH3D solves the dimensionless equations in a vertically-stretched and horizontally boundary-fitted system. These equations were derived by Sheng (1986) and are shown in Appendix A.
4.2 Review of Advection Schemes
Numerical simulation of circulation and transport processes is very important for many engineering and scientific studies in coastal waters, lakes and estuaries. In particular, when quantifying the transport of pollutants and water quality parameters, the accuracy of




the advection scheme is very crucial.
In the past decades, many numerical schemes have been developed and applied to the study of tracer advection. Probably the simplest scheme is the first order upwind scheme. However, fewer people are using it now because of its high numerical diffusion. Another popular scheme is the Leap Frog scheme which is used in the princeton Ocean Model (POM). This scheme does not have as much dissipation as the up vind scheme but suffers a strong dispersive error. Numerical filter has to be used to suppress the strong oscillation, thus increasing the diffusion in the combined scheme. More recently QUICKEST scheme (Leonard 1979) has been applied to CH3D (Sheng 1989), because of its simplicity, efficiency, third-order accuracy with little dissipation and dispersion.
Many numerical schemes have been developed recently to eliminate spurious oscillations while maintain high order accuracy. One of the schemes is the TVD (Total Variation Diminution) scheme which used so-called "flux limiting" technique, e.g., Roe (1985). Although TVD scheme is totally free from numerical overshoots and undershoots near discontinuities, some studies showed that the scheme also has many disadvantages such as "clipping" of local extrema and "steepening" of (what should be) gentle gradients (Leonard 1991). Leonard (1991) proposed the concept of "universal limiter" and developed an algorithm, ULTIMATE, which can be used to arbitrary high order explicit conservative finite difference schemes. By comparing the most popular 36 advection schemes (including TVD and ULTIMATE), Cahyono (1993) concluded that the ULTIMATE scheme seemed to be the most "attractive." He applied these advection schemes to six numerical test cases for both one-dimensional and two-dimensional flows and found that the ULTIMATE scheme




27
was more general than the other schemes and was easier to apply. Lin et al. (1997) also applied this scheme in his two-dimensional transport model and simulated the water quality and sediment transport in Humber Estuary, U. K..
The next section will describe how the ULTIMATE algorithm is incorporated in CH3D.
4.3 ULTIMATE-QUICKEST Scheme
The ULTIMATE algorithm can be used to explicit conservative finite difference schemes of any order accuracy. The higher order of accuracy, the more interpolation points will be involved and the less efficiency of computation. For most large scale and relatively smooth flows the cost-effective third-order upwind (QUICKEST) ULTVIATE scheme gives excellent practical results (Leonard 1991). The following sections will only present the formulation of the ULTIMATE-QUICKEST scheme for the three-dimensional nonorthogonal curvilinear-coordinate numerical model CH3D.
4.3.1 Third Order Transient Interpolation Model
Considering unsteady one-dimensional pure advection of a scalar 4 at constant velocity u, the equation is
.___.+u__ =0. (4.9)
at ax
Applying center-space explicit method, equation (4.9) can be interpreted into the following




finite difference form
in l= C_-cli,- ? (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*, +p )-SGN(c)( -4)n)], (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 I on all involved center values,
P) r(i- 1). (4.13)
If we consider second order accuracy, one of the well-known schemes is the Lax-Wendroff scheme. The right face value for the L-W scheme can be written as




r 2=[( 1)-( 1-44)]. (4.14)
For QUICKEST, the third order upwind interpolation scheme (Leonard 1979), the right face value can be written in terms of the second order Lax-Wendroff value and an upwind biased curvature (Leonard 1991),
(3)( =(2)( (1-c2) 1 + SGN(c)(dn 2n+_)+ 1 SGN(c)(d 24 ++q)] (4.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
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: cn+l= n [
=Ir C [cr4 r- t] i' (4.16)
where cr= A._t and At are the Courant numbers at the right and left faces,
L X AX




respectively. Note that cp(i= 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 4. (upstream), 4d(downstream) and 4 ,(centered between the two) are defined for the convenience of considering both positive and negative velocities, as shown in figure 4.2.




i-1 i i+1 i-1 i i+1
C.V. Uf C.V. uu
(a) u > 0 (b) u < 0
(a) u1>O (bC) u < O

Figure 4.2 Local monotonicity criteria and universal limiters. If the right face velocity is positive, i.e., u>0, then define = q= n n n u<0, then define U =j+2, n+1 n+ 7. Considering positive right face velocity, i.e., u>0 (Figure 4.2 a), for 4n in the locally monotonic range n ,n n ,n,, we get the monotonic condition for control volume face value: S &,. (4.17)
In order to maintain monotonicity, the updated center value must also be constrained by n+1 tn+l n+1
u n+q 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 value4,u at n+1 time step will not be greater than the value at n time step, i.e., 4,,1n .




Hence, the left-hand side inequality in equation (4.18) is less restrictive than cIt' n+l. (4.19)
For the right-hand side of equation (4.18), the monotonic criteria require the right face

value 41f to be limited by

n+
C :g %"

(4.20)

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

(4.21)

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

+1 =C- _(cAfc Cl ), then inequality (4.19) and equation (4.22) can be combined into the following form: wf
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:

n c
nmin(ci nc c)+ n-n
4"a 4pmin[d_.

) + ,,. ,
min(c4 u,c) c) + n
The upper limiter for f is superr min[" ], the lower limiter is
lower= 4n. 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 cfis the control
volume right face Courant number, c, is the next left face Courant number, cr is the next
right face Courant number.
Table 4.1 One-dimensional universal limiters Right face Local Right face limiters for
velocity, monotonicity
Uf upper limiter, lower limiter,
I_ upper lower
positive, nmin(c,# ct c n nn
positive, n nncm)+ax ctnc
>0 4 4 d c
CC
positive, nn[ n ma(c4 +4 n n
>0 4u 4 d c c d)
C
negative, nmnn+m.
<0 (n, :>n> (gn a[d, min(Cruc) c u
<0 u) C) _) d) max[) &U]
Cf
negative, c)nc+).
<0 Cu< n n n MaX(Cr~uCr)<0 min [ Cdr&rC
I _cf

(4.25)




4.3.3 ULTIMATE-QUICKEST Algorithm
Based on the knowledge of the third order transient interpolation model and the universal limiters, the ULTIMATE-QUICKEST algorithm can be established as following:
* Designate upstream(u), downstream(d), and central(c) nodes on the basis of SGN(uf)
for each control volume face.
* Compute DEL= -g1-, ADEL=IDELI and ACURV=I d-2 +$I for each
face. If ACURV>ADEL(non-monotonic), set = and proceed to the next face.
* Calculate the limiters( 4,per, #tower) for each face based on Table 4.1 and set up the
QUICKEST face value 4f based on equation (4.15).
* Limit (4j by limiters upper and lower"
* Update the control volume center value 4 c according to equation (4.22).
4.4 CH3D with ULTIMATE-QUICKEST Scheme
4.4.1 Dimensionless Advection Equation in CH3D
Considering the advection terms only in equation (A. 16), the dimensionless threedimensional scalar advection equation in CH3D (without diffusion terms) is:
I ~t +R R1(
1 +( ) Rb [( )+. r(go + )]+ (H 4~)=0, (4.26)
ft ai KtIa a-9 aa
where: f= H .COR, XREF XREF REF ar
wh =tCOR, =u v=v' U g, a -e
ZREF UREF UREF UREF (XREFf




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

(_)1(, )n
(04) r' (
Ai

R- )x
[1 +

-RbZ
AUF
With some mathematical manipulation, we get:
n+ H c -Rb Aix
Al f V~ AAt
-Rb-[9y AU n1 Z Z Defining:

~ n
RAHC- c
n*I
IC
- H (Jgo)
-At xfy Ox TXF=Rbu X
Hf (Jg)y
A At YV0f
TYF=RbvY,
TZF=Rb 6 At
z AG

TH, = t-Hx,f g-O)x,
At
TXL = Rb x ,
TYL=Rb
TZL=Rb At
b so

(4.29)

(4.27)

(4.28)

r-- )71m y,
f]




equation (4.28) can be simplified as:
+1 nJ41C)4
= (RAHQ- C)(4.3n
- [(TX )x,- (rx.n)- l(m), -( )-4, -[(rz.)-, -(z (4.30)
4.4.2 Universal Limiters for CH3D
According to Leonard & Niknafs (1990), the one-dimensional limiter method described above can be directly applied to three-dimensional control volume by taking three direction sweeps (g, rl, a) separately.
First considering t sweep only, we take the first and second terms on the right-hand side of equation (4.30):
CC+ =(RAHC)* -[(T *x~-(TXL).)" (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 CO Cm,,
positive, {n n
u?0 cn n n XdC
min[(TX)-4x.,(MX)-4x] + (RAHC)-4,
(TXF)
positive, x n
u?0 # n g qn max Xd
max[(ma).4 x,,(M ). ,x] +(RAHC) n (TXF
negative, nmax n
u, min[(TXR),, x,(TXR).,] -(RAHC)4, +
(TXF)
negative, n n
Uf<0 n x xd )' + .
max[(X)~,(T)4x ]-(RAHC)4n
_________ I_ (TXF)

Then, considering I1 and a sweeps respectively,
n=(RH n [(Ty. .-(TY). ,,

(4.32)

and

tn+l= (RAHQ*4n -[( F) -(TZL), ZCo zc). ( ZF.I- (TZ),l j

(4.33)




the limiters on faces normal to i 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
face monotonicity velocity in ri
in TI direction upper lower
direction, limiter, limiter,
Vf
positive, m n{
v + Y n on dd
min[(TYL)- .,(TYL)*y] +(R,,AHC)*.n- ;
(TMF) positive, n
max[(TYL).- ,(TYL).,4 ] + (RAHC) n.
negative, n{
v<0 n no g >a nYd
min[(TYR)~ ,(TYR)4,]- (RAHC)' + n
negative, n. n
vY<0 nn (u 4Yd n
max[(TER)- ,(TYR)- o]- (RAHC)* Yn + n, I(TMF I




Table 4.4 Universal limiters on faces normal to a direction
Right Local Right face limiters for
face monotonicity velocity in a
in a direction upper lower
direction, limiter, limiter,
CJf IOw.
positive, n
S 0 n positive, a ,
)>O0 n n n n ma Zd
max[(TZL)-z,(TZL)* ]+(RAHC).4
(Tze))
min[(TZR)-#z,(TZR)*z] (RAHC)- n~+' (TZF)
negative, "
,Wf ~f<0 '~~Zc4 Zd CCd
max[(TZR), ,(TZ R) z-(RAHC)
(ZF) _n

Where: in Tables 4.2, 4.3 and 4.4,

TXR=Rbix At(r)x TYR=Rby Al Y(r (. 4

TZR=RbG zt.
to

(4.34)




40
4.4.3 Third Order Transient Interpolation Method (QUICKEST) for CH3D
The one-dimensional QUICKEST interpolation formula (equation 4.15) can be directly used in CH3D to obtain the control volume face values if considering three direction sweeps separately. The only difference is that the Courant numbers have to involve the dimensionless coefficients because the equations in CH3D are dimensionless.
The QUICKEST right face value of control volume in three directions for CH3D can be expressed as:
-n (n n
(1-C2- 1+SGN(c ,.-2 1-SGN(c,
6 2 2 (
6 2 n+1-2 }+#1-1)+ 2 n+2-20+ 1)] (4.35)
in 4 direction,
W (4j+1;
I- 1cy ISN( n I-- S,+4N;c)+ 2 n+n
6 2 N(C,) (4.36)
in il direction,
and
Ik n n C,, n n
2~ 2
(1-c 1 +SGN(c) 1 S GN(c,)
6 2 2n n ng-+ n n(
- 2 2.(4 .37-2)k + 2-2,+,+(*1
6 2 2 (4.37)




in a direction, where:
At I (UREF) At 1 R At
c =u- u RbuS (REF)(COR) A '
C v A (URM At I =R V At
YfA l ( -)(OR i Anl
At (UREF) At 1 At
Ao (XREF) (COR) Ac AG
9, A o (AREF) (COR) Ao Ao (4.38)
are control volume face Courant numbers in ll, and a directions respectively.
4.4.4 Applying ULTIMATE Algorithm to CH3D
Splitting method. This is the most direct method to extend one-dimensional ULTIMATE-QUICKEST scheme to the three-dimensional CH3D. Simply implement three sweeps in E, lI, and a directions consecutively while updating the central value within each sweep. This method solves only one-dimensional advection equation at every intermediate time step, and can make better use of the one-dimensional limiters. The advantage of this method is its accuracy and its applicability on relative large Courant number conditions (C < 1).
Non-splitting method. The basic procedure is similar but with only one sweep. 1) Apply one-dimensional ULTIMATE strategy described above in E, T", ando directions respectively to find the limited control volume face values; 2) Update the central value according to equation (4.30). This method requires solving the three-dimensional advection equation at a time using the one-dimensional limiters. The Courant number condition for this method is more restrictive than that for splitting method:




42
At 1
u + Ivl+ Iw.L (4.39)
Ax Ay AZ
Both of the methods were tested and the results are included in the next chapter.




CHAPTER 5
NUMERICAL TESTS OF ADVECTION SCHEMES
Leonard (1991) performed many one-dimensional pure advection tests which showed the advantages of the ULTIMATE algorithm against other schemes. This chapter will focus on two- and three-dimensional tests in a non-orthogonal curvilinear coordinate system. A test grid was set up and five bench mark tests (advection of a cubic hill, a cone hill, and a Gaussian hill in two-dimensional uniform flow, advection of a column in two-dimensional circular flow, and advection of a cubic block in three-dimensional uniform flow) were conducted. Results of the splitting and non-splitting ULTIMATE-QUICKEST methods were also compared with different Courant numbers.
5.1 Test Grid System
The test domain is a rectangular basin with a uniform depth of 6000 m, and is divided into 100x100 grid cells with i=1 to 101 in X direction andj=1 to 101 in Y direction (Figure 5.1). Instead of a uniform grid spacing for the entire area, a variable grid spacing and nonorthogonal grid-lines are established in order to test the new scheme in the non-orthogonal curvilinear grid system. The grid points on the four boundaries are defined as follows:




(1) Along the southern boundary (Y=0): xi=,x(i- 1), 1< i< 31;
xi=x1_+ x[0.5+(51-i)/20], 321i 51;
xi=x51+^x(i-51), 52gik 101, (5.1)
(2) Along the northern boundary (Y=Y ): xi=tx(i-1), 1: is 31;
x,=xt_1+Ax[0.5+(i-32)/20], 32 i 51;
xi=x51'+x(i-51), 52 (3) Along the western boundary (X= 0): yj =AY( 1), 1 yj=yj _+Ay[0.5+(51-j)/20], 32 j 51;
Yj =Y51 +ty(f- 51), 52 j 101, (5.3)
(4) Along the eastern boundary (X=X ): yJ=Ay(f- 1), 1 y=y_+Ay[0.5+(f-32)/20], 32 yj=y51 +Ay(f- 51), 52 where Ax-1200m and y=-1000m.
Thus, the grid is orthogonal and uniform with Ax=-1200m and &y=1000m in the following four regions: (1) 1 Four two-dimensional and one three-dimensional bench mark tests are selected to verify the new scheme.




0
25000 50 00 75 00 100000 X (m)
Figure 5.1 The 100x100 non-orthogonal curvilinear test grid system.
5.2 Two-dimensional Tests
5.2.1 Test 1: Advection of A Cubic Hill in Uniform Flow

The initial concentration distribution of this test is defined as:
If =1.0, Ix-xol0jdx and ly-yoj:sdy;
1 =0.0, elsewhere,

(5.5)




46
is located on the lower left comer of the grid domain shown in Figure 5.1, which is close to the origin (x=-0 and y=O).
The flow field is uniform: u=O.18m/s along the x direction, and v=0.15m/s along the y direction. Three numerical schemes (the upwind, QUICKEST and non-splitting ULTIMATE-QUICKEST) were tested. Figures 5.2 5.7 give the results after 1800 time steps with At=-200s (100 hours).
Theoretically, the pure advection of any type of profile will give the exact same shape of the profile as the initial one. As seen in Figure 5.2, the cubic hill should keep its original shape after 100 hour advection in uniform flow. Numerically it is hard to obtain the exact same profile. Different numerical schemes give different accuracy. In this test the upwind scheme gives the lowest accuracy due to its severe numerical diffusion. Although the QUICKEST scheme keeps the higher order accuracy, it cannot avoid the spurious oscillations at sharp gradient regions (shown in figure 5.4). Only the ULTIMATE-QUICKEST scheme achieves the higher order accuracy without any undershooting and overshooting (Figure 5.5).




x 10

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 x0=1.26x104m and yo=0.95x104m, which is close to the origin.




48
100
80
60
e40
0
0
01
8
6 10
X 104 4 6 8
2 4 x 10,
2
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
100
80
60
40,
20
01
8
6 10
x lO1 4 6
2 4 xlO
2
Y (m) X (m)
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 x0=1.26x104m and yo=0.95x104m, which is close to the origin.




x 10

x 10

Y (m)

X (m)

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




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

Figure 5.7 Comparison of numerical and analytical solutions for the advection of a cubic hill in uniform flow. 2-D projection along the y axis through the center of the analytical solution. The ULTIMATEQUICKEST 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.0 (x -xo)2 + (y-o) /Ro, (x -xo)2 + y-yo)2 R; {= 0.0, elsewhere,

(5.6)

where x0=12600m, yo0=9500m and Ro=8000m.
The flow field is the same as in Test 1: u=0.18m/s along the x direction, and v=0.15m/s along the y direction. Three numerical schemes (the upwind, QUICKEST and non-splitting ULTIMATE-QUICKEST) were tested. Figures 5.8 5.13 give the results after 1920 time steps with At=--200s (about 106.7 hours).




52
This test was designed to show the performance of the schemes on local extremum. Due to the high diffusion, the profile is severely damped and spreads around if the upwind scheme is used (Figure 5.9). The QUICKEST scheme maintains the third order accuracy, but unfortunately it cannot avoid the undershooting at the base of the profile (Figure 5.10). The ULTIMATE-QUICKEST scheme seems to have the best performance. It removes the spurious oscillation while maintaining the same order of accuracy as QUICKEST scheme (Figure 5.11).




80
c 60,
40
20,
0
8
6 10
x10. 4 6 8
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.26x10'm and yo=0.95x10'm, which is close to the origin.




100
80
c 60
SD 40o
0 -"
20, 01
88
x 10 4 6
" 4 xlo
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 yo=0.95x104m, which is close to the origin.




1I II

x 10

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 x0=1.26x104m and yo=0.95x104m, which is close to the origin.




56
100.
80, 60,.
0
o
20 0
8
6 10
x 04 6 8
2 4 x 10
2
Y (m) X (m)
Figure 5.11 Numerical solution for the advection of a cone hill in uniform flow using the ULTIMATE-QUICKEST scheme after 106.7 hours. The center of the cone hill is initially located at xo=1.26x104 m and yo=0.95x104m, which is close to the origin.




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

Figure 5.13 Comparison of numerical and analytical solutions for the advection of a cone hill in uniform flow. 2-D projection along the y axis through the center of the analytical solution. The ULTIMATEQUICKEST 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:
4 =exp[- ], (5.7)
2R2
-0
where, xo=16200m, yo=13500m and Ro=4000m.
The flow field is the same as in Tests 1 and 2: u=0.18m/s along the x direction, and v=-0.15m/s along the y direction. Three numerical schemes (the upwind, QUICKEST and non-splitting ULTIMATE-QUICKEST) were tested. Figures 5.14 5.19 give the results after 1800 time steps with At=200s (100 hours).




58
Still, in this test the upwind scheme is proved to be with high damping. It seems that for the smoothly changed profile, both the ULTIMATE-QUICKEST and the QUICKEST schemes have similarly good performance. Sometimes for a very steep peak value like in this test case, the "clipping of the peak value" is difficult to avoid when using the ULTIMATEQUICKEST scheme.




c 60
cc 40-...
20 0
8
6 10
x 10 4 6
2 4 x 104
22 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 x0= 1.62x 104m and yo=1.35x 104m, which is close to the origin.




60
100
80
c 60.2
40
C
20
0
8
6 10
x 104 4 6
2 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 x0=1.62x104m and y0=1.35x104m, which is close to the origin.




61
100
80 60 40
20-m
0
8
6 10
X 104 4 6
2 2 4 x 10
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 x0=1.62x104m and yo=1.35x104m, which is close to the origin.




100 .
80 .
c 60
0
,D 40C.)
20
O 0
8
6 10
8
x 10, 4 6
2 4 x 10,
2
Y (m) X (m)
Figure 5.17 Numerical solution for the advection of a Gaussian hill in uniform flow using the ULTIMATE-QUICKEST scheme after 100 hours. The center of the Gaussian hill is initially located at xo=1.62x104m and yo=1.35x104m, which is close to the origin.




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

Figure 5.19 Comparison of numerical and analytical solutions for the advection of a Gaussian hill in uniform flow. 2-D projection along the y axis through the center of the analytical solution. The ULTIMATEQUICKEST 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:

4= 1.0, (x-xo)2+(Y-yo0)2R; 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 (o=1/59400 1/s, and the Cartesian velocity can be described as: u= -(o(y-yc), v=W(x-xc). Three numerical schemes (the upwind, QUICKEST and non-splitting ULTIMATE-QUICKEST) were tested. Figures 5.20 5.25 give the results after 1800 time




steps with At=200s (100 hours).
Since all the tests described above are advections in uniform flow. This test was designed to see the performance of different schemes in non-uniform flow. Again the test results show that the ULTIMATE-QUICKEST scheme appears to be the best scheme.

100
80 S60
S40
0 C
20 0

x 10'

x 10

Y (m) X (m)

I
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.




x 10

x 10

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 104m and y0=6.9x104m.




66
oo
100
0
= 42-4 x -10
20
0
6 10
x 104 24 x 10 4
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.08x10m and yo=6.9x104m.




10080.
0 60. a 40.
0
20
x 104

2 4 x10
2
Y (m) X (m)
Figure 5.23 Numerical solution for the advection of a column in circular flow using the ULTIMATE-QUICKEST scheme after 100 hours. The center of the column is initially located at x0=7.08x104m and y0=6.9x 104m.

8
6 10
4 6

A4




100
80
- Nnorkf olusn
20
1 2 3 4 5 6 7 8 9 10 11 X(m) xa
Figure 5.24 Comparison of numerical and analytical solutions for the advection of a column in circular flow. 2-D projection along the x axis through the center of the analytical solution. The ULTIMATEQUICKEST scheme was used.

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

5.3 Three-dimensional Test
5.3.1 Test 5: Advection of a Cubic Block in 3-D Uniform Flow
The initial concentration distribution for this test is:

=1.0, Ix-xol dx, ly-y01o dy, Iz-z0o dz;
1=0.0, elsewhere,

(5.9)

where, x0o=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,




QUICKEST and non-splitting ULTIMATE-QUICKEST) were tested. Figures 5.26 5.28 give the results after 1560 time steps with at-200s (about 86.7 hours).

100000

x (m)

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 corner is the initial condition.




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 comer is the initial condition.

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

-4000
-6000

100000

25000 50000 75000
X (M)




71
The figures show that there is little difference between the QUICKEST and the ULTIMATE-QUICKEST schemes. A detailed analysis, however, does show significant advantage of the ULTIMATE-QUICKEST scheme. Table 5.1 shows the analysis of the maximum percentage of over and under prediction and the error ratio of mass conservation. There is no overshooting or undershooting for the ULTIMATE-QUICKEST scheme. Table 5.1 Comparisons of numerical schemes for Test 5: The advection of a cubic block in uniform flow
Numerical Scheme Erroe, Err under Err..
1 upwind -44.10% 0.0% 2.85 x10-6 %
2 QUICKEST 24.8% -7.7% 2.25 x10-6 %
3 ULTIMATE-QUICKEST 0.0% 0.0% 2.41 x105 %
Where, in the table, the Percentage of Overshooting ErrVer is defined as: Errver (4max) 1- (a)
Err axl00%. (5.10)
The Percentage of Undershooting Err wider is defined as:
Err'wider = (mid)_ (ki n)aX100%. (5.11)
(4 ax a EL
The Percentage of Mass Conservation Error Errmm is defined as:
Err., = Total Final Mass-Total Initial Mass x100%. (5.12)
Total Initial Mass
Where, ((mx) a and (ki.)a are the maximum and minimum values of analytical solution, (




72
and (4i) are the maximum and minimum values of numerical solution.
5.4 ULTIMATE-QUICKEST Scheme With Splitting
The bench mark tests using the non-splitting ULTIMATE-QUICKEST method seem to give us very encouraging results. Further study shows that they are only applicable to very small Courant numbers (C0.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 nonsplitting and splitting methods at Cmax=0.45.




100
80
60
40
0
C),
20
01
8
6 10
x 11 4 4 8
2 4 x 1d
2
Y (m) X (m)
Figure 5.29 Numerical solution for the advection of a cubic hill in uniform flow (Test 1 case) using the non-splitting ULTIMATE-QUICKEST method. (Cx)max=(Cy)max=0.45.




S40.
20-.
0
8
6 10
x 4 6 8
2 4 xli
2
Y (m) X (m)
Figure 5.30 Numerical solution for the advection of a cubic hill in uniform flow (Test 1 case) using the splitting ULTIMATE-QUICKEST method. (Cx)max=(Cy)max=0.45.
Based on these tests, It seems that there is little difference between the results of splitting and non-splitting schemes when the Courant number is very small (C<0.1). As the




75
Courant number becomes larger and approaches unity, the results of splitting method are much better than that of the non-splitting method.
5.5 Summary of Tests
This chapter has described several numerical test problems: (1) Advection of a cubic hill in 2-D uniform flow; (2) Advection of a cone hill in 2-D uniform flow; (3) Advection of a Gaussian hill in 2-D uniform flow; (4) Advection of a column in 2-D circular flow; and
(5) Advection of a cubic block in 3-D uniform flow. These test problems represent different conditions of the real world. For each test problem, three advection schemes: (1) upwind, (2) QUICKEST, and (3) ULTIMATE-QUICKEST schemes were compared. The test grid is specially designed for verification of the schemes based on the non-orthogonal curvilinear coordinate system. For all the five test problems, the upwind scheme has too much diffusion. The QUICKEST scheme cannot avoid spurious oscillation at sharp gradient regions. Only the ULTIMATE-QUICKEST scheme seems to have the best performance. For the ULTIMATE-QUICKEST scheme, both the splitting and non-splitting method were implemented. Further tests with different Courant numbers show that both the Splitting and non-splitting methods work well if the Courant number is less than 0.1. However, when the Courant number is greater than 0.1 and less than unity, the non-splitting method is proved to be unacceptable due to the severe deformation of the original test profile, while the splitting method can still keep the same good result no matter what Courant number it is. Thus the splitting ULTIMATE-QUICKEST scheme will be used in this study.




CHAPTER 6
INDIAN RIVER LAGOON FLUSHING CALCULATION
6.1 Flushingz 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' in3. The tidal prism is about 0. 15x109 m3(The estuarine maximum volume minus the estuarine minimum volume). Assuming semi-diurnal tides with T= 12.4 hours, the flushing time can be obtained as, t = 5.65 days, thus, the R50 value is less than 2.83 days.
The second method used here is the Fraction of Fresh Water Method. The total river discharges into the MR (the runoff discharge is neglected) is estimated as 53.17 m 3/s based on the summation of averaged one year historical data at the major 17 rivers (will be discussed later). Assuming the average estuarine salinity is 20 ppt and the ocean salinity is 35 ppt, the flushing time can be determined by Equation (2.3): t = 157.3 days. The R50 value is 78.7 days.
The large difference between results of the Tidal Prism Method and the Fraction of




77
Fresh Water Method indicates that these simple methods are not very useful when studying such a complex system like Indian River Lagoon. This points to the need of the hydrodynamic method which uses a hydrodynamic model to calculate the flushing time.
6.2 IRL Flushing Calculation Using CH3D
6.2.1 IRL Grid System
Before we can apply CH3D with the third order ULTIMATE-QUICKEST advection scheme to the simulation of the flushing characteristics of Indian River Lagoon, an accurate grid system has to be set up. In order to represent the Indian River Lagoon's complex geometry and bathymetry, a high resolution boundary-fitted grid is necessary. By applying Thompson's (1983) grid generation method, Sheng et al (1995) generated a non-orthogonal curvilinear grid for Indian River Lagoon (Figure 6.1). The grid is of very high resolution with horizontally 477x43 cells (I=1 to 477 in the north-south direction and J=1 to 43 in the westeast direction) and vertically 4 layers (K=1 to 4).




3.2E+06 3.15E+06
E
I
3.1 E+06 3.05E+06 3E+06

500000

550000
X UTM (m)

Figure 6.1 The Indian River Lagoon grid system.

600000




6.2.2 Forcing Terms / Boundga Conditions
The major forcing functions for LRL circulation need to be determined once the grid is set up. Four major forcing functions are considered in this study: surface elevations, winds, river discharges and salinity distributions. The first three forcing functions (surface elevation, wind and river discharge) are given by time series of measured data at selected locations. More than one year continuous data were prepared (September 1997 to December 1998) in order to represent the seasonal variation. The fourth forcing function, salinity distribution, is taken into account by calculating the baroclinic circulation term in the numerical model. Other less dominant forcing functions, e.g., precipitation, evaporation, surface runoff (except for Merritt Island) and ground water seepage, etc., are not considered.
At the open boundary, a simplified one-dimensional advection equation is used. At tidal boundaries, tidal back flushing was considered by calculating inflow and outflow advection fluxes based on the direction of velocities. At river boundaries, only inflow condition is considered, because river discharge is always positive. Details of inflow and out flow conditions can be found in Appendix B.
Water surface elevation data. The Indian River Lagoon estuarine system has four tidal inlets connecting the system to the Atlantic Ocean. As one boundary condition, the surface elevation at an inlet is given by a time series of measured data. These data were collected at four EDEP tide stations at the inlets and leveled relative to NAVD88. The data used in this study contains water surface elevations every 30 minutes over a period of more than one year (September 1997 to December 1998). Figure 6.2 shows the time series of surface elevation




80
data at the four inlets over a one-year period (September 1997 December 1998). It is interesting to see that there is a declination of mean sea level from north to south along the shore. According to the data, the one year averaged mean sea level at the Ponce de Leon Inlet is about 11 cm higher than that at the Ft. Pierce Inlet.
Wind data. The free surface boundary condition requires the wind stress at every grid c6ll. In this study, the measured wind speed and direction at selected locations are interpolated over the entire grid. Five FDEP and USGS wind stations are selected and their locations can be seen in Figure 6.3. Figures 6.4 6.8 show time series of hourly-averaged wind speed at X (eastern is positive) and Y (northern is positive) directions respectively at the five stations over a period of more than one year (September 1997 December 1998).
River discharge data. A total of 17 rivers are considered for the system. The locations of these rivers can be seen in Figure 6.1 and Table 7. L 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 (1, J) locations of all the river and runoff sites. Figures 6.9 6.12 show the time series river discharge plots of some of the major rivers. Smaller rivers and runoff from Merritt Island are not shown in these figures.




Table 6.1 (I, J) locations of fresh water input (I, J) location
I J
1 Tumbull Creek 101 13
2 Big Flounder 109 13
3 Sykes Creek 216 24
4 Horse Creek 258 7
5 Eau Gallie River-1 269 3
6 Eau Gallie River-2 273 9
7 Crane Creek 280 7
8 Turkey Creek 289 5
9 Goat Creek 303 7
10 Trout Creek 309 7
11 Fellsmere Canal 317 9
12 South Prong River 325 1
13 North Canal 365 12
14 Main Canal 377 12
15 South Canal 388 12
16 North Fork 439 11

South Fork

450

1 169-173 21
2 Northern Indian River 176 199 20
3 203-212 20
4 215-234 20
5 BananaRiver 177 199 28
6 204-213 28




265

200

300

Sebastian Inlet
I .1.

265 0 100 200 300

|Ft. Pierce Inlet)

265 0 100

200

265
Julian Day 1997

200
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.

50
0
> -50
z E -100
-S

-100

3UU

I tLuc
, L kj AL A It
, t , 2.2_lo ll ,L




FDEP 872-1147

3.2E+06
USGS 02248380
FDEP
872-1456
3.15E+06
FDEP 872-1789
I-
3.1 E+06
3.05E+06
FDEP 872-2213
3E+06 I I I I
500000 550000 600000
X UTM (m)

Figure 6.3 Location map of five wind stations.




10
5
a/)
CL
-5
-10
265 0 100 200 300
Julian Day 1997 Julian Day 1998
Figure 6.4 Wind data at FDEP station 872-1147 (Ponce Inlet) from September 1 (Julian Day
244), 1997 December 31 (Julian Day 365), 1998.

Station 872-1147x (East-West) I

, I I IFDEP E