Modeling the circulation and water quality in Charlotte Harbor Estuarine System, Florida

MISSING IMAGE

Material Information

Title:
Modeling the circulation and water quality in Charlotte Harbor Estuarine System, Florida
Physical Description:
xix, 325 leaves : ill. ; 29 cm.
Language:
English
Creator:
Park, Kijin
Publication Date:

Subjects

Subjects / Keywords:
Coastal and Oceanographic Engineering thesis, Ph. D   ( lcsh )
Dissertations, Academic -- Coastal and Oceanographic Engineering -- UF   ( lcsh )
Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 2004.
Bibliography:
Includes bibliographical references.
Statement of Responsibility:
by Kijin Park.
General Note:
Printout.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 003100805
System ID:
AA00008994:00001

Full Text










MODELING THE CIRCULATION AND WATER QUALITY
IN CHARLOTTE HARBOR ESTUARINE SYSTEM, FLORIDA












By

KIJIN PARK


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


2004
















ACKNOWLEDGMENTS

I would like to thank my advisor, Dr. Y. Peter Sheng, for his guidance, support and

financial assistance throughout my study. In addition, much appreciation is owned to my

other committee members, Dr. Robert G. Dean, Dr. Robert J. Thieke, Dr. Louis H. Motz and

Dr. K. Ramesh Reddy, for their review of my dissertation.

I would like to thank the South Florida Water Management District and Southwest

florida Water Management District for sponsoring research project and providing data for

Charlotte Harbor estuarine system. I would like to express my thanks to Justin, Jeffery,

Yangfeng, Taeyoon, Jun, Vadim, and Vladimir whose help with class, research and writing

this dissertation. Many thanks go to Becky, Lucy, Sonna, Sidney, Kim and Halen for making

life easier.

I would like to dedicate this dissertation to my parents whose love and support made

this degree possible. Last, but not least, I would like to thank my wife, Kyung-Mi, who have

been praying for me to be faithful, kind, and honest all the time.














TABLE OF CONTENTS

Page

ACKNOWLEDGMENTS ........... ................................... ii

LIST OF TABLES ........................ ......... ................ vi

LIST OF FIGURES .................................................... ix

ABSTRACT ......................................................... xvii

CHAPTER
1 INTRODUCTION .......................................... ............ 1

2 CHARLOTTE HARBOR CHARACTERIZATION .......................... 9

2.1 Clim ate ....................................... .............. 11
2.2 Hydrodynamics ............................... ................ 11
2.2.1 Tidal Stage, Discharge at Inlet, and Tidal Circulation .......... 13
2.2.2 Freshwater Flow ....................................... 14
2.2.3 Salinity .............................................. 14
2.3 W ater Quality ................................... .............. 15
2.3.1 Nutrients ............................................. 16
2.3.2 Dissolved Oxygen ................................... ... 18
2.3.3 Phytoplankton ......................................... 19
2.4 Sedim ent ............................... ..................... 20
2.5 Light Environment ............................................. 21

3 HYDRODYNAMIC AND SEDIMENT TRANSPORT MODEL ................ 23

3.1 Governing Equation .............................. .............. 24
3.1.1 Hydrodynamic M odel ................................... 24
3.1.2 Sediment Transport Model ............................... 26
3.2 Boundary and Initial Conditions ................................. 27
3.2.1 Boundary Conditions ................................... 27
3.2.2 Initial Conditions ...................................... 30
3.3 Heat-Flux at Air-Sea Interface ................................... 30
3.3.1 Short-W ave Solar Radiation .............................. 31
3.3.2 Long-W ave Solar Radiation .............................. 32

iii









3.3.3 Sensible and Latent Heat Fluxes .................... .... 33

4 WATER QUALITY MODEL ............................................ 36

4.1 Mathematical Formulae ......................................... 38
4.2 Phytoplankton Dynamics ........................................ 40
4.2.1 M odeling Approach ..................... .............. 40
4.2.2 Relationship between Phytoplankton and Nutrients ............ 42
4.3 Nutrient Dynamics ............................................. 43
4.4 Oxygen Balance ......................... ...................... 45
4.5 Effects of Temperature and Light Intensity on Water Quality Model ...... 59
4.5.1 Tem perature .........................................59
4.5.2 Light intensity ........................... .............. 61
4.6 Light Attenuation M odel ........................................ 61
4.7 Model Parameters and Calibration Procedures ....................... 69

5 APPLICATION OF CIRCULATION AND TRANSPORT MODEL ............. 81

5.1 A High-Resolution Curvilinear Grid for Charlotte Harbor Estuarine System 82
5.2 Forcing Mechanism and Boundary Condition ........................ 86
5.3 Simulations for 1986 Hydrodynamics .............................. 94
5.3.1 Sensitivity and Calibration Simulations ................ .... 99
5.3.2 Results of the 1986 simulation ...................... ... 101
5.4 Simulations for 2000 Hydrodynamic .............................. 118
5.4.1 Sensitivity and Calibration Simulations ............... ... 118
5.4.2 Results of the 2000 Simulation .................... .... 130
5.4.3 Application of the 2000 hydrodynamic Simulations .......... 146

6 APPLICATION OF WATER QUALITY MODEL ......................... 167

6.1 Forcing Mechanism and Boundary Condition for Circulation .......... 167
6.2 Initial and Boundary Condition for the Water Quality Model ........... 172
6.3 Simulations of W ater Quality in 1996 ............................. 180
6.3.1 Calibration .............................. ............. 180
6.3.2 Results for 1996 Water Quality Simulations ................ 191
6.4 Simulations of W ater Quality in 2000 ...................... .... 212
6.4.1 Verification ......................... ................. 212
6.4.2 Results of 2000 Water Quality Simulations ................. 215
6.4.3 Application of 2000 Water Quality Simulations .............. 235

7 CONCLUSION AND DISCUSSION .................................... 258

APPENDIX
A FLOW CHARTS FOR CH3D_IMS .................................... 263









B DIMENSIONLESS EQUATIONS IN CURVILINEAR BOUNDARY-FITTED AND
SIGM A GRID .................................................. 278

C COMPARISON OF WATER QUALITY MODELS ........................ 284

D NUMERICAL SOLUTION TECHNIQUES FOR WATER QUALITY MODEL.. 289

E NUTRIENT DYNAMICS ............................................. 294

F SEDIMENT FLUX MODEL ......................................... 306

G MODEL PERFORMANCE TEST WITH PARALLEL CH3D_IMS ........... 309

H TIDAL BENCH MARKS FOR CHARLOTTE HARBOR ................... 313

REFERENCES ....................................................... 314

BIOGRAPHICAL SKETCH ............................................. 325














LIST OF TABLES


Table page

1.1 Component models of the CH3D-IMS. ................................... 6

2.1 Ratios of nitrogen and phosphorous constituents ........................ 18

3.1 Mean latitudinal values of the coefficient .............................. 33

4.1 Average values of oxygen uptake rates of bottom ........................... 54

4.2 The spectrum of incident sunlight data. .............................. 64

4.3 Coefficient ranges to use in stand along light model ........................ 67

4.4 The best fit light model coefficients for Charlotte Harbor estuarine system ..... 67

4.5 Dixon and Gray's model coefficients for the Charlotte Harbor estuarine system.. 67

4.6 Temperature adjustment functions for water quality parameters. ............... 70

4.7 Water quality parameters related to conversion rate ......................... 71

4.8 Water quality parameters related to phytoplankton and zooplankton. ........... 71

4.9 Water quality parameters in the nitrogen dynamics .......................... 72

4.10 Water quality parameters in the phosphorous dynamics. .................... 73

4.11 Water quality parameters in the oxygen balance ....................... .. 74

4.12 The relationship between water quality parameters and model constituents ...... 79

5.1 Descriptions of 1986 and 2000 river boundary conditions for Charlotte Harbor
estuarine system................................................ 87

5.2 Locations of tidal stage and velocity-salinity measured stations by USGS.. ...... 96









5.3 The effect of removing selected boundary conditions on the accuracy of simulated
water level, velocity and salinity in July 1986. Value shown are average RMS
differences vs. baseline simulation at all data stations. .................... 99

5.4 The effect of varying bottom roughness, z0, on the accuracy of simulated water level,
velocity and salinity in July 1986. Value shown are average RMS errors at all
data stations.. ............................................. 100

5.5 The effect of varying horizontal diffusion, AH, on the accuracy of simulated water
level, velocity and salinity in July 1986. Value shown average RMS errors at all
data stations.. ................................................... 101

5.6 A summary of boundary conditions and model parameters used in the 1986
simulation. ....................................... ............ 102

5.7 Calculated RMS errors between simulated and measured water level for 1986
simulation. ........... ... ................................... 102

5.8 Calculated RMS errors between simulated and measured current velocity for 1986
sim ulation........................................ ............... 104

5.9 Calculated RMS errors between simulated and measured salinity for 1986
simulation. .................................. ................... 105

5.10 The effect of horizontal grid resolution, on the accuracy of simulated water level
and salinity. Values shown are average RMS errors for 2000 calibration at all
available stations. Values shown in parenthesis are % RMS error normalized by
maximum values.. ............................................... 121

5.11 The effect of vertical grid resolution, on the accuracy of simulated water level and
salinity. Values shown are average RMS errors for 2000 calibration at all
available stations. Values shown in parenthesis are % RMS error normalized by
maximum values.. ............................................... 122

5.12 The effect of varying bottom roughness, z0, on the accuracy of simulated water level
and salinity in 2000. Values shown are average RMS errors all data stations. 123

5.13 The effect of varying salinity advection scheme on the accuracy of simulated water
level and salinity in 2000. Values shown are average RMS errors at all data
stations. ...................................... ................. 128

5.14 The effect of modifying bathymetry on the accuracy of simulated water level and
salinity in 2000. Values shown are average RMS errors for all data stations. 129









5.15 A summary of boundary conditions and model parameters used in the 2000
sim ulation. ............ ...... ................................ 130

5.16 Calculated RMS errors between simulated and measured water level for 2000
sim ulation ...................................................... 131

5.17 Calculated RMS errors between simulated and measured salinity for 2000
sim ulation. ..................................................... 133

5.18 The effect of hydrologic alternations on 2000 water level. Value shown are
average RMS differences with baseline simulation for all selected stations. .. 154

5.19 The effect of hydrologic alternations on 2000 salinity. Value shown are average
RMS differences with baseline simulation for all selected stations........... 155

6.1 Locations of water quality measured stations ............................ 173

6.2 Sediment types for Charlotte Harbor water quality simulations. .............. 174

6.3 Water quality parameters, baseline values, and ranges used in the sensitivity
analysis. ....................................................... 183

6.4 Sensitivity analysis results in RMS difference w.r.t. baseline for 1996 water quality
calibration simulation. ............................................ 184

6.5 The water quality model coefficients used for the Charlotte Harbor simulation. .. 185

6.6 The temporally averaged water quality species concentrations for baseline 2000
sim ulation. ..................................................... 240

6.7 Normalized RMS differences of water quality species concentrations at 11 stations
during April to June 2000, showing the effect of no causeway islands. ...... 240

6.8 Normalized RMS differences of water quality species concentrations at 11 stations
during April to June 2000, showing the effect of no ICW. ................ 241

6.9 Normalized RMS differences of water quality species concentrations at 11 stations
during April to June 2000, showing the effect of no causeway islands and ICW241

C. 1 Comparison of water quality models. ................................ 287

G. 1 CPU time for the parallel, shared memory, CH3D procedures on SGI origin
platform ...................................................... 311

H.1 Tidal datum referred to Mean Low Low Water (MLLW), in meter.. .......... 313

viii

















LIST OF FIGURES


Figure Page

1.1 Charlotte Harbor estuarine system and its subarea boundaries. ................. 2

2.1 Drainage basins of Charlotte Harbor estuarine system ....................... 10

2.2 Seasonal wind pattern in Florida. ....................................... 12

2.3 Average monthly concentration of dissolved oxygen in upper Charlotte Harbor, site
CH-006, 1976-84. ................................................ 19

4.1 The connection between nitrogen, phosphorous and carbon cycle .............. 44

4.2 CBOD cycle and DO cycle ............................................ 46

4.3 Comparison of wind-dependent re-aeration formula ........................ 49

4.4 The relationship between POM flux and SOD flux related in the oxidation and
reduction of organic matter in sediment column ........................ 53

4.5 Effect of dissolved oxygen on sediment consumption and SOD release ......... 57

4.6 The scatter plots for kd(PAR) during calibration period with best fit coefficients 68

4.7 Systematic calibration procedure ....................................... 80

5.1 Boundary-fitted grid (92 x 129) used for numerical simulation for Charlotte Harbor
estuarine system .................................................. 84

5.2 Bathymetry in the boundary-fitted grid for Charlotte Harbor estuarine system (92 x
129) ....................................................... .... 85

5.3 Tidal forcing and river discharges for 1986 simulations of Charlotte Harbor
circulation .. ................................................... 89

5.4 Wind velocity for 1986 simulations for Charlotte Harbor circulation ........... 90

ix









5.5 Tidal forcing and river discharges for 2000 simulations of Charlotte Harbor
circulation ...................................................... 9 1

5.6 Wind velocity for 2000 simulation of Charlotte Harbor circulation ............. 92

5.7 Air temperature for 2000 simulation of Charlotte Harbor circulation ........... 93

5.8 Locations of the available 1986 water level and discharge measurement stations of
U SG S .......................................................... 97

5.9 Locations of available 1986 velocity and salinity measurement stations of USGS 98

5.10 Comparison between simulated and measured water level in July 1986 ....... 108

5.11 Comparison between simulated and measure spectra of water level in July 1986 110

5.12 Comparison between simulated and measured current velocity in July 1986 ... 111

5.13 Comparison between simulated and measured salinity in July 1986 .......... 115

5.14 Typical flow pattern of Charlotte Harbor estuarine system during one tidal cycle for
A ugust 6, 1986............... ................................... 116

5.15 The 29-day residual flow and salinity for Charlotte Harbor estuarine system during
July 2 to July 30, 1986. ........................................... 117

5.16 Locations of the available 2000 water level and salinity measured stations at
Caloosahatchee River operated by SFWMD. .......................... 119

5.17 The comparison of the coarse grid (71x92) and the fine grid (92x129) for Charlotte
Harbor estuarine system. .......................................... 120

5.18 A comparison between simulated and measured salinity at Shell Point using
ultimate QUICKEST, QUICKEST, and upwind advection schemes.. ....... 124

5.19 A comparison between simulated and measured salinity at Fort Myers using
ultimate QUICKEST, QUICKEST, and upwind advection schemes.. ....... 125

5.20 A comparison between simulated and measured salinity at BR31 using ultimate
QUICKEST, QUICKEST, and upwind advection schemes................ 126

5.21 Simulated longitudinal-vertical salinity along the Caloosahatchee River at slack-
water before flood tide on September 7, 2000........................... 127

5.22 Comparison between simulated and measured water level in 2000. .......... 136

x










5.23 Comparison between simulated and measured salinity at S79 in 2000. ........ 137

5.24 Comparison between simulated and measured salinity at BR31 in 2000. ....... 138

5.25 Comparison between simulated and measured salinity at Fort Myers in 2000. .. 139

5.26 Comparison between simulated and measured salinity at Shell Point in 2000. .. 140

5.27 Comparison between simulated and measured salinity near Sanibel Causeway in
2000 ........................................................... 14 1

5.28 Comparison between simulated and measured temperature at Fort Myers in 2000.
.................... ...................................... .. 142

5.29 Typical flow pattern of San Carlos Bay during ebb and flood tide for August, 7 on
2000.. ......................................................... 143

5.30 One-year residual flow in San Carlos Bay in 2000......................... 144

5.31 One-year residual salinity distribution in San Carlos Bay in 2000............. 145

5.32 The locations of Sanibel Causeway and IntraCoastal Waterway and stations for
comparing the effect of hydrologic alterations ......................... 148

5.33 The comparison of bathymetry and shoreline for each hydrologic alteration case
scenarios which are Baseline, the absence of IntraCoastal Waterway, and the
absence of causeway ............................................. 149

5.34 The comparisons of water level for three cases at three selected stations: ST05 (Pine
Island Sound), ST08 (San Carlos Bay), and ST10 (Caloosahatchee River mouth).
........................................................... 150

5.35 The comparisons of surface and bottom salinity for three cases at three selected
stations: ST05 (Pine Island Sound), ST08 (San Carlos Bay), and ST10
(Caloosahatchee River mouth) ...................................... 151

5.36 The comparisons of surface residual flow and salinity fields for three cases. ... 152

5.37 The comparisons of bottom residual flow and salinity fields for three cases. ... 153

5.38 The vertical-longitudinal salinity profiles along the axis of the Caloosahatchee River
during wet season in 2000. ........................................ 159









5.39 The vertical-longitudinal salinity profiles along the axis of the Caloosahatchee River
during dry season in 2000. .................................... 160

5.40 Time histories of river discharge at S-79 and the locations of 1, 10, and 20 ppt
surface salinity along the Caloosahatchee River during 2000 simulation period.161

5.41 Time histories of river discharge at S-79 and the 1 ppt salinity location along
Caloosahatchee River ......................................... 162

5.42 The 1-day averaged 20 ppt surface salinity location and 30-day averaged 10 ppt
surface salinity location during 2000 simulation period. .................. 163

5.43 The locations of 10 ppt surface salinity due to river discharge rate at S-79 during
2000 baseline simulation. ..................................... 164

5.44 The relationship between location of specific salinity value vs. river discharge at
S-79. ......................................................... 165

