UFL/COEL 2000/007
USE OF THREEDIMENSIONAL HYDRODYNAMIC MODEL AND ACCURATE ADVECTION SCHEME FOR MODELING FLUSHING DYNAMICS AND DEVELOPING SEGMENTATION SCHEMES IN INDIAN RIVER LAGOON
by
Haifeng Du
Thesis
2000
USE OF THREEDIMENSIONAL HYDRODYNAMIC MODEL AND
ACCURATE ADVECTION SCHEME FOR MODELING FLUSHING DYNAMICS AND
DEVELOPING SEGMENTATION SCHEMES IN INDIAN RIVER LAGOON
By
HAIFENG DU
A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
UNIVERSITY OF FLORIDA
2000
Copyright 2000
by
Haifeng Du
ACKNOWLEDGMENTS
First, I would like to express my gratitude to my advisor, Dr. Sheng, for his direction and financial support in making this study a success. I would also want to thank the other two members of my committee, Dr. Dean and Dr. Thieke, for carefully reviewing this thesis.
My appreciation is extended to the St. Johns River Water Management District. As part of the Indian River Lagoon Pollutant Load Reduction Model Development Project, this study could not be fulfilled without the funding from them.
I would like to thank all the colleagues in room 429: Justin, Detong and Chenxia, for their help and assistance on the model application, developing and debugging; Dave, for being a good buddy on field work; Vadim and Jack, for much advice on making the slides; Kijin and Jun, for being good studying and working partners.
Special thanks go to Sidney, Vik, Vernon, J. J., and Chuck, all the coastal lab specialists, for their guidance and assistance with field work. Grateful thanks are also extended to Becky, Lucy, Sandra, Helen and Doretha, for the administrative work and help.
Finally I would like to thank my parents, my wife and my daughter, for their love and support throughout my graduate research.
TABLE OF CONTENTS
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 ULTIM4ATEQUICKEST 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 ULTIMATEQUICKEST Algorithm ....................... 34
4.4 CH3D with ULTIMATEQUICKEST Scheme ....................... 34
4.4.1 Dimensionless Advection Equation in CH3D ................ 34
4.4.2 Universal Limiters for CH3D ............................. 36
4.4.3 Third Order Transient Interpolation Method (QUICKEST) for CH3D ................................................. 40
4.4.4 Applying ULTIMATE Algorithm to CH3D .................. 41
5 NUMERICAL TESTS OF ADVECTION SCHEMES ....................... 43
5.1 Test Grid System .............................................. 43
5.2 Twodimensional Tests ......................................... 45
5.2.1 Test 1: Advection of A Cubic Hill in Uniform Flow ........... 45
5.2.2 Test 2: Advection of A Cone Hill in Uniform Flow ............ 51
5.2.3 Test 3: Advection of A Gaussian Hill in Uniform Flow ......... 57
5.2.4 Test 4: Advection of A Column in Circular Flow ............. 63
5.3 Threedimensional Test ......................................... 68
5.3.1 Test 5: Advection of a Cubic Block in 3D Uniform Flow ...... 68
5.4 ULTIMATEQUICKEST Scheme With Splitting ..................... 72
5.5 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 9Segment Scheme ...................... 105
7.2.1 Hydrodynamic Calibration .............................. 105
7.2.2 Flushing Rate Calculation ............................... 107
7.2.3 50% Renewal Time .................................... 112
7.2.4 Flushing 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 BoundaryFitted Coordinates .........126 Numerical Algorithms of CH3D)............................ 130
B OPEN BOUNDARY INFLOW/OUTFLOW CONDITIONS............... 131
LIST OF REFERENCES............................................. 132
BIOGRAPHICAL SKETCH.......................................... 136
LIST OF TABLES
Table page
2.1 Comparison between different estuarine flushing calculation methods .........13 3.1 Outline of the general information about Indian River Lagoon ................17
4.1 Onedimensional universal limiters................................... 33
4.2 Universal limiters on faces normal to direction......................... 37
4.3 Universal limiters on faces normal to 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 eightsegment scheme by Sheng et al .............. 16
3.2 The Indian River Lagoon 18segment scheme by SJRWMD ................ 22
3.3 Changes of seagrass coverage and nutrient loadings in each of the 18 segments by
SJRWMD over the past 50 years (1943 1992). The # are segment numbers for the ninesegm 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 nonorthogonal curvilinear test grid system .................. 45
5.2 Analytical solution for the advection of a cubic hill in uniform flow after 100 hours.
The center of the cubic hill is initially located at 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
ULTIMATEQUICKEST 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. 2D projection along the x axis through the center of the analytical
solution. The ULTIMATEQUICKEST scheme was used.................. 51
5.7 Comparison of numerical and analytical solutions for the advection of a cubic hill in
uniform flow. 2D projection along the y axis through the center of the analytical solution. The ULTIMATEQUICKEST scheme was used.................. 51
5.8 Analytical solution for the advection of a cone hill in uniform flow after 106.7 hours.
The center of the cone hill is initially located at xo=1.26x104m and 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
ULTIMATEQUICKEST 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. 2D projection along the x axis through the center of the analytical solution. The ULTIMATEQUICKEST scheme was used.................. 57
5.13 Comparison of numerical and analytical solutions for the advection of a cone hill in
uniform flow. 2D projection along the y axis through the center of the analytical solution. The ULTIMATEQUICKEST scheme was used.................. 57
5.14 Analytical solution for the advection of a Gaussian hill in uniform flow after 100
hours. The center of the Gaussian hill is initially located at xo=1.62x104m and yo=1.35x104m, which is close to the origin............................. 59
5.15 Numerical solution for the advection of a Gaussian hill in uniform flow using the
upwind scheme after 100 hours. The center of the Gaussian hill is initially located at x0=1.62x104m and y0=1.35x104m, which is close to the origin............. 60
5.16 Numerical solution for the advection of a Gaussian hill in uniform flow using the
QUICKEST scheme after 100 hours. The center of the Gaussian hill is initially
located at xo=1.62x104m and 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
ULTIMATEQUICKEST 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. 2D projection along the x axis through the center of the analytical solution. The ULTIMATEQUICKEST scheme was used .................. 63
5.19 Comparison of numerical and analytical solutions for the advection of a Gaussian hill
in uniform flow. 2D projection along the y axis through the center of the analytical solution. The ULTIMATEQUICKEST scheme was used .................. 63
5.20 Analytical solution for the advection of a column in circular flow after 100 hours.
The center of the column is initially located at x0=7.08x104m and y0=6.9x104m. 64
5.21 Numerical solution for the advection of a column in circular flow using the upwind
scheme after 100 hours. The center of the column is initially located at 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
ULTIMATEQUICKEST 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. 2D projection along the x axis through the center of the analytical solution. The ULTIMATEQUICKEST scheme was used .................. 68
5.25 Comparison of numerical and analytical solutions for the advection of a column in
circular flow. 2D projection along the y axis through the center of the analytical solution. The ULTIMATEQUICKEST scheme was used .................. 68
5.26 Numerical solution after 86.7 hours for the advection of a cubic block in uniform
flow using the upwind scheme. The cubic block on the left front 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 ULTIMATEQUICKEST 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 nonsplitting ULTIMATEQUICKEST 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 ULTIMATEQUICKEST 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 8721147 (Ponce Inlet) from September 1 (Julian Day
244), 1997 December 31 (Julian Day 365), 1998 ........................ 84
6.5 Wind data at USGS station 02248380 (Haulover Canal) from September 1 (Julian
Day 244), 1997 December 31 (Julian Day 365), 1998 .................... 85
6.6 Wind data at FDEP station 8721456 (Titusville Brewer Causeway) from September
1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998 ........... 86
6.7 Wind data at FDEP station 8721789 (Banana River, south end) from September I
(Julian Day 244), 1997 December 31 (Julian Day 365), 1998 .............. 87
6.8 Wind data at FDEP station 8722213 (Ft. Pierce Inlet) from September 1 (Julian Day
244), 1997 December 31 (Julian Day 365), 1998 ........................ 88
6.9 River discharge data at Eau Gallie River, Crane Creek and Turkey Creek from
September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 89
6.10 River discharge data at Goat Creek, Trout Creek and Fellsmere Canal from
September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 90
6.11 River discharge data at South Prong river, North Canal and Main Canal from
September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 91
6.12 River discharge data at South Canal, North Fork 1 2 and South Fork 1 2 from
September 1 (Julian Day 244), 1997 December 31 (Julian Day 365), 1998. .. 92
6.13 Relative remaining tracer mass (in percentage of the initial mass) in each of the eight
segments of IRL, obtained using ULTIMATEQUICKEST advection scheme with 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 ULTIMATEQUICKEST 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 ULTIMATEQUICKEST advection scheme with higher diffusion of AH=5m2/s .......................................98
6.17 Comparison of 50% renewal time at each of the eight segments of IRL:
Test 1: ULTIMATEQUICKEST advection scheme with no diffusion, A,=0,
Test 2: ULTIMATEQUICKEST advection scheme with moderate diffusion, A=1m2/s,
Test 3: Upwind advection scheme with moderate diffusion, A=lm2/s,
Test 4: ULTIMATEQUICKEST advection scheme with higher diffusion, A=5m2/s.
............ .................................................. 99
7.1 Maps of Scheme A, an 8segment scheme and Scheme B, an 18segment scheme.
........... .................................................. 102
7.2 Indian River Lagoon 9segment scheme............................... 104
7.3 One year measured and simulated water surface elevation at FDEP Station 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 oneyear simulation using splitting ULTIMATEQUICKEST advection scheme and normal diffusion coefficient (AH =lm2/s). The lower figure is a closer look during the first 100 days ...............................110
7.6 Flushing results of one tracer method: Relative remaining mass in each of the nine
segments and the entire Indian River Lagoon based on oneyear simulation using ULTIMATEQUICKEST advection scheme and normal diffusion coefficient (AH =lm2/s). The upper figure includes results in segment 1, 2, 3, 4 and 5; the lower figure includes results in segment 6, 7, 8, 9 and the entire Indian River Lagoon (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 (19431992) within each of the
18 segments according to SJRWMD ................................. 117
7.10 Contour map of R50 value within each of the nine segments (ninetracer method).
..... .. ... .. ...... ............. ........ .. ... ... ... ........... 118
7.11 Interpolated map of the Indian River Lagoon bottom sediment size (D50).
.. ... . .. ........ ...... ........ .......... ....... .............. 119
A. 1 Transformation from z coordinate to 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 THREEDIMENSIONAL HYDRODYNAMIC MODEL AND
ACCURATE ADVECTION SCHEME FOR MODELING FLUSHING DYNAMICS AND
DEVELOPING SEGMENTATION SCHEMES IN INDIAN RIVER LAGOON By
Haifeng Du
August 2000
Chairperson: Dr. Y. Peter Sheng
Major Department: Civil and Coastal Engineering
Recent studies show that many Florida estuaries are facing the problems of excessive nutrient loading and consequent effects of degenerated water quality and seagrass loss. For example, in the Indian River Lagoon (IRL), the average phosphorus and nitrogen loadings have been doubled, and the total seagrass coverage has decreased about 50% during the past 50 years. In order to improve the IRL water quality, restore the seagrass beds and protect the estuarine ecosystem, a goal of maximum loadings that IRL can sustain (Pollutant Load Reduction Goal) needs to be established.
As part of the IRLPLR (Pollutant Load Reduction) Model Development Project under the supervision of Dr. Y. P. Sheng, this study focuses on two major components. The first component is the description and the application of an accurate numerical scheme
(ULTIMATEQUICKEST) for the modeling of tracer advection process. The second component is the application of the threedimensional hydrodynamic numerical model CH3D with the ULTIMATEQUICKEST scheme to the Indian River Lagoon flushing and segmentation study. A boundaryfitted nonorthogonal curvilinear grid has been deployed in order to best represent the complex geometry of Indian River Lagoon. Vertically arstretching technique is employed. Various numerical tests and comparisons have been. conducted to verify the superior accuracy of the new advection. scheme in nonorthogonal curvilinear grid. Numerical sensitivity studies were also performed to see how different advection schemes and model parameters could influence the flushing results. Previous IL segmentation schemes have been reviewed and compared and a ninesegment scheme has been proposed. Flushing simulations based on the entire Indian River Lagoon and each of the nine segments were carried out separately using 1tracer method and 9tracer methods. The simulation results show that the Banana River, which is almost an enclosed water body, has the lowest flushing rate with R50 value (50% renewal time) equal to about 200 days in the northern section of the Banana River, and R50 value equal to about 100 days in the southern section of the Banana River. However, some segments, which have inlets directly connected to the ocean, have very high flushing rates, e.g., Segments 7 (including Sebastian Inlet) and 9 (including Ft. Pierce Inlet and St. Lucie Inlet). The R50 values of these two segments are about 2 days and 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
Seagrass beds provide habitats 'and nursery areas for estuarine animals, enhance water quality by removing nutrients and stabilizing sediments, and are one of the most productive ecosystems on the planet (Zieman 1982). Recent studies showed that the penetration of photosynthetically active radiation (PAR) through the water column is the primary determinant of the depth limitation of seagrass system (Kenworthy et al. 1991, Morris and Vimstein 1993). With increasing population and development, many Florida estuaries such as Indian River Lagoon, are facing the problems of declining seagrass beds due to excessive nutrient/pollutant loadings, degrading water quality and increasing light attenuation. While some of the nutrients/pollutants can be quickly flushed out of the estuary through inlets, some of them will remain in the estuary for a long time. If the nutrient/pollutant loading rate is faster than the estuarine natural flushing rate, the estuarine water quality will degenerate gradually. Hence, in order to improve estuarine water quality, we need to first understand estuarine hydrodynamics and flushing characteristics.
Generally speaking, flushing is the capability of an estuary or a part of an estuary to replace its old water mass with new ocean water mass. It is conceptually very simple but
2
complex to calculate in reality. There are many methods that have been used to quantify estuarine flushing characteristics. Traditionally volumetric driven methods were used to calculate flushing characteristics (Ketchum 1951, Stommel and Arons 1951, Dyer and Taylor 1973). Despite their simplicity, these methods have several drawbacks and are not reliable to represent the real world situation. First, most of these methods assume complete mixing. Second, they consider a single principal tidal constituent but exclude wateI: movements caused by other tidal constituents and nontidal forces. Third, estuarine hydrodynamics are not incorporated. More recently the hydrodynamic driven method (Edinger et al. 1998, Sheng et al. 1995), which considers the estuarine hydrodynamics, is being used in favor of the traditional volumetric driven methods.
1.2 Scope of The Study
As part of the Indian River Lagoon Pollutant Load Reduction (IRLPLR) model project sponsoredbythe St. John's River Water Management District (SJRWMD), this study focuses on the quantification of the Indian River Lagoon flushing characteristics, using a 3D hydrodynamic model CH3D (Sheng 1986a, 1986b, 1987, 1989). Here the Indian River Lagoon (IRL) is a linked system of three water bodies: Mosquito Lagoon, Banana River Lagoon and Indian River Lagoon, not only the Indian River itself. A relatively fine boundaryfitted curvilinear grid (477x43x4) for the Indian River Lagoon was used in order to fit the complex geometry and bathymetry. A thirdorder accurate and oscillationfree advection scheme has also been developed for CH3D to calculate the flushing time of a passive tracer.
3
Because of IRL's large size and complex geometry, flushing characteristics are very different from place to place. Some regions in the lagoon may have good flushing and some may have very poor flushing. A reasonable segmentation scheme, which reflects the spatial variation in flushing characteristics in IRL, is necessary. Sheng et al. (1995) calculated the IRL flushing characteristics by using the verticallyintegrated version ofCH3D and provided a preliminary segmentation scheme with eight segments. However, Sheng et al. (1995) only considered the flushing due to tidal circulation forced by simple harmonic tides at the tidal inlets. This segmentation scheme will be compared with that produced from the seagrass perspective (Vimstein et al. 1994). A ninesegment scheme has been provided and more accurate flushing calculations for the entire lagoon and each individual segment will be conducted by using the fully threedimensional version of CH3D with a highly accurate advection scheme and considering realistic circulation driven by tide, wind and density over a one year period.
In summary, the scope of the study includes:
* Apply the ULTIvLATEQUICKEST advection scheme (Leonard 1991) to the threedimensional nonorthogonal curvilinear coordinate system.
0 Verify the scheme in nonorthogonal curvilinear coordinate system by conducting
various numerical tests.
0 Develop CH3D flushing model with the ULTIMATEQUICKEST advection scheme.
* Apply CH3D flushing model to the study of Indian River Lagoon flushing
characteristics.
* Analyze flushing results and compare different Indian River Lagoon segmentation
4
schemes and propose an optimal segmentation scheme.
Chapters two and three present the background knowledge for this study. The definition of flushing, different flushing methods, and previous flushing and segmentation studies of Indian River Lagoon will be reviewed in these two chapters. Chapters four and five cover the numerical model and test results, including CH3D hydrodynamic model, flushing model, the ULTIMATEQUICKEST scheme with CH3D, and scheme verification. Chapters six and seven present various IRL segmentation schemes and flushing results. Conclusions and recommendations are presented in chapter eight.
CHAPTER 2
FLUSHING IN ESTUARIES
2.1 Definitions and Classifications
An estuary is where the river meets the sea (Fischer et al. 1979). Pritchard (1967), based on his study on James River in Virginia, defined an estuary as "a semienclosed coastal body of water which has a free connection with the open sea and within which sea water is measurably diluted with fresh water derived from land drainage." This definition, however, does not work very well with all kinds of estuaries in the world. In Florida, for example, most of the estuaries are in the form of coastal lagoons, e.g., Indian River Lagoon, Sarasota Bay, and Biscayne Bay, etc. Kjerfve (1994) defined the coastal lagoon as, "a shallow coastal water body separated from the ocean by a barrier, connected at least intermittently to the ocean by one or more restricted inlets, and usually oriented shoreparallel."
Since almost half of the world's population are now living around estuaries, growing human activities, (e.g., effluent disposal, dredging of navigational channels, reclamation of wetlands or open water, agriculture runoff, etc.) are causing adverse impacts on the estuarine environments.
One of the major concerns nowadays is the impact of pollutant loadings on estuarine water quality, habitat, and living resources. Understanding estuarine flushing is essential to
6
understanding the impact of pollutant loadings on the estuarine ecosystem. Flushing
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 hixing.
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=(0S (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 onedimensionally from the head of the estuary down to the mouth of the estuary, and assumes that the high tide and low tide are in phase for all segments. The inner segment (segment 0) is defined based on the assumption that the intertidal volume (tidal prism) in segment 0, 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 + [(1a)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 EquationBased Estuarine Flushing Model" and applied it to a few Gulf of Mexico estuaries. Focusing on "tidal flushing" rather than "river flushing," this method incorporates the functionality of the "tidal prism method," though it still uses fresh water as the tracer. By tracking a single input of fresh water over several time steps (tidal periods) until it is flushed out of the estuary, the flushed fresh water volume is defined as:
+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"=VnI 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 backflushing from the mouth based on the following equation,
Vt,=P(1 )SO (2.8)
35
where, Vb is the fresh water volume reentered to the estuary by backflushing, P is the tidal prism, So is the salinity at the mouth, the ocean water salinity is assumed as 35 ppt. This method is simple and can apply to the flushing calculations of coastal lagoons. However, it does not include the estuarine hydrodynamic process. The narrow focus on tidal flushing also excludes many other nontidal flushing features.
2.2.5 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 crosssection averaged velocities to drive the model; (3) Besides tide and/or river discharge (fresh water input), none of other major driving forces, such as wind and density gradients, is considered in these methods; (4) Most of these methods assume steady and uniform estuarine conditions, and instant and complete mixing in estuarine sections; (5) Some of the older methods ignore the effect of tidal backflushing; (5) These models are onedimensional models along the head of the estuary to the mouth of the estuary;
(6) Estuarine hydrodynamics are not incorporated in either of these methods. Because of the narrow focus and unrealistic assumptions, the results based on some of theses methods become very unreliable. For example, according to the "freshwater flushing rate method," it takes 715 days to replace Tampa Bay's fresh water and in contrast, only 5 days are required to flush Tampa Bay according to the "tidal prism flushing method" (Biggs 1986).
More recently, with the development of more powerful computers, the hydrodynamic
12
method is becoming more popular (Sheng et al. 1995, Edinger et al. 1998). This method requires a quantitative understanding of the estuarine hydrodynamics under the influence of multiple forcings. Usually the continuity equation and the Reynold's averaged NavierStokes equations of motion need to be solved in order to obtain the surface elevation and the flow field. Then the estuarine flushing time (e.g., 50% renewal time) can be determined by tracking the time history of relative remaining tracer mass contained in any segment of the estuary (or the entire estuary) which is initially filled with a uniform tracer concentration (e.g., 100). The conservative advectiondiffusion equation needs to be solved along with the equations of motion until the total tracer mass in the segment (or the entire estuary) drops to below 50% of its initial total mass. Rather than considering one tidal cycle as the time step for estimating the flushing time, this method solves the transport equation, the equations of motion, and the continuity equation simultaneously at a much smaller time step (on the order of seconds or minutes). At each time step the tracer concentration will be calculated according to the updated new surface elevation and flow field. The onedimensional, twodimensional or threedimensional model can be used. The flushing results are very sensitive to the numerical schemes used for solving the advectiondiffusion equation. Traditional first order upwind scheme usually overestimates the flushing rate due to excessive numerical diffusion. Certain higher order schemes may give negative concentration because of numerical dispersion. An accurate and oscillationfree advection scheme is needed for flushing calculation, especially in those regions where the flows are advectiondominated.
Table 2.1 summarizes all of the methods for calculating flushing time in estuaries.
13
Table 2.1 Comparison between different estuarine flushing calculation methods
Methods Driving Forces Time Dimension Mixing Tidal Back Estuarine
step Type flushing hydrodynamics
Tidal Prism Harmonic tide, I tidal OD instant and not not
Estuarine cycle ("box" complete considered considered
volumes, Tidal model) in estuary
Prism
Fraction of Harmonic tide, I tidal OD 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 ID 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 OD considering considered not
Equation Estuarine cycle ("box" vertical considered
volumes, Tidal model) stratification
prism, Fresh by simple
water input adjustment
of
coefficient
Hydro Real tide, seconds 1, 2 or 3D Eddy considered considered
dynamic wind, river or viscosity
discharge, minutes and
density diffusivity
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 24 km wide) and shallow (typically 13 m deep) estuary with four inlets connecting it to the Atlantic Ocean and numerous rivers and canals discharging fresh water into it. The northern end of Mosquito Lagoon is connected to the ocean through Ponce de Leon Inlet. The southern part of Indian River Lagoon is connected to the ocean by Ft. Pierce Inlet and St. Lucie Inlet. Sebastian Inlet connects the middle part of Indian River Lagoon to the ocean. The southern end of Mosquito Lagoon is connected to Indian River Lagoon through Haulover Canal. Banana River is an almost enclosed water body with only two narrow connections with the Indian River Lagoon. Indian River Lagoon is a micro tidal system with the semidiurnal M2 tide as the principal constituent (Smith 1987). Tide is the major force causing flushing for most of the Indian River Lagoon, except in the northern IRL and Banana River, where tidal action is weaker. The geometry of the Indian River Lagoon is very complicated because of highly irregular boundaries, barrier islands, and numerous causeways. Indian River Lagoon has a mean depth
15
of about two meters, with a maximum depth of nearly ten meters and a minimum depth less than half meter. Dredging of the intercoastal waterway has added to the complexity of the bathymetry. Due to the complex geometry and bathymetry, the flushing characteristics are very different within different segments of the Indian River Lagoon. Segments with strong tidal action and less constriction usually have higher flushing rates, while segments with weak tidal action and more constriction have slower flushing rates. Based on the natural boundaries and the preliminary flushing study (Sheng et al. 1995), Indian River Lagoon was divided into eight segments (Figure 3.1). Table 3.1 gives an outline of some general information about Indian River Lagoon and its eight segments.
3.2E+06
3.15E+06
 \ (3)
II
3.1E+63.05E+06
>.' :(8)
3E+06'
500000 525000 550000 575000 600000
X UTM (m) Figure 3.1 The Indian River Lagoon eightsegment scheme by Sheng et al..
17
Table 3.1 Outline of the general information about Indian River Lagoon
Segment 1 2 3 4 5 6 7 8 Total
North 405 516 56 free 67 free 78 free Boundary shore shore shore highway highway connection connection connection
South 405 3 516 56 free 67 free 78 free Boundary shore highway highway highway connection connection connection shore
West
Boundary shore shore shore shore shore shore shore shore
East
Boundary shore shore shore shore shore shore shore shore
Ponce de Sebastian Ft Pierce
Tidal Inlet Leon Inlet     Inlet  Inlet
St Lucie
Inlet
Annually 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 (IRLPLR) model is being developed by the University of Florida (Dr. Y. P. Sheng is the Principal Investigator), with funding from the St. Johns River Water Management District (SJRWMD). The project aims to develop a Pollution Load Reduction (PLR) model which incorporates a circulation model, a sediment transport model, a water quality model, a light model and a sea grass model. An important first step is to develop a segmentation scheme and an understanding of the flushing characteristics. As part of the PLR model development effort, this study will provide an accurate calculation of flushing characteristics within the entire lagoon and each of the nine segments determined in this study. The flushing model sensitivity study was also conducted.
3.3 Previous Indian River Lagoon Segmentation Studies
Smith (1993, 1996) used a onedimensional lagoon wide model to calculate the flushing time in Indian River Lagoon. He divided Indian River Lagoon into 16 segments (excluding Mosquito Lagoon and Banana River) and Banana River into 3 segments. He then quantified the intertidal volume by using harmonic constituents (principal semidiurnal and diurnal) recorded at 31 NOS stations. He multiplied the tidal amplitudes by the surface areas, and incorporated the phases to represent the tidal volumes and the tidal movements. The model first estimated the averaged onedimensional flow between the segments and then solved a onedimensional advectiondiffusion equation. In his 1996 version model, he added "Flexible Grid" technique which allows control of the extent to which segment boundaries move along the longitudinal axis of the lagoon with the flow. If the boundaries can move freely with the flow, the transport equation is solved with diffusion only, and if the boundaries are fixed, then the transport equation is solved with advection only. Based on one year simulation he found that 50% renewal time is about 140 days if transport by advection (fixed grid), and 50% renewal time is not reached until the end of the simulation if transport by diffusion (fully flexible grid). The results suggest that it tends to underestimate flushing if using fully flexible grid and overestimate flushing if using fixed grid. Another series of simulations using completely flexible grid show that in the southern and central subbasins of Indian River Lagoon 50% renewal times were about 5 and 12 days respectively while the northern subbasin of Indian River Lagoon and the Banana River reached only about 30%
20
renewal by the end of one year simulation. Although the method he used is comparable to the traditional modified tidal prism method, his results still can be used as a rough indication of the IRL flushing characteristics.
Sheng et al. (1993a, 1993b) used a onedimensional lagoonwide model and a threedimensional lagoon wide model to calculate the flushing in the Indian River Lagoon. They found that the onedimensional model tended to overestimate the flushing rate. Further, as part of the preliminary study for IRLPLR model, Sheng et al. (1995) conducted a preliminary study of the Indian River Lagoon segmentation scheme. A vertically integrated twodimensional version of the threedimensional curvilineargrid hydrodynamic model CH3D was used to calculate the IRL flushing. Because of lack of data, they assumed semidiurnal M2 tides with amplitudes of 5 cm, 15 cm and 30 cm at Ponce, Sebastian and Ft. Pierce inlets respectively. Without considering any other forcing functions, they simulated the flushing characteristics corresponding to four segmentation schemes. Based on the flushing results and the hydrodynamic and water quality data from 24 FDEP or USGS monitoring stations, they found that the eightsegment scheme (Figure 3.1) appeared to be the best scheme.
In 1994, St. Johns River Water Management District (SJRWMD) conducted a spatial analysis of seagrass changes versus changes in nutrient loadings over the past 50 years (Virnstein et al. 1994). Based on changes of seagrass coverage and drainage subbasin boundaries, they divided Indian River Lagoon into 18 segments (Figure 3.2). Then they calculated change of sea grass coverage and changes of nutrient loadings for each segment from 1943 to 1992, and regressed percentage change of seagrass coverage on percentage changes of nutrient loadings. The results as shown in Figure 3.3 indicate that most of the
21
segments have significant increases in nutrient loadings and corresponding decreases in seagrass coverage during the past 50 years except segment 15 where the fast flushing rate appears to have flushed increased nutrient loadings into the ocean.
In this study, simulation of flushing characteristics will be first conducted with eight segments developed by Sheng et al. (1995), and finally with the nine segments developed based on the result of comparison between eightsegment (Sheng et al. 1995) and eighteensegment schemes (SJRWMVD 1994). The ninesegment scheme will be described in following chapters of this thesis.
3.2E+06
550000
X UTM (m)
Figure 3.2 The Indian River Lagoon 18segment 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 Ploading,
by segrhent
0 100 200 300 400 500 12
3
% Change in TSSloadln
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 threedimensional momentum equations can be derived from Newton's second law of motion which is based on the concept of conservation of momentum. The continuity and scalar transport equations can be derived from the mass conservation criteria. With the assumptions of incompressible flow, vertically hydrostatic water column, Boussinesq approximation, and turbulent eddy viscosity concept, these equations were expressed in the following forms (Sheng 1989):
au+ 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 verticallystretched and horizontally boundaryfitted system. These equations were derived by Sheng (1986) and are shown in Appendix A.
4.2 Review of Advection Schemes
Numerical simulation of circulation and transport processes is very important for many engineering and scientific studies in coastal waters, lakes and estuaries. In particular, when quantifying the transport of pollutants and water quality parameters, the accuracy of
the advection scheme is very crucial.
In the past decades, many numerical schemes have been developed and applied to the study of tracer advection. Probably the simplest scheme is the first order upwind scheme. However, fewer people are using it now because of its high numerical diffusion. Another popular scheme is the Leap Frog scheme which is used in the princeton Ocean Model (POM). This scheme does not have as much dissipation as the up 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, thirdorder accuracy with little dissipation and dispersion.
Many numerical schemes have been developed recently to eliminate spurious oscillations while maintain high order accuracy. One of the schemes is the TVD (Total Variation Diminution) scheme which used socalled "flux limiting" technique, e.g., Roe (1985). Although TVD scheme is totally free from numerical overshoots and undershoots near discontinuities, some studies showed that the scheme also has many disadvantages such as "clipping" of local extrema and "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 onedimensional and twodimensional flows and found that the ULTIMATE scheme
27
was more general than the other schemes and was easier to apply. Lin et al. (1997) also applied this scheme in his twodimensional transport model and simulated the water quality and sediment transport in Humber Estuary, U. K..
The next section will describe how the ULTIMATE algorithm is incorporated in CH3D.
4.3 ULTIMATEQUICKEST Scheme
The ULTIMATE algorithm can be used to explicit conservative finite difference schemes of any order accuracy. The higher order of accuracy, the more interpolation points will be involved and the less efficiency of computation. For most large scale and relatively smooth flows the costeffective thirdorder upwind (QUICKEST) ULTVIATE scheme gives excellent practical results (Leonard 1991). The following sections will only present the formulation of the ULTIMATEQUICKEST scheme for the threedimensional nonorthogonal curvilinearcoordinate numerical model CH3D.
4.3.1 Third Order Transient Interpolation Model
Considering unsteady onedimensional pure advection of a scalar 4 at constant velocity u, the equation is
.___.+u__ =0. (4.9)
at ax
Applying centerspace 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 wellknown schemes is the LaxWendroff scheme. The right face value for the LW scheme can be written as
r 2=[( 1)( 144)]. (4.14)
For QUICKEST, the third order upwind interpolation scheme (Leonard 1979), the right face value can be written in terms of the second order LaxWendroff value and an upwind biased curvature (Leonard 1991),
(3)( =(2)( (1c2) 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.
i1 i i+1 i1 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 lefthand side inequality in equation (4.18) is less restrictive than cIt' n+l. (4.19)
For the righthand 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 nonuniform 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)+ nn
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 Onedimensional 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 ULTIMATEQUICKEST Algorithm
Based on the knowledge of the third order transient interpolation model and the universal limiters, the ULTIMATEQUICKEST algorithm can be established as following:
* Designate upstream(u), downstream(d), and central(c) nodes on the basis of SGN(uf)
for each control volume face.
* Compute DEL= g1, ADEL=IDELI and ACURV=I d2 +$I for each
face. If ACURV>ADEL(nonmonotonic), 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 ULTIMATEQUICKEST 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 a9 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
nondimensionalized 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, = tHx,f gO)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 onedimensional limiter method described above can be directly applied to threedimensional 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 righthand 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 onedimensional QUICKEST interpolation formula (equation 4.15) can be directly used in CH3D to obtain the control volume face values if considering three direction sweeps separately. The only difference is that the Courant numbers have to involve the dimensionless coefficients because the equations in CH3D are dimensionless.
The QUICKEST right face value of control volume in three directions for CH3D can be expressed as:
n (n n
(1C2 1+SGN(c ,.2 1SGN(c,
6 2 2 (
6 2 n+12 }+#11)+ 2 n+220+ 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
(1c 1 +SGN(c) 1 S GN(c,)
6 2 2n n ng+ n n(
 2 2.(4 .372)k + 22,+,+(*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 onedimensional ULTIMATEQUICKEST scheme to the threedimensional CH3D. Simply implement three sweeps in E, lI, and a directions consecutively while updating the central value within each sweep. This method solves only onedimensional advection equation at every intermediate time step, and can make better use of the onedimensional limiters. The advantage of this method is its accuracy and its applicability on relative large Courant number conditions (C < 1).
Nonsplitting method. The basic procedure is similar but with only one sweep. 1) Apply onedimensional ULTIMATE strategy described above in E, 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 threedimensional advection equation at a time using the onedimensional limiters. The Courant number condition for this method is more restrictive than that for splitting method:
42
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 onedimensional pure advection tests which showed the advantages of the ULTIMATE algorithm against other schemes. This chapter will focus on two and threedimensional tests in a nonorthogonal curvilinear coordinate system. A test grid was set up and five bench mark tests (advection of a cubic hill, a cone hill, and a Gaussian hill in twodimensional uniform flow, advection of a column in twodimensional circular flow, and advection of a cubic block in threedimensional uniform flow) were conducted. Results of the splitting and nonsplitting ULTIMATEQUICKEST methods were also compared with different Courant numbers.
5.1 Test Grid System
The test domain is a rectangular basin with a uniform depth of 6000 m, and is divided into 100x100 grid cells with i=1 to 101 in X direction andj=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 gridlines are established in order to test the new scheme in the nonorthogonal curvilinear grid system. The grid points on the four boundaries are defined as follows:
(1) Along the southern boundary (Y=0): xi=,x(i 1), 1< i< 31;
xi=x1_+ x[0.5+(51i)/20], 321i 51;
xi=x51+^x(i51), 52gik 101, (5.1)
(2) Along the northern boundary (Y=Y ): xi=tx(i1), 1: is 31;
x,=xt_1+Ax[0.5+(i32)/20], 32 i 51;
xi=x51'+x(i51), 52
(3) Along the western boundary (X= 0): yj =AY( 1), 1
yj=yj _+Ay[0.5+(51j)/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+(f32)/20], 32
yj=y51 +Ay(f 51), 52
where Ax1200m and y=1000m.
Thus, the grid is orthogonal and uniform with Ax=1200m and &y=1000m in the following four regions: (1) 1
Four twodimensional and one threedimensional bench mark tests are selected to verify the new scheme.
0
25000 50 00 75 00 100000 X (m)
Figure 5.1 The 100x100 nonorthogonal curvilinear test grid system.
5.2 Twodimensional Tests
5.2.1 Test 1: Advection of A Cubic Hill in Uniform Flow
The initial concentration distribution of this test is defined as:
If =1.0, Ixxol0jdx and lyyoj: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 nonsplitting ULTIMATEQUICKEST) were tested. Figures 5.2 5.7 give the results after 1800 time steps with At=200s (100 hours).
Theoretically, the pure advection of any type of profile will give the exact same shape of the profile as the initial one. As seen in Figure 5.2, the cubic hill should keep its original shape after 100 hour advection in uniform flow. Numerically it is hard to obtain the exact same profile. Different numerical schemes give different accuracy. In this test the upwind scheme gives the lowest accuracy due to its severe numerical diffusion. Although the QUICKEST scheme keeps the higher order accuracy, it cannot avoid the spurious oscillations at sharp gradient regions (shown in figure 5.4). Only the ULTIMATEQUICKEST scheme achieves the higher order accuracy without any undershooting and overshooting (Figure 5.5).
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 ULTIMATEQUICKEST 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. 2D 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. 2D 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 + (yo) /Ro, (x xo)2 + yyo)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 nonsplitting ULTIMATEQUICKEST) were tested. Figures 5.8 5.13 give the results after 1920 time steps with At=200s (about 106.7 hours).
52
This test was designed to show the performance of the schemes on local extremum. Due to the high diffusion, the profile is severely damped and spreads around if the upwind scheme is used (Figure 5.9). The QUICKEST scheme maintains the third order accuracy, but unfortunately it cannot avoid the undershooting at the base of the profile (Figure 5.10). The ULTIMATEQUICKEST scheme seems to have the best performance. It removes the spurious oscillation while maintaining the same order of accuracy as QUICKEST scheme (Figure 5.11).
80
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 ULTIMATEQUICKEST 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. 2D 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. 2D 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 nonsplitting ULTIMATEQUICKEST) were tested. Figures 5.14 5.19 give the results after 1800 time steps with At=200s (100 hours).
58
Still, in this test the upwind scheme is proved to be with high damping. It seems that for the smoothly changed profile, both the ULTIMATEQUICKEST and the QUICKEST schemes have similarly good performance. Sometimes for a very steep peak value like in this test case, the "clipping of the peak value" is difficult to avoid when using the 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
20m
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 ULTIMATEQUICKEST scheme after 100 hours. The center of the Gaussian hill is initially located at xo=1.62x104m and yo=1.35x104m, which is close to the origin.
Figure 5.18 Comparison of numerical and analytical solutions for the advection of a Gaussian hill in uniform flow. 2D projection along the x axis through the center of the analytical solution. The ULTIMATEQUICKEST scheme was used.
Figure 5.19 Comparison of numerical and analytical solutions for the advection of a Gaussian hill in uniform flow. 2D projection along the y axis through the center of the analytical solution. The 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, (xxo)2+(Yyo0)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(yyc), v=W(xxc). Three numerical schemes (the upwind, QUICKEST and nonsplitting ULTIMATEQUICKEST) were tested. Figures 5.20 5.25 give the results after 1800 time
steps with At=200s (100 hours).
Since all the tests described above are advections in uniform flow. This test was designed to see the performance of different schemes in nonuniform flow. Again the test results show that the ULTIMATEQUICKEST scheme appears to be the best scheme.
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
= 424 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 ULTIMATEQUICKEST 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. 2D 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. 2D projection along the y axis through the center of the analytical solution. The ULTIMATEQUICKEST scheme was used.
5.3 Threedimensional Test
5.3.1 Test 5: Advection of a Cubic Block in 3D Uniform Flow
The initial concentration distribution for this test is:
=1.0, Ixxol dx, lyy01o dy, Izz0o 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 nonsplitting ULTIMATEQUICKEST) were tested. Figures 5.26 5.28 give the results after 1560 time steps with at200s (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 ULTIMATEQUICKEST 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 ULTIMATEQUICKEST schemes. A detailed analysis, however, does show significant advantage of the ULTIMATEQUICKEST scheme. Table 5.1 shows the analysis of the maximum percentage of over and under prediction and the error ratio of mass conservation. There is no overshooting or undershooting for the ULTIMATEQUICKEST scheme. Table 5.1 Comparisons of numerical schemes for Test 5: The advection of a cubic block in uniform flow
Numerical Scheme Erroe, Err under Err..
1 upwind 44.10% 0.0% 2.85 x106 %
2 QUICKEST 24.8% 7.7% 2.25 x106 %
3 ULTIMATEQUICKEST 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 MassTotal 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 ULTIMATEQUICKEST Scheme With Splitting
The bench mark tests using the nonsplitting ULTIMATEQUICKEST method seem to give us very encouraging results. Further study shows that they are only applicable to very small Courant numbers (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 nonsplitting ULTIMATEQUICKEST 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 ULTIMATEQUICKEST method. (Cx)max=(Cy)max=0.45.
Based on these tests, It seems that there is little difference between the results of splitting and nonsplitting schemes when the Courant number is very small (C<0.1). As the
75
Courant number becomes larger and approaches unity, the results of splitting method are much better than that of the nonsplitting method.
5.5 Summary of Tests
This chapter has described several numerical test problems: (1) Advection of a cubic hill in 2D uniform flow; (2) Advection of a cone hill in 2D uniform flow; (3) Advection of a Gaussian hill in 2D uniform flow; (4) Advection of a column in 2D circular flow; and
(5) Advection of a cubic block in 3D uniform flow. These test problems represent different conditions of the real world. For each test problem, three advection schemes: (1) upwind, (2) QUICKEST, and (3) ULTIMATEQUICKEST schemes were compared. The test grid is specially designed for verification of the schemes based on the nonorthogonal curvilinear coordinate system. For all the five test problems, the upwind scheme has too much diffusion. The QUICKEST scheme cannot avoid spurious oscillation at sharp gradient regions. Only the ULTIMATEQUICKEST scheme seems to have the best performance. For the ULTIMATEQUICKEST scheme, both the splitting and nonsplitting method were implemented. Further tests with different Courant numbers show that both the Splitting and nonsplitting methods work well if the Courant number is less than 0.1. However, when the Courant number is greater than 0.1 and less than unity, the nonsplitting method is proved to be unacceptable due to the severe deformation of the original test profile, while the splitting method can still keep the same good result no matter what Courant number it is. Thus the splitting ULTIMATEQUICKEST scheme will be used in this study.
CHAPTER 6
INDIAN RIVER LAGOON FLUSHING CALCULATION
6.1 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 semidiurnal tides with T= 12.4 hours, the flushing time can be obtained as, t = 5.65 days, thus, the R50 value is less than 2.83 days.
The second method used here is the Fraction of Fresh Water Method. The total river discharges into the 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 ULTIMATEQUICKEST advection scheme to the simulation of the flushing characteristics of Indian River Lagoon, an accurate grid system has to be set up. In order to represent the Indian River Lagoon's complex geometry and bathymetry, a high resolution boundaryfitted grid is necessary. By applying Thompson's (1983) grid generation method, Sheng et al (1995) generated a nonorthogonal curvilinear grid for Indian River Lagoon (Figure 6.1). The grid is of very high resolution with horizontally 477x43 cells (I=1 to 477 in the northsouth direction and J=1 to 43 in the 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 onedimensional advection equation is used. At tidal boundaries, tidal back flushing was considered by calculating inflow and outflow advection fluxes based on the direction of velocities. At river boundaries, only inflow condition is considered, because river discharge is always positive. Details of inflow and out flow conditions can be found in Appendix B.
Water surface elevation data. The Indian River Lagoon estuarine system has four tidal inlets connecting the system to the Atlantic Ocean. As one boundary condition, the surface elevation at an inlet is given by a time series of measured data. These data were collected at four 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 oneyear period (September 1997 December 1998). It is interesting to see that there is a declination of mean sea level from north to south along the shore. According to the data, the one year averaged mean sea level at the Ponce de Leon Inlet is about 11 cm higher than that at the Ft. Pierce Inlet.
Wind data. The free surface boundary condition requires the wind stress at every grid 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 hourlyaveraged wind speed at X (eastern is positive) and Y (northern is positive) directions respectively at the five stations over a period of more than one year (September 1997 December 1998).
River discharge data. A total of 17 rivers are considered for the system. The locations of these rivers can be seen in Figure 6.1 and Table 7. 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 River1 269 3
6 Eau Gallie River2 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 169173 21
2 Northern Indian River 176 199 20
3 203212 20
4 215234 20
5 BananaRiver 177 199 28
6 204213 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 8721147
3.2E+06
USGS 02248380
FDEP
8721456
3.15E+06
FDEP 8721789
I
3.1 E+06
3.05E+06
FDEP 8722213
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 8721147 (Ponce Inlet) from September 1 (Julian Day
244), 1997 December 31 (Julian Day 365), 1998.
Station 8721147x (EastWest) I
, I I IFDEP E