5.45 The relationship between salinity at Fort Myers station vs. river discharge at S-79.
........................................... .............. 166

6.1 Tidal forcing and river discharges for 1996 simulations of Charlotte Harbor.. 169

6.2 Wind velocity for 1996 simulations of Charlotte Harbor. ................... 170

6.3 Air temperature for 1996 simulations of Charlotte Harbor ................... 171

6.4 Locations of 1996 water quality measurement stations operated by EPA ....... 176

6.5 Locations of 2000 water quality measurement stations operated by SFWMD and
SW FW M D ..................................................... 177

6.6 Light intensity at water surface for 1996 and 2000 simulations ............... 178

6.7 Segments for Charlotte Harbor estuarine system .......................... 179

6.8 The scatter plots for water quality constituents during calibration period ....... 189

6.9 Temporal water quality variations at CH002 station in 1996. ................ 194

6.10 Temporal water quality variations at CH004 station in 1996. ............... 195

6.11 Temporal water quality variations at CH005 station in 1996. ............... 196

6.12 Temporal water quality variations at CH006 station in 1996. ............... 197

xii










6.13 Temporal water quality variations at CH007 station in 1996. ............... 198

6.14 Temporal water quality variations at CH09B station in 1996. ............... 199

6.15 Temporal water quality variations at CH009 station in 1996. ............... 200

6.16 Temporal water quality variations at CH010 station in 1996. ............... 201

6.17 Temporal water quality variations at HB002 station in 1996. ............... 202

6.18 Temporal water quality variations at HB006 station in 1996. ............... 203

6.19 Temporal water quality variations at HB007 station in 1996. ............... 204

6.20 Temporal water quality variations at CH013 station in 1996. ............... 205

6.21 Simulated dissolved oxygen concentration in Charlotte Harbor estuarine system for
August 21, 1996. ................................................ 206

6.22 Simulated chlorophyll a concentration in Charlotte Harbor estuarine system for
August 21, 1996. ................................................ 207

6.23 Simulated dissolved ammonium nitrogen concentration in Charlotte Harbor
estuarine system for August 21, 1996................................. 208

6.24 Simulated soluble organic nitrogen concentration in Charlotte Harbor estuarine
system for August 21, 1996. ................................... 209

6.25 Simulated soluble reactive phosphorous concentration in Charlotte Harbor estuarine
system for August 21, 1996. ................................... 210

6.26 Simulated soluble organic phosphorous concentration in Charlotte Harbor estuarine
system for August 21, 1996. ................................... 211

6.27 The scatter plots for water quality constituents in 2000. .................. 213

6.28 Temporal water quality variations at CH002 station in 2000. ............... 218

6.29 Temporal water quality variations at CH004 station in 2000. ............... 219

6.30 Temporal water quality variations at CH005 station in 2000. ............... 220

6.31 Temporal water quality variations at CH006 station in 2000. ............... 221









6.32 Temporal water quality variations at CH007 station in 2000. ............... 222

6.33 Temporal water quality variations at CH09B station in 2000. ............... 223

6.34 Temporal water quality variations at CH009 station in 2000. ............... 224

6.35 Temporal water quality variations at CHO10 station in 2000. ............... 225

6.36 Temporal water quality variations at CES02 station in 2000. ............... 226

6.37 Temporal water quality variations at CES03 station in 2000. ............... 227

6.38 Temporal water quality variations at CES08 station in 2000. ............... 228

6.39 Temporal water quality variations at CH013 station in 2000. ............... 229

6.40 The comparison between simulated dissolved oxygen concentration and the
possible causes hypoxia: river discharge, salinity, temperature, and re-aeration and
SOD fluxes at CH006 water quality measured station. ................... 230

6.41 Simulated longitudinal-vertical salinity and dissolved oxygen concentration along
the Peace River at 1 pm on June 18 (Julian Day 170), 2000................ 231

6.42 Simulated longitudinal-vertical salinity and dissolved oxygen concentration along
the Peace River at 1 pm on October 6 (Julian Day 280), 2000.............. 232

6.43 Simulated near-surface chlorophyll a concentration in Charlotte Harbor estuarine
system for February 9, May 9, August 7, November 5, 2000. .............. 233

6.44 Simulated near-bottom dissolved oxygen concentration in Charlotte Harbor
estuarine system for February 9, May 9, August 7, November 5, 2000. ...... 234

6.45 Comparisons of simulated surface chlorophyll a concentration for three cases at
three selected stations: ST05 (Pine Island Sound), ST08 (San Carlos Bay), and
ST10 (Caloosahatchee River mouth) ................................. 237

6.46 Comparisons of simulated surface chlorophyll a concentration fields in San Carlos
Bay after 90 days of simulation for three cases ................... ..... 238

6.47 Comparisons of simulated surface dissolved ammonium nitrogen (NH4)
concentration fields in San Carlos Bay after 90 days of simulation for three cases.
. .. . . . . . 2 3 9

6.48 The water quality species at CH004 water quality measured station before and after
100 % nitrogen load reduction from Peace River ........................ 247

xiv










6.49 The water quality species at CH006 water quality measured station before and after
100 % nitrogen load reduction from Peace River ........................ 248

6.50 The water quality species at CH004 water quality measured station before and after
100 % phosphorous load reduction from Peace River. ................... 249

6.51 The water quality species at CH006 water quality measured station before and after
100 % phosphorous load reduction from Peace River. ................... 250

6.52 The water quality species at CES02 water quality measured station before and after
100 % nitrogen load reduction from Caloosahatchee River. ............... 251

6.53 The water quality species at CES08 water quality measured station before and after
100 % nitrogen load reduction from Caloosahatchee River. ............... 252

6.54 The water quality species at CES02 water quality measured station before and after
100 % phosphorous load reduction from Caloosahatchee River. ........... 253

6.55 The water quality species at CES08 water quality measured station before and after
100 % phosphorous load reduction from Caloosahatchee River. ........... 254

6.56 Dissolved oxygen and Chlorophyll a concentrations at CH006 water quality
measured station before and after 100 % organic matter load reduction from Peace
River using DiToro's sediment flux model ................... .. 255

6.57 Dissolved oxygen concentrations at CH004 and CH006 water quality measured
stations before and after 50 % SOD reduction and 75% SOD reduction ...... 256

6.58 The comparison of hypoxia area at Upper Charlotte Harbor according to varying
SOD constant rate at 20C. ........................................ 257

A. 1 Flow chart for main CH3D program. ................................ 264

A.2 Flow chart for driving subroutine for the time stepping of the solution. ........ 266

A.3 Flow chart for initializing sediment transport model. ...................... 271

A.4 Flow chart for main sediment transport model ............................ 272

A.5 Flow chart for initializing water quality model ........................... 277

A.6 Flow chart for main water quality model. ....................... ..... 275









A.7 Flow chart for computing temperature and light attenuation functions for water
quality model .......................... .......................... 277

D. 1 The vertical one-dimensional z-grid ............................ .... 292

E.1 Nitrogen cycle. ......................... ........................... 296

E.2 Phosphorous cycle ................................................. 302

G.1 Parallel speedup gained in performing the simulation on the SGI Origin platform.312














Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

MODELING THE CIRCULATION AND WATER QUALITY
IN CHARLOTTE HARBOR ESTUARINE SYSTEM, FLORIDA

By

Kijin Park

August 2004

Chair: Y. Peter Sheng
Department: Civil and Coastal Engineering

This study aims to develop an enhanced version of a three-dimensional curvilinear-

grid modeling system, CH3D-IMS, which include a 3-D hydrodynamics model, a 3-D

sediment transport model, and a 3-D water quality model, to simulate circulation and water

quality of the Charlotte Harbor Estuarine System and to provide quantitative assessment of

various management practices.

In the past decade, the upper Charlotte Harbor system has been suffering summer

hypoxia in bottom water. Field study indicated that hypoxia in the upper Charlotte Harbor

is related to a strong stratification caused by high freshwater flows and dissolved oxygen

fluxes at the air-sea and sediment-water interfaces. To simulate the hypoxia event, models

of oxygen balance and oxygen fluxes at the air-sea and sediment-water interfaces in previous

versions of CH3D-IMS are enhanced. The three dimensional temperature model and

physics-based light model in CH3D-IMS are also enhanced to enable better understanding


xvii








of the temporal and spatial variations of temperature and light and their effect on water

quality processes.

The hydrodynamics component of the integrated model CH3D-IMS for Charlotte

Harbor has been successfully calibrated, using hydrodynamic data gathered by the National

Oceanic and Atmospheric Administration (NOAA) and the United States Geological Survey

(USGS) in 1986 and 2000. This calibrated model was applied to assess the impact of the

removal of the Sanibel Causeway and the IntraCoastal Waterway on the circulation in the

San Carlos Bay area and to provide a quantitative evaluation of the minimum flow and level

(MFL) for Caloosahatchee River. The results show that the hydrologic alterations would not

noticeably affect the circulation, salinity, and water quality in San Carlos Bay except near the

mouth of Caloosahatchee River. The minimum flow required to produce a salinity of no

more than 10 ppt at Fort Myers is about 18 m3/s.

The water quality model of the CH3D-IMS was calibrated with the systematic

calibration procedure and validated using hydrodynamic, sediment, and water quality data

provided by the USGS and the United States Environmental Protection Agency (USEPA)

in 1996, and by the USGS, the South Florida Water Management District (SFWMD) and the

Southwest Florida Water Management District (SWFWMD) in 2000. This validated model

was used to examine the temporal and spatial dynamics of factors which can affect hypoxia

in upper Charlotte Harbor, such as freshwater inflow, tidal variation, sediment oxygen

demand (SOD), water column oxygen consumption, and dissolved oxygen (DO) re-aeration.

The model results suggest that hypoxia in the upper Charlotte Harbor System is primarily

caused by a combination of vertical salinity stratification and SOD, while water quality also

affected by the oxygen re-aeration and water column oxygen consumption. This model was


xviii









applied to assess the effects of hydrologic alterations and to provide a preliminary evaluation

of pollutant load reduction goal (PLRG). Due to lack of detailed data and insufficient

understanding on the causes of SOD and how SOD is related to the loading, the present

model cannot simulate the complete effect of nutrient load reduction on the DO

concentration in the estuary. Further research and more complete data are needed

A systematic calibration procedure has been developed for a more efficient and more

objective calibration of the water quality model. A consistent framework for systematic

calibration is formulated, which include in the following steps: model parameterization,

selection of calibration parameters, and formulation of calibration criteria.














CHAPTER 1
INTRODUCTION

Charlotte Harbor (Figure 1.1) is a shallow estuarine system in southwest Florida. The

estuary receives freshwater from the Caloosahatchee, Peace, and Myakka Rivers; is

connected to the Gulf of Mexico through the Boca Grande Pass, Gasparilla Pass, Captiva

Pass, Blind Pass, and San Carlos Bay; and provides water resources for several counties in

Southwest Florida. Charlotte Harbor is dominated by rivers that flow into the coastal area.

While most estuaries in Southwest Florida are influenced by the Gulf of Mexico, the

characteristics of the Charlotte Harbor system are strongly influenced by large rivers such as

the Caloosahatchee and Peace Rivers. Large fluctuations of river flow between wet and dry

seasons strongly affect the salinity and water characteristics in Charlotte Harbor (Estevez,

1998).

Population growth and development in the surrounding areas during the past few

decades have led to concerns over human impacts on the quality of the estuarine system.

Industrial and agriculture development also increase environment pollution. Growth and

development will cause an increased demand for fresh water and a corresponding increase

in urban, agricultural, and industrial waste. The inflow of freshwater is essential to the

integrity and health of the estuarine system. Increased freshwater withdrawal or diversion,

or increased wastewater discharges to the rivers and streams that flow into the estuary will

create environmental stress in the estuary (McPherson et al., 1996).












I enice


yakka River


\ Deer Praine Creek


Peace River


Big Slough Canal I r
i "V ... ,. SarasotaCo. ; DeSotoCo.
.. -ortliorr .. .. ..-. 1k--- .-----------
Chariotte Co.
'?.''" .\ oIJbean .
S .. .-, ./ Shell Creek
S ,. ------ -.,,-
ti44"Piirutni Gorda
IV.' tt" "" '' 4, ----
'' Gaspanila ,p
Sound .Charott,
t s ailla Harbor
.. -. Pass A 6' ,, *" :


L: !'' 'Carl -

o ** .-*. .* ...:i. i:;', ..j. o '

:,,,c .:- :;'* ..-.. .:. i ..^ ,. .^..3 ,.
A i.:'


* ,U !, 1 .i .,
EXPLANATION

-- SubareaDivide
* ,-.*' : -.,.,<^ ; :* .*^


Figure 1.1 Charlotte Harbor estuarine system and its subarea boundaries


--s
S0




Study Area >
1tv' 7


Charlotte Co.
Lee Co.


N




Kilometers


10 20


,H H H
10 0






3

Water quality in the southern Charlotte Harbor estuarine system (the Caloosahatchee

Estuary, San Carlos Bay, and Pine Island Sound) appear to be influenced by freshwater

discharge from the Franklin Lock and Dam, also known as United States Army Corps of

Engineers Structure Number 79 (S-79). Waste load allocation studies conducted by the

Florida Department of Environmental Regulation (Degrove, 1981; Degrove, and Nearhoof,

1987; Baker, 1990) concluded that the Caloosahatchee Estuary had reached its nutrient

loading limits as indicated by elevated chlorophyll-a and depressed dissolved oxygen (DO)

levels. Similarly, McPherson and Miller (1990) concluded that increased nitrogen loading

would result in undesirable increases in phytoplankton and benthic algae. In order to

effectively manage the loading of pollutants such as nutrients in a shallow estuary such as

Charlotte Harbor, and control eutrophication, it is necessary to have a quantitative

understanding of the transport and transformation processes of nutrients in system.

Increased nutrient loading from tributaries, the atmosphere, and bottom sediments

have been known to cause eutrophication in estuaries. Eutrophication varies both spatially

and temporally. Eutrophication does not necessarily occur in regions with high nutrient

loading because the circulation and sediment dynamics affect nutrient fate and transport. In

the Upper Charlotte Harbor where a pycnocline develops frequently during high river flow

conditions, hypoxia usually develops within two days of the formation of the pycnocline,

when coupled with significant sediment oxygen demand (SOD) (CDM, 1998). In fact,

hypoxia can occasionally extend all the way to the Boca Grande area. This strong linkage

between hydrodynamics and water quality dynamics suggests the importance of having a

comprehensive understanding of hydrodynamics.

Nutrient dynamics in estuaries are not only determined by biological and chemical






4

processes, but is also strongly affected by weather and climate (wind, air temperature, etc.),

hydrodynamics (wave, tide, current, turbulence mixing, etc.), and sediment transport

processes resuspensionn, deposition, flocculation, etc.). In this study, integrated modeling

of the system is conducted to understand the complex water quality processes, which are

strongly linked to hydrodynamics and sediment processes. The primary objective of this

study is to use models and field data to produce a detailed characterization of

hydrodynamics, sediment, and water quality dynamics within the Charlotte Harbor estuarine

system.

The Chesapeake Bay model (Cerco and Cole, 1994) and Moreton Bay model

(McEwan and Garbric, 1998) are coupled hydrodynamics-water quality models. These

coupled models are more suitable for managing nutrient loads and predicting eutrophication-

related problems than uncoupled models. Because these models were not coupled with a

dynamic model for sediment transport, however, they could not accurately consider

sediment-process effects such as resuspension, deposition, flocculation, and settling on

nutrient dynamics in estuaries. Therefore, these loosely coupled models cannot account for

nutrient release by sediments in episodic events (Chen and Sheng, 1994)

Chen and Sheng (1994) developed a coupled hydrodynamic-sediment-water quality

model and applied it to Lake Okeechobee. Their 3-D model includes a hydrodynamic model,

a sediment transport model, and a water quality model with a nitrogen and a phosphorous

cycle on a rectangular grid. Measured nitrogen and phosphorous dynamics during episodic

events in 1993 and monthly sampling events in 1989 were accurately simulated with

absorption-desorption reactions and the exchange of nutrients between the sediment and

water column.






5

While Chen and Sheng (1994) developed a rectangular grid model, Yassuda and

Sheng (1996) developed an integrated model in curvilinear grids and applied it to Tampa

Bay. Their water quality model, based on the Chen and Sheng (1994) model, includes a

nitrogen cycle and an oxygen cycle, but not a phosphorous cycle. Zooplankton distribution

is allowed to influence the phytoplankton dynamics. Their integrated model also includes

a wave model, a light model, and a seagrass model. Although Yassuda and Sheng (1996)

simulated the observed hypoxia in Tampa Bay during 1991, their model did not include

dynamic fluxes of oxygen at the air-sea interface (re-aeration) and the sediment-water

interface (sediment oxygen demand). Their model could not simulate the daily fluctuation

and vertical stratification of DO. Although the phytoplankton growth rate in their model is

controlled by temperature and light, it is difficult to reproduce the vertical distribution of

phytoplankton because the temperature is not simulated. The light model in Yassuda and

Sheng (1996) has not been sufficiently validated with Tampa Bay data, thus limiting the

predictability of their model for vertical distribution of phytoplankton and dissolved oxygen.

Their model could not produce the low dissolved oxygen phenomena in upper Charlotte

Harbor because there are no dissolved oxygen fluxes at the air-water and sediment-water

interfaces, such as reaeration and sediment oxygen demand (SOD).

Sheng (2000) developed the framework for an integrated modeling system: CH3D-

IMS (http://ch3d.coastal.ufl.edu). The integrated modeling system (Table 1.1) is based on

the curvilinear-grid hydrodynamic model CH3D (Sheng, 1987, 1989; Sheng et al., 2002) and

also includes a wave model, a sediment transport model (Sheng et al., 2002a), a water quality

model, (Sheng et al., 2001b), a light attenuation model (Sheng et al., 2002c, Christian and

Sheng, 2003), and a seagrass model (Sheng et al., 2002d). Sheng et al. (2002) applied the






6

CH3D-IMS to the Indian River Lagoon (IRL), Florida. The water quality model (Sheng et

al., 2001b) includes a nitrogen cycle, a phosphorous cycle, a phytoplankton cycle, and a DO

cycle. The CH3D-IMS was validated with comprehensive data collected from the IRL.

Table 1.1 Component Models of the CH3D-IMS
Component Model Model Name

Hydrodynamic Model CH3D
Flow Model
Salinity Transport Model

Wave Model SMB

Sediment Transport Model CH3D-SED3D

Water Quality Model CH3D-WQ3D
Dissolved Oxygen Model
Phytoplankton and zooplankton
Model
Nitrogen Model
Phosphorous Model

Light Attenuation Model CH3D-LA

Seagrass Model CH3D-SAV

The CH3D-IMS was able to successfully simulate the observed circulation, wave,

sediment transport, water quality, light attenuation, and seagrass biomass in IRL during 1998.

However, several aspects of the model were identified for further improvement. For

example, although a temperature model was developed during the later phase of the IRL

study, it was not sufficiently validated with data. The water quality model, while capable of

simulating the annual variation in DO, could not simulate the diurnal DO variation.

In this study, the CH3D-IMS is enhanced and used to study the circulation and water

quality dynamics in the Charlotte Harbor estuarine system. The validated model will be used

to determine minimum flow and level and pollutant load reduction goal by management

agencies. Specially, the CH3D-IMS used for this study includes the simulation of






7

temperature in addition to the flow and salinity simulations. An air-sea heat flux model is

implemented along with the temperature equation. The water quality model is enhanced to

include a re-aeration model and a sediment oxygen demand (SOD) model, such that hypoxia

and daily fluctuation and vertical stratification of DO can be simulated. Moreover, the same

basic light attenuation model developed for the IRL is used for this study.

Although the factors that cause hypoxia have not been positively identified, it is

generally accepted that hypoxia is related to strong stratification caused by high freshwater

inflow. The modified CH3D-IMS is used in this study to examine the major factors of the

observed hypoxia in 2000, including freshwater inflow, tidal variation, SOD, water column

oxygen consumption, and DO reaeration.

This study includes the following objectives:

* To validate the CH3D-IMS with extensive hydrodynamic and water quality data from
Charlotte Harbor,

* To improve the simulation of dissolved oxygen processes, including the observed
daily dissolved oxygen fluctuation, and hypoxia in bottom water during high flow
event with vertical stratification,

To validate the ability of temperature model and heat flux model to simulate seasonal
and spatial temperature distribution in Charlotte Harbor,

To apply the physics-based light attenuation model to improve Charlotte Harbor
simulations, and

To demonstrate the feasibility of the validated modeling system for simulating the
effects of various anthropogenic impacts on the Charlotte Harbor estuarine system.

More specifically, the objectives of my research are to:

Reproduce observed circulation and transport dynamics in 1986 and 2000,

Reproduce observed water quality dynamics in 1996 and 2000,

Investigate the effect of the Sanibel causeway islands and IntraCoastal Waterway on
the flow, salinity, and water quality distribution in San Carlos Bay and Pine Island









Sound,

* Reproduce the main factors for bottom water hypoxia and habitat loss in Charlotte
Harbor estuarine system,

* Provide a preliminary determination of pollutant load reduction goal (PLRG) and
minimum flow and levels (MFL),

* Perform sensitivity tests to quantify the influence of each model parameter and to
determine the most sensitive parameters and most plausible parameter set, and

* Establish a systematic calibration procedure.

To achieve these goals, the characteristics of Charlotte Harbor estuarine system are

presented in Chapter 2. Chapter 3 discusses the general hydrodynamic circulation and

transport of temperature, salinity, and sediment in the estuarine system. This chapter will

include initial and boundary conditions for the hydrodynamic and transport models. As an

air/sea interface boundary condition, a heat flux model is also presented in this chapter.

Chapter 4 describes the water quality processes and models of phytoplankton, zooplankton,

nitrogen, phosphorous, dissolved oxygen, and light attenuation processes in estuarine system.

The effects of temperature and light intensity, numerical solution technique, and model

parameters and calibration procedure are also discussed in this chapter. Chapter 5 presents

hydrodynamic field data and model simulations of hydrodynamics in Charlotte Harbor

estuarine system. These model simulations include model calibration and verification of

1986 and 2000 data, simulation of the effect of causeway islands and navigation channel, and

determination of minimum flow and levels (MFL). Chapter 6 presents water quality and

sediment field data and model simulations of circulation, sediment transport and water

quality processes in the study area. These model simulations include calibration and

verification of 1996 and 2000 data and determination of pollutant load reduction goal

(PLRG). Discussions and conclusions are presented in Chapter 7.














CHAPTER 2
CHARLOTTE HARBOR CHARACTERIZATION

Charlotte Harbor, a coastal-plain estuarine system is one of the largest estuarine

systems on the southwest florida coast and is an important part of the Gulf of Mexico

watershed. It is located on the southwest comer of the Florida peninsula, between 260 20'

and 270 10'N and 81 40' and 82 30' W. As shown in Figure 1.1, the Charlotte Harbor

estuarine system is sub-divided into Upper Charlotte Harbor, Lower Charlotte Harbor, Pine

Island Sound, Matlacha Pass, San Carlos Bay, Gasparilla Bay, Peace River, Myakka River

and Caloosahatchee River. The drainage area (Figure 2.1) consists of the Peace, Myakka,

and Caloosahatchee River Watersheds and the coastal area and islands that drain directly into

the harbor. The estuary has a surface area of 648 km2, a drainage area of a more than 10000

km2, and a average depth of 2.1 m. The upper harbor has an average depth of 2.6 m, and the

lower harbor has an average depth of 1.6 m (Stoker, 1986). The estuary is separated from

the Gulf of Mexico by barrier islands and is connected to the Gulf through two major inlets:

Boca Grande and San Carlos; and through several smaller passes: Gasparilla, Captiva,

Redfish (McPherson et al., 1996).

In 1995, the Charlotte Harbor National Estuary Program (CHNEP) was established

by U.S. Environmental Protection Agency (USEPA) and the State of Florida. In its planning

documents, CHNEP identified three draft major problems: (1) hydrologic alterations, (2)

nutrient enrichment, and (3) habitat loss. In order to address these problems, it is useful to






10

have a comprehensive understanding of the climate, hydrodynamics, sediment properties, and

biological and chemical characteristics of the system.


Figure 2.1 Drainage basins of Charlotte Harbor estuarine system (Stoker, 1992)








2.1 Climate

The climate of the Charlotte Harbor estuarine system is subtropical and humid. The

mean annually average temperature is 72 F, with low mean of 600F in December and

January and a high mean of 80 F during the summer (McPherson et al., 1996).

Annual rainfall averages 132 cm, of which more than half occurs from June through

September, during local thundershowers and squalls. Rain during the fall, winter, and spring

seasons is usually the result of large frontal systems and tends to be more broadly distributed

than rain associated with local thunderstorms and squalls. The period from October through

February is characteristically dry, with November usually being the driest month. The

months of April and May also are characteristically dry. Low rainfall in April and May

coincides with high evaporation and generally results in the lowest streamflow, lake stage,

and ground-water levels of the year (Hammett, 1990).

The annual average wind speed is 3.9 m/s from the east. Four typical seasonal wind-

field patterns are shown in Figure 2.2. In the Winter months, the easterly trade winds

dominate the region south of latitude 270 N, while the westerlies dominate the area north of

latitude 290N. Spring and Summer generally exhibit more southerly winds, and Fall is

characterized by easterly or northeasterly winds. Wind speed can exceed 10 m/s during the

passage of Winter storms or during Summer squalls, hurricanes and tornadoes (Wolfe and

Drew, 1990).

2.2 Hydrodynamics

The hydrodynamic behavior of the Charlotte Harbor estuarine system is affected by

physical characteristics (bathymetry and geometry) as well as weather, climate, and

oceanographic and hydrologic characteristics, such as wind, atmosphere heating and cooling,






12

evaporation and precipitation, tidal-stage oscillations, discharge through tidal inlets, tidal

velocity, and freshwater inflow.


Figure 2.2 Seasonal wind pattern in Florida (Echternacht, 1975)








2.2.1 Tidal Stage, Discharge at Inlet, and Tidal Circulation

Tides along the Gulf coast of West-Central Florida in the vicinity of Charlotte Harbor

have a range of 30 to 140 cm and are of the mixed type with both diurnal and semidiurnal

characteristics (Goodwin and Michaelis, 1976). Spring tides, which have the largest range,

sometimes have only a diurnal fluctuation, whereas neap tides, which have the smallest

range, approach semidiurnal conditions of two nearly equal high and low water levels per

day.

The boundary between the Charlotte Harbor estuarine system and the Gulf of Mexico

extends about 40 mi from Gasparilla Pass on the north to San Carlos Pass on the south.

Tidal characteristics in the Gulf of Mexico are nearly uniform from Gasparilla Island to the

western face of Sanibel Island, but are of a larger range off the southern shore of Sanibel

Island (Goodwin, 1996).

Circulation in the system is primarily driven by Gulf tides, entering the system

through San Carlos Pass, Boca Grande Pass, Captita Pass, Redfish Pass, Gasparilla Pass and

Blind Pass. San Carlos Pass has a maximum depth of 5.3 m and a width of 3.3 km; Boca

Grande Pass has a maximum depth of a 18.3 m and a width of 1.28 km. The discharge

through Boca Grande Pass is about twice the discharge through San Carlos Bay and three to

four times the discharge through Captiva and Redfish Passes (Goodwin, 1996). The

geometric narrowing that occurs at passes focuses tidal energy, resulting in high velocities

associated with the large volume of water moving through the pass. This tidal energy is

dispersed inside the harbor and influences the harbor and tidal rivers as much as 40 to 44 km

upstream from the Peace River mouth. There is about a 2-hour lag between tide phases at

Boca Grande and tide phase in the upper harbor near the mouth of the Peace River (Stoker,






14

1992). Tide and wind keep the water column well mixed in the southern part of the estuary.

In the northern part of the estuary, vertical stratification can develop during moderate to high

fresh water inflows and can persist for weeks after a high freshwater inflow event (Sheng,

1998).

2.2.2 Freshwater Inflow

The majority of the freshwater that enters Charlotte Harbor come from the Myakka

River, Peace River and Caloosahatchee River. Average flows are 17.8 m3/s, 56.9m3/s and

56.7 m3/s, respectively. Flows in the Myakka and Peace Rivers are largely unregulated,

while flow in the Caloosahatchee is controlled by operation of the Franklin Lock about 43

km upstream from the mouth. Discharge in the Peace and Myakka Rivers tend to peak in

August and September when rainfall totals are generally the greatest. Discharges are usually

lowest in April and May (Stoker, 1992). The Caloosahatchee River discharge does not

always correspond to rainfall patterns in the basin, since it is controlled by S-79.

Analyses of long-term streamflow trends in the Charlotte Harbor area have indicated

statistically significant decreases in streamflow at several gages in the Peace River basin

from 1931 to 1984 (Hammett, 1990). The long-term decrease in streamflow of the Peace

River is probably related to the increased use of ground water and subsequent decline of the

potentiometric surface of the upper Florida aquafer (Hammett, 1990). A sustained

significant reduction in streamflow could result in an increase of salinity in Upper Charlotte

Harbor, possibly approaching the Gulf of Mexico salinity (Sheng, 1998). The impact of

freshwater reduction is of major ecological and economic significance.

2.2.3 Salinity

As in any other typical estuarine system, Charlotte Harbor generally exhibits






15

significant horizontal gradients in salinity. The higher salinity values in the adjacent Gulf

of Mexico fluctuate around 36 ppt, whereas the lower salinity levels nearly zero in wet

season occur near the mouth of creeks and rivers. Seasonal changes in salinity occur

primarily in response to changes in freshwater inflow from the Peace, Myakka, and

Caloosahatchee River basins. Other sources of freshwater, including direct rainfall, runoff

form coastal areas, ground-water seepage, and domestic influence, have smaller and usually

more local effects on salinity in the estuary (McPherson et al., 1996).

Stoker (1992) described salinity characteristics in the system based on data collected

from June 1982 to May 1987. Salinity generally is the lowest during the wet season between

July and September, and is the highest during the dry season from January through March.

Salinity also varies daily in response to tidal fluctuation. Peak salinity is near the flood-tide

stage, and lowest salinity is near the ebb-tide stage.

Due to the shallow depth and because of significant vertical mixing, salinity is

generally not stratifed in the southern part of the estuary. Vertical salinity stratification in

the Upper Charlotte Harbor is a common seasonal occurrence (Environmental Quality

Laboratory, Inc., 1979). In high river inflow events, a stable vertical salinity gradient is

created which suppresses vertical mixing unless there are sufficient mixing by wind or tide.

2.3 Water Quality

Water quality refers to the condition of water (e.g., dissolved oxygen concentration,

chlorophyll concentration, light attenuation, etc.) relative to legal standard, social

expectations or ecological health. Overall, water quality in the Charlotte Harbor estuarine

system is fair or good, but some areas have poor water quality or declining trends (Estevez,

1998). The water quality of an estuary is strongly influenced by hydrodynamic processes






16

(e.g., circulation and flushing), and chemical and biological processes (e.g., ammonification,

mineralization, decomposition, algae uptake, excretion, and mortality, etc.) in the estuary and

basin.

Color levels and concentrations of nitrogen, phosphorous and chlorophyll-a in

Charlotte Harbor exhibit pronounced salinity-related gradients which extend from the head

of the estuary to its mouth (Morrison, 1997). The water quality in the Caloosahatchee system

is more degraded than the water quality in the Myakka or Peace systems. Oxygen depletion

is common upstream of Franklin Lock. Nutrient and chlorophyll levels are high, and algal

blooms occur regularly in the tidal river (Estevez, 1998).

2.3.1 Nutrients

Nutrient availability is a key factor in the regulation of primary productivity in

estuarine and coastal water (Ketchum, 1967). Increased nutrient loads related to the urban

development of coastal basins have been implicated in estuarine nutrient enrichment,

increased phytoplankton productivity, and increased phytoplankton biomass (Jaworski, 1981),

and in declines of seagrass communities (Orth and Moore, 1983).

The distribution of nutrients in the system is mainly the result of nutrient input from

rivers, freshwater and tidal flushing, and recycling processes in the estuary (McPherson and

Miller, 1990). The major factor that influences estuarine nutrient distribution is freshwater

inflow from rivers, which contributes substantial nutrient loads and flushes nutrients

seaward. NOAA estimates that Charlotte Harbor receives about 2,500 tons of nitrogen as

total Kjeldahl nitrogen, or TKN, and 1,000 tons of phosphorous per year. Relative to its

dimensions and flushing characteristics, phosphorous loads are high, signifying a nitrogen

limited system (Estevez, 1986).






17

Molar ratios of dissolved inorganic nitrogen to dissolved inorganic phosphorous were

below the Redfield ratio of 16, averaging 5.7 at Caloosahatchee River during two sampling

periods, which are from 1985 to 1989 and from 1994 to 1995 (Doering and Chamberlain,

1997). By contrast, total nitrogen to total phosphorous ratios were above 16, averaging 35

in Caloosahatchee River. With these results, dissolved inorganic N:P ratios (<16) suggest

that nitrogen could limit phytoplankton productivity in agreement with McPherson and

Miller (1990), while total N:P ration (>16) indicated that phosphorous could limit

productivity (Doering, 1996). In fact, nutrient addition studies conducted by FDER found

N and P to be co-limiting (Degrove, 1981)

The distributions of phosphorous in the system were nearly conservative and a

function of river phosphorous concentration, flow, and physical mixing (McPherson and

Millor, 1990). A large amount of phosphorous from the watershed is carried by freshwater

discharge into the tidal reaches of the Peace River, with subsequent dilution by the lower-

nutrient seawater entering the estuary from the Gulf of Mexico. Concentrations of total

phosphorous averaged about 0.08 mg/L in Pine Island Sound, 0.15 mg/L in the tidal

Caloosahatchee River, and 0.62 mg/L in the tidal Peace River according to USGS nutrient

data from 1982 to 1989 (McPherson and Miller, 1990).

Most of the nitrogen in the rivers and estuary is organic nitrogen (McPherson and

Miller, 1990). Organic nitrogen concentrations decreased over the salinity gradient,

indicating river input as a source. The relatively low concentrations of inorganic nitrogen

could limit plant growth in the estuary (McPherson and Miller, 1990). Ratios of nitrogen and

phosphorous constituents at the Peace River, Myakka River, and Upper Charlotte Harbor

determined by USGS nutrient data during 1975 through 1990 are shown in Table 2.1 (Pribble









et al., 1998)

Table 2.1 Ratios of nitrogen and phosphorous constituents (Pribble et al., 1998)
RATIO Peace River Myakka River Upper Harbor

NH3: TN 0.0327 0.0349 0.0338

N03 : TN 0.4378 0.0505 0.2442

ON: TN 0.5614 0.9096 0.7355

P04: TP 0.9384 0.9371 0.9361

OP : TP 0.0616 0.0683 0.0650

Concentrations of ammonia were highly variable along the salinity gradient and were

in about the same range as concentrations in the rivers (McPherson and Miller, 1990).

Ammonia concentrations increased in the deeper water of Charlotte Harbor during summer

(Fraser, 1986). Ammonia enrichment probably was related to density stratification and to

low concentrations of dissolved oxygen in bottom waters (McPherson et al, 1996).

Concentrations of nitrate and nitrite nitrogen were nonconservative and decreased sharply

along the salinity gradient (McPherson and Miller, 1990). The sharp decline of the nitrate

and nitrite nitrogen in the low salinity regions indicates a substantial removal of nitrogen

from the water column due to biological uptake.

2.3.2 Dissolved Oxygen

Dissolved oxygen is critical for survival of plants and animals in fresh and salt water,

and a major constituent of interest in water quality study. Dissolved oxygen concentrations

in the near surface water of the system ranged from about to 6 to 8 mg/L during daylight

sampling in 1982 84 (Stoker, 1986). Dissolved oxygen concentrations of nearbottom water

of the estuary generally are lower than near surface concentrations.

Bottom water hypoxia (dissolved oxygen <2 mg/L) in Charlotte Harbor has been

reported periodically by the Environmental Quality Laboratory (EQL) since 1975 (Heyl,






19

1996). Dissolved oxygen concentrations were measured monthly from 1976 to 1984 in the

river mouth of Peace River near Punta Gorda (CH-006) (Fraser, 1986). The average monthly

near-surface concentrations declined from 8.5 to 6.7 mg/L from January to July and then

began to rise (Figure 2.3). Near-bottom average monthly concentrations at this area were

highest in February, declined slowly through May, and then declined more rapidly until July

(Fraser, 1986). The hypoxia conditions during summer are attributed to strong stratification,

which cause restricted reaeration, and to SOD. After breakup of the stratification, the

concentration increased from October to December.



9 II I I-I-I\I-I-I--

8 NEAR SURFACE



0 6
7.


X 5




0 NEAR BOTTOM
C0 2


J F M A M J J A S O N D

Figure 2.3 Average monthly concentration of dissolved oxygen in upper Charlotte
Harbor, site CH-006, 1976-84 (Fraser, 1986)

2.3.3 Phytoplankton

Phytoplankton is an important component of water quality processes and a major

primary producer in coastal and estuarine waters. The temporal and spatial variability of






20

phytoplankton productivity and biomass in Charlotte Harbor have been investigated by

Environmental Quality Laboratory, Inc (EQL) (1987). Phytoplankton productivity and

biomass (as chlorophyll_a) in the system are relatively low most of time. Productivity

ranged from 5 to 343 (mgC/m3)/h and averaged 59 (mgC/m3)/h from 1985 to 1986

(McPherson et al., 1990). Chlorophyll_a concentrations ranged from 1 to 46 mg/m3 and

averaged 8.5 mg/m3. Both productivity and biomass were greater during summer near the

mouth of tidal rivers which has middle range salinity of 6 to 12 ppt (McPherson et al., 1990).

Phytoplankton productivity and biomass in the system are affected by freshwater

inflow that lowers salinity, increases nutrient availability, and reduces light penetration in the

water column. The nutrient rich colored water is diluted by seawater at middle range salinity

of 10 to 20 ppt, so that availability of light increases and sufficient nutrient concentrations

remain available from runoff, to stimulate productivity and growth of phytoplankton in these

areas (McPherson et al., 1990). Of the major nutrients for phytoplankton productivity,

inorganic nitrogen is in lowest supply and most critical in limiting phytoplankton

productivity and growth in the system (Fraser and Wilcox, 1981).

The composition of the phytoplankton in the system varied with location and season

(McPherson et al., 1990). Diatoms were dominant in 55 % of 289 phytoplankton samples

collected in the system between 1983 and 1984, cryptophytes in 35 %, cyanophytes (blue-

green algae) in about 6 %, dinoflagellates in about 4 %, and other classes in 1 % (McPherson

et al., 1996).

2.4 Sediment

Sediment quality data from Charlotte Harbor are scarce. Organic carbon, nitrogen

and phosphorous data exist in only two studies. Hwang (1966) obtained data in 1965 from






21

119 stations throughout the study area. The FDEP Coastal Sediment Contaminant Survey

obtained carbon, nitrogen phosphorous, and metal data in the sediment column from 33

stations during 1985 through 1989, in all areas except the Gulf of Mexico (Sloane, 1994).

Mean carbon/nitrogen ratios in the six different areas of Charlotte Harbor ranged from 8.9

at Lower Charlotte Harbor to 16.4 at Gasparilla Sound (Schropp, 1998). Despite data

limitation, the available data indicate that Charlotte Harbor is relatively free of sediment

contamination in comparison with Tampa Bay and Biscayne Bay, where sediment

contaminations are widespread and are present in higher concentrations (Schropp, 1998).

2.5 Light Environment

The amount of photosynthetically active radiation (PAR) in natural water is

fundamentally important in determining the growth and vigor of aquatic plants (McPherson

et al., 1996). Light is reflected, absorbed, and refracted by dissolved and suspended

substances in the water column, and by water itself. The controlling factors are chlorophyll-

a, dissolved substances, and non-algal particulate matter (Christian and Sheng, 2003).

Dissolved and suspended matter are the major causes of light attenuation in the

system: phytoplankton and chlorophyll-a are generally minor causes of attenuation

(McPherson and Miller, 1990). On average, non-chlorophyll suspended matter accounted

for 72% of light attenuation, dissolved organic matter (color) accounted for 21%,

phytoplankton chlorophyll for 4%, and water for the remaining 3% (McPherson et al., 1996).

Water color can cause light attenuation at tidal rivers because the source of the water color

is terrestrial dissolved organic matters. Dissolved organic matter has little effect on light

attenuation in much of the southern part of the estuary, where suspended matter is the major

cause of light attenuation (McPherson and Miller, 1996). The source of suspended matter






22

includes the bottom of the estuary, which consisted of very fine to fine sand (Hwang, 1966)

and organic detrital-material (McPherson and Miller, 1990), and the major rivers.














CHAPTER 3
HYDRODYNAMICS AND SEDIMENT TRANSPORT MODEL

Biological processes in an estuarine system are strongly affected by hydrodynamics

and sediment transport processes. Therefore, as the first step to study nutrient cycling in an

estuary, hydrodynamics and sediment transport processes must be investigated and

understood.

The three-dimensional hydrodynamics model CH3D (Sheng et al., 2001) and the

associated sediment transport model CH3D-SED3D (Sheng et al., 2000a) are used for

numerical simulations in this study. The model framework has been improved and modified

from earlier versions in order to develop an integrated model that couples hydrodynamics,

sediment and water quality dynamics. The code of the enhanced integrated model is

optimized to achieve more efficient coupling and more systematic structure. The detail

structure of enhanced model is explained in Appendix A as flow chart. The application of

the circulation and transport model to produce a detailed characterization of hydrodynamics

within system is the first step in the development of the integrated model of the system.

The CH3D model for Tampa Bay and Indian River Lagoon did not include the

simulation of temperature distribution. In this study, the temperature field is simulated by

solving the temperature equation with an air-sea interface boundary condition. Temperature

is an important factor for baroclinic circulation and water quality processes, in which almost

all reaction parameters in the nutrient cycle are a functions of temperature. Therefore,






24

improved simulation of temperature is expected to produce more accurate baroclinc and

water quality simulations.

3.1 Governing Equation

Hydrodynamics and sediment transport processes in estuaries are complicated, three-

dimensional and time-dependent. Mathematical descriptions of these processes generally

require simplifying assumptions.

3.1.1 Hydrodynamic Model

The governing equations that describe the velocity and surface elevation fields in

shallow water are derived from the Navier-Stokes equations. In general, four simplifying

approximations are applied: 1) the flow is incompressible, which results in a simplified

continuity equation, 2) the horizontal scale is much larger than the vertical scale such that

the hydrostatic pressure distribution is valid, 3) the Boussinesq approximation can be used

to simplify the treatment of the baroclinic terms, 4) the eddy-viscosity concept, which

assumes that the turbulent Reynolds stresses are the product of mean velocity gradients and

"eddy viscosity", can be employed. With the above assumptions, the continuity equation and

x- and y- momentum equations have the following form (Sheng, 1983):

au av w
-+-+ -= 0 (3.1)
ax ay az

au a uv auw H U (2 a2u u ( f ua
-+ + +- =-g + fv +A+ +- (3.2)
at ax ay az ax +x2 y2

av avu avv avw (a2v a a 2 2
+ + =-g --fu+A + +- Av (3.3)
at ax ay az ay x y2 ) z az
where u(x,y,z,t), v(x,y,z,t), w(x,y,z,t) are the velocity components in the

horizontal x- and y-directions, and vertical z-direction; t is time; ;(x, y, t) is the free surface






25

elevation; g is the gravitational acceleration; and AH and Av are the horizontal and vertical

turbulent eddy coefficients, respectively.

In Cartesian coordinates, the conservation of salt and temperature can be written as

aS auS avS awS =a a (D a Ds a ( 3S4
-+-+-+-=- D +- D- +- DV- (3.4)
at ax ay az ax axr ay a az r a )
aT auT +vT awT = a T a a K aT)l a T a(
-+- -+- + K- +- K- +- Kv- (3.5)
at ax ay az ax ax) ax y ay az az
where S is salinity; T is temperature; DH and KH are the horizontal turbulent eddy

diffusivity coefficients for salinity and temperature, respectively; and Dv and Kv are the

vertical turbulent eddy diffusivity coefficients for salinity and temperature, respectively.

Since the length scales of horizontal motion in estuarine systems are much greater

than those of vertical motion, it is common to treat vertical turbulence and horizontal

turbulence separately. In shallow estuaries, the effect of the horizontal eddy viscosities on

circulation are much smaller than the effect of the vertical eddy viscosity, although the

horizontal eddy viscosity is typically 2-3 orders of magnitude larger than the vertical eddy

viscosity.

Vertical turbulent mixing is an important process, which can significantly affect

circulation and transport in an estuary. Since turbulence is a property of the flow instead of

the fluid, it is essential to use a robust turbulence model to parameterize the vertical turbulent

mixing. In this study, the vertical eddy coefficients ( A, D, and K ) are computed from

a simplified second-order closure model developed by Sheng and Chiu (1986), and Sheng

and Villaret (1989).






26

Various form forms of the equation of state can be used. The present model uses the

equation given by (Eckart, 1958)

P
(a + 0.698P)
P = 5890 + 38T 0.357T2 + 3S (3.6)
a = 1779.5 +11.25T 0.0475T2 (3.8 + 0.01T) S
where T is in C S is in ppt and p is in g /cm3.

The complete details of model equations in the curvilinear boundary-fitted and sigma

coordinates have been derived (Sheng, 1989), and are presented in Appendix B.

3.1.2 Sediment Transport Model

The suspended sediment model includes the advection-diffusion processes, which are

computed by the hydrodynamics model, as well as such processes as erosion, deposition,

flocculation, settling, consolidation, and entrainment (Sheng, 1986; Metha, 1986).

The governing equation that represents the transport of suspended sediments is given

by:

ac auc Fvc a(w+w,)c 3 ( ac a ac ( ac
+-+-+ = B +- B, +- Bv (3.7)
at ax ay az ax ax ay ay az az
where c is the suspended sediment concentration, w is the settling velocity of suspended

sediment particles (negative downward), BH is the horizontal turbulent eddy diffusivity, and By

is the vertical turbulent eddy diffusivity.

Four simplifying approximations are implied in the above equation: 1) the concept

of eddy diffusivity is valid for the turbulent mixing of suspended sediments, 2) the suspended

sediment dynamics are represented by the concentrations of two particle size groups (Sheng

et al., 2002a), 3) the suspended sediment concentrations are sufficiently small that all






27

particles follow turbulent eddy motions, 4) the SSC is sufficiently low that non-Newtonian

behavior can be neglected.

In this study, the determination of settling, flocculation, deposition, erosion,

flocculation, and consolidation processes are based on the previous work of Sheng and Lick

(1979), Sheng(1986), Metha (1989), Sheng et al. (1990), Chen and Sheng (1994) and Sun

and Sheng (2001).

3.2 Boundary and Initial Conditions

Obtaining solutions for Equations (3.1) Through (3.7) requires the specification of

appropriate boundary and initial conditions.

3.2.1 Boundary Conditions

The boundary conditions at the free surface (o=0) in a non-dimensional, vertically

stretched, boundary-fitted coordinate system are:

au H
For Hydrodynamics: A- = Tw

av H


For Salinity: =- -0 (3.8)

ST HPrv
For Temperature: -t = v qT
da Ey

For Sediment: sici + = 0
H au"
where Ev is a vertical Ekman number and Prv is a vertical Prandtl number; The Cartesian

wind stress, e", is calculated using


r, = paCds Vw w
<= PCV 2+2
^=/yuw+vw


(3.9)






28

where p, is the air density (0.0012 g/cm3), uw and v, are the components of wind speed

measured at some height above the sea level. Cd., the drag coefficient, is given as a function

of the wind speed measured at 10 meters above the water surface by (Garrat, 1977).

Cds = 0.001(0.75 + 0.067W,) (3.10)

Boundary conditions, if specified in a Cartesian coordinate system, such as wind

stress, must be transformed before being used in the boundary-fitted equations. For example,

the surface stress in the transformed system is given by

=-- +-2.


a a sx a (3.11)
The second surface boundary condition is the kinematic free surface boundary

condition which states

w=- -+u-- v (3.12)
at ax ay
The boundary conditions at the bottom in a non-dimensional, vertically stretched (a

= -1), boundary-fitted coordinate system are

For Hydrodynamics:
Av ,= = HrZCd [gub + 2gUbVb + g22]ub
da Eb A,
av H U2
Av Tb r HrZCd gIlu2 2gl2Ubb 1 V 22 Vb
u "b A,,y
as
For Salinity: -- = 0 (3.13)

For Temperature: aT- 0

For Sediment: w sc; + =D E,
H aoy






29

where ub and vb are the contravarient velocity components at the first grid point above the

bottom. E; is the erosion rate and Di is the deposition rate for sediment group i. Cd is drag

coefficient which is a function of the size of bottom roughness elements, Zo, and the height

at which ub is measured, so long as z, is within the constant flux layer above the bottom. The

drag coefficient is given by (Sheng, 1983)


Cd= (3.14)
Cd 10 / Z0
where c=-0.4 is the von Karman constant, zo=kk/30 and k, is the bottom roughness.

Along the shoreline where river inflow may occur, the conditions are generally

u = u(x, y, a, t)
v = v(x, y, or, t)
w=0
S =S(x,y, a, t) (3.15)
T = T(x, y, a, t)
c, = Ci (X, y, a, t)
Contravarient velocity components provide lateral boundary conditions similar to

those in Cartesian systems (x,y). Along solid boundaries, the normal velocity component

must be zero to satisfy the no-slip condition. In addition, normal derivatives of salinity,

temperature and suspended sediment concentration are assumed to be zero. When flow is

specified at a boundary, the normal velocity component is prescribed. Tidal boundary

conditions are specified using water level, (, directly When tidal boundary conditions are

given in terms of C, the normal velocity component is assumed to be of zero slope, while

tangential velocity component may be either zero, of zero slope, or computed from the

momentum equations. During an ebb tide, the concentrations of salinity, temperature, and

suspended sediment flowing out are calculated using a 1-D advection equation while during






30

a flood tide the offshore concentrations are generally prescribed as either fixed or time

varying.

3.2.2 Initial Conditions

To initiate a simulation, the initial spatial distribution of C, u, v, w, S, T, and ci must

be specified. When these values are unknown, "zero" initial fields can be used. When these

values are known at a limited number of locations, an initial field can be generated by spatial

interpolation. In the principle, the interpolated field should satisfy the conservation equation

for a particular variable. For practical simulations, a "spin-up" period is required to damp

out transients caused by the initial condition, which does not satisfied conservation. The

length of a spin up period is variable and depends on such factors as basin size, flushing

time, and current velocity.

3.3 Heat Flux at the Air-Sea Interface

The ocean receives energy through the air-sea interface by exchange of momentum

and heat. The temperature of ocean waters varies from place to place and from time to time.

Such variations are indications of heat transfer by currents, absorption of solar energy, and

loss by evaporation. The size and character of the temperature variations depends on the net

rate of heat flow into or out of water body (Pickard and Emery, 1990). The transfer of heat

across the air-sea interface determines the distribution of temperature in the ocean as a

surface boundary condition for the temperature equation, as follow:

r)T
poKv = q, (3.16)

where qT is the net heat flux across the surface water.

To estimate the net heat flux at the air-sea interface, it is necessary to consider the

following processes: net heat flux across the surface, qT; heating by incoming solar radiation






31

insolationn), q, ; warming by conductive heat exchange across the surface (sensible heat

flux), q,,h ; cooling by net outgoing long wave radiation, q. ; heat loss as water evaporated

(latent heat flux), q, ; heat flux from precipitation. Other source of heat flow, such as that

from the earth's interior, change of kinetic energy of waves into heat at the surf; heat from

chemical or nuclear reactions, are small and can be neglected. The net heat exchange at the

surface can be divided into four terms:

q = q qb qh q (3.17)

Direct measurement of these fluxes is the best way to provide the net heat flux at the

air-sea interface. However, directly measured air-sea fluxes are only available at very few

stations to allow calculation of air-sea interface over a large area. Instead, directly measured

air-sea fluxes are used for developing, calibrating, and verifying the parametric formulae for

estimating the fluxes from the primary variables including wind speed, air temperature, water

temperature, and cloud cover (Taylor et al., 2000). These required parametric formulae are

then used to compute heat fluxes over a large area.

3.3.1 Short-Wave Solar Radiation

The main source of heat flux through the air-sea interface is short-wave solar

radiation, q, (incoming solar radiation), received either directly or by reflection and scattering

from clouds and the atmosphere. The incoming solar radiation is based on the empirical

formula of Reed (1977):

q, = Ssin y(A + Bsin r)(1-0.62n -0.0019r)(1- a) (3.18)

where S is solar constant =1353 w/m2

y is solar elevation angle

n is cloud cover = 0 1








a is the albedo = 0.06

A and B are empirically determined coefficients for each category of reported total

cloud amount and the solar elevation

The solar elevation angle was computed to the nearest 0.10 from date and time for the

latitude and longitude of the study area, using the following equations (Miller and

McPherson, 1995):

360
365.242
6= 12 + 0.1236sin (p 0.0043cos ( + 0.1538sin 2(p + 0.0608cos 2(p
a = 279.9348 + (p + 1.9148sin ( 0.0795cos (p + 0.00199 sin 2(p+ 0.0016cos 2(p
Y =15 (r- -)-A (3.19)
i = arcsin (0.39785077 x sin ea)
sin fl = sin ysin Kt+ cos ycos Kccos Y

where (p is the angular fraction of the year, in degree; d is Julian date; 6 is true solar noon,

in hours; T is the solar hour angle, in degrees; T is the Greenwich Mean Time, in hours; ,

is the longitude, in degrees; a is an estimate of the true longitude of the sun, in degrees; K is

the solar declination, in degrees; y is the latitude, in degrees; and P3 is the solar elevation

angle, in degrees.

3.3.2 Long-Wave Solar Radiation

The back radiation term, qb, is the net amount of energy lost by the sea as long-wave

radiation. The amount of long wave flux is dependent on surface water temperature,

atmospheric temperature, humidity and cloud cover (Clark et al, 1974)

qb = EUoT4 (0.39 0.05e5) )(1- An2) + 4eoT 3 (T a T) (3.20)

where e is the emittance of the sea surface = 0.98

o is Stefan-Boltzman constant = 5.673 x 108 W/m2









e is the water vapor pressure

Ta and Ts are the air and sea temperatures in K.

A is a cloud cover coefficient which varies with latitude

The Antoine constants (http://www.owlnet.rice.edu/~ceng301/31 .html) give the vapor

pressure as a function of temperature so that the approximate value for water vapor pressure

is:


e = exp 16.5362 -- (3.21)
T, 38.997
where T, is the surface water temperature (K).

A theoretical calculation of the mean values of the coefficient A for different latitudes

has been made by M. E. Berliand (1952). In this calculation, the mean frequency of clouds

of different layers at each latitude is taken into account. The values obtained for the

coefficient I are given in Table 3.1 (Budyko, 1974). The reduction of values of this

coefficient in low latitudes is explained mainly by a higher mean altitude of clouds in these

regions.

Table 3.1 Mean latitudinal values of the coefficient 1
750 700 650 600 550 500 450 400

1 0.82 0.80 0.78 0.76 0.74 0.72 0.70 0.68

350 300 250 200 150 100 50 00

A 0.65 0.63 0.61 0.59 0.57 0.55 0.52 0.50

3.3.3 Sensible and Latent Heat Fluxes

The most significant part in terms of heat transfer from the sea to the atmosphere is

latent heat flux. The rate of heat loss is equal to the rate of vaporization times the latent heat

of evaporation. Sensible heat flux is due to temperature gradient in the air above the sea.

The rate of loss or gain of heat is proportional to the temperature gradient, heat conductivity






34

and the specific heat of air at constant pressure. Sensible and latent heat fluxes are estimated

by using the bulk aerodynamic equations (Mellor, 1996):

q= pc,CHUo (T. Ta) (3.22)

q, = pCELUl (h, h, ) (3.23)
where p is the density of air

CH and CE are transfer coefficients for sensible and latent heat respectively

Uo is the wind speed at 10 m height above surface

T, and T, are the sea surface and air temperatures

cp is the specific heat of air, 1.0048 x 103 J/kgK

L is the latent heat of evaporation, 2.4 x 106 J/kg

h, and h,, are specific humidities at the surface and at a 10 m height above the surface

The air density is determined from the ideal gas equation in the form,

p = P,/(RTa), where P, is air pressure and R is universal gas constant, 287.04 J/kg K.


The traditional estimates of CR and CE over the ocean tend to support a fairly constant value

over a wide range of wind speed. Smith (1989) recommended a constant "consensus" value

CE = (1.2 0.1 )x 10' for winds between 4 and 14 m/s. DeCosmo et al. (1996) also suggested

a near constant value with CE = (1.12 0.24)x10-3 for winds up to 18 m/s. For the Stanton

number, CH, Friehe and Schmitt (1976) obtained slightly different values for unstable and

stable conditions, 0.97 x 103 and 0.86 x103 respectively. Smith (1989) suggested CH = 1.0

x 103. CE= 1.2 xl03and CH = 1.0 x 103 were used for this study.

The formulae used for air specific heat, cp and latent heat for evaporation, L, have

been taken from Stull (1988),


c, = 1004.67 + 0.84q,


(3.24)






35

L = 106 [2.501- 0.00237 (T, 273.16)] (3.25)

An empirical relation (based on the Clausius-Clapeyron theory) provides the

saturation, water vapor partial pressure, e, = 10(0.7859+003477T)/(I+O.00412T) mb. A general


relation between specific humidity and partial pressure is q = (0.622e/P)/(l- 0.378e/P)


, where P is the atmospheric pressure (Mellor, 1996). Specific humidities at surface and at

a 10 m height above the surface, h, and ha, can be calculated using these empirical formulae

with T, and Ta respectively.














CHAPTER 4
WATER QUALITY MODEL

In this chapter, a water quality model for simulating water quality processes occurring

in the Charlotte Harbor estuarine system in both the water and sediment columns is

presented. This is an extension of an earlier model developed to simulate water quality for

Lake Okeechobee (Sheng et al., 1993; Chen and Sheng, 1994), Tampa Bay (Yassuda and

Sheng, 1996) and the Indian River Lagoon (Sheng et al., 2001 and 2002). These models

include the effect of sediment transport on nutrient dynamics through the explicit use of a

sediment transport model and the incorporation of a nutrient resuspension flux.

The water quality model incorporates on the interactions between oxygen balance,

nutrient dynamics, light attenuation, temperature, salinity, phytoplankton and zooplankton

dynamics. To develop the water quality model, the mass conservation principle can be

applied to each water quality parameter, as it relates to phytoplankton and zooplankton

dynamics, nitrogen and phosphorous cycles, and oxygen balance.

The nitrogen and phosphorous cycles in an estuarine system, are modeled through a

series of first order kinetics. Nutrient concentrations in the estuarine system are constantly

changing in time and space due to loading from rivers, exchanges with the ocean, seasonal

climatic changes, biogeochemical transformations, hydrodynamics, and sediment dynamics.

Nitrogen species include ammonia nitrogen (NH3), soluble ammonium nitrogen (NH4),

nitrate and nitrite nitrogen (N03), particulate ammonium nitrogen (PIN), soluble organic






37

nitrogen (SON), particulate organic nitrogen (PON), phytoplankton nitrogen (PhyN), and

zooplankton nitrogen (ZooN). Similar to the nitrogen species, phosphorous species include

soluble reactive phosphorous (SRP), particulate inorganic phosphorous (PIP), soluble organic

phosphorous (SOP), particulate organic phosphorous (POP), phytoplankton phosphorous

(PhyP), and zooplankton phosphorous (ZooP).

Phytoplankton kinetics are the central part of this water quality model, since the

primary water quality issue in the estuarine system is eutrophication (Boler et al., 1991).

Phytoplankton population is a complex variable to measure in the field. However, the lack

of data on each specific species prevented a more detail characterization, the entire

phytoplankton community is represented in this study by a single state variable and

quantified as carbonaceous biomass. Chlorophyll_a concentrations, for comparison with

observations, are obtained through division of computed carbonaceous biomass by the

carbon-to-chlorophyll_a ratio.

The oxygen balance couples dissolved oxygen to other state variables. Reaeration

through the air-sea interface, and phytoplankton production during photosynthesis are the

main source for oxygen. Oxidation, nitrification, respiration, mortality and SOD reduce

oxygen in the system. Oxidation of organic matter and carbonaceous material, respiration by

zooplankton and phytoplankton, and oxygen consumption during the nitrification process are

collectively grouped into the CBOD (Carbonaceous-Biogeochamical Oxygen Demand)

variable, which is a sink for dissolved oxygen (Ambrose et al., 1994).

The methods of coupling with hydrodynamic and sediment transport models, the

simulated parameters, the assumptions, the chemical/biological processes of the CH3D water

quality model (CH3D-WQ3D) was compared with those for existing water quality models,






38

specially Water Quality Analysis ans Simulation Program (WASP), the integrated

Compartment water quality model developed by the US Army Corps (CE-QUAL-ICM) in

Appendix C.

Temperature, salinity, and light are important parameters that effect the rate of

biogeochemical reactions. Most of the transformation processes in the nutrient cycle are

affected by temperature. Salinity also influences DO saturation concentration and is used in

the determination of kinetics constants that differ in saline and fresh water. Light intensity

affects the photosynthesis process and thus algae growth rate.

The water quality model is enhanced to include a reaeration model and a sediment

oxygen demand (SOD) model and newly coupled with temperature model and physics-based

light attenuation model. Applying enhanced water quality model and a more accurate light

model and temperature model could improve the vertical distribution and daily fluctuation

of both phytoplankton and dissolved oxygen. Furthermore, the bottom water hypoxia at

upper Charlotte Harbor which is caused by SOD and vertical stratification can be reproduced

and analyzed in this study.

4.1 Mathematical Formulae

The water quality equations are derived from a Eulerian approach, using a control

volume formation. In this method, the time rate of the concentration of any substance within

this control volume is the net result of (i) concentration fluxes through the sides of the

control volume, and (ii) production and sink terms inside the control volume. The

conservation equation for each of the water quality parameters is given by:

V.(i )=V.[DV(u)] + Q (4.1)
at i (4.1)
Wi (ii) (iii) (iv)






39

where (I) is the evolution term (rate of change of concentration in the control volume), (ii)

is the advection term (fluxes into/out of the control volume due to the advection of the flow

field), (iii) is the diffusion term (fluxes into/out of the control volume due to turbulent

diffusion of the flow field), and (iv) is the sink and source term, representing the kinetics and

transformations due to sorption/desorption, oxidation, excretion, decay, growth, and bio-

degradation. In the finite difference solution of the water quality model, the advection and

horizontal diffusion terms are treated explicitly, whereas the vertical diffusion and

biogeochemical transformations are treated implicitly. The detail descriptions of numerical

solution technique are described in Appendix D. The water quality equations in the

curvilinear non-orthogonal boundary fitted system (E, ir, o) are given by:

a-=q E, a Dv a w-Ra(Hw)
at H -Scv v o) (a

R H Fonu(pj)+ a( ogv(pi)1

a ( 'i (4.2)
+ Ev g Hg" (+ ,v Hg12 +Fi}(.


+ E_ 8 Hg_+ H22 ?
SCH[o a ]
where (po represents any water quality parameter, 90 is the Jacobian of horizontal


transformation, ( g" g2, g22 ) are the metric coefficients of coordinate transformation,

and Qj represent biogeochemical processes. Equation (4.2) is in dimensionless form and the

dimensionless constants are defined in Appendix B.

In the following sections, the biogeochemical processes controlling the sink/source

term of Equation (4.2) will be discussed in detail for nutrient dynamics, zooplankton and









phytoplankton dynamics, and oxygen balance in system.

4.2 Phytoplankton Dynamics

The overall water quality in the system is markedly influenced by the dynamics of

zooplankton and phytoplankton communities (Boler et al., 1991). Phytoplankton dynamics

and nutrient dynamics are closely linked, since nutrient uptake during phytoplankton growth

is the main process to remove dissolved nutrients from the water, and phytoplankton and

zooplankton respiration and mortality are major components of nutrient recycling. The

diurnal variation of dissolved oxygen is related with photosynthetic oxygen production

during the day and oxygen consumption due to phytoplankton respiration during the night.

4.2.1 Modeling Approach

Phytoplankton kinetics are represented by growth, respiration, non-predator-mortality,

grazing by zooplankton, and a settling term. The phytoplankton sources and sinks in the

conservation equation can be written as:

SaPhy -K K. + -K+ WS g. PhyC u ZooC (4.3)
at az

Where PhyC is phytoplankton biomass, expressed as carbon (gCm3); /J, is phytoplankton

growth rate (1/d); Kas is respiration rate (1/d); K. is non-predator mortality (1/d); WSgae is

the phytoplankton settling velocity(m/d); ZooC is zooplankton biomass (gCm-3); and /z is

zooplankton growth rate (1/d).

Phytoplankton growth is determined by the intensity of light, by the availability of

nutrients, and by the ambient temperature. Light limitation is formulated according to the

photo inhibition relationship (Steel, 1965). The quantity of the growth limitation factor for

nutrients is related with a half-saturation constant. The half-saturation constant refers to the

concentration of the nutrient at which the growth rate is one half its maximum value. This






41

results in a hyperbolic growth curve. An exponential increasing function is applied for

temperature limitation..

The minimum formation approach has been used to combine the limiting factor of

light and each limiting nutrient. The minimum formation is based on "Liebig's law of the

minimum" which states that the factor in shortest supply will control the growth of algae

(Bowie et al., 1980)

,u = (. )max f (T)" f (1) f(N,P)

r/T-20 min I exp -NH4 +NO SRP (4.4)
=({a) *<0 rmin -exp 1-- 4-- -, j
max ( I s H +NHH +NO, H, +SRP


where (f)nmax is the phytoplankton maximum growth rate (1/day); 0 is temperature

adjustment coefficient; T is temperature (oC); I is the light intensity, calculated by the light

attenuation model; Is is the optimum light intensity for algae growth; H,, is half saturation

concentration for nitrogen uptake (gN -3); Hp is half saturation concentration for

phosphorous uptake (gP -3); NH4 is ammonium concentration (gN -3); NO3 is nitrate

concentration (gN -3); and SRP is soluble reactive phosphorous concentration (gP -3).

Respiration and mortality are considered to be an exponentially increasing functions

of temperature:

Kas =(K,)rT "
K = O( (4 .5 )
K. = (K.), O)r -T
where (K.,)Tr and (K,,)T are respiration and mortality rate at Tr (1/day); and Tr is reference

temperature of respiration and mortality.

For phytoplankton, literature values of algae settling velocity, which account for the

limited vertical motion of these organisms will be used.

Zooplankton are included in water quality models primarily because of their effects






42

on algae and nutrients. Phytoplankton and zooplankton dynamics are closely tied through

predator-prey interaction. Phytoplankton dynamics are of major concern in this study while

no attempt is made to investigate zooplankton dynamics due to lack of zooplankton data.

Zooplankton is only considered as the predators of phytoplankton, utilizing their available

biomass as food supply. Zooplankton kinetics, influenced by growth, respiration and

mortality, are represented in a source and sink term as (Bowie, 1985):

ZooC4.6)
at = (/z KZ Kz") ZooC (4.6)

where /z is zooplankton growth rate (1/day); K,. is respiration rate (1/day); and K, is

mortality (1/d)

Zooplankton growth is represented by a temperature-dependent maximum growth

rate, which is limited by phytoplankton availability:

T-2 PhyC TrsPh
"/ ( Z) HT 20 phy(
Smax hy+ PhyC Trs phy


where (Pz)max is the zooplankton maximum growth rate (1/day); Ois temperature adjustment

coefficient; HPhy is half saturation concentration for phytoplankton uptake (gC -3); and Trsphy

is threshold phytoplankton concentration for zooplankton uptake (p.g/1).

4.2.2 Relationship between Phytoplankton and Nutrients

Phytoplankton biomass is quantified in units of carbon. In order to express the effects

of phytoplankton on nitrogen and phosphorous, the ratio of nitrogen-to-carbon and

phosphorous-to-carbon in phytoplankton biomass must be specified. Global mean values of

these ratios are well known (Redfield et al., 1966). The amounts of nitrogen and

phosphorous incorporated in algae biomass is quantified through a stoichiometric ratio.

Thus, total nitrogen and total phosphorous in the model are expressed as:









TotN = NH4 + NH3 + N03 + SON + PON + PIN + Anc PhyC + Anc ZooC (4.8)
TotP = SRP + SOP + POP + PIP + Apc PhyC + Apc ZooC

where TotN is total nitrogen (gN -3); NH4 is dissolved ammonium nitrogen (gN -3); NH3

is ammonia nitrogen (gN -3); N03 is nitrate and nitrite nitrogen (gN -3); SON is soluble

organic nitrogen (gNm-3); PON is particulate organic nitrogen (gN -3); PIN is particulate

inorganic nitrogen (gN -3); Anc is Algae nitrogen-to-carbon ratio (gN/gC) ; TotP is total

phosphorous (gP -3); SRP is soluble reactive phosphorous (dissolved phosphate) (gPm-3);

SOP is soluble organic phosphorous (gP -3); POP is particulate organic phosphorous (gP -3);

PIP is particulate inorganic phosphorous (gP -3); and Apc is algae phosphorous-to carbon

ratio (gP/gC).

The connection between the carbon, nitrogen and phosphorous cycle is shown in

Figure 4.1. Phytoplankton uptakes dissolved ammonium nitrogen, nitrate and nitrite

nitrogen, and soluble reactive phosphorous during production and releases dissolved

ammonium nitrogen, soluble reactive phosphorous and organic nitrogen, organic

phosphorous during respiration and mortality processes. Zooplankton has similar kinetic

processes as phytoplankton. The measured phytoplankton as algae mass per volume, was

converted to phytoplankton carbon with a algae to carbon ratio. The amounts of nitrogen and

phosphorous from phytoplankton can be converted with a nitrogen to carbon ratio and a

phosphorous carbon ratio, respectively.

4.3 Nutrient Dynamics

Nutrients are essential elements for life processes of aquatic organisms. Nutrients

of concern include carbon, nitrogen, phosphorous, silica and sulfur. Among these nutrients,

the first three elements are utilized most heavily by zooplankton and phytoplankton. Since

carbon is usually available in excess, nitrogen and phosphorous are the major nutrients






44

regulating the ecological balance in an estuarine system. Nutrients are important in water

quality modeling for several reasons. For example, nutrient dynamics are critical

components of eutrophication models since nutrient availability is usually the main factor

controlling algae bloom. Algae growth is typically limited by either phosphorous or nitrogen

(Bowie et al., 1980). Details on the nutrient dynamics including all the equations used by

the water quality model to calculate the nitrogen and phosphorous, can be found in

Appendix-E.


C) up Zooplankton (ZooC)


Figure 4.1 The connection between nitrogen, phosphorous and carbon cycle








4.4 Oxygen Balance

Dissolved oxygen (DO) refers to the volume of oxygen contained in water. Five state

variables participate in the dissolved oxygen balance: phytoplankton carbon (PhyC),

ammonia (NH4), nitrate (NO3), carbonaceous biochemical oxygen demand (CBOD), and

dissolved oxygen (DO). A summary is illustrated in Figure 4.2. The methodology for the

analysis of dissolved oxygen dynamics in natural water, particularly in streams, rivers, and

estuaries is reasonably well-developed (O'Connor and Thomann, 1972). Dissolved oxygen

evolution depends on the balance between production from photosynthesis, consuming from

respiration and mortality, and exchanges with the atmosphere and sediment.

The main physical mechanisms influencing DO concentration are horizontal and

vertical dispersion and diffusion. Vertical diffusion occurs across the air-sea interface as a

function of wind, waves, currents, and DO saturation rate. Laboratory experiments show that

the bottom shear stress controls the dissolved oxygen diffusive layer thickness and the flux

at the sediment-water interface (Steinberger and Hondzo, 1999). Slow oxygen diffusion rates

and high oxygen demand by sediment results in a thin aerobic layer. Two distinct sediment

zones are created in the sediment column: an aerobic layer and an anaerobic layer. The

thickness of the aerobic soil zone is influenced by oxygen concentration in the overlying

water column, and concentration of the reduced compounds in the anaerobic soil zone. The

model dissolved oxygen cycle includes the following processes

1) Reaeration

2) Carbonaceous oxygen demand (CBOD)

3) Nitrification

4) Sediment oxygen demand (SOD)









5) Photosynthesis and respiration



Carbonaceous Oxygen Demand


Dissolved Oxygen


Sediment Oxygen Demand
SOD 1.047 r-r, DO
H K02+DO
Water column


K CBOD Ku H 3 NO 3 6K Do NH 4
H. +DO H,, + DO 14 H,+DO

Aerobic Layer


Figure 4.2 DO and CBOD cycles

Reaeration

Reaeration is the process of oxygen change between the atmosphere and sea surface.

Typically, dissolved oxygen diffuses into surface waters because dissolved oxygen levels in

most natural waters are below saturation. However, when water is super-saturated as a result

of photosynthesis, dissolved oxygen returns the atmosphere. Dissolved oxygen saturation

in seawater is determined as function of temperature and salinity (APHA, 1985)









1.575701 105 6.642308x107 1.2438x10' 8.621949 x10"
In DO, =-139.34411+ + 3 4+

Sa 31929X10_ 1.9428x10 3.8673x10
1.80655 T T2 2


where DO, is equilibrium oxygen concentration, mg/1, at standard pressure

T is temperature, K, oK = OC + 273.150

Sa is salinity, ppt

The reaeration process is modeled as the flux of dissolved oxygen across the water

surface:


V d KO AEA (DO, DO) (4.10)
dt

where V and As are volume and surface area of the water body

In case where the air-sea interface is not constricted, the volume is V = A, Az The


equation for reaeration can be expressed as

dDO K
= AE (DO -DO) (4.11)
dt Az

where DO is dissolved oxygen concentration (mg/1); KAE is reaeration coefficient (m/day).

Many empirical formulas have been suggested for estimating reaeration rate

coefficient specially in the river. Bowie et al. (1985) have reviewed thirty-one reaeration

formulas, and have tried to evaluate the performance of the each formulas. Most formulae

have been developed based on hydraulic parameters, most often depth and velocity. This

review of stream reaeration has shown that no one formula is best under all conditions, and

depending on the data set used, the range of the reaeration coefficients in the data set, and

error measurement selected, the best formula may change. Among these formulas, the most






48

common method of simulating reaeration in rivers is the O'Connor-Dobbins formula. This

method has the widest applicability being appropriate for moderate to deep streams with

moderate low velocities. With approximately 2.09x105 cm2/s diffusivity of oxygen in

natural waters, the O'Connor-Dobbins formula can be expressed as


K = 3.93 H05 (4.12)

For standing water, such as lake, impoundments, and wide estuaries, wind becomes

the predominant factor in causing reaeration. The oxygen-transfer coefficient itself can be

estimated as a function of wind speed by a number of formulas. Chapra (1997) compared

four common wind-dependent reaeration formulas: Broeckeret al. (1978), Banks and Herrera

(1977), O'Connor (1983), and Wanninkhof et al. (1991). The comparison shown in Figure

4.3 show all these methods except Broecker et al.'s have similar reaeration coefficients when

wind speed is less than 5 m/s. When wind speed is greater than 5 m/s, Bank and Herrera's

formula produce the middle range of reaeration coefficient among these three formulas.

This formula uses various wind dependencies to attempt to characterize the difference

regimes that result at air-water interface as wind velocity increase (Banks 1975; Banks and

Herrera, 1977).

K= =0.728UOw -0.317U +0.0372U2 (4.13)

Since estuary gas transfer can be affected by both water and wind velocity, effort to

determine reaeration in estuaries combines elements of current and wind-driven approaches.

Thomann and Fitzpatrick (1982) combined the two approaches for estuaries affected by both

tidal velocity and wind,


KAE =3.93 + 0.728Uo5 -0.317Uw -0.0372Uw (4.14)
HL / -(






49

where U0 is depth averaged velocity (m/s); H is a depth (m); U, is wind speed (m/s).

For the Charlotte Harbor estuarine system, The aeration coefficient is assumed to be

proportional to the water velocity, depth, and wind speed following Thomman and

Fitzpatrick (1982).




10-

9

8

7-

6 -
Broecker et al. (1978) ,'
1 5-

4 "
Wanninkhof et al (1991), -
3-
Banks and
Herrera (1977)


1 O'Connor( 1983)


0 2 4 6 8 10
Uw (m/s)

Figure 4.3 Comparison of wind-dependent reaeration formulas.

Carbonaceous Oxygen Demand

The use of carbonaceous oxygen demand (CBOD) as a measure of the oxygen-

demanding processes simplified modeling efforts by aggregating their potential efforts

(Ambrose et al., 1994). Oxidation organic matter, nitrification, non-predatory mortality and

respiration by zooplankton and phytoplankton are nitrogenous-carbonaceous-oxygen-

demand, collectively combined into the state variable CBOD.






50

The kinetic pathway of CBOD is represented in the source term of the equation as

(Ambrose et al., 1994):

For water column:
kCBOD = WSCBOD (1- fdCBOD) CBOD] K DO CBOD
at az CBOD + D Hcoo. +DO
(4.15)
-D_ K -NO3 +Aoc (K, PhyC +K ZooC)
4 14 Hno3 + DO
For sediment column:
a a DO
-CBOD = +-[wsCBOD (1- fdBOD ) CBOD K, CBOD
at az HCBOD + DO
5 32 H (4.16)
KDN no3 NO
4 14 Hno3 +DO

where fdcBoo corresponds to the fraction of the dissolved CBOD; WScBOD is the settling

velocity for the particulate fraction of CBOD ; K. is oxidation coefficient which is a

temperature function; HCBOD is a half saturation constant for denitrification; KDN is

denitrification constant; HNo3 is a half saturation rate for denitrification; and Aoc is oxygen-

carbon ratio (gO2/gC).

The consumption of oxygen, as a function of water column CBOD decay, can be

expressed as CBOD oxidation;

aDO= K DO CBOD (4.17)
at HCBOD +DO
Note that non-oxidative processes such as settling, denitrification, and mortality do

not contribute to dissolved oxygen depletion, and are not included in the expression.

Nitrification

The transformation of reduced forms of nitrogen to more oxidized forms

(nitrification) consumes oxygen. Although nitrification is also a nutrient transformation, this

section addresses oxygen consumption. First order kinetics is the most popular approach for








simulating nitrification in natural system:

a 64 DO
SDO =-- KN NH4 (4.18)
at 14 Hnit + DO

where KNN is nitrification rate which is a function of temperature; and H,,i is the half-

saturation constant for the bacteria growth.

Sediment Oxygen Demand

Sediment oxygen demand is a dissolved oxygen flux at water/sediment interface due

to the oxidation of organic matter in bottom sediments. The particulate organic matters are

from a source outside the system such as wastewater particulate or leaf litter materials, and

generated inside system as occurs with plant growth in highly productive environments.

According to accumulate these particulate organic matter in the bottom sediment, the

sediment oxygen demand will increase due to oxidation of the accumulated organic matter

at aerobic sediment layer. Figure 4.4 shows the mechanism of SOD flux with particulate

organic matters in the sediment column. Particulate organic matter (POM) is delivered to

the sediments by settling. Within the anaerobic sediment layer, the organic carbon, sulfate,

nitrate undergo reduction reactions to yield dissolved methane, sulfide, dissolved ammonium

nitrogen. These reduced species (CH4, H2S, NH4) diffuses upward to the aerobic layer

where these are oxidized. During these oxidation reactions, SOD is generated. Any residual

reduced species that are not oxidized in the aerobic layer is diffused back into the water

column where additional oxygen is consumed by oxidation. Developing an SOD model is

quite complicated task since many aerobic and anaerobic reactions in the sediments are

involved and many chemical species are included in mass balance equation. Therefore, these

redox chemistry with POM are usually treated as a composite characteristics of the particular

system. Recently, techniques have been developed for investigating these factors.






52

DiToro et al. (1990) developed a model of the SOD process in a mechanistic fashion

using the square-root relationship of SOD to sediment oxygen carbon content. Using similar

analysis as applied to carbon, they also evaluate the effect of nitrification on SOD. In this

model, carbon and nitrogen diagenesis are assumed to occur at uniform rates in a

homogeneous layer of the sediment of constant depth (active layer). The sediment oxygen

demand and sediment fluxes are calculated by the concentrations of particulate organic

carbonaceous material and of particulate organic nitrogenous material in this active layer.

The more detail description of sediment flux model developed by DiToro and Fitzpatrick

(1993) is presented in Appendix F. Their framework has been applied to the Chesapeake

Bay (Cerco and Cole, 1994; DiToro and Fitzpatrick, 1993) with sediment flux model which

include ammonia and nitrate flux and the sulfide, oxygen, phosphorous and silica flux.

A calculation of the detailed redox chemistry of the sediment interstitial water is

required for a detail understanding of the situation for each reduced species and chemical

parameters for redox reaction such as pH, Eh in the system. In the absence of data for CH4,

H2S, Iron, Methane, and C02, and detail understanding of redox chemistry of the system,

applying this model could create more uncertainty than the simple empirical formula based

on the measurement technique.

The process of oxygen demand in the sediment/water interface is usually referred to

as sediment oxygen demand (SOD) because of the typical mode of measurement: enclosing

the sediments in the chamber and measuring the change in the dissolved oxygen

concentration at several time increments. The major factors affecting SOD are: temperature,

oxygen concentration at the sediment water interface, organic and physical characteristics of

the sediment, and current velocity over the sediments (Bowie, 1985).











Water Column


Anaerobic Layer


OM, S04 Reduction ,

- Production of
CH4. NH4. H2S


;H4, NH4, H2S


POM flux


SOD flux


Figure 4.4 The relationship between POM flux and SOD flux related in the oxidation and
reduction of organic matter in sediment column.


CH4 NH4 1-12S
9 9






54

Typical values of SOD are listed in Table 4.1 (Chapra, 1997). In general, values from

about 1 to 10 gO/m2-day are considered indicative of enriched sediments. According to this

table, SOD at a mud bottom is higher than SOD at a sandy bottom. In Charlotte Harbor

estuarine system, the sediment type near upper Charlotte Harbor is finer and includes more

clay and mud than those for the other study area. The finer sediment which has clay, silt, and

mud would have higher SOD as more organically rich sediments. From this assumption,

SOD can be applied as a function of sediment bottom type (sediment size).

Table 4.1 Average Values of Oxygen Uptake Rates of Bottom (Chapra, 1997)
SOD rate at 200C (gO2/m-day)
Bottom type and Location V
Average Value Range

Sphaerolitus (10 g-dry wt -2) 7
Municipal sewage sludge:
Outfall vicinity 4 2-10
Downstream of outfall 1.5 1-2
Estuarine mud 1.5 1-2
sandy bottom 0.5 0.2-1
Mineral soil 0.07 0.05-1

The effect of temperature and sediment type on SOD can be represented by

SB(T)= SB,20 ST. 0T-20 (4.19)

where S,20o is an areal SOD rate at 20 C (gO2/m2-day ) which is user defined value; ST is

a fractional coefficient for sediment type. (When sediment type 2, ST=1, when sediment

type >3, ST=0.5) ; 0 is temperature coefficient. Zison et al. (1978) have reported a range of

1.04 to 1.13 for 0. A value 1.065 is commonly employed and is used in this study.

Oxygen is another factor that affects sediment oxygen demand. Sediment oxygen

consumption is reduced as oxygen concentration in the overlying water decrease. Lam et al.

(1984) use a Michealis-Menten relationship to represent the dependence, by a saturation

relationship,








DO
S,(DO) = SB (T) (4.20)
KSOD + DO

where Ksoo is half saturation rate for SOD. Lam et al. (1984) have suggested a value for Ksoo

of 1.4 mg/1.

The decay of substrate is assumed to balance continued settling resulting in a steady-

state sediment concentration of oxygen demand substrate. According to this assumption, the

kinetic equation for sediment oxygen demand is (Ambrose et al., 1994):

aDO SOD (4.21)
at H

where H is water depth (m); and SOD is sediment oxygen demand (as measured), gOzm2-

day. SOD can be calculated as a function of temperature, dissolved oxygen at the water-

sediment interface, and sediment bottom type based on characteristics of SOD measurement.

Exchange of material between the water column and benthic sediment is an important

component of the eutrophication process. Sediment oxygen demand may comprise a

substantial fraction of total system oxygen consumption (Cerco and Cole, 1995).

Oxygen consumption in the sediments depends upon water-column temperature and

oxygen availability, and sediment type. As temperature increases, respiration in the sediment

increases. Sediment oxygen consumption is reduced as oxygen concentration in the

overlying water decreases. Therefore, the kinetic equation for sediment oxygen consumption

(SOC) in sediment column can be represented as (Cerco and Cole, 1995):

aDOsed SOC 1 20 T-20 DOw,,,er (422)
=_ SB,20" .- (4.22)
at H H KSOD + DOw,,e

where SB.2 is a function of sediment bottom type (from table 4.1).

The processes that create sediment oxygen demand are little affected by the






56

concentration of oxygen in the overlying water. When oxygen is unavailable to fulfill

sediment oxygen demand, the demand is exported to the water column. The exported

demand may be in the form of reduced iron, manganese, methane, or sulfide, which are

represented in the model as sediment oxygen demand and provides a function which

computes additional release as oxygen consumption in the sediment is restrained (Cerco and

Cole, 1995).

aDO SOD 1 BOT-20 K (4.23)
S -- (4.23)
at H H B20 Kso + DO
The relationship between sediment oxygen consumption and sediment oxygen

demand is represented in Figure 4.5. The SOD is negligible when DO much higher than

KsOD. When dissolved oxygen is absent from the water column, the maximum oxygen

demand is released to the water as sediment oxygen demand.










--------- SOC (Sediment Oxygen Consumption) at sediment column
2 SOD (Sediment Oxygen Demand) at water column
where S, (T= 200C) = 2 g/m2/day
1.5 KSOD =2 g/m3


1.5



1-




-2 Dissolved Oxygen (g/m3)
--0.5 -\




-1.5 ....--


-2 -


Figure 4.5 Effect of dissolved oxygen on sediment oxygen consumption and SOD release








Photosynthesis and Respiration

The photosynthesis and respiration of phytoplankton can add and deplete significant

quantities of oxygen from natural systems. The produced oxygen concentration by

photosynthesis depends on the form of the nitrogen species accessed for phytoplankton

growth. One mole carbon dioxide can produce one mole oxygen when ammonium is the

nitrogen source, while one mole carbon dioxide produces 1.3 moles oxygen when nitrate is

the nitrogen source, according to Morel's equation (1983)

106CO2 + 16NH4 + H2PO +106H20 -> protoplasm + 10602 +15H+
(4.24)
106CO2 +16NOQ + H2PO4 +122H20 +17H+ -> protoplasm +13802

The simple representation of the respiration process can be used to determine how

much oxygen would be consumed in the decomposition of a unit mass of organic carbon,

6CO2 + 6H2 _0 C6HI206 + 602 (4.25)

The dissolved oxygen-to-carbon ratio in respiration can be calculate from this

equation.


Aoc 6(32) =2.67 gO/gC (4.26)
6(12)

The equation that describes photosynthesis and respiration on dissolved oxygen is:

aDO
t =[(1.3-0.3. PJ ) a Kas Kx Aoc PhyC (4.27)

where P, is nitrogen preference coefficient; Pa is phytoplankton growth rate, which is a

function of the intensity of light, the availability of nutrients, and the ambient temperature;

and K.a and K. are respiration and mortality, which are functions of temperature.

The mass balance equation for dissolved oxygen is written by combining all oxygen

transformation processes.








For water column:
aDO K DO
= + AE (DO DO) K, D .CBOD
at A HCBOD +DO
64 DO SOD (4.28)
--KNN NH4
14 Hurr +DO Az
+ [(1.3 0.3. P, K K,. Aoc PhyC
which include reaeration, oxidation by CBOD, nitrification, sediment oxygen demand, and

photosynthesis and respiration terms.

Fore sediment column:

aDO K DO CBOD 64 KNN DO NH SOC (4.29)
=*-K. D .CBOD--K 4NH4 (4.29)
at HCBOD +DO 14 HuNr +DO Az

which include oxidation by CBOD, nitrification, and sediment oxygen consumption.

4.5 Effects of Temperature and Light Intensity on Water Quality Processes

Temperature and light intensity are the most important parameters for transformation

processes. To achieve a better understanding of water quality processes, it is necessary to

improve the spatial and temporal variation of these parameters.

4.5.1 Temperature

In the nutrient cycle, almost all the reaction parameters are affected by temperature,

such as zooplankton and phytoplankton growth, respiration and mortality, nitrification,

denitrification, NH3 stability, mineralization, oxidation, sediment oxygen demand, and

sorption/desorption reactions. The effect of temperature on reaction rates can be explained

by the Van't Hoff-Arrhenius equation, as follows:

K)- AH (4.30)
dt RT2

where K is reaction rate at temperature T, AH is the amount of heat required to bring the

molecules of the reactant to the energy state required for the reaction, and R is universal gas








constant.

Integrating Equation (4.44) from temperature T, to T2 gives:


K2 exp AH( T ) (4.31)
Ki R-T 2 ) (

where K, and K2 are reaction rate at temperature T, and T2 respectively. The temperature


adjustment function (9 = exp ) is almost constant at the temperature range of
RTIT2



interest (0 ~ 300C), ranging from 1.01 to 1.2. This equation can be rearranged into a more

useful form as:


K(T) = K(Tref)0( (4.32)

where K(T) is the reaction coefficient at temperature T, such as Pa (phytoplankton growth

rate), K,,, (phytoplankton respiration rate), Ka, (phytoplankton mortality rate),p, (zooplankton

growth rate), K, (zooplankton respiration rate), Ku (zooplankton mortality rate), KM

(ammonia instability), KouN (ammonification rate), KNN (nitrification rate), do,

(sorption/desorption rate for organic nitrogen), d,, (sorption/desorption rate for inorganic

nitrogen), Kopm mineralizationn rate), doP (sorption/desorption rate for organic phosphorous),

d,gp (sorption/desorption rate for inorganic phosphorous), KD (oxidation rate),or SOD

(sediment oxygen demand). Each rate term has a unique temperature adjustment function.

Most models, which use exponential temperature functions, assume a reference

temperature of 20 C (Chen and Orlob, 1975; Thomann and Frizpatrick, 1982). Eppley

(1972) showed that an exponential relationship describes the envelope curve of growth rate

versus temperature data. The determination was made with a large number of studies, with








many different species.

4.5.1 Light intensity

Light intensity affects the photosynthesis process and thus the phytoplankton growth

rate. The effects of light intensity on nutrient cycling is often modeled by a light intensity

limiting function as follow (Steele, 1974)


f(I) =-exp 1- (4.33)


where I is the light intensity, and I, is the optimum light intensity for phytoplankton growth

According to the Lambert-Beer equation, the light intensity over the water depth is:

I(z) = I0o exp [-Kd (PAR) z] (4.34)

where I(z) is the light intensity at depth z, and Io is the light intensity at the water surface.

K/PAR) is a function of suspended sediment concentration, algae concentration, and color.

This value was calculated by the light attenuation model, which will be discussed in the next

section.

4.6 Light Attenuation Model

One of the most important variables controlling phytoplankton photosynthesis is

"photosynthetically active radiation" (PAR), or light, in the range of wavelengths from 400-

700nm, which provides the predominant source of energy for autotrophic organisms (Day

et al., 1989). Absorption and scattering of light by water and dissolved and suspended

matter determine the quantity and spectral quality of light at a given depth (Jerlov, 1976;

Prieur and Sathyendranath, 1981), which in tern affect the photosynthesis of aquatic plants.

One way to develop a light attenuation model is to find a simple regression

relationship between Kd (PAR), calculated from light measurements, and water quality






62

measurements collected at the same instant as the Kd (PAR) vales. This type of empirical

model is a simple way to relate light attenuation to water quality, at a certain time. This

method was used by Mcpherson and Miller (1994) in Tampa Bay and Charlotte Harbor. A

physics-based light attenuation model was developed as part of the CH3D-IMS (Sheng et al.

2001c, Christian and Sheng, 2003) and successfully applied to the Indian River Lagoon. This

light model is adopted for the Charlotte Harbor study.

In the light model, the vertical light attenuation coefficient (Kd) is a function of solar

zenith angle (/uo), scattering (b) and total absorption (a,) (Kirk, 1984)

1 [ + G \ a l1/2
Kd =- 2 +G(o)-a,-b]
A) (4.35)
and G( o) = io 2

with g, =0.473 and g2=0.218 determined for the mid point of euphotic zone (Kirk, 1984).

The scattering coefficient, b, is determined only using particles since scattering due

to particles is much greater than scattering due to water (Gallegos, 1994). Scattering can be

described as a function of turbidity (Morel and Gentili, 1991) as follow:

b(A) = (550/2) [Turbidity] (4.36)

Total absorption (a,) is partitioned into attributes of water (a,), phytoplankton (aph),

dissolved color (adr), and detritus (ad).

a, = a + aph+ adc + ad (4.37)

The absorption of water (am) can be determined from literature values (Smith and

Baker, 1981) with 1 nm linearly interpolated from 5 nm found in the literature. Chlorophyll-

specific absorption (aph) is calculated for the model from the linear relationship between the

maximum absorption and analytical chlorophyll_a concentration (Dixon and Gary, 1999):








aph () = (aph )MAX Normalized Spectra(A)
(4.38)
(aph MAX = 0.0209. (Chlorophyll a, corrected)

where (aph)MAX is the maximum absorption of chlorophyll-a. To calculate normalized_spectra

(X), individual spectra are normalized to the maximum absorption (437 440 nm), and

averaged for all samples for an overall normalized spectra.

The absorption by dissolved color, adc, for each wavelength in the visible spectrum

can be found using a negative exponential function (Bricaud et al., 1981)

ad () = 440 exp [-Sdc (A 440)] (4.39)

where g440 is the absorption by dissolved color at 440 nm and sd, is spectral slope. Dixon

and Gary (1999) calculated empirical absorption at 440 nm and spectral slope as a function

of color (in PCU) at Charlotte Harbor in the form:

g440= 0.0667 [color], n = 129, r2 = 0.9329
(4.40)
sd = 0.00003. [color]- 0.0178, n = 129, r2 = 0.5111

Absorption due to organic and mineral detritus is represented as a function of

turbidity (Gallegos, 1994):

ad = d (A) [Turbidity] (4.41)

where turbidity is in NTUs and ad(X) is the wavelength specific absorption cross section of

turbidity as calculated in:

d (A) = b + 400 exp [-sd (A 400)] (4.42)

in which aOb is the longwave absorption cross section, 400oo is the maximum detritus

absorption at 400 nm, and sd is exponential slope.

The total absorption and scattering are used in Equation (4.49) to calculate the








vertical attenuation coefficient, KjA), for each wavelength depending on color,

chlorophylla and turbidity. This KjA) value and incident irradiation E(/A) can be used in

the equation for calculating irradiance Ez(A) at the reference depth, Zr:


Ez(2) = Eo(A)-exp[-Kd (A). Z


(4.43)


The spectrum of incident sunlight data from table F-200 in Weast (1977) is used for

incident spectral information, Eo(), as in Gallegos's work (1994). These data are shown

in Table 4.2.

Table 4.2 The spectrum of incident sunlight data (Gallegos, 1994)
X (nm) Eo A (nm) Eo A (nm) Eo

400 4.780
405 5.568 505 8.108 505 8.332
410 6.003 510 8.026 510 8.340
415 6.052 515 7.894 515 8.323
420 6.135 520 7.970 520 8.305
425 6.017 525 8.130 525 8.289
430 5.893 530 8.163 530 8.271
435 5.940 535 8.133 535 8.267
440 6.659 540 8.051 540 8.263
445 7.152 545 7.993 545 8.239
450 7.548 550 7.993 550 8.213
455 7.826 555 7.982 555 8.207
460 7.947 560 7.937 560 8.201
465 7.963 565 8.055 565 8.180
470 7.990 570 8.160 570 8.157
475 8.119 575 8.265 575 8.136
480 8.324 580 8.318 580 8.114
485 8.014 585 8.375 585 8.106
490 7.990 590 8.387 590 8.089
495 8.113 595 8.369 595 8.052
500 8.119 600 8.359 600 8.013

Eo(A) and Ez(A) are integrated over the visible spectrum to get PARo for the incident


PAR and PARz for the calculated PAR at the reference depth. The spectrally sensitive






65

attenuation coefficient Kd (PAR) is calculated from these integrated values by using the


rearranged form of the Lambert-Beer equation:

1 PAR
Kd(PAR)= ln PA (4.44)
zr PARo

This Kd (PAR) will be used in the model to calculate light levels throughout the


water column as a function of measured incident light intensity. Dixon and Gray (1999)

calibrated the light attenuation model with measured data and applied the model to determine

the light requirements for seagrasses of the Charlotte Harbor estuarine system. The results

of the optical model are very good, with mean percentage agreement between modeled and

observed Kd of 110%. Using water quality data from all stations as inputs, Dixon and Gray

(1999) determined that chlorophyll, color, and turbidity account for 4%, 66%, and 31% of

water column light attenuation, respectively. The maximum annual average chlorophyll

contribution is 6%, with an individual maximum of 18% during a phytoplankton bloom.

Color dominated water column attenuation, ranging from 40% to 78%. Stations in the lower

Harbor and southern sites showed increased attenuation due to turbidity, compared to the

upper Harbor site (Dixon and Gary, 1999). For the southern sites, a much larger portion of

light attenuation is produced by turbidity, up to 55% for the station near Captiva pass.

To provide a light attenuation component which can be coupled with water quality,

hydrodynamic, and sediment model components in the Charlotte Harbor estuarine system.

A stand alone light model needs to be developed and calibrated with measured data. The

stand alone light model has been calibrated by finding best fit between simulated and

measured light attenuation with various coefficient sets.

For calibration of the light model, orb, o40o, Sd, and sy are allowed to vary within






66

certain ranges. The predicted k/PAR) values for the stand alone light model are compared

to the corresponding k/PAR) values from data provided by the SFWMD. If the predicted

values do not match data well, then attempts are made to try to find model coefficients to

produce better fit. In this case, the model is run and the coefficients are varied between

maximum and minimum literature values specified in Table4.3 (Christian and Sheng, 2003).

More than one thousand runs of light model simulations were conducted with this data set.

The RMS errors were calculated for each run to test model performance of each set of light

model coefficients. The root mean square error is an indication of the average discrepancy

between observations and model results. In addition to root mean square error, model

calibration was assessed via plots of model output and observation. Scatter plots of model

output and observed data provide an indication of overall model performance. The best

RMS errors for each of the tests is 0.6111 -1 with best fit coefficients shown in Table 4.4

Dixon and Gray (1999) found the coefficients in Table 4.5 provided the best fit to all the data

they examined. They applied uroo and sy as a function of color. The RMS error with this

coefficient is 0.6511 m'. Even though, they reduced adjustable coefficients with relationship

of 4oo and sy with color, the RMS error is greater than that for adjusting all four coefficients.

In this study, the best fit coefficients were used for simulating light attenuation. Figure 4.6

contains calibration period scattering plots with best fit of coefficient sets. The location of

circles indicate the correlation between model predictions and observed data. A perfect

match between model and observed data is indicated by the diagonal line on each graph.

Circle above the line shows over prediction, while circles below the line indicate that under

prediction.









Table 4.3 Coefficient ranges for use in stand along light model.
Coefficients Minimum value Maximum value

Obi 0.004 0.090

0400 0.255 0.600
Sd 0.011 0.019

S_ 0.011 0.019

Table 4.4 Best fit light model coefficients for Charlotte Harbor estuarine system.
Coefficient Value

O.bl 0.01 m-'NTU-1

a0 0.59 m'NTU-'

sd 0.014 nm'

S 0.013 nm-'

Table 4.5 Dixon and Gray's model coefficients for the Charlotte Harbor estuarine system.
Coefficient Value

abi 0.064 m-'NTU-'

0"400 0.0014xcolor+0.0731 m-'NTU'
sd 0.01125 nm-

Sy 0.00003xcolor+0.0178 nm-







68






PARPS model
......... perfect model

5 -



46
"


*





0 .lrl"





1 2 3 4 5 6

KpAR observed

Figure 4.6 The scatter plots for Kd(PAR) during calibration period with best fit coefficients.






69

The water quality and sediment model outputs needed for the light model are

turbidity in NTUs, chlorophylla concentration in l.g/L, and color in Pt. Units. To calculate

solar angle, this model needs the latitude and simulation day and time. Because the sediment

model simulates the total suspended solids (TSS), it is necessary to develop a regression

between TSS (mg/L) and turbidity (NTU) as follow:

Turbidity = 0.2677 x TSS + 0.9665 for 1996 data
Turbidity = -0.0153 x TSS + 2.602 for 2000 data

The stand alone light model was tested, calibrated, and then coupled with the models

of hydrodynamics, sediment and water quality.

4.7 Model Parameters and Calibration Procedures

Reaction terms described in the previous section contain many model parameters

which must be determined before using the model. The determination of model parameters

is generally very difficult because they depend on many physical and biochemical factors,

such as the location of the estuary, temperature and tidal variation, and point or non-point

source loadings of nutrients and other chemical materials. In practice, parameters are

selected from a range of feasible values, tested in the model, and adjusted until an optimal

agreement between simulated and measured values is obtained. Two ways to determine

feasible ranges of model parameters are field observation and laboratory experimentation.

When field observations and laboratory experimentation are not available, the feasible ranges

are obtained from literature or previous modeling studies. A list of the important kinetic

parameters and their literature values are given in Table 4.6 4.11.








Table 4.6 Temperature adjustment functions for water quality parameters
Paramete description Unit literature Reference
r I__ value
(phy) temperature function for 1.01-1.2 Di Toro et al. (1980)
phytoplankton growth 1.09 Pribble et al.(1997)
1.08 Sheng et al. (2001)
(. ) temperature function for 1.045 Ambrose (1991)
kas phytoplankton 1.05 Pribble et al.(1997)
respiration and mortality 1.08 Sheng et al. (2001)
(8 ) temperature function for 1.01-1.2 Di Toro et al. (1980)
OO zooplankton growth, 1.04 Sheng et al. (2001)
respiration and mortality

(QONM temperature function for 1.0 ~ 1.04 Bowie et al. (1985)
ammonification 1.07 Pribble et al.(1997)
1.02 Sheng et al. (2001)
(QNN) temperature function for 1.02 ~ 1.08 Bowie et al. (1985)
/ nitrification 1.08 Pribble et al.(1997)
1.08 Sheng et al. (2001)
(ON temperature function for 1.02 1.09 Bowie et al. (1985)
/ denitrification 1.04 Pribble et al.(1997)
1.045 Sheng et al. (2001)

(OAI temperature function for 1.08 Sheng et al. (2001)
ammonia instability

(OPM temperature function for 1.08 Sheng et al. (2001)
mineralization

(OSD) temperature function for 1.08 Sheng et al. (2001)
sorption/desorption
(oD temperature function for 1.02-1.15 Bowie et al. (1985)
oxidation 1.08 Sheng et al. (2001)

(Oso temperature function for 1.045 Bowie et al. (1985)
sediment oxygen 1.08 Sheng et al. (2001)
demand








Table 4.7 Water quality parameters related to conversion rate
Paramete description Unit literature Reference
r __ value__

Aoc Phytoplankton Carbon gO2/gC 2.67 Ambrose(1991)
/ Oxygen rate 2.67 Cerco and Thomas
2.67 (1995)
Sheng et al. (2001)
ChlaC Phytoplankton Carbon gC/gChl 10-~ 112 Bowie et al. (1985)
/ Chlorophyll_a rate a 100 Pribble et al.(1997)
60 Cerco and Thomas
50 (1995)
Sheng et al. (2001)

ACN Phytoplankton Carbon gN/gC 0.05 -0.43 Jorgensen (1976)
/ Nitrogen rate 0.15 Pribble et al.(1997)
0.167 Cerco and Thomas
0.15 (1995)
Sheng et al. (2001)

Ac Phytoplankton Carbon gP/gC 0.005-0.03 Jorgensen (1976)
/ Phosphorous rate 0.027 Cerco and Thomas
0.025 (1995)
Sheng et al. (2001)

Table 4.8 Water quality parameters related to phytoplankton and zooplankton
Paramete description Unit literature Reference
r ___ value

(lp ) maximum phytoplankton I/day 0.2 ~ 8 Bowie et al. (1985)
max growth rate 2.25 -2.5 Cerco and Thomas
1.06-~ 2.68 (1995)
Sheng et al. (2001)
H, Nitrogen half saturation pg/l 1.5 ~ 400 Bowie et al. (1985)
rate for phytoplankton 0.5 Pribble et al.(1997)
uptake 1 Cerco and Thomas
10 (1995)
Sheng et al. (2001)
H Phosphorous half P.g/l 1. 105 Bowie et al. (1985)
saturation rate for 1 Pribble et al.(1997)
phytoplankton uptake 1 Cerco and Thomas
2-4 (1995)
Sheng et al. (2001)









p optimum light intensity ly/day 225-600 Canale et al. (1976)
for phytoplankton 300 Sheng et al. (2001)
growth

K, phytoplankton I/day 0.02- 0.24 Jorgenson (1976)
respiration rate 0.03- 0.09 Cerco and Thomas
0.03- 0.05 (1995)
Sheng et al. (2001)

Kas phytoplankton non- l/day 0.01-~ 0.22 Jorgenson (1976)
predator mortality 0.03- 0.09 Cerco and Thomas
0.02- 0.06 (1995)
Sheng et al. (2001)

ws phytoplankton settling m/da 0. 3. Bowie et al. (1985)
rate y 0. ~ 0.25 Cerco and Thomas
0.05 0.1 (1995)
Sheng et al. (2001)

Hphy phytoplankton half pg/l 200-2000 Bowie et al. (1985)
saturation rate for 800-1200 Sheng et al. (2001)
zooplankton uptake

TrSphy phytoplankton threshold pg/l 1-200 Bowie et al. (1985)
for zooplankton uptake 200 Sheng et al. (2001)

()max maximum zooplankton I/day 0.15 -0.5 Bowie et al. (1985)
max growth rate 0.18 -0.2 Sheng et al. (2001)

K zooplankton respiration I/day 0.003-0.07 Bowie et al. (1985)
rate 5 Sheng et al. (2001)
0.01

Ks zooplankton non- I/day 0.001-0.36 Jorgensen (1976)
predator mortality 0.015-0.05 Sheng et al. (2001)
5

Table 4.9 Water quality parameters in the nitrogen dynamics
Paramete description Unit literature Reference
r __value

KONM Ammonification rate I/day 0.001-0.4 Bowie et al. (1985)
0.1 Pribble et al. (1997)
0.015 Cerco and Thomas
0.01 (1995)
Sheng et al. (2001)









KNN nitrification rate I/day 0.004-0.11 Bowie et al. (1985)
0.08 Pribble et al.(1997)
0.07 Cerco and Thomas
0.01-0.02 (1995)
Sheng et al. (2001)

H,0t DO saturation rate for mg/1 0.1-2.0 Ambrose (1994)
nitrification 2 Pribble et al.(1997)
2 Sheng et al. (2001)

KDN denitrification rate I/day 0.02- 1.0 Bowie et al. (1985)
0.09 Pribble et al.(1997)
0.09 Sheng et al. (2001)

Hno3 DO saturation rate for mg/1 0.-~ 1.0 Bowie et al. (1985)
denitrification 0.1 Sheng et al. (2001)

p partition coefficient of I.E-5 Simon (1989)
PON/SON 1.E-6-9E-6 Sheng et al. (2001)

P partition coefficient of 5.E-6-1.E- Simon (1989)
PIN/NH4 5 Sheng et al. (2001)
3.E-5-4.E-
3

don sorption/desorption rate I/day 0.02 Bowie et al. (1985)
for SON/PON 0.08 Cerco and Thomas
0.01-0.02 (1995)
Sheng et al. (2001)

dan sorption/desorption rate l/day 0.01 Sheng et al. (2001)
for PIN/NH4

KpDN preference partition 0.5 Cerco and Thomas
coefficient of mortality 0.5 (1995)
for SON / PON Sheng et al. (2001)

Table 4.10 Water quality parameters in the phosphorous dynamics
Paramete description Unit literature Reference
r __ value __

KopM mineralization rate l/day 0.001 0.6 Bowie et al. (1985)
2.27 Pribble et al. (1997)
0.1 Cerco and Thomas
0.1 (1995)
Sheng et al. (2001)

p partition coefficient of 8.E-6 -1.E- Sheng et al. (2001)
POP/SOP 4









p partition coefficient of 1.E-6 6E- Sheng et al. (2001)
PIP/SRP 4

dop sorption/desorption rate I/day 0.08 Cerco and Thomas
for SOP/POP 0.01 (1995)
Sheng et al. (2001)

dap sorption/desorption rate l/day 0.01 Sheng et al. (2001)
for PIP/SRP

PPDP preference partition 0.5 Cerco and Thomas
coefficient of mortality 0.5 (1995)
for SOP/ POP Sheng et al. (2001)

Table 4.11 Water quality parameters in the oxygen balance
Paramete description Unit literature Reference
r value

KD oxidation rate I/day 0.02 0.6 Bowie et al. (1985)
0.05 Sheng et al. (2001)

HcBOD DO half saturation rate mg/1 1.5 -400 Bowie et al. (1985)
for oxidation 0.5 Pribble et al.(1997)
1 Cerco and Thomas
10 (1995)
Sheng et al. (2001)

fdcBoo partition coefficient of 0.5 Bowie et al. (1985)
particular/dissolved 0.3-0.5 Sheng et al. (2001)
CBOD
SOD Sediment oxygen gO2/m 0.02-10. Thomann (1972)
demand 2-day 0.0 10.7 Bowie et al. (1985)

HSOD DO half saturation rate mg/1 0.01- 0.22 Jorgenson (1979)
for sediment oxygen 0.03- 0.09 Cerco and Thomas
demand 0.02- 0.06 (1995)
Sheng et al. (2001)
KAE reaeration rate I/day 0. 3. Bowie et al. (1985)
0. 0.25 Cerco and Thomas
0.05 0.1 (1995)
S_ Sheng et al. (2001)







75

Model calibration is the first stage testing or tuning of the model to a field data not

used in the original construction of the model. Such tuning is to include consistent and

rational set of theoretically defensible parameters and inputs (Thomann, 1992). Proper

calibration of the water quality model requires having accurate representation of the inflow

and loads of nutrients into the water body and selecting appropriate model parameters.

The reaction equation shown in previous section is a function of its concentration and

the water quality parameters connected to it by the indicated processes. Within each reaction

equation, there are numerous kinetic parameters and additional parameters. Water quality

model in CH3D-IMS consists of 13 state variable equations with over 40 interrelated

parameters. The interactions of water quality model parameters and the constituent equations

as shown in Table 4.12 clearly indicate the complexity of the calibrating these types of

models. The column of Table 4.12 show that each modeled constituent equation contains

between 2 and 15 different water quality parameters. According to each row, each parameter

can be found in up to 10 different constituent equations. Therefore, changing one parameter

to improve the calibration of one constituent will simultaneously affect many other

constituents. Traditionally, calibration of water quality models has been performed manually

using a trial-and-error parameter adjustment procedure. The process of manual calibration

depending on the number of model parameters and the degree of parameter interaction. With

complexity of water quality model, the traditional process is a very tedious and time

consuming task. It is necessary to apply more systematic and efficient calibration procedure

for reducing calibration time and effort.

Based on the cascading effect of adjusting interrelated parameters, the efficient

calibration of the water quality model should begin with the parameters that affect the more







76

constituents and the more sensitive parameters. First of all, each water quality parameter

can be ranked with these sensitivity and relativity as parameterization. According to this

order, high ranked parameters will be adjusted to reproduce major pattern of all constituents,

and then lower ranked parameters will be calibrated for detail characteristics of local

constituents. This procedure will reduce numerous calibration efforts and increase efficiency

and effectiveness.

The calibration procedure involves optimization of numerical measures (objective

functions) that compare observations of the state of the system with corresponding simulated

values. The most commonly used objective function adopted in calibration is the root mean

square errors between the observed and simulated model response. The root mean square

error is an indication of the average discrepancy between observations and model results.

It is computed as follow:


RMS= (O- P)2 (4.45)
n
where RMS is root mean square error

O is observation

P is model prediction

n is number of observation

In addition with root mean square error, model calibration was assessed via plots of

model output and observation with correlation coefficient R2. Scatter plots of model output

and observed data provide an indication of overall model performance. The correlation

coefficient is defined by









2 _S (XO-P-n .-P)2
SS* SSY,, ( 2-nOP2-n 2)
SS = (0,-)2 = 02 -n 2 (4.46)

S S (p = ( -P)2 p2P2-nP 2
SS= (O, )(P P)= O.-P- n .P

The process of model calibration is illustrated in Figure 4.7. In systematic calibration

procedure, parameters are adjusted according to order of sensitivity and relativity according

to parameterization for optimization of certain criteria (objective functions) that measure the

goodness-of-fit of the simulation model. The process is repeated until a specified stopping

criterion is satisfied.

Formation of a proper framework for systematic calibration involves the following

key elements:

* Sensitivity analysis
* Model parameterization and choice of calibration parameters
* Specification of calibration criteria

The best way to calibrate water quality parameters is the automatic calibration, in

which parameters are adjusted automatically according to a specific search scheme for

optimization of certain calibration criteria. The process is repeated until a specified stopping

criterion is satisfied. In this study, automatic calibration procedures have been developed

and tested, using a Gauss-Newton method.

To test the automatic calibration procedure of the CH3D water quality model, a 1996

baseline run was adopted as representative of the true and valid field condition. All

parameters from baseline run were considered representative of conditions that could

hypothetically exist in the field. A small number of these parameters were perturbed slightly,

and the automatic calibration procedure was employed. If the procedure is correct, these






78

perturbed parameters should converge on pre-perturbed values. The results of the test show

that the changed parameters did converge on pre-perturbed baseline values, with the Gauss-

Newton method. The method is therefore considered valid and accurate.

In real case, a trial run was conducted by using all the measured data to find best fit

parameters. The calibration procedure continuously failed as calibration parameters were

automatically adjusted outside the upper and lower bounds. The accuracy of the model

calibration generated by this procedure, relies heavily on the quality and quantity of field

observation, the model structure, and the nature of the system. Limited data, and

uncertainties in water quality processes caused the procedure to fail to generated a valid

calibration condition. To apply an automatic calibration procedure to this, or other water

quality models, further investigation is required. These investigations should be focus on the

optimization algorithm and refinement of the water quality processes.










Table 4.12 The relationship
p between water quality parameters and model constituent 4


CA CZ PIN NH4 N03 SON PON PIP SRP SOP POP DO CBOD #of
Eqns


I: ~i


. 1 1


# of 15 8 2 10 7 8 7 2 8 8 7 14 7 103
Parameters


AGRM
HALFN
HALFP
KAEX
KAS
WAS
HALFA
TRESHA
ZGRM
KZEX
KZS
SONM
NITR
DENR
PCON
DRON
KPDN
PCAN
DRAN
SOPM
PCOP
DROP
KPDP
PCIP
DRIP
SODM
FD5
AKD
AOC
CHLAC
CAN
ACP
TOPT
OPTL
AKNIT
AKDEN
AKBOD
AKSOD
AKAIR


. t .. .. ...... .t .... .... ... ....... .................. ........


......... ........... .. .......... ...........? .......... ?........... ... .......... ........... .......... ...........t .......... ,.......... .............
.......... .......... ...... .......... .......... .......... .......... .......... .......... .......... ........... .............

......".......... ....... ........L ... ...........i ...........! .......... i.......... i....!...... .......... i.......... i....!.....i... .......

......... i .. ...... .......... ".".......... '..'......... .......... f...........f .......... "..'........ .. f.......... .......... ." .......... .......
a........ ... ........ 4....


...i ." f.......... f.......... "..'......... ;. .......... .. .......... .: .......... ". .......... ". .......... .'.......... .......... .............
.. .. ... .. .. ... ... ..4... .....
.a ~ ~ a .j. .a.. a... ...~a a
.. .. .. .. .. .. . .. .. .


. ........ :....... ...........f ...........f .......... .......... ...........f .......... .. .......... '........... .......... ." .......... ............

................. : .......... ".......... .......... : .......... .......... : .......... : .......... i. ..... ".......... .'.......... ............

........... .......... ........... .... ......... .......... :...........: .......... .......... ........... .......... :.......... ............
1 i 1 1 i i 1 1 1 1



........... ........ .... .......... ......... ........... .......... ........... .......... ........... .......... ............


........... 4.......... 4. ........... .. :. ... ...".......... .......... :.......... -".......... :.......... ...........: .......... i.............
........... ....... .. ......... ...... ... .... .......... ...........i ........... ...........i .......... i........... .....L ... ............
........... .......... -........... .......... .......... ..,......... ,........... ........... ...........i .......... ..... '.......... :............
i i"1 i 1 i i i










........... ...... ............ ....f ...... ....... L.:. ... .....1 ......... .."L ........... : .......... .......... .......... .......
........... 7:........... ........... ...........7 .......... ". T ... T .......... ...........f .......... 7:.......... .: .......... .............
........... .. .. .......... ... ..........! .. ." .f : . .i . .
......... 7 ........ i F T ." ........7. .... ... .......... .......... f.......... f.......... .......... .,.......... ..'...........
........... ........... .......... .......... ........... .......... .. ..........
1. :. ....... .......... .................






........... ,.. ......... .......... ........... ... .........., ........... ..'.......... .......... ,....!....... ......L .. .......... i...........i .............
.......... .." .......... .i.......... .......... I ..... .i.......... ..........." .......... .......... ..1.....i .......i .......... i.............
..~~~~~~~.......... ............ ........ .......... ...........,"............ ...........,".......... t .......... .... i.. .i ..... .... i ............

........... ......... ........... .......... ........... ........... .......... :...........: ...........: .......... ". ........... z .......... :.............
1 1 1 : : : :









........... f.......... f.......... f.......... f......... ......f......... ..... .......... }........... .......... ".. .............
... ....... ........... ........... .......... .' .......... -'........... .......... .......... .......... ...........4 ........... ...... -- .... .....
........... i...........: ...........- .......... .......... i.......... ............ ...........i .......... ...........i .......... i..........1 .............
...........: .......... .......... :.......... .......... ...........: .......... 4.......... .......... 4...........i .......... .......................
1 1 1 1 1 1 i 11 i 1 1











.... *.. i .. ....... ......... .. .......... .......... ,. .......... .......... f......... .. ......... .. .......... *........... ............
.. .. .... ... .. .. . .... .......... ....... .......... f ...........'. ... ..f ......... ...........} ............
........... --.......... ...........- .......... .......... ........... ..........: .......... ...........i .......... .. .......... .. .......... ............
1 i 1 !1 1 : 1











...... .... .......... .... .... ..... ..... .......... .......... ........... .......... ............

...'" f...........f ......... .......... f. .. .... ." ......... "..".......... .......... f...........f .......... }.......... }...... .... ............
i i 1 1 : : 1 :

i 1 1 1











........... .......... .......... .......... ........... .......... .......... .......... .......... .......... .......... 4..... ......
..........i1 .... .... 1.. .............
I I"" "" "
I I"'{ "'






I I"
1 11
I I




111



1 1 1 1 1 1 1 1 1
1 1

1 1 1 1 1 1
1 1 1 1 1

.441A .A 4 4
.. i


5
1
1
5
7
1
2
2
2
3
6
2
3
2
2
2
2
2
2
2
2
2
2
2
2
1
2
2
10
2
6
5
1
1
3
2
2
1
1









INITIAL VALUE for WO parameter


Execute WQ model


Adjust WQ river boundary Condition


Execute WO model


Adjustcheck river BNDparam dt


change i parameter for sensitivity test


Execute WQ model


Errolast parameter?


Sensitivity analysis:
check RMS differneces between base run and each sensitivity test

parameterization :
make order with sensitivity result and relativity table


Adjust i ranked parameter


Execute WO model

Error Check --last ranked parameter.


Calculate statistics : scatter plot, RMS, R2
Plot time series for each species with measured datat
Figure 4.7 Systematic calibration procedure















CHAPTER 5
APPLICATION OF CIRCULATION AND TRANSPORT MODEL

Circulation in Charlotte Harbor is driven by tide, wind, and density gradient. The

numerical model of circulation and transport described in Chapter 3 was applied to the

Charlotte Harbor estuarine system, a large shallow estuary in southwest Florida. Model

applications include 1) the simulation of three-dimensional circulation during summer 1986

2)one year simulation of flow, salinity, temperature, and sediment transport in 2000 and 3)

simulation of the impact of removal of Sanibel Causeway on the circulation in San Carlos

Bay and Pine Island Sound

A major purpose of the first two simulations is the calibration and validation of the

Charlotte Harbor circulation and transport model using field data obtained by USGS, NOAA,

SWFWMD, and SFWMD in Summer 1986 and January to December in 2000. During the

calibration process, a few model coefficients and/or boundary/initial conditions are adjusted

to produce more accurate overall simulation of observed data collected in July 1986 by

USGS (Hammet, 1992; Stoker, 1992; Goodwin, 1996) covered the entire Charlotte Harbor

estuarine system and include water level, horizontal currents and salinity. Hence July 1986

data were used for short-term calibration. To supplement the short-term calibration of July

1986, the long-term model calibration was conducted using one year water level and salinity

data in Caloosahatchee River during 2000.