UFDC Home  UF Institutional Repository  UF Theses & Dissertations  Internet Archive   Help 
Material Information
Record Information

Table of Contents 
Title Page
Page i Acknowledgement Page ii Table of Contents Page iii Page iv Page v List of Figures Page vi Page vii Page viii Page ix Page x Page xi Page xii Page xiii List of Tables Page xiv Page xv Abstract Page xvi Page xvii Chapter 1. Introduction Page 1 Page 2 Page 3 Page 4 Page 5 Page 6 Page 7 Page 8 Page 9 Page 10 Chapter 2. Nutrient dynamics in lakes and estuaries Page 11 Page 12 Page 13 Page 14 Page 15 Page 16 Page 17 Page 18 Page 19 Page 20 Page 21 Page 22 Page 23 Page 24 Page 25 Page 26 Page 27 Page 28 Page 29 Page 30 Page 31 Page 32 Page 33 Page 34 Page 35 Page 36 Page 37 Page 38 Page 39 Page 40 Page 41 Page 42 Page 43 Chapter 3. Hydrodynamics and sediment dynamics in lakes and estuaries Page 44 Page 45 Page 46 Page 47 Page 48 Page 49 Page 50 Page 51 Page 52 Page 53 Page 54 Page 55 Page 56 Page 57 Page 58 Page 59 Page 60 Page 61 Page 62 Page 63 Page 64 Page 65 Page 66 Page 67 Page 68 Page 69 Page 70 Page 71 Page 72 Page 73 Page 74 Page 75 Page 76 Page 77 Page 78 Page 79 Page 80 Page 81 Page 82 Page 83 Chapter 4. Modeling the hydrodynamics and sediment Page 84 Page 85 Page 86 Page 87 Page 88 Page 89 Page 90 Page 91 Page 92 Page 93 Page 94 Page 95 Page 96 Page 97 Page 98 Chapter 5. Modeling nutrient dynamics Page 99 Page 100 Page 101 Page 102 Page 103 Page 104 Page 105 Page 106 Page 107 Page 108 Page 109 Page 110 Page 111 Page 112 Page 113 Page 114 Page 115 Page 116 Page 117 Page 118 Page 119 Page 120 Page 121 Page 122 Page 123 Page 124 Page 125 Chapter 6. Field data collected in Lake Okeechobee Page 126 Page 127 Page 128 Page 129 Page 130 Page 131 Page 132 Page 133 Page 134 Page 135 Page 136 Page 137 Page 138 Page 139 Page 140 Page 141 Page 142 Chapter 7. Model applications to Lake Okeechobee Page 143 Page 144 Page 145 Page 146 Page 147 Page 148 Page 149 Page 150 Page 151 Page 152 Page 153 Page 154 Page 155 Page 156 Page 157 Page 158 Page 159 Page 160 Page 161 Page 162 Page 163 Page 164 Page 165 Page 166 Page 167 Page 168 Page 169 Page 170 Page 171 Page 172 Page 173 Page 174 Page 175 Page 176 Page 177 Page 178 Page 179 Page 180 Page 181 Page 182 Page 183 Page 184 Page 185 Page 186 Page 187 Page 188 Page 189 Page 190 Page 191 Page 192 Page 193 Page 194 Page 195 Page 196 Page 197 Page 198 Page 199 Page 200 Page 201 Page 202 Page 203 Chapter 8. Field data collected in Tampa Bay Page 204 Page 205 Page 206 Page 207 Page 208 Page 209 Page 210 Page 211 Chapter 9. Model applications to Tampa Bay Page 212 Page 213 Page 214 Page 215 Page 216 Page 217 Page 218 Page 219 Page 220 Page 221 Page 222 Page 223 Page 224 Chapter 10. Conclusions and recommendations Page 225 Page 226 Page 227 Page 228 Bibliography Page 229 Page 230 Page 231 Page 232 Page 233 Page 234 Page 235 Page 236 Page 237 Page 238 Page 239 Page 240 Page 241 Page 242 Page 243 Page 244 Page 245 Appendix A. Computing pressure gradients from current data Page 246 Page 247 Appendix B. Flow charts of the 3D model system Page 248 Page 249 Appendix C. Numerical solution algorithms in hydrodynamic and sediment transport models Page 250 Page 251 Page 252 Page 253 Page 254 Page 255 Appendix D. Field data of 1989 spring survey in Lake Okeechobee Page 256 Page 257 Page 258 Page 259 Page 260 Page 261 Page 262 Page 263 Page 264 Page 265 Page 266 Page 267 Page 268 Page 269 Page 270 Page 271 Page 272 Page 273 Page 274 Page 275 Page 276 Page 277 Page 278 Page 279 Page 280 Page 281 Appendix E. Field data of 1993 spring storm event in Lake Okeechobee Page 282 Page 283 Page 284 Page 285 Appendix F. 3 runs for 1993 spring storm event in Lake Okeechobee Page 286 Page 287 Page 288 Page 289 Page 290 Page 291 Page 292 Page 293 Page 294 Appendix G. Field data collected in Tampa Bay Page 295 Page 296 Page 297 Page 298 Page 299 Page 300 Page 301 Page 302 Page 303 Page 304 Page 305 Page 306 Appendix H. Determining model coefficients Page 307 Page 308 Page 309 Biographical sketch Page 310 Page 311 Page 312 Page 313 
Full Text 
EFFECTS OF HYDRODYNAMICS AND SEDIMENT TRANSPORT PROCESSES ON NUTRIENT DYNAMICS IN SHALLOW LAKES AND ESTUARIES By Xinjian Chen 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 ACKNOWLEDGEMENTS First of all, I would like to express my very deep appreciation to the chairman of my advisory committee, Professor Peter Y. Sheng, for his advice, encouragement, and financial support. Special thanks are extended to Professors Ashish J. Mehta, D. Max Sheppard, Robert J. Thieke, and K. Ramesh Reddy for serving as the members of my doctoral advisory committee, revising the dissertation and attending the final exam. Thanks are also due to Professors Robert G. Dean, Hsiang Wang, Michel K. Ochi, and Daniel M. Hanes for their teaching efforts during my study in Gainesville. Countless pieces of help provided by Mr. Subarna Malakar and Mr. Sidney Schofield in solving my computer problems are appreciated. I am also grateful to Professors Werner Zielke and Mark Markofsky of the Institut fuir Str6mungsmechanik der Universitit Hannover for introducing me to the field of sediment transport and water quality modeling and providing financial support during my stay in Germany. Finally, I would like to thank my wife, parents, and parentsinlaw for their love, support and patience during this study. TABLE OF CONTENTS ACKNOWLEDGEMENTS ............................ ii LIST OF FIGURES ................................ vi LIST OF TABLES ................................. xiv ABSTRACT .................................... xvi CHAPTERS 1 INTRODUCTION ............................... 1 1.1 Study Background ............................ 1 1.2 Study Objectives ...... ... ........... ........... 4 1.3 Outline of Presentation .......................... 8 2 NUTRIENT DYNAMICS IN LAKES AND ESTUARIES ......... 11 2.1 Introduction ................ ......... .... ... 11 2.2 Phosphorus Dynamics in Lakes and Estuaries ............... 13 2.3 Nitrogen Dynamics in Lakes and Estuaries ................. 20 2.4 Algal and Zooplankton Growth ......................... 29 2.5 Effects of Temperature and Light Intensity on Nutrient Dynamics . 35 2.6 Effects of Hydrodynamics and Sediment Processes on Nutrient Dynamics 37 2.7 Conclusions ......................... .. ..... 42 3 HYDRODYNAMICS AND SEDIMENT DYNAMICS IN LAKES AND ESTUARIES .................................. 44 3.1 Circulation .. .. . .. .. . .. . .. . .. . .. . 44 3.2 Wave Boundary Layer .......................... 49 3.3 Turbulent Mixing ............................. 55 3.4 Sediment Dynamics ............................ 64 4 MODELING THE HYDRODYNAMICS AND SEDIMENT TRANSPORT PROCESSES ......................... 84 4.1 Governing Equations ........................... 84 4.2 Boundary and Initial Conditions ........................ 88 4.3 Model Parameters ................................. .. 91 4.4 Numerical Solution Algorithm ...................... 98 5 MODELING NUTRIENT DYNAMICS ...................... 99 5.1 Governing Equations ........................... 99 5.2 Boundary and Initial Conditions ....................... 104 5.3 Modeling Nutrient Transformation Processes ............... 107 5.4 Model Parameters. ...... ............................. 113 5.5 Numerical Solution Algorithm ...................... 121 6 FIELD DATA COLLECTED IN LAKE OKEECHOBEE .......... 126 6.1 Physical Setting .............................. 126 6.2 Field Measurements ............................ 127 6.3 1989 Spring Survey ............................ 134 6.4 1993 Spring ShortTerm Survey ....................... 139 7 MODEL APPLICATIONS TO LAKE OKEECHOBEE .......... 143 7.1 Modeling of the 1993 Storm Event .................... 143 7.2 1D Simulations of 1989 Spring Survey .................. 150 7.3 3D Simulations of 1989 Spring Survey .................. 164 7.4 Sum m ary ................................. 200 8 FIELD DATA COLLECTED IN TAMPA BAY ................ 204 8.1 Physical Setting .............................. 204 8.2 Field Measurements ............................ 207 8.3 A Storm Event in February, 1992 .................... 208 8.4 A Storm Event in March, 1993 ...................... 210 9 MODEL APPLICATIONS TO TAMPA BAY .................. 212 9.1 Modeling the 1992 Storm Event ..................... 212 9.2 Modeling the 1993 Storm Event ..................... 216 9.3 Conclusions ................................ 219 10 CONCLUSIONS AND RECOMMENDATIONS ................ 225 10.1 Conclusions ................................ 225 10.2 Recommendations ............................. 228 BIBLIOGRAPHY ................................. 229 APPENDICES A COMPUTING PRESSURE GRADIENTS FROM CURRENT DATA . 246 B FLOW CHARTS OF THE 3D MODEL SYSTEM .............. 248 C NUMERICAL SOLUTION ALGORITHMS IN HYDRODYNAMIC AND SEDIMENT TRANSPORT MODELS ...................... 250 C.1 3D Hydrodynamics Model ........................ 250 C.2 1D Hydrodynamics Model ........................ 252 C.3 3D Sediment Transport Model ....................... 253 C.4 Numerical Stability ........................... 254 D FIELD DATA OF 1989 SPRING SURVEY IN LAKE OKEECHOBEE . 256 E FIELD DATA OF 1993 SPRING STORM EVENT IN LAKE OKEECHOBEE 282 F 3 RUNS FOR 1993 SPRING STORM EVENT IN LAKE OKEECHOBEE 286 G FIELD DATA COLLECTED IN TAMPA BAY ................ 295 H DETERMINING MODEL COEFFICIENTS .................. 307 H.1 1D Simulation of Lake Okeechobee Phosphorus Dynamics ...... 307 H.2 3D Simulation of Lake Okeechobee Phosphorus Dynamics ...... 308 H.3 1D Simulation of Tampa Bay Nitrogen Dynamics . . . . ... ..308 BIOGRAPHICAL SKETCH ........................... 310 LIST OF FIGURES 2.1 Nutrient dynamics and the relevant hydrodynamics and sediment transport processes in lakes and estuaries . . . . . ... ..12 2.2 The twofilm model for the volatilization . . . . . .... .. 29 2.3 Relative growth rate of algae in a phosphorus limiting aquatic system . . . . . . . . . . . . . . . .. .. 31 3.1 Critical shear stresses for erosion with and without wave action (from M aa, 1986) . . . . . . . . . . . .... .. 50 3.2 Settling velocity as a function of suspended sediment concentra tion, with flocculation occurring at high concentrations (after Hwang and Mehta, 1989) . . . . . . . . . . . .... .. 71 3.3 Erosion rate versus bottom shear stress (after Lavelle et al., 1984; the numbers are for different sediments and are explained in Lavelle et al., 1984) . . . . . . . . . . . . . . .. .. 78 4.1 A typical bed density profile from Lake Okeechobee in the vicinity of Station C (after Hwang and Mehta, 1989) . . . . ... .. 96 6.1 Major features of Lake Okeechobee (after Sheng et al., 1989a). 128 6.2 Bathymetry of Lake Okeechobee (after Sheng et al., 1989a). . 129 6.3 Bottom sediment types of Lake Okeechobee (after Reddy and Graetz, 1989) . . . . . . . . . . . . .... .. 130 6.4 Mud Thickness in Lake Okeechobee (after Hwang and Mehta, 1989).131 6.5 A portable platform for field measurements . . . . .... .. 133 6.6 Locations of synoptic stations in Lake Okeechobee during the 1989 Spring Survey . . . . . . . . . . . . ... .. 136 6.7 Spectra of measured current, wind speed, and SSC at Station C in Lake Okeechobee . . . . . . . . . . . .... .. 140 7.1 Model results of uvelocities, vvelocities, and suspended sediment concentrations at the middle layer (Layer 2) and the top layer (Layer 3) at Station L002 in Lake Okeechobee during the storm event in February, 1993 . . . . . . . . . . ... .. 145 7.2 Correlation coefficients between DO and SRP, DO and TDP, and DO and TP measured at Station L002 in Lake Okeechobee during the storm event in February, 1993 . . . . . . . ... .. 148 7.3 Measured and simulated SRP of a run with adjusted pH at the top layer at Station L002 in Lake Okeechobee during the storm event in February, 1993 . . . . . . . . . . ... .. 151 7.4 Measured and simulated TDP of a run with adjusted pH at the top layer at Station L002 in Lake Okeechobee during the storm event in February, 1993 . . . . . . . . . . ... .. 152 7.5 Measured and simulated TP of a run with adjusted pH at the top layer at Station L002 in Lake Okeechobee during the storm event in February, 1993 . . . . . . . . . . . .... .. 153 7.6 Measured and simulated SSC at Station C during week 2 of the 1989 Spring Survey in Lake Okeechobee . . . . . ... .. 155 7.7 Measured and simulated SSC at Station D during week 2 of the 1989 Spring Survey in Lake Okeechobee . . . . . ... .. 156 7.8 Measured and simulated SSC at Station E during week 2 of the 1989 Spring Survey in Lake Okeechobee . . . . . ... .. 157 7.9 Settling velocity obtained from laboratory by Hwang and Mehta (1989) and the modified settling velocity considering the shear stress effect . . . . . . . . . . . . . ... .. 160 7.10 Comparison of model results of suspended sediment concentration using two different settling velocity formula . . . . ... .. 161 7.11 Results of Monte Carlo simulations of SSC at Station C during Week 2 of the 1989 Spring Survey in Lake Okeechobee.. . . 163 7.12 Contours of the rootmeansquare error of the Monte Carlo sim ulations of SSC at Station C during week 2 of the 1989 Spring Survey in Lake Okeechobee . . . . . . . . . ... .. 164 7.13 Results of the lutocline simulation at Station C during a threeday event in the 1989 Spring Survey in Lake Okeechobee . . ... .. 165 7.14 SSC profiles in the lutocline simulation at Station C during a 3day event in Week 2 of the 1989 Spring Survey in Lake Okeechobee. 166 7.15 A Cartesian grid system for the 3D simulation of Lake Okeechobee hydrodynamics, sediment, and nutrient dynamics . . ... .. 167 7.16 Measured and predicted significant wave height and period at Sta tion C in Lake Okeechobee. Model parameters in SMB model has been adjusted . . . . . . . . . . . . ... .. 170 7.17 Contours of waveinduced bottom shear stress from the 1D TKE m odel . . . . . . . . . . . . . . . . .. 172 7.18 Comparison of model (1D TKE) results of bottom shear stresses induced by a single wave to those induced by group waves. . 174 7.19 Simulated uvelocities, vvelodcities, and suspended sediment con centrations at Station C of Lake Okeechobee form May 21 to June 16, 1989 . . . . . . . . . . . . ...... .. 176 7.20 Simulated uvelocities, vvelocities, and suspended sediment con centrations at Station D of Lake Okeechobee form May 21 to June 16, 1989 . . . . . . . . . . . . . . . .. 177 7.21 Simulated uvelocities, vvelocities, and suspended sediment con centrations at Station E of Lake Okeechobee form May 21 to June 16, 1989 . . ........... ................... .. 178 7.22 Spatial distribution of simulated suspended sediment concentra tion at the top and the bottom layers of Lake Okeechobee at noon of May 27, 1989 . . . . . . . . . . . . .... .. 179 7.23 Spatial distribution of simulated suspended sediment concentra tion at the top and the bottom layers of Lake Okeechobee at noon of June 3, 1989 . . . . . . . . . . . . ... .. 180 7.24 Spatial distribution of simulated suspended sediment concentra tion at the top and the bottom layers of Lake Okeechobee at noon of June 10, 1989 . . . . . . . . . . . . .... .. 181 7.25 Spatial distribution of simulated suspended sediment concentra tion at the top and the bottom layers of Lake Okeechobee at noon of June 16, 1989 . . . . . . . . . . . . .... .. 182 7.26 Spatial distribution of simulated SRP concentration at the top and the bottom layers of Lake Okeechobee at noon of May 27, 1989. 183 7.27 Spatial distribution of simulated SRP concentration at the top and the bottom layers of Lake Okeechobee at noon of June 3, 1989. 184 7.28 Spatial distribution of simulated SRP concentration at the top and the bottom layers of Lake Okeechobee at noon of June 10, 1989. 185 7.29 Spatial distribution of simulated SRP concentration at the top and the bottom layers of Lake Okeechobee at noon of June 16, 1989. 186 7.30 Spatial distribution of simulated TP concentration at the top and the bottom layers of Lake Okeechobee at noon of May 27, 1989. 187 7.31 Spatial distribution of simulated TP concentration at the top and the bottom layers of Lake Okeechobee at noon of June 3, 1989. 188 7.32 Spatial distribution of simulated TP concentration at the top and the bottom layers of Lake Okeechobee at noon of June 10, 1989. 189 7.33 Spatial distribution of simulated TP concentration at the top and the bottom layers of Lake Okeechobee at noon of June 16, 1989. 190 7.34 Correlations among SSC, SRP, and TP in the synoptic data and the 3D model results for the 1989 synoptic survey in Lake Okee chobee . . . . . . . . . . . . . . . .. .. 193 7.35 Time series of simulated SSC, SRP, and TP at Station C of Lake Okeechobee during the 1989 Spring Survey . . . . ... .. 195 7.36 Model results of time series of sediment and PIP resuspension/deposition fluxes and SRP mass in Lake Okeechobee during the 1989 Spring Survey . . . . . . . . . . . . . . . ... .. 196 7.37 Spatial distribution of simulated (with the WASP transforma tions) SRP concentration at the top and the bottom layers of Lake Okeechobee at noon of May 27, 1989 . . . . . . ... .. 199 7.38 Spatial distribution of simulated SRP concentration (by the 3D model developed by Dickinson et al, 1992) at the top and the bottom layers of Lake Okeechobee at noon of May 27, 1989. . 201 8.1 Location of Tampa Bay and the experiment stations (after Sheng et al., 1993b) . . . . . . . . . . . . . ... .. 205 9.1 Significant wave heights calculated from SMB model for the time period of the 1992 storm event . . . . . . . . .... .. 213 9.2 Model results of uvelocity, vvelocity, and sediment concentration at Station A during the storm event of 1992 . . . . ... .. 215 9.3 Comparison of model results of NH+N and measured data during the storm event of 1992 at Station A . . . . . . ... .. 217 9.4 Model results of resuspension fluxes of (a) sediments and (b) NH+ N during the storm event of 1992 at Station A . . . ... .. 218 9.5 Model results of uvelocity, vvelocity, and sediment concentration at Station A during the storm event of 1993 . . . . ... .. 220 9.6 Model results of uvelocity, vvelocity, and sediment concentration at Station B during the storm event of 1993 . . . . ... .. 221 9.7 Comparison of model results of NH+N and measured data during the storm event of 1993 at Station A . . . . . . ... .. 222 9.8 Model results of resuspension fluxes of (a) sediments and (b) NH+ N during the storm event of 1993 at Station A ... . . .. 223 9.9 A test run with zero desorption rate of ammonium N at Station A during the storm event of 1992 . . . . . . . ... .. 224 B.1 Flow chart of the 3D model system of hydrodynamics, sediment transport, and nutrient dynamics in lakes and estuaries . . 248 B.2 Flow chart of the 3D submodel for nutrient dynamics . ... .. 249 D.1 Measured wind at Station A in Lake Okeechobee from May 21 to June 16, 1989 . . . . . . . . . . . . .... .. 256 D.2 Measured wind at Station B in Lake Okeechobee from May 21 to June 16, 1989 . . . . . . . . . . . . .... .. 257 D.3 Measured wind at Station C in Lake Okeechobee from May 21 to June 16, 1989 . . . . . . . . . . . . .... .. 258 D.4 Measured wind at Station D in Lake Okeechobee from May 21 to June 16, 1989 . . . . . . . . . . . . .... .. 259 D.5 Measured wind at Station E in Lake Okeechobee from May 21 to June 16, 1989 . . . . . . . . . . . . .... .. 260 D.6 Measured uvelocities, vvelocities, and suspended sediment con centrations at Station A from May 21 to June 16, 1989 . . 261 D.7 Measured uvelocities, vvelocities, and suspended sediment con centrations at Station B from May 21 to June 16, 1989 . . 262 D.8 Measured uvelocities, vvelocities, and suspended sediment con centrations at Station C from May 21 to June 16, 1989 . . 263 D.9 Measured uvelocities, vvelocities, and suspended sediment con centrations at Station D from May 21 to June 16, 1989 . . 264 D.10 Measured uvelocities, vvelocities, and suspended sediment con centrations at Station E from May 21 to June 16, 1989 . . 265 D. 1 Measured uvelocities, vvelocities, and suspended sediment con centrations at Station F from May 21 to June 16, 1989 . . 266 D.12 Contours of synoptic data of suspended sediment concentrations at the top and the bottom layers of Lake Okeechobee at noon of M ay 21, 1989 . . . . . . . . . . . . ... .. 267 D.13 Contours of synoptic data of suspended sediment concentrations at the top and the bottom layers of Lake Okeechobee at noon of M ay 27, 1989 . . . . . . . . . . . . ... .. 268 D.14 Contours of synoptic data of suspended sediment concentrations at the top and the bottom layers of Lake Okeechobee at noon of June 3, 1989 . . . . . . . . . . . . . .... .. 269 D.15 Contours of synoptic data of suspended sediment concentrations at the top and the bottom layers of Lake Okeechobee at noon of June 10, 1989 . . . . . . . . . . . . .... .. 270 D.16 Contours of synoptic data of suspended sediment concentrations at the top and the bottom layers of Lake Okeechobee at noon of June 16, 1989 . . . . . . . . . . . . . ... .. 271 D.17 Contours of synoptic data of SRP concentrations at the top and the bottom layers of Lake Okeechobee at noon of May 21, 1989. 272 D.18 Contours of synoptic data of SRP concentrations at the top and the bottom layers of Lake Okeechobee at noon of May 27, 1989. 273 D.19 Contours of synoptic data of SRP concentrations at the top and the bottom layers of Lake Okeechobee at noon of June 3, 1989. 274 D.20 Contours of synoptic data of SRP concentrations at the top and the bottom layers of Lake Okeechobee at noon of June 10, 1989. 275 D.21 Contours of synoptic data of SRP concentrations at the top and the bottom layers of Lake Okeechobee at noon of June 16, 1989. 276 D.22 Contours of synoptic data of TP concentrations at the top and the bottom layers of Lake Okeechobee at noon of May 21, 1989. . 277 D.23 Contours of synoptic data of TP concentrations at the top and the bottom layers of Lake Okeechobee at noon of May 27, 1989. . 278 D.24 Contours of synoptic data of TP concentrations at the top and the bottom layers of Lake Okeechobee at noon of June 3, 1989. . 279 D.25 Contours of synoptic data of TP concentrations at the top and the bottom layers of Lake Okeechobee at noon of June 10, 1989. . 280 D.26 Contours of synoptic data of TP concentrations at the top and the bottom layers of Lake Okeechobee at noon of June 16, 1989. . 281 E.1 Measured wind during a storm event at Station L002 in Lake Okee chobee in February, 1993 . . . . . . . . . ... .. 282 E.2 Measured uvelocities, vvelocities, and suspended sediment con centrations at the middle and the bottom arms at Station L002 in Lake Okeechobee in February, 1993 . . . . . . ... .. 283 E.3 Measured total phosphorus (TP), total dissolved phosphorus (TDP), and soluble reactive phosphorus (SRP) at Station L002 in Lake Okeechobee in February, 1993 . . . . . . . . ... .. 284 E.4 Measured DO and pH at the bottom layer at Station L002 in Lake Okeechobee in February, 1993 . . . . . . . . .... .. 285 F.1 Measured and simulated SRP of the first run at Station L002 in Lake Okeechobee during the storm event in February, 1993. . 286 F.2 Measured and simulated TDP of the first run at Station L002 in Lake Okeechobee during the storm event in February, 1993. . 287 F.3 Measured and simulated TP of the first run at Station L002 in Lake Okeechobee during the storm event in February, 1993. . 288 F.4 Measured and simulated SRP of the second run at Station L002 in Lake Okeechobee during the storm event in February, 1993. . 289 F.5 Measured and simulated TDP of the second run at Station L002 in Lake Okeechobee during the storm event in February, 1993. . 290 F.6 Measured and simulated TP of the second run at Station L002 in Lake Okeechobee during the storm event in February, 1993. . 291 F.7 Measured and simulated SRP of the third run at Station L002 in Lake Okeechobee during the storm event in February, 1993. . 292 F.8 Measured and simulated TDP of the third run at Station L002 in Lake Okeechobee during the storm event in February, 1993. . 293 F.9 Measured and simulated TP of the third run at Station L002 in Lake Okeechobee during the storm event in February, 1993. . 294 G.1 Measured uvelocites and vvelocities at Station A in Tampa Bay during Julian days 30 to 40, 1992 . . . . . . . ... .. 295 G.2 Measured wind at Station A in Tampa Bay during Julian days 30 to 40, 1992 . . . . . . . . . . . . . . .. .. 296 G.3 Measured free surface variation about the mean water level, salin ity, and temperature at Station A in Tampa Bay during Julian days 30 to 40, 1992 . . . . . . . . . . . .... .. 297 G.4 Measured uvelocity, vvelocity, and suspended sediment on the storm day at Station A of Tampa Bay during Julian day 35.6 to 36.6, 1992 . . . . . . . . . .. ........... .. 298 G.5 Measured SRP, TP, and ammonium N concentrations at Station A in Tampa Bay during Julian days 35.6 to 36.6, 1992 . ... .. 299 G.6 Measured SRP, NH+N, and redox profiles in Tampa Bay. . 300 G.7 Measured salinity, temperature, pH, and dissolved oxygen concen tration at Station A in Tampa Bay during Julian days 70.59 to 71.1, 1993 . . . . . . . . . . . . . . .. .. 301 G.8 Measured wind at Station A in Tampa Bay during Julian days 70.59 to 71.1, 1993 . . . . . . . . . . . ... .. 302 G.9 Measured uvelocity, vvelocity, and suspended sediment concen tration at Station A in Tampa Bay during Julian days 70.59 to 71.1, 1993 . . . . . . . . . . . . . . . .. 303 G.10 Measured wind at Station A in Tampa Bay during Julian days 70.5 to 76.3, 1993 . . . . . . . . . . . .... .. 304 G.11 Measured uvelocity, vvelocity, and suspended sediment concen tration at Station B in Tampa Bay during Julian days 70.5 to 76.3, 1993 . . . . . . . . . . . . . . . . .. .. 305 G.12 Measured ammonium N and SRP concentrations at Station A in Tampa Bay during Julian days 70.5 to 71.1, 1993 . . ... .. 306 LIST OF TABLES 3.1 Measured exponents m in Krone's Equation . . . . .... .. 69 4.1 Values of PB, eM, and rT determined by Hwang and Mehta (1989). indicates that these values were obtained by combining the ero sion rate data resulting from tests 1, 2 and 3 . . . . ... .. 95 5.1 Literature values for Hp and H. . . . . . . . ... .. 116 5.2 Literature values for H . . . . . . . . . ... . 116 5.3 Literature values for w . . . . . . . . . . . .. 117 5.4 Literature values for DOP mineralization rate . . . ... .. 118 5.5 Literature values for K2, K3, and K4 . . . . . . ... .. 120 6.1 The synoptic surveys conducted in the spring of 1989 in Lake Okeechobee . . . . . . . . . . . . . ... .. 135 6.2 Locations and depths of the stations in Lake Okeechobee ..... 135 6.3 Elevations (above lake bottom) of the instruments arms at the six stations . . . . . . . . . . . . . . . .. 137 7.1 rc, and Eo values determined by the 1D sediment model . . 154 7.2 Correlation coefficients between measured and simulated SSC, SRP, and TP. The effect of water depth variation with time on wave induced bottom shear stress was considered . . . . ... .. 192 7.3 Correlation coefficients between measured and simulated SSC, SRP, and TP. The effect of water depth variation with time on wave induced bottom shear stress was not considered . . . ... .. 192 7.4 Results of some sensitivity runs of the 3D nutrient model . .197 H.I Results of some sensitivity tests for Lake Okeechobee phosphorus dynamics . . . . . . . . . . . . . . ... .. 308 H.2 Model coefficients used in the 3D simulations of Lake Okeechobee phosphorus dynamics . . . . . . . . . . .... .. 309 H.3 Model coefficients used in the 1D simulations of Tampa Bay ni trogen dynamics . . . . . . . . . . . . .... .. 309 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 EFFECTS OF HYDRODYNAMICS AND SEDIMENT TRANSPORT PROCESSES ON NUTRIENT DYNAMICS IN SHALLOW LAKES AND ESTUARIES By Xinjian Chen April 1994 Chairman: Dr. Y. Peter Sheng Major Department: Coastal and Oceanographic Engineering Hydrodynamic processes (turbulent mixing, wave action, circulation, etc.) and sediment transport processes (settling, resuspension, deposition, flocculation, etc.) can significantly affect the nutrient dynamics in lakes and estuaries. Previous water quality models (e.g., WASP model) did not use sophisticated models for hydrodynam ics and sediment dynamics and therefore lack the ability of quantitatively modeling the complex nutrient dynamics in lakes and estuaries. In addition, most of the previ ous water quality models used very large grid (box) size and time step (much larger than the time and length scales of hydrodynamics and sediment transport processes) which could result in totally incorrect model predictions. Because of the dependence of nutrient dynamics on hydrodynamics and sediment dynamics, it is essential that a water quality model should be able to quantify hydro dynamic and sediment transport quantities in lakes and estuaries and that the grid size and the time step in model applications should be smaller than the length and time scales of hydrodynamics and sediment transport. This study developed 1D and 3D nutrient models and aggregated these nutrient models with models for hydrody namics and sediment transport originally developed by Sheng et al. to'study effects of hydrodynamics and sediment transport processes on nutrient cycles in shallow lakes and estuaries. The integrated 1D and 3D models for hydrodynamics, sediment transport, and nutrient dynamics were applied to Lake Okeechobee and Tampa Bay. From field data and model applications, this study obtained four major conclusions: (1) resuspensions of sediments and nutrients in Lake Okeechobee and Tampa Bay were mainly caused by windinduced waves, (2) the resuspension flux of nutrients is about 2 to 3 orders of magnitude larger than the diffusive flux, (3) the release of nutri ents from suspended sediments is affected by pH and dissolved oxygen levels in Lake Okeechobee and Tampa Bay, and (4) because the time step of a realtime simulation (1015 min.) is generally much smaller than the time scale of desorptionadsorption reactions of nutrients, equilibrium models for desorptionadsorption reactions often overpredict the release of nutrients from resuspended sediments. CHAPTER 1 INTRODUCTION 1.1 Study Background Excess loadings of nutrients into lakes and estuaries around the world has re sulted in augmented growth of various autotrophs (an autotroph is an organism which can supply its own food without requiring a specified exogenous factor for normal metabolism), of which some may in turn reach high trophic levels. Both the growth and the mortality of autotrophs (such as algae) require oxygen for respiration and for decomposition. In a stratified system, due to weak mixing process in the water column, oxygen is hardly transported downward and the dissolved oxygen (DO) level near the bottom of the water column can be very low. In this case, if the biological oxygen demand (BOD) at the bed is still high, a very stressed ecosystem (noxious condition) may be the consequence. For example, about 66 million kilograms of nitro gen and 6.5 million kilograms of phosphorus enter into Chesapeake Bay in an average year. As a result, the bay has become eutrophic with increasing anoxia, decline in submerged vegetation, severe reductions in harvests of rockfish and oysters, and con taminated sediments where marine organisms feed. In Tampa Bay, increased nutrient loading has resulted in hypoxia and decline in seagrass and fisheries, particularly in Hillsborough Bay where sediment oxygen demand becomes higher (Spaulding et al., 1989). Lewis and Estevez (1988) report that dissolved oxygen (DO) levels in Hillsbor ough Bay are below Florida state water quality criteria for DO levels for about 6090 days per year. In Lake Okeechobee, the annual external loading of phosphorus varied more than fivefold in less than three years (Schelske, 1989). In 1980, the phosphorus loading into the lake was about 0.22 million kilograms. It increased to 1.16 million 2 kilograms in 1982. The annual average of orthophosphorus was as high as 45 /ug/l, indicating that orthophosphorus supplies could be greater than could be utilized for algal growth. In order to effectively manage the loading of pollutants such as nutrients in a lake or an estuary and control the eutrophication, it is necessary to have a quanti tative understanding of the transport and transformation processes of nutrients such as nitrogen and phosphorus in the aquatic system. The nutrient dynamics in lakes and estuaries are not only determined by biological and chemical processes, but are also strongly affected by climate (wind, air temperature, etc.), hydrodynamics (wave, tide, current, turbulence mixing, etc.) and sediment transport processes (resuspen sion, deposition, flocculation, etc.). It is apparent that there exist interactions among all these processes. For example, wind can generate waves and create circulation in lakes or estuaries. Waves and currents can not only transport pollutants through ad vective/diffusive processes, but can also cause resuspension of the bottom sediments, which is strongly coupled with the food chain of organisms and nutrient dynamics. The major connections among hydrodynamics, sediment transport dynamics, and nutrient dynamics in lakes and estuaries are shown schematically in Figure 2.1. Because a physical process usually occurs more quickly than a biochemical one, hydrodynamics and sediment transport processes are often the dominant controlling factors of a marine ecosystem. It is now widely recognized (Kemp et al., 1982; Simon, 1989; Sheng et al., 1989a) that a thorough understanding of the ecology of lakes and estuaries requires a comprehensive knowledge on the interactions among biochemical and physical processes affecting the system. Numerous water quality models (e.g., Chen and Orlob, 1975; Di Toro et al., 1977; Harleman, 1977; Voinov and Akhremenkov, 1990; Ambrose et al., 1991) have been developed since the primary work of Streeter and Phelps (1925). However, most of these models lacked the coupling of hydrodynamics and sediment transport processes 3 with biochemical processes. Although some of the recent models (e.g., Lam et al., 1983; Ambrose et al., 1991) did consider physical processes such as the benthic flux of nutrients due to the sediment resuspension, they can not be used to simulate realtime nutrient variations caused by a shortterm physical process such as an episodic storm event because these models are unable to resolve such a shortterm event. Most of the previous water quality models (e.g., Chen and Orlob, 1975; Ambrose et al., 1991) have been used to simulate longterm variations of nutrient dynamics in lakes and estuaries. In these applications, because the time step (usually a few hours to one day) was larger than the time scale of hydrodynamics and sediment transport (a few minutes to a few hours), a lot of hydrodynamic processes (e.g., circulations, wave actions etc.) and sediment transport processes (e.g., deposition, resuspension, etc.) which can have significant effects on nutrient cycles have not been accurately considered. For example, the internal loading of nutrients from bottom sediments in Lake Okeechobee can not be accurately calculated in a water quality model if a time step of 6 hours is used to simulate nutrient dynamics in the lake because the resuspension event can be shorter than 6 hours and the net flux of sediments at the watersediment interface can alter direction. During one time step (6 hours), the bottom sediments can be eroded in the first 3 hours due to an increase of wind speed and the suspended sediments can be deposited to the bottom in the second 3 hours when the wind becomes calm. In this case, the internal loading of nutrients from bottom sediments calculated by a conventional water quality model may be zero because the average of net flux of bottom sediments in 6 hours may be zero. Apparently, the internal loading of nutrients from bottom sediment so calculated is not correct, because it is not different from that of no resuspension/deposition in 6 hours. In reality, even the net loading of sediments from the lake bottom is zero in 6 hours, the net loading of nutrients due to resuspension may not be zero because of the desorption of nutrients from sediments after the resuspension. 4 Problems with these previous water quality models are beyond the time step and grid size used by them. Even with small time step and grid size, most of these models still can not be used to simulate shortterm nutrient dynamics in shallow lakes and estuaries because of the following reasons: 1. Wave boundary layer dynamics and wavecurrent interactions are simulated in none of the previous water quality models. 2. It is assumed that the water column is well mixed, and a spatial and time independent eddy diffusivity is used in most previous water quality models (e.g., Chen and Orlob, 1975; Ambrose et al., 1991). 3. Because most previous water quality models were not coupled with a 3D time dependent free surface model of hydrodynamics, they could not accurately sim ulate the horizontal and vertical transport of different water quality components due to advective and diffusive motions. 4. Because most previous water quality models were not coupled with a sophisti cated model for sediment dynamics (e.g., constant settling velocities were often used in previous water quality models even for fine sediments), they could not accurately consider effects of sediment processes such as resuspension, deposi tion, flocculation, and settling on nutrient dynamics in lakes and estuaries. 5. If a small time step (e.g., 15 minutes) is used, the assumption that the desorp tion/adsorption reactions reach equilibrium instantaneous may be not correct because the time scale of desorption/adsorption reactions can be much larger than the time step used in the model. However, most previous water quality models used the instantaneous model for desorption/adsorption reactions. 1.2 Study Objectives Because nutrient concentrations in the sediment column are usually 2 to 3 or ders of magnitude higher than these in the water column (Simon, 1989; Sheng et al., 5 1989a), the resuspension of bottom sediments and the desorption/adsorption of nu trients from/onto sediments are two important processes controlling nutrient budget in the water column. The resuspension of nutrients is especially significant in shallow lakes and estuaries because it is mainly caused by windinduced waves which can be very large under severe storm conditions. However, the resuspension of nutrients during a storm event does not necessarily mean an increase of nutrient concentra tion in the water column after the storm is gone because the deposition of sediments can bring the resuspended nutrient back to the bottom. An increase of nutrient concentration occurs when the desorption process takes place during the sediment resuspension. Recently, some researchers have studied the adsorption/desorption ki netics of chemical species in the presence of solid particles (e.g., Gschwend et al., 1987; Gibbs, 1983). However, little has been done to quantify the influence of various sediment transport processes (erosion, deposition, turbulent mixing, settling) on the distribution and cycling of pollutants or nutrients, particularly in the nearbed region. Hydrodynamics, sediment dynamics, and nutrient dynamics in shallow lakes and estuaries may be different from those in deep ones. For example, because of the shallowness and the turbulent mixing induced by wind, a shallow lake is often well mixed (Harleman and Shanahan, 1983), while a deep lake can have a significant density stratification. Another important hydrodynamic characteristic in a shallow lake is the water surface setup wherein the free surface at the downwind side of the lake is superelevated by the wind force above the elevation of the upwind side. Compared to deep lakes, setup in shallow lakes is enhanced (Harleman and Shanahan, 1983). Because of the setup, a decrease in wind speed often causes a seiche motion in a shallow lake. Sheng et al. (1989a) recorded strong seiche motion in Lake Okeechobee (a shallow lake in South Florida with a mean depth of about 3 meters). Their data show that the peak value of water current was associated with seiche motion in the lake. 6 Because of the shallowness, waveinduced bottom shear stress in shallow lakes and estuaries is usually much higher than in deep ones (for the same wave height and wave period). As a result, wave action in shallow lakes and estuaries plays an important role in sediment dynamics. Therefore, it is essential that modeling sediment dynamics in shallow lakes and estuaries requires a sophisticated model for wave boundary layer dynamics and wavecurrent interaction. Biochemical transformations in nutrient cycles in shallow lakes and estuaries can also be different from deep ones. In shallow lakes and estuaries, the water column is usually wellaerated because of the shallowness and the turbulent mixing induced by wind (Somlyody, 1983). Because the depthintergal light intensity in shallow lakes and estuaries is higher than in deep ones, photosynthesis in the water column of shallow lakes and estuaries is quicker than that of deep ones. On the other hand, wave induced sediment resuspension in shallow lakes and estuaries can affect nutrient cycles by causing the reduction of light intensity and loading nutrient from the sediment column to the water column. Because of the importance of hydrodynamics and sediment transport processes to nutrient dynamics in shallow lakes and estuaries, this study intended to develop a coupled numerical model of hydrodynamics, sediment transport, and nutrient dynam ics. By using the numerical model and some field data collected in Lake Okeechobee and in the shallow region of Tampa Bay (Hillsborough Bay), this study investigated effects of hydrodynamics and sediment transport processes on nutrient dynamics in shallow lakes and estuaries. In detail, this study would like to answer the following questions: 1. What is the relationship among wind, current, and suspended sediment concen tration in Lake Okeechobee? 2. How important is horizontal transport in affecting the distributions of suspended sediment and phosphorus components in Lake Okeechobee? 3. How important are wave actions in affecting sediment resuspension process in Lake Okeechobee and Tampa Bay? 4. How lutocline, if there is, is built up in the muddy zones of Lake Okeechobee and Tampa Bay? 5. How important is sediment resuspension to nutrient dynamics in the water col umn? 6. Is the release of SRP from sediments affected by pH, temperature, and dissolved oxygen concentration in Lake Okeechobee? 7. Is the instantaneous equilibrium model for adsorption/desorption reactions ap plicable in a realtime simulation? 8. How important are adsorption/desorption kinetics in affecting nutrient dynam ics in lakes and estuaries? 9. How important is settling of particulate nutrients in affecting shortterm and longterm nutrient dynamics in lakes and estuaries? 10. How large can resuspension fluxes of nutrients be in Lake Okeechobee and Tampa Bay during episodic events? In order to accomplish these objectives, the following steps have been taken: 1. Field data on wind, waves, suspended sediment concentration, current, tide, dissolved oxygen (DO), and various nitrogen and phosphorus species of nutri ents were collected in Tampa Bay and Lake Okeechobee, 2. Field data were analyzed, 3. Numerical models of hydrodynamics and sediment transport originally devel oped by Sheng (e.g. Sheng and Chen, 1992) were used to compute currents and sediment concentrations in Lake Okeechobee and Tampa Bay. 8 4. Numerical models of nutrient dynamics were developed and used to study phos phorus dynamics in Lake Okeechobee and nitrogen dynamics in Tampa Bay. 5. Some important factors controlling nutrient concentrations in the water column (e.g. advections, wave actions, sediment resuspension, the adsorption/desorption kinetics of phosphorus and nitrogen, algal uptake of nutrients) were studied us ing the nutrient models. 6. Model results of current, suspended sediment concentration, and nutrient con centrations in Lake Okeechobee and Tampa Bay were analyzed and compared to the data. 7. Based on model results and field data on currents, suspended sediment concen trations, and nutrient concentrations, some conclusions were drawn. 1.3 Outline of Presentation Chapter 2 discusses phosphorus and nitrogen dynamics in lakes and estuaries. Ef fects of temperature, light intensity, hydrodynamics and sediment transport processes on nutrient dynamics are also discussed in this chapter. Chapter 3 discusses the general hydrodynamics of circulation, waves, and wave current interaction in shallow lakes and estuaries, followed by a discussion of sediment transport processes in these water bodies. Previous studies on the relevant subjects are reviewed. Chapter 4 describes a 3D model and a 1D model for hydrodynamics and sed iment transport. This chapter contains (1) assumptions used in the 3D and 1D models for hydrodynamics and sediment dynamics in shallow lakes and estuaries, (2) governing equations for current and suspended sediment concentration, (3) specifi cation of boundary and initial conditions for current and suspended sediment con centration in the models, (4) determination of model parameters, and (5) solution algorithms used in the models. 9 Chapter 5 describes a 3D model and a 1D model for nutrient dynamics in shallow lakes and estuaries. Similar to Chapter 4, this chapter contains (1) assumptions used in the models for phosphorus dynamics and nitrogen dynamics, (2) governing equations for phosphorus and nitrogen species, (3) specification of boundary and initial conditions for different phosphorus and nitrogen species, (4) models for kinetic pathways in phosphorus and nitrogen cycles, (5) determination of model parameters, and (6) solution algorithms. Chapter 6 presents field data collected in Lake Okeechobee. These data include a fourweek survey in the spring of 1989 and a threeday storm event survey in Febru ary, 1993. Measured wind, current, suspended sediment concentration, and nutrient concentrations are analyzed in this chapter. Chapter 7 presents model simulations of hydrodynamics, sediment dynamics, and nutrient dynamics presented in Lake Okeechobee. These model simulations include the following: 1D simulations of sediment and phosphorus dynamics at L002 (a station) during a threeday storm event in February 1993, 1D simulations of sediment dynamics at several stations during the 1989 Spring Survey, and 3D simulations of hydrodynamics, sediment transport, and phosphorus dynamics of the fourweek survey in spring of 1989. Chapter 8 presents field data collected in Tampa Bay. Measured wind, suspended sediment concentration, nutrient concentrations during two storm events in February 1992 and March 1993 are analyzed in this chapter. Chapter 9 presents applications of the 1D model of hydrodynamics, sediment transport, and nutrient dynamics to Tampa Bay. These applications include (1) a oneday storm event observed in February 1992, (2) a oneday storm event observed in March 1993, and (3) a oneweek event observed in March 1993. Chapter 10 summarizes this study. Some conclusions are drawn from the appli cations of the numerical models to Lake Okeechobee and Tampa Bay. Limitations 10 and shortcomings of the numerical models for hydrodynamics, sediment dynamics, and nutrient dynamics used or developed for this study have been discussed. Recom mendations are provided for further modifications of the models and further research on modeling hydrodynamics, sediment transport processes, and nutrient dynamics in shallow lakes and estuaries. CHAPTER 2 NUTRIENT DYNAMICS IN LAKES AND ESTUARIES 2.1 Introduction Nutrients serve as raw materials for primary production of organic matter in lakes and estuaries. The importance of nutrient dynamics in controlling the primary and secondary production in lakes and estuaries is becoming increasingly apparent. As more and more information on nutrient dynamics in lakes and estuaries is developed, our understanding of the roles of various transformation pathways continues to be refined. Excess nutrient loading into a lake or an estuary can result in increased growth of various autotrophs. Algae and other autotrophs require numerous nutri ents, including nitrogen, phosphorus, carbon, silicon, and sulfur, etc. Among these nutrients, the first three are utilized most heavily by algae. Since carbon is generally abundant in lakes and estuaries, nitrogen and phosphorus are the two major nutrients regulating the ecological balance. Nutrient concentrations in lakes and estuaries are constantly changing in time and space due to loadings from rivers, exchanges with ocean, seasonal climatic changes, biochemical transformations, as well as hydrodynamics and sediment dynamics. Fig ure 2.1 presents a conceptional picture of the major phosphorus and nitrogen path ways and the relevant hydrodynamic and sediment transport processes in lakes and estuaries. Although somewhat simplified, this figure shows how various phosphorus and nitrogen components are transformed and how phosphorus and nitrogen cycles in lakes and estuaries are connected with hydrodynamics and sediment transport processes. Wind (NH3). NH "f(NH41.d Alae SP PIP 15i ,,SON:Z, I pop '"9 9o t 9 11" ^Sl~Dr\ "Jr 9z.: t .... __(N HV e N3 ~~   "s13 __.__ Nf "N20 NO3 pop : Pathways: 1: Uptake of Algae by Zooplankton 2: Uptake of SRP. ammonium, and nitrate by Algae 3: Release of SRP and ammonium due to algal excretion and respiration 4: Release of SRP and ammonium due to zooplankton excretion and respiration 5: Release of DOP and SON during mortality of algae 6: Release of DOP and SON during mortality of zooplankton 7: Mortality of algae SMuAerobic Layer OP 10 SRP "PIP Anaerobic Lay., 8: Mortaljity of zooplankton 9: Adsorption eowption 10: Mineralization of DOP 11: Ammonlflcation 12: Nitrification 13: Denitiriflcatloni 14: Volatilization of ammonia 15: Instability of ammonium 16: Settling of algae Figure 2.1: Nutrient dynamics and the relevant hydrodynamics and sediment trans port processes in lakes and estuaries 13 2.2 Phosphorus Dynamics in Lakes and Estuaries Phosphorus enters lakes or estuaries from weathering of soil and rock and subse quent runoff, point source discharge such as sewage treatment plants, and from dairy farms and agricultural fertilizers. The magnitude of these pathways and their con stituent makeups ultimately determine phosphorus availability. However, phosphorus concentration in the water column of a lake or an estuary is not merely determined by these external sources. Internal sources such as resuspended bottom sediments containing inorganic and organic phosphorus are also very important. Phosphorus can be classified into two groups: dissolved phosphorus and particu late phosphorus. The criterion for this classification of phosphorus species is to use a 0.45 im membrane filter to separate the two fractions. Typically, dissolved (or sol uble) phosphorus in lakes and estuaries includes soluble reactive phosphorus (SRP), dissolved organic phosphorus (DOP), and dissolved inorganic phosphorus (DIP). Par ticulate phosphorus includes alga particulate phosphorus (green, bluegreen, and di atom), bacteria particulate phosphorus, zooplankton particulate phosphorus, particu late inorganic phosphorus (PIP), and particulate organic phosphorus (POP). Individ ual lakes or estuaries can have very different compositions of particulate phosphorus. For example, a lake with high external loading of organic matter can have a bac terial biomass greater than the phytoplankton biomass (Grobbelaar, 1985). Algae associated particulate phosphorus seldom dominates in lakes unless the lake system is very eutrophic or an algal bloom is occurring. Bacterial phosphorus is usually not simulated in a water quality model because once bacteria are saturated with phospho rus their excretion of internal phosphorus molecules will be equal to the phosphorus molecules they have actively taken up from the external pool (Jansson, 1988). The cycle of phosphorus in lakes and estuaries is in dynamic flux with continu ous movement among the different forms of phosphorus mentioned above. In most lakes and estuaries particulate phosphorus (PP) concentration is much higher than 14 dissolved organic phosphorus concentration, which in turn is higher than dissolved in organic phosphorus concentration (DIP). Generally, phosphorus in lakes is partitioned into about 60 to 70 percent PP, 20 to 30 percent DOP, and 5 to 10 percent DIP (Va liela, 1984). However, the relative proportions of the different forms of phosphorus vary from one aquatic system to another. The Organization for Economic Coopera tion and Development (OECD, 1982) found that the percentage of orthophosphorus ranges from 20 percent to 45 percent in many of the European and U.S. lakes. Schnoor and O'Connor (1980) found dissolved phosphorus to range from 35 to 75 percent in 81 northern U.S. lakes, and phytoplankton phosphorus to account for 10 to 40 percent of the total phosphorus. In seawater, because algal concentrations are usually lower than those in lakes, dissolved organic phosphorus concentration may be lower than dissolved inorganic phosphorus concentration. In the English Channel, most of the phosphorus present is DIP in winter and DOP in summer when algae are abundant (Valiela, 1984). As shown in Figure 2.1, biochemical transformation processes in the phosphorus cycle include (1) mineralization of organic phosphorus, (2) uptake of soluble reactive phosphorus by algae, (3) conversion of algal particulate phosphorus to zooplankton particulate phosphorus due to the uptake of algae by zooplankton, (4) excretion of SRP by algae and zooplankton, and (5) sorptiondesorption reaction of inorganic and organic phosphorus. 2.2.1 Mineralization of DOP The mineralization process of DOP is a biological decomposition process, which is mediated by bacteria. Since dissolved organic phosphorus comprises a major portion of phosphorus in many lakes and estuaries, the speed of mineralization of DOP directly influences the SRP level in lakes and estuaries and further influences the growth rate of algae. The mineralization of DOP is a relatively fast process, it can take place in a 15 few hours (compared to the mineralization of carbon and nitrogen, which takes place in a few days. See, for example, Golterman, 1973). The decomposition of dissolved organic phosphorus has been studied extensively. Thompson et al. (1954) found that the mineralization of dissolved organic phospho rus is influenced by pH. An increase in pH causes a temporary increase in the rate of mineralization of dissolved organic phosphorus. Temperature can also affect the speed of the mineralization of DOP. High temperature can stimulate the mineral ization process of DOP and thereby liberate organicbound phosphorus (Jensen and Andersen, 1992). The mineralization rate of DOP is usually modeled by a firstorder equation as follows (J0rgensen, 1983b): S= KDP2 (2.1) wherein, P2 is the DOP concentration and KD is a rate constant for the mineralization of DOP. KD is a function of pH and temperature. 2.2.2 Uptake of SRP by Algae Soluble reactive phosphorus (SRP), sometimes also called dissolved reactive phos phorus (DRP), can be taken up by algae as a nutrient. In a phosphorus limiting system, algal growth is dependent on the available SRP in the water column. Soluble reactive phosphorus is a combination of orthophosphorus and low molec ular weight organic phosphorus (Barica and Allen, 1988). Since the low molecular weight organic phosphorus can be rapidly cycled and ultimately taken up by algae, SRP is an acceptable measure of dissolved bioavailable phosphorus (Barica and Allen, 1988). Bioavailable phosphorus is commonly defined as the sum of immediately available phosphorus and that which will become available by a naturally occurring process (Bostrom et al., 1988). NaOHExtractable P, commonly assumed to correspond to aluminum and ironbound phosphorus, is the sediment form most readily available 16 to algae (Bostrom et al., 1988). HClExtractable P, corresponding to calciumbound phosphorus, and particulate organic P (POP) were found not to support algal growth. Through the uptake of SRP by algae, SRP is converted to algal particulate phos phorus as shown in Figure 2.1. The uptake rate of SRP by algae is proportional to the growth rate of algae. This can be seen from the mass conservation of SRP (excluding other processes): ap Oa(p3c) (2.2) or a49 A 9c. 19P3 t Oc3 t 3t (2.3) where P1 is the SRP concentration, p3 is phosphorus per unit biomass of algae, and c, is the algal biomass concentration. Because p3 varies only slowly within a narrow range (0.01 to 0.02), the second term on the righthand side of Equation (2.3) is negligible compared to the first term, and the uptake rate 0Pj/Ot becomes ^= P3 1 P3pA (2.4) where P3(= p3c.) is the concentration of algal particulate phosphorus and /i,' is the growth rate of algae due to the uptake of SRP (see Section 2.4). 2.2.3 Conversion of Algal P to Zooplankton P While the uptake of SRP by algae converts the dissolved inorganic phosphorus to algal particulate phosphorus, the uptake of algae by zooplankton converts the algal particulate phosphorus to zooplankton particulate phosphorus. As will be seen in Section 2.4, different zooplankton species use different mechanisms to eat algae at different speeds. The conversion rate of algal P to zooplankton P is related to the uptake rate of algae by zooplankton. Let p4 be the phosphorus per unit biomass of zooplankton, P4 the concentra tion of zooplankton particulate phosphorus, and c, the concentration of zooplankton biomass, then OP4 Oc9C 4 (2O5) t = P4 5 + C  (2.5) Because p4 is almost constant during the grazing of algae, the second term on the right hand side can be neglected. Therefore, the rate of increase of zooplankton particulate phosphorus (P4) during the grazing of algae by zooplankton becomes O P4 Oc, 4P= p =p4Cc, = P4 (2.6) wherein p, is the growth rate of zooplankton due to the grazing of algae (see Section 2.4). 2.2.4 Excretion of Algae and Zooplankton The excretion process of dissolved phosphorus by algae and zooplankton is a regeneration process of phosphorus from organic phosphorus. Although there are little or no data on the excretion of algae and zooplankton under field conditions, there are some excretion data from laboratory studies and in situ field measurements within enclosed systems. The excretion process is dependent of the feeding rate, density of algae or zooplankton, density and quality of the food, temperature, etc. (Hooper, 1972). For example, it was found (Martin, 1968) that the excretion of phosphorus by zooplankton is minimal when phytoplankton are abundant and vice versa. The reason for it is because phospholipids are being stored and/or used in egg production when food is abundant. When food is scarce, the stored lipids are used as energy source and phosphorus is then excreted. If excretion rates are constant, the SRP increase due to excretion of algae and zooplankton will be proportional to the biomass of algae and zooplankton (J0rgensen, 1983a): OP,  = K,,c. + K.,c, (2.7) wherein Ka_ and K,_ are rate constants of excretion of SRP by algae and zooplankton, respectively. 2.2.5 SorptionDesorption Reactions Sorptiondesorption processes play important roles in phosphorus dynamics in lakes and estuaries. Due to the fact that adsorbed phosphorus in sediments is gen erally two to three orders of magnitude higher than the phosphorus concentration in the water column, desorption of adsorbed phosphorus during a resuspension event can cause a significant increase of phosphorus concentration in the water column. On the other hand, sorption of phosphorus ions by sediment particles implies the removal of phosphorus from solution, an often used method in wastewater treatment. The term sorption means either the process of adsorption, or the process of ab sorption, or both. If the conversion of phosphorus from soluble phase to solid phase is restricted to the surface, it is regarded as an adsorption reaction. An absorption reaction, on the other hand is the penetration of phosphorus into the solid phase. Since it is generally very difficult to distinguish these two reactions, the less specific term sorption is frequently used. Sorptiondesorption reactions of inorganic and organic phosphorus are influenced by temperature, pH value, and the concentration of dissolved oxygen (Berkheiser et al., 1980). Sorption isotherms have been used to described the relationship between the amount of P sorbed and that remaining in water at a constant temperature in equilibrium conditions. The Freundlich and Langmuir equations (Berkheiser et al., 1980) below are frequently used to represent sorption isotherms. The Freundlich equation reads as Pad = KP, (2.8) where Pd is the phosphorus sorbed per unit sediment particles, Pdi, is the dissolved phosphorus concentration, and K and n are constants. This equation shows a linear relationship between adsorbed phosphorus (Pad) and dissolved phosphorus (Pdi,) in a loglog plot, log Pd,, = n log Pd. + log K (2.9) 19 On the other hand, the Langmuir equation shows a linear relationship between the reciprocal values of adsorbed and dissolved phosphorus: 1 a 1 1 +a (2.10) Pad a d P or P, Pd , P.ad P, Pd. (2.11) a P. + Pdj where a is a constant related to the sorption energy and P. is the maximum value of adsorbed phosphorus when Pdi, goes to positive infinity. P, is dependent on temper ature and pH value. Most water quality models assume that sorptiondesorption reactions take place so quickly that the adsorbed and the dissolved phosphorus reach equilibrium instan taneously (e.g. Ambrose et al., 1991; Dickinson and Huber, 1992). This may be a reasonable assumption when the time step of the water quality simulation is large compared to the time it takes to reach equilibrium. However, sometimes a water quality model may use a small time step during which the adsorbed and the dissolved phosphorus may not achieve equilibrium. In this case a kinetics model which describes the process rate is needed. As pointed out by Witkowski and Jaffe (1987), analyt ical chemists have repeatedly found that complete recovery of contaminants from soil/sediments frequently requires lengthy extraction periods, abrasive mixing, and strong solvents. These observations are in contradiction with the equilibrium models which assume that sorptiondesorption reactions are accomplished instantaneously. The frequently used kinetics models for sorptiondesorption reactions are the first order models (e.g., Berkheiser et al., 1980):  = D,Pd + SrPdi, (2.12) where Dr is the desorption rate of adsorbed phosphorus and S, is the sorption rate of dissolved phosphorus. 20 Because at equilibrium 9Pad/lt = 0, we have S = P'd P (2.13) D, P2 . where P'd and Pd, are the adsorbed and dissolved phosphorus, respectively, at the equilibrium condition and Pc is the partition coefficient. Therefore, Equation (2.12) becomes aPod 1 D,(P d PcPd.) (2.14) Wu and Gschwend (1988) developed a retarded intraparticlediffusion model to simulate sorptiondesorption kinetics of organic compounds. In their model, sediment particles making up of natural soils or sediments were viewed as porous spherical aggregates of finer grains in which macroscopically sorbed chemicals occur micro scopically either sorbed to the natural organic matter on the solid matrices or freely dissolved in the intraparticle pore water. The diffusive penetration of hydrophobic compounds into or out of such particles is retarded by microscopic scale linear re versible sorption exchange between intraparticle pore water and adjacent solids. Wu and Geschwend's model seems to have a more solid theoretical background than the widely used firstorder models. However, Wu and Geschwend's model for sorption desorption is complicated to use in a water quality model. Furthermore, their model has no validation at all. They compared their model results to the results of the first order model. If the size distribution is not very wide, the retarded intraparti cle diffusion model gives almost the same results as the firstorder model (Wu and Geschwend, 1988). 2.3 Nitrogen Dynamics in Lakes and Estuaries Nitrogen enters into lakes or estuaries from point and diffuse sources on the land, atmospheric diffusion, upwelling deep ocean, and biological fixation. The gaseous form of elemental nitrogen N2 comprises about 78% of the atmosphere. Nitrogen is the major part of the atmospheric air. Although the gaseous form of nitrogen is 21 relatively insoluble in water, it dissolves sufficiently so that nitrogen deficiencies would never limit photosynthesis if N2 were available for photoplankton assimilation. Similar to phosphorus, nitrogen components in lakes or estuaries also include ad sorbed organic nitrogen, soluble organic nitrogen (SON), ammonium nitrogen NH+, adsorbed ammonium nitrogen, ammonia nitrogen NH3, nitrate NO;, and nitrite NO2. As shown in Figure 2.1, kinetic pathways in the nitrogen cycle in lakes and estuaries are also shown, together with the relevant processes of hydrodynamics and sediment dynamics. This figure indicates that nitrogen in lakes and estuaries can undergo the following biological, chemical, and physical processes: 1. mineralization of organic nitrogen, i.e., ammonification of soluble organic nitro gen to ammoniumnitrogen, 2. nitrification of ammoniumnitrogen in the presence of oxygen; in this process, ammonium is first converted to nitrite and then converted to nitrate, 3. volatilization of ammonianitrogen, 4. denitrification of nitrate to gaseous nitrogen N2 in the anaerobic sediment layer, 5. fixation of nitrogen N2 by aquatic phytoplankton, 6. uptake of nitrate and ammonia by algae, and 7. conversion of algal particulate nitrogen to zooplankton particulate nitrogen due to grazing of algae by zooplankton. Like the phosphorus cycle, the nitrogen cycle is influenced by the hydrodynamics and sediment dynamics in lakes and estuaries. For example, the resuspension of sed iments from the bottom and the desorption process afterwards provide a significant pathway for nitrogen nutrient to enter into the water column, while the adsorption process and the deposition of sediment particles are major sinks for nitrogen in the water column. The vertical turbulent mixing not only transports the resuspended nitrogen from the near bottom layer to the top layer of the water column, but also 22 influences the desorption/adsorption reactions of nutrients from/onto suspended sedi ments because of the flocculation and breaking of sediments caused by the turbulence. The following is a brief review of nitrogen dynamics in lakes and estuaries. Because of the similarity to those in the phosphorus cycle, nitrogen kinetic pathways associated with algae and zooplankton are omitted in the discussion. 2.3.1 Ammonification Ammonification is the biological process of formation of ammonium from soluble organic nitrogen. It is the first step of nitrogen mineralization, in which organic nitrogen is converted to the more mobile, inorganic state. The mineralization process is a pathway of regenerating the nutrient in a form usable by plants. Ammonification is actually the decomposition process of soluble and particulate nitrogen compounds of dead organisms and those excreted by plants and animals. This decomposition process is brought about by the various species of proteolytic bacteria. Generally, the decomposition is quite efficient, and at no time does any sub stantial concentration of free or combined dissolved amino acids build up. However, some particularly stable organic nitrogen compounds resist attack by bacteria and eventually sink to the bottom where they are incorporated into the sediments. The rate of ammonification is often expressed as a firstorder reaction (Rao et al., 1984): aN, S= K4N, (2.15) wherein, K4 is the rate constant of ammonification which is a function of water temperature, pH, and C/N ratio of the residue (Reddy and Patrick, 1984). 2.3.2 Nitrification The second step of mineralization of organic nitrogen is nitrification, in which nitrite is formed by the oxidation of ammonium and is further oxidized to nitrate: NHs + 3/202 * HNO2 + H1120 (2.16) 23 HN02 + 1/202 * HN03 (2.17) Nitrification requires oxygen as the electron acceptor and thus can only occur in an aerobic condition, i.e. in the water column and in the aerobic layer of the sediment column. This pathway is typically associated with the metabolism of certain bacteria. Two groups of bacteria are distinguished: one derives its energy for cell synthesis by the oxidation of ammonium, the other by the oxidation of nitrite. The dominant species of the former group is Nitrosomonas, while the latter group is Nitrobacter. In most cases, species of Nitrosomonas and Nitrobacter are found together; otherwise nitrite might accumulate to phototoxic levels. Several models for the nitrification have been proposed: (1) zeroorder equation, (2) firstorder equation, and (3) Monod equation of population dynamics. Some researchers (e.g. Srinath et al., 1976) found that the nitrification rate is independent of the substrate concentration N,: ON. K3 (2.18) at where Ni is the concentration of soluble NH+N and K3 is the nitrification rate coefficient. Others (e.g. Broadbent and Clark, 1965) showed that the nitrification rate is indeed a linear function of the substrate concentration:  = K3N,, (2.19) 0t Because the rate of nitrification is determined by the growth rate of nitrifying bacteria, ON,  = Pn (2.20) where a. is a constant and p, is the growth rate of nitrifying bacteria, Stratton and McCarty (1969) and Garland (1978) used the Monod expression for the rate of nitrifying bacterial growth in natural streams such as the Clinton River and Trent River, N,,i t" =/ ....Hn + N,"' (2.21) 24 where u.,. is the maximum growth rate of nitrifying bacteria and H., is the half saturation constant for the bacterial growth. The nitrification rate is therefore ONT i N. (2.22) '9t H. + N," The Monod equation gives a continuous transition between zeroorder kinetics and the firstorder kinetics based on the substrate concentration of ammonium. For N, much lower than H,, it approaches the firstorder equation; and for N, much higher than H,2, it approaches the zeroorder equation. Since the H, value is typically smaller than N, in natural aquatic systems, zeroorder kinetics were often used in the previous studies (e.g. Hall and Murphy, 1980; Watanabe et at., 1980). 2.3.3 Denitrification The microbial reduction of nitrite and nitrate with the liberation of molecular nitrogen and nitrous oxide is called denitrification. Denitrification is essentially a respiratory mechanism in which nitrate replaces molecular oxygen, i.e., nitrate respi ration. The end product of denitrification is the gaseous form of elemental nitrogen N2, which can escape to the atmosphere. Therefore, this process is a sink for nitrate. However, denitrification is not the sole sink of nitrate. Another sink is the utilization of nitrate as a nutrient source of plants (this may be termed nitrate assimilation). Both transformations involve reductive pathways, but the end products of nitrate respiration are volatilized while the products of nitrate assimilation are incorporated into the cell material. Denitrification process in lakes and estuaries is associated with denitrifying bac teria such as Pseudomonas denitrificans, which use NO0 as an electron acceptor to oxidize organic matter anaerobically, releasing N2 gas by the following reaction (Day et al., 1989): 5C6H1206 + 24HN03 30C02 + 42H20 + 12N2 (2.23) 25 In fact, this reaction consists two steps. In the first step of denitrification, nitrous oxide NO20 is produced from NO by denitrifers: CsH1206 + 6HN03  6C02 + 911H20 + 3N20 (2.24) Nitrous oxide, which is produced as an intermediate product in denitrification, is unavailable for phytoplankton assimilation. Anthropogenic acceleration of N20 production could result in a general heating of earth's atmosphere because N20 reacts with and breaks down atmospheric ozone 03 which normally helps maintain the earth's heat balance (McElroy et al., 1978). While nitrification occurs in the water column and the aerobic layer of the sedi ment column, denitrification occurs only in the anaerobic layer of the sediment col umn. The nitrate used in the denitrification process comes from the water column and the aerobic layer by diffusion between the aerobic and anaerobic layers. Seitzinger (1988) reviewed various studies on the denitrification process in freshwater and coastal marine ecosystems. He found that the major source of nitrate for denitrification in most river, lake, and marine sediments comes from the nitrification process in the sediments, not the diffusion of nitrate from overlying water into the sediments. Den itrification accounts for a major portion of loss of mineralized nitrogen during the mineralization of SON in sediments. In rivers and lakes, 76100% of nitrogen flux across the sedimentwater interface is N2, but only 570% in estuarine and coastal marine sediments because of the presence of salinity (Heath, 1992). Denitrification can be measured directly as total N2 production from degassed sediments. However, because N2 is very abundant in the atmosphere, long incubations (about one week) are needed to obtain significant N2 concentration changes. The process can also be measured indirectly as N20 production in the presence of acetylene which inhibits the reduction of N20 to N2. The denitrification process can be described by a zerothorder rate equation, a firstorder rate equation, or the standard MichaelisMenten equation. Dawson and 26 Murphy (1972) found that the denitrification is of zerothorder aN3 K2 (2.25) wherein N3 is the concentration of nitrate and K2 is the denitrification rate constant. Stanford et al. (1975) found the denitrification to be firstorder at low nitrate nitrogen concentration. Reddy et al. (1978) also showed that under carbonlimiting conditions the denitrification followed the firstorder kinetics: aN3  = K2N3 (2.26) Bowman and Focht (1974) showed that the denitrification followed the Michaelis Menten kinetics: _N_ N3 __ K2 H + N3 (2.27) where K2_., is the maximum denitrification rate and H,3 is the halfsaturation con stant for denitrification. 2.3.4 Nitrogen Fixation Nitrogen fixation is a biological process by which organisms reduce N2 gas to ammonia. This process is brought about by freeliving bacteria or bluegreen algae. It is an energyintensive anaerobic reaction: N2 + 3H2 2NH3 (2.28) Quantitative estimation of the flux rate of nitrogen fixation in the nitrogen cycling is quite uncertain, because this process is governed by a number of physical and chemical factors. For example, if ammonium is abundant in the system fixation can be inhibited because bacteria and algae use the nitrogen salt rather than N2. pH value can affect nitrogen fixation because the abundance of certain organisms such as Azotobacter is characteristically sensitive to pH value. Temperature also has a profound influence on N2 fixation. There is little activity at low temperature and warming promotes the microbial uptake of N2. 27 In lakes and other fresh water systems, N2 fixation may represent a major com ponent in the overall nitrogen budget and this pathway commonly contributes 30 to 80 percent of the total annual nitrogen income for an entire water body (Homrne, 1978). However, in estuarine waters, N2 fixation tends to be relatively unimportant for overall nitrogen budgets except in specialized environments. The reason why the N2 fixation in estuaries is commonly less important than in lakes is uncertain. One possible reason may be because of the relative high ammonium concentration. However, coastal waters containing low ammonium concentration still have little N2 fixation. There has been some suggestion that turbulent mixing may cause sufficient stress to destroy the structure of delicate heterocysts (Fogg, 1978). Another possibil ity suggested by Howarth and Cole (1985) is that the high sulfate concentration in seawater compared to freshwater may inhibit the ability of nitrogenfixers to obtain sufficient quantities of molybdenum needed for synthesis of nitrogenase. 2.3.5 Nitrogen Volatilization The dissolved form of ammonium in water is generally not stable and can exist in its gaseous form or ammonia (NH3). Because the ammonia concentration in the atmosphere is very low, ammonia in the water column can escape to the air. This is the volatilization process of ammonia. The volatilization of ammonia is a sink for nitrogen in an aquatic system. This process is dependent on the pH value, as can be seen from the following reaction: NH+ * NH3 + H+ (2.29) [NH3] K (230 [NH+] = [H+] (2.30) or log1o( [ ) =pIH pK, (2.31) where K. is the dissociation constant for ammonium ion and pK, = loglo K. is 9.3. 28 It can be seen from Equation (2.31) that ammonia only consists a small portion of the total ammonia and ammonium at neutral pH and that ammonia concentration equals ammonium concentration at pH = pK.=9.3. During the summer months, a pH value higher than 9.0 can be found in some very eutrophic lakes. For example, a pH of 11.0 has been reported in the shallow Sollerod Lake of Denmark (J0rgensen, 1989). In this case, the volatilization of ammonia can be very significant. A first order rate equation can be used to describe the kinetics of the ammonia volatilization process: VN4 = vn(hN4 N4) (2.32) where N4 is the concentration of ammonia nitrogen at the top layer of the water column, N4 is the ammonia concentration in the air, v., is the rate constant of volatilization, and h, is Henry's constant. Equation (2.32) can be derived from the socalled twofilm model (J0rgensen, 1983c), which is often used in chemical engineering. The basic idea of a twofilm model assumes that a stagnant liquid film of thickness dL and a stagnant gas film of thickness dG exist at the waterair interface (Figure 2.2). The transport of a volatile component through these films is due to the diffusion Vw = K.(N4 N4) (2.33) V9 = (N N) (2.34) wherein V, and V. are the rates at which ammonia is transported across the water film and the air film, respectively, KI, and K. are transfer coefficients of water film and air film, respectively, T is the absolute air temperature, and R is the gas constant. V, and V. should be equal because of the mass conservation across the two films. By using Henry's law, N4i = h ,N4 (2.35) z N: Air AirWater N dG Interface .  N; ..   Water N4 Figure 2.2: The twofilm model for the volatilization. we get VN, =V RTK.= ,,K (hN4' N) (2.36) or Vn K.K, (2.37) RTK, + hK( in Equation (2.32). 2.4 Algal and Zooplankton Growth 2.4.1 Algal Growth For an aquatic system with steady nutrient, temperature, and light intensity, the growth of algae is exponential: one single cell gives two, and two cells give four cells, etc. In other words, the rate of growth is proportional to the number of cells at any time: dM dM =..M (2.38) where M is the number of algal cells and y, is the growth rate constant, which is a function of temperature, light intensity, and nutrient concentration. In a phosphorus limiting aquatic system, the growth rate constant j, of algae may be described by an 30 empirical equation as e.g. of Monod (1949) P1 S= '.... H (2.39) p+ P, where Hp is the halfsaturation concentration of SRP for algal growth and p.... is the maximum growth rate of algae. /ta,. is dependent on temperature and light intensity. An example of Equation (2.39) is shown in Figure 2.3 (Golterman, 1975). If nitrogen is also limiting in the system, Equation (2.39) becomes P5 N3+N4 P A N3 + N4 ,(2.40) = H, + P, H + N3 + N4 where H, is the halfsaturation concentration of nitrogen for algal growth. The growth rate of freshwater algae (e.g., Cyanobacteria) can be affected by salin ity toxicity. Cerco and Cole (1992) used an empirical equation to represent the effect of salinity on freshwater Cyanobacteria growth: f(S) (2.41) S? + S2 where S is salinity (ppt) and St is the salinity at which Microcystis growth is halved. Some researchers (e.g., Jorgensen, 1976) argued that the growth rate of algae is more likely dependent on the internal nutrient content of the cell as defined by the cell quota concept (Equation (2.42)). According to the cell quota model, external nutrients is taken up by the cell and stored; subsequent cell growth is dependent on the internal nutrient concentrations in the cell. The advantage of this model over Monod kinetics is that the uptake of nutrients and cell growth are decoupled. Here, = q (2.42) where Qe is the cell quota and q, is the minimum cell quota, or the subsistence cell quota at zero growth rate. The external nutrients and cell quota formulations are equivalent if the internal cell concentration is assumed to be in dynamic equilibrium with the external concen tration (Di Toro, 1980). Di Toro (1980) estimated that the uptake of internal cell o0 100  0 80 60  Calculated 40 *Y 40 Experimental 20 0 100 200 300 400 500 600 700 800 900 1000 ug P/I Figure 2.3: Relative growth rate of algae in a phosphorus limiting aquatic system. storage occurs in 0.125 0.5 of the time required for the uptake of external nutrients. The internal cell concentrations are thus in dynamic equilibrium with the external concentration. Equation (2.42) is difficult to use directly since the cell quotas and minimum cell quotas are species specific. However, this difficulty can be avoided by expressing p. as a function of internal phosphorus, nitrogen, and chlorophyll a concentrations (Riley and Stefan, 1988): Pc C N, CaN (2.43)W S= .... P  N (2.43) where Pc and N. are internal concentrations of phosphorus and nitrogen, respectively; Pi. and Nm,, are the minimum ratios of phosphorus to chlorophyll a and nitrogen to chlorophyll a, respectively, for algal growth; and Ca is the concentration of chlorophyll a. Algal biomass is typically represented by chlorophyll a, which can be measured fluorometrically. The ratio of chlorophyll a to carbon in algal cell is about 0.01 to 0.02, which is almost the same as the ratio of phosphorus to carbon in algal cell 32 (OECD, 1982). The final concentration of algal biomass co is dependent on the algal growth rate u, the respiration rate r., the mortality of algae K_,, the grazing rate g, by zooplankton, and the settling rate of algae wa: = ,c, r.c. g.c g+ 0 (2.44) 2.4.2 Growth of Zooplankton Zooplankton play an important role in the recycling of phosphorus and nitrogen in lakes and estuaries by uptaking algae and excreting phosphorus and nitrogen after ingesting algae. Zooplankton grazing of algae enhances the cycling speed of phos phorus and nitrogen. Devol (1979) found that zooplankton can be a major supply of SRP to algae during periods of low SRP concentrations. Additionally, the selective zooplankton grazing of green algal cells shifts the algal population towards bluegreen algal domination in the summer. The growth of zooplankton in lakes and estuaries is determined by its present abundance, its grazing rate, and the efficiency with which it converts algal carbon into zooplankton carbon. The eating mechanism of zooplankton differs greatly among zooplankton species nonselectivee filters, selective filters, and raptors). Nonselective filters eat algae by grazing with a constant rate, regardless of algae concentration. Se lective filters eat algae by grazing, but have developed an ability to vary their filtering rate. The filtering rate increases as the concentration of algal biomass decreases, but approaches a fixed filtering rate which is determined by the physical characteristics of selective filters. If the algal concentration is high, selective filters can alter their particlesize selectivity and lower their filtering rate, and thus make more efficient use of algae and conserve energy. A formulation for the filtering rate us has been proposed by Canal et al. (1976): S,,c + Ah (2.45) S = #So. c+a (2.45) 33 where us_,. is the maximum filtering rate of selective filters (dependent on tempera ture), Sm is the minimum filtering rate multiplier, Ah is the algal concentration when the multiplier is (1 + S.,)/2, and c. is the concentration of algal biomass. The growth rate of selective filters is then OCs t = isCs (2.46) where Cs is the concentration of selective filters. Raptorial species eat algae by selecting, snatching, and devouring. The eating rate of raptorial species increase with the concentration of algal biomass but reaches a saturation level. The mathematical formulation for the eating rate UR of raptors can be expressed in a Monod form (Canal et al., 1976): AR = IR.. Ca(2.47) R a + Ha where Ha is the halfsaturation algal level for raptors and R., is the maximum snatching rate of raptors. uR_., is a function of temperature. 2.4.3 Algal NonPredatory Mortality Algal mortality involves the conversion of algal phosphorus and nitrogen to bio logically available nutrients. This can be conceived as a two step process (DePinto et al., 1986). First, algal death and lysis of the cell membrane releases any stored excess inorganic phosphorus and nitrogen. DePinto et al. (1986) found that cellu lar phosphorus exceeding the minimum cell quota was released in an available form rapidly after cell death and lysis. Lysed material is usually 2050 percent of the cell's organic content (Cole, 1982), and this excess inorganic phosphorus is immediately available for algal uptake. Secondly, bacteria mineralize the remainder of the cell organic phosphorus to available phosphorus. The nonpredatory mortality of algae can be modeled by a firstorder equation (Gargas, 1976): OCa  = a ,C, (2.48) where K., is the mortality rate of algae and is temperaturedependent. 2.4.4 Algal Settling Nutrient loss to sediments is a significant fraction of the total incoming nutrient load from rivers or oceans in lakes and estuaries. The loss to sediments is via two main mechanisms: (1) algal sedimentation and (2) nonalgal nutrient sedimentation. Some of the nutrients are regenerated and later mixed with the water column due to biologic movement and sediment resuspension. Different algal classes may have different settling velocities. Takamura and Ya suno (1988) measured algal sinking rates in a shallow hypereutrophic lake and found bluegreens to settle much slower than green and diatom algae. The sinking rate of senescent algae was much higher than algae in exponential growth. Diatoms sank as living cells and bluegreens sank as detritus. The cell structure of the bluegreen algae Microcystis collapsed rapidly after losing its function of carbon fixation (Takamura and Yasuno, 1988). Since algal radius is generally in the order of ym's, algae can undergo floccu lation, just like normal sediment particles (Section 3.4). Alldredge and Gotschalk (1988) reported the formation of flocs, known as marine snow or aggregates, after di atom blooms in marine water. Riebesell's (1991) studies in the North Sea show that flocculation is an important mechanism for determining the fate of organic matter produced after the algal bloom. Recently, Jackson and Lochmann (1992) studied the effect of coagulation on nutrient and light limitation of an algal bloom. Their major conclusions are (1) loss of algal cells to coagulating particles can occur when algal cells are growing at a fairly constant rate, placing a cap on the concentrations that algae can achieve and (2) the vertical flux of algae from the top layer is enhanced when flocculation occurs. 35 In contrast to the settling, algae can also have buoyancy. The mechanisms an algal cell uses to control buoyancy are very complicated. Bluegreen algae use three mechanisms to control buoyancy: (1) cell pressure increases at high light intensities causing the collapse of cell gas vacuoles and decreased buoyancy; (2) algae alter the rate of gas vacuole synthesis in response to the cell C/N internal ratio (a low C/N ratio increases gas vacuole synthesis); and (3) the cell may accumulate highdensity storage products that act as ballast and affect cell buoyancy (Spencer and King, 1989). Some bluegreen algae that use one or more of these mechanisms are Anabaena sp., Aphanizomenon flosaque, Microcystis, Oscillatoria, and Nostoc muscorum. 2.5 Effects of Temperature and Light Intensity on Nutrient Dynamics 2.5.1 Temperature In the nutrient cycles shown in Figure 2.1, almost all of the transformation pro cesses are affected by temperature. Generally speaking, temperature will accelerate the processes. For example, the desorption rate of inorganic phosphorus from sedi ment particles increases as temperature increases. However, exception also exist. For example, many investigations have confirmed that the nitrification rate below 5C and above 40C is very slow. The maximum rate of nitrification lies between 30 and 35C. The reason for this is presumably because of physiological differences in dominant bacterial strains (Alexander, 1965). The effect of temperature on reaction rates can be explained by the Van't Hoff Arrhenius equation as follows: d(_nK) AH (2.49) dT RT2 (4) where K is the 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 the universal gas constant. 36 Integrating Equation (2.49) from temperature T1 to T2 gives: K2 __H K = exp[RlT (T2 Ti)] (2.50) Since exp(AT ) is almost constant at the temperature range of interest (0 30 C), we have (Chen and Orlob, 1975) K(T) = K20OT20 (2.51) wherein K(T) is the rate constant at temperature T, K20 is the rate constant at 20C, and 0 is a unitless temperature coefficient ranging from 1.0 to 1.1. For some biological processes, Equation (2.51) can not appropriately describe the temperaturedependency because the growth of organisms may stop below or above certain temperature but reaches maximum at an optimum temperature. This phenomenon has been studied by some researchers. For example, Lassiter and Kearns (1974) proposed a formula as follows: K(T) = Kopiexp[a(T TT,)]( T T T (2.52) T*max T'opt wherein Kopt is the rate constant at optimum temperature, Tt is the optimum tem perature, T,., is the maximum temperature at which the rate constant is zero, and a is a coefficient with the unit C1. 2.5.2 Light Intensity Light Intensity affects the photosynthesis process and thus the algal growth rate. As will be discussed later, one of the effects of sediment resuspension on nutrient cycling is that resuspended sediment will reduce the light penetration and weaken the photosynthesis process in the water column. The effects of light intensity on nutrient cycling is often modeled by a light inten sity limiting function as follows (Steele, 1974): 1 1_1 f2(I) = Te *T (2.53) where I is the light intensity and I, is the optimum light intensity for algal growth. According to BeerLambert's law light adsorbed as it passes through a solute or solution is proportional to the light intensity and the optical path length: AI = eAzl (2.54) where the extinction coefficient e is a function of suspended sediment concentration and algal concentration in the water column: S= f(c,c0) (2.55) If e is constant over the water depth, then the light intensity over the water depth is I(z) = Ie (2.56) where I(z) is the light intensity at depth z and I. is the light intensity at the water surface. 2.6 Effects of Hydrodynamics and Sediment Processes on Nutrient Dynamics Water and sediment particles are carriers of nutrient species in lakes and estuaries. Thus the transport and transformation processes of nutrients as discussed in Section 2.2 and 2.3 are thus closely linked with the dynamics of the water movement and the sediment transport processes. For example, as shown in Figure 2.1, nitrogen cycling is connected to hydrodynamics and sediment dynamics as follows: Organic nitrogen contained in the soil particle is brought from the soil surface to the water column due to the excess of bottom shear stress over the critical shear stress of the bed. In the water column, part of the organic nitrogen on particles can change into dissolved form due to desorption because the resuspended sediment particles 38 contain a large amount of organic nitrogen. This process will not stop until the concentration of organic nitrogen on sediment particles and the soluble organic nitrogen (SON) in water reaches equilibrium. Ammonium NH+ can be adsorbed onto the clay particle due to electrical at traction and be nitrified to nitrate NO0 due to the presence of 02 in the water column. In the water column, nitrogen species are mixed due to the turbulence mixing process. As a result, the vertical distributions of dissolved nitrogen components become more uniform. If adsorbed nitrogen concentrations are lower than their equilibrium values, sed iment particles will adsorb nitrogen species in dissolved form. Through the settling of sediment particles, nitrogen will fall back to the bottom. In this case sediment particles act like a sink. Another important effect of hydrodynamics and sediment transport on nutrient budget in lakes and estuaries is the horizontal transport due to advection and diffu sion. Adsorbed and dissolved forms of nutrient species are transported by currents from one part of a lake or an estuary to another part. During this process, dissolved nutrients can generally follow the the current circulation, while nutrients in adsorbed (particulate) form can be transported and deposited to the bottom with sediment. The resuspension of bottom sediments is an important mechanism controlling the nutrient regeneration in lakes and estuaries. It is often found that the nutrient con centration in the bottom sediments is one or two orders of magnitude higher than that in the water column (Simon, 1988 and 1989). McClelland (1984) produced a nutrient box model for Tampa Bay and found that all point and nonpoint sources accounted for only onethird of the nitrogen released from the sediments. Adding inputs from the point and nonpoint sources to the diffusive benthic flux still only accounted for 39 about 50% of the nitrogen needed to support the observed primary production. Mc Clelland suggested that other sources supply the nitrogen via sediment resuspension and insitu water column regeneration. Johansson and Squires (1989) computed a preliminary nutrient budget for Tampa Bay and concluded that internal loading of nitrogen associated with sediment resuspension events can be quite significant. In an episodic storm event, large amount of nutrients can be brought into the water column with sediment due to large current or wave induced bottom shear stress. In the water column, because of the relatively higher redox potential, stronger mixing, and sufficient light, some transformation (e.g., nitrification) processes of nutrients can occur faster than those in the bed. Therefore, the resuspension of sediment from bottom not only increases the nutrient concentration in the water column, but may also accelerate the nutrient cycling. Because the time scale of sediment dynamics is much smaller than that of the transformation processes, the effect of resuspension of bottom sediment on nutrient cycling is very significant. On the other hand, since the food chain in lakes and estuaries is built by the constant nutrient cycling between water column and benthic sediment, it is very important to clarify the role of benthic sediment fluxes in the cycling of nutrients. It is problematic to quantify the benthic flux due to the competing influences of molecular diffusion, resuspension, and ground water seepage, and difficulties in quan tifying each one of them. Some past water quality modeling studies (e.g., HydroQual, 1987) treated the net benthic flux as an adjustable parameter in the model, which was adjusted to achieve a reasonable fit between model results and limited data. Thus, the conventional empirical treatment does not lead to a robust estimate of the ben thic flux and introduces major uncertainty in the use of the water quality model as a predictive and management tool. A definitive model of benthic flux is needed to significantly reduce the uncertainty of water quality models. Harleman (1977) sim ply ignored the exchange of nitrogen between sediment and water. Lam and Jaquet 40 (1976) proposed an empirical formula for computing the upward flux of phosphorus due to the windwave induced sediment resuspension in Lake Erie, J = kp, P' w, (2.57) p. pPw where k is a constant, p., and p, are water density and dry sediment density, respec tively, and w, is the excess wave orbital velocity. Beck and Somlyody (1982) assumed a simple relationship between wind speed W and we, and proposed J = kipw P. W_ (2.58) P. p for Lake Balaton, where k, and m are constants. Lam et al. (1983) used Equation (2.57) in a simplified Lake Erie phosphorus model for the Western Basin. There have been a limited number of studies of nutrient cycling between ben thic sediments and water column in estuaries (Simon, 1988 and 1989; Kemp et al., 1982) and lakes (Sheng et al., 1989a; Sheng, 1993; Reddy and Graetz,1989). Simon's studies strongly suggested that during periods of sediment resuspension, desorption of ammonium from sediment solids can be the major pathway for enriching the water column in a shallow estuary. However, Simon's resuspension experiment was con ducted in laboratory rather than in the field. The field studies of Sheng et al. (1989a and 1989b). and Reddy and Graetz (1989) confirmed that particulate phosphorus in water column increased significantly during sediment resuspension events in a shallow lake. Recent studies by Simon (1989) and Sheng (1993) show that the resuspension flux of nutrient in shallow water is typically 2 to 3 orders of magnitude larger than the diffusive flux from the bottom. Simon's study in the Potomac estuary found that particulate forms of nutrient deposited onto the bottom sediments in the upper reaches became resuspended and released into the water column a season later. In addition to generating benthic flux of nutrients, sediment resuspension can also affect algal populations by adsorbing the incoming light from the water surface and reducing light penetration through the water column. High sediment resuspen 41 sion may result in dramatic reduction of light intensity and thus affect the rate of photosynthesis. Because hydrodynamics and sediment dynamics can directly contribute to nutri ent transport and nutrient regeneration in lakes and estuaries, the study of nutrient dynamics should be coupled with the study of hydrodynamics and sediment dynam ics. The U.S. E.P.A. developed a water quality model called WASP5 (Ambrose et al., 1991), which assumed that lakes or estuaries consist of a few boxes and solved the nutrient equations for each box. However, the coupling of nutrient dynamics with hydrodynamics and sediment dynamics in this model is rather poor because a few box can not give a good resolution to hydrodynamics and sediment dynamics in lakes and estuaries. In addition, WASP5 generally uses a large time step (a few hours to a day), which can be larger than the time scale of hydrodynamics and sediment transport processes. Because of the poor resolution of hydrodynamics and sediment dynam ics, the prediction ability of WASP5 model is questionable. Harleman (1977) used a onedimensional model to simulate ammonia concentrations in an idealized estuary. Two cases were studied: one with the realtime tidal current, and the other with throughflow velocity due to freshwater inflow. He found that ammonia concentra tions predicted in the latter case (throughflow) differ by as much as 100 percent from the realtime tidal averages even though the same biochemical model is used. This indicates that the consideration of shortterm hydrodynamics and sediment transport processes is very important in a water quality model. To study phosphorus dynamics in Lake Okeechobee with realtime hydrodynamic and sediment transport processes, Dickinson et al. (1992) developed a 3D phosphorus model, which coupled with a 3D model of hydrodynamics and sediment dynamics developed by Sheng et al. (1991). The time step used in the 3D phosphorus model was the same as that in the 3D model of hydrodynamics and sediment transport processes. The transformation processes considered in their model were similar to 42 the WASP5 model except for that Dickinson et al. differentiated algae into three groups (green, bluegreen, and diatom). However, since no data exists for each of the algal groups thus differentiation has not been validated yet. Same as WASP5, Dickinson et al. assumed instantaneous equilibrium between dissolved and adsorbed nutrients, which is not necessary true in many cases (Witkowski and Jaffe, 1987). 2.7 Conclusions It is very difficult to study the transport/transformation processes of nutrients in lakes or estuaries for the following reasons: The measurement and analysis of nutrient species generally takes a very long time and sometimes can only be done in laboratories. To date, there is still no method available to measure the instantaneous variations of different nutrient concentrations in lakes or estuaries. As a result, the measured data may lack dynamics. Nutrient cycles in lakes or estuaries are influenced by hydrodynamics and sedi ment transport processes. Hydrodynamics and sediment dynamics in lakes and estuaries are very complicated. Transformation processes are affected by algal activities, aquatic plants, mi crobial communities and other biochemical factors. Quantitative relationships between many of these factors and the nutrient cycles are still unknown. Nutrient cycles are also affected by geological factors such as soil texture and cation exchange capacity (CEC) of soil. For example, CEC affects not only par ticle aggregation, but also the kinetics of NH+ and the volatilization process of NH3. Nutrient transformations in lakes and estuaries are generally mediated by dif ferent bacteria. The activities of these bacteria are often determined by the partial pressures of other species such as oxygen and carbon. 43 Meteorological factors (wind, humidity, temperature, sunshine etc.) can also alter nitrogen cycling either directly or indirectly. There is not enough data to quantify these complicated effects. Because of the these difficulties, numerical modeling for nutrient dynamics in lakes and estuaries is still in its infancy, although a lot of models have been developed (see, for example, the book edited by Jorgensen and Gromiec in 1989). None of the existing numerical models for nutrient dynamics has the real ability of predicting an event, particularly a shortterm event. Some models may fit field data well, but are generally sitespecific because of the assumptions contained in these models may not be correct for other aquatic systems. Moreover, because of the lack of complete understanding of all the processes involved in nutrient cycles, many model parameters have to be tuned to fit field data. CHAPTER 3 HYDRODYNAMICS AND SEDIMENT DYNAMICS IN LAKES AND ESTUARIES As mentioned previously, biological processes in an ecosystem such as a lake or an estuary are strongly affected by hydrodynamics and sediment transport processes. For example, the amount of nutrient resuspended during an episodic event could be as high as the total amount of dissolved nutrient diffused from the bottom due to molec ular diffusion over a period of a few months or even a few years (Simon, 1989). The resuspension of sediments was found to be a dominant factor controlling the nutrient budget in an estuary or a lake. Therefore, as the first step to study nutrient cy cling in an estuary or a lake, hydrodynamics and sediment transport processes should be investigated and understood. Important hydrodynamics and sediment transport processes affecting nutrient dynamics in lakes and estuaries include circulation, wave action, wavecurrent interaction, erosion, deposition, settling, flocculation, and con solidation. This chapter reviews all these processes. 3.1 Circulation Circulation is the physical process of water movement in a lake or an estuary. Circulation in a lake is generally winddriven. In a large lake, differential heating over the lake can sometimes result in horizontal densitydriven circulation. Circulation in an estuary could be driven by tide, wind, and density gradients. Estuarine circulation induced by the horizontal density gradient between fresh water and salt water is often referred to as gravitational circulation (Hansen and Rattray, 1966). Fresh water runoff, which is less dense than salt water in the ocean side of the estuary, tends to stay at the surface layer of the estuary and flows towards the ocean. The salt water remains at the bottom layer of the estuary and has the tendency 45 of flowing upstream because of the horizontal density gradient. The stratification in an estuary can be described by Harleman's estuary number which is defined as: PFr2/RT, where P = tidal prism = amount of ocean and river water entered into the estuary in one tidal cycle, Fr2 is the Froude number defined as U R is the river discharge rate, T is the tidal period, U, is the maximum flood velocity, and h is the mean water depth. For wellmixed estuary, tidal range is generally large and the river discharge is generally low. Therefore, a wellmixed estuary has large Harleman's estuary number. Detailed discussion and study of this "gravitational circulation" can be found in Pritchard (1956), Hansen and Rattray (1966), and Dyer (1973). Tidal circulation is the water movement driven by tides at the mouth of the estu ary. If the tidal range is large and the river flow is low, tide is the most important force driving the current, especially for a shallow estuary such as Tampa Bay, which has an averaged water depth of 3.7 meters. Tidal force generally consists of more than one astronomic component. In Tampa Bay, the predominant astronomic constituents are the lunar semidiurnal M2 and solar diurnal 01 (Goodwin, 1987). Because of the com plexity of the geometry and bathymetry of an estuary, tidal circulations are often very complex. For example, the interaction between tidal currents and estuary boundary can sometimes cause a residual circulation pattern in which the timeaveraged tidal currents are ebbdirected on one side of an estuary cross section and flooddirected on the other side (Kjerfv and Proehl, 1979). Tee (1976) studied the tidal circulation in the Minas Basin of the Bay of Fundy. Results of his numerical model showed that the timeaveraged tidal circulation manifested itself as large horizontal eddies with diameters on the order of the width of the estuary. The rotation of these eddies indi cated that they were not caused by the Coriolis force but by the interaction between tidal force and coastal boundary. Winddriven circulation can be found in both lakes and estuaries. It is the current induced by the wind action on the free surface of the lake or the estuary. Winddriven 46 circulation is particularly pronounced in shallow lakes and lagoons. Depending on the wind and the water body, winddriven circulation can be very important even when gravitational and tidal circulations are important in an estuary. For example, in the event of an oilspill, the movement of oil at the water surface is primarily driven by the wind. Since oil is generally less dense than water, the velocity of the thin oil layer on the water surface can be very large. Because wind and other meteorological factors can be highly variable both tem porally and spatially, winddriven circulation is very complicated. Sheng and Lick (1979) developed a threedimensional model for simulating winddriven circulation in Lake Erie. Sheng et. al (1991) developed Cartesiangird and curvilinear grid three dimensional models to study winddriven circulation in Lake Okeechobee. Their study showed that winddriven circulation in Lake Okeechobee is often associated with se iche oscillations due to the sudden increase of wind over the lake in the afternoons. A recent study by Lee and Sheng (1993) indicated that winddriven current in Lake Okeechobee is affected by the vertical variation of lake temperature. In order to successfully simulate strong near surface current in the lake during seiches, it was necessary to include temperature in their model simulation. The development of numerical models for circulations in lakes and estuaries has advanced significantly during the last two decades due to increased demand for quan titative assessment of the human impact on water quality and the rapid improvement of the computer resources and numerical methods. The fast development of high speed computers has facilitated remarkable progress in threedimensional freesurface timedependent circulation models. Earlier freesurface threedimensional models (e.g., Leenderste and Liu, 1975) used explicit finitedifference methods to solve the threedimensional equations (Equa tions (4.1) (4.4)). The explicit algorithm is simple and has relatively little numerical diffusion but requires a very small time step to resolve the gravity wave propagation, 47 hence is very timeconsuming. Leenderste and Liu (1975) used a time step of 10 seconds for modeling circulation in Chesapeake Bay and San Francisco Bay. The threedimensional mode] system developed by Markofsky et al. (1986) was basically the same as Leenderste and Liu's model. They also used a very small time step (15 seconds) for Weser Estuary, Germany. To overcome the shortcoming of small time step in explicit 3D freesurface time dependent models, Simons (1974) and Sheng et al. (1978) applied the timesplitting technique to models of Lake Ontario and Lake Erie, respectively. The timesplitting technique is based on the fact that while the surface gravity wave propagates rapidly, the rate of change of internal flow (vertical flow structure) is much slower so that for the computation of the vertical flow structure it is possible to use a much larger time step than that for the external flow (water levels and the verticallyintegrated flow). Using the timesplitting technique, Simons was able to use a time step of 100 seconds for the external mode (flow) and a time step of 15 minutes for the internal mode with a horizontal grid of 5 km. However, using different At's for external and internal steps could lead to inaccurate results, due to inconsistency between the external and internal solutions. Thus, Sheng and Butler (1982) developed a 3D model which generally uses the same large steps for the external and internal modes. For the external mode, they treated all terms in the continuity equation and the surface slopes in the verticallyintegrated momentum equations implicitly and all other terms explicitly so that a larger time step is allowed for the external computation. By using the factorized implicit method for the external computation, the computing time for each time step was significantly reduced. Two types of vertical grid systems have been used in threedimensional free surface models. One is the zgrid and the other is the agrid. In zgrid 3D models, the vertical grid system is composed of fixed straight lines (layers), while in agrid 48 3D models, grids are vertically stretched through the following transformation: 0, = ((, Y,(3.1) zH(x,y,t) where a is the new vertical coordinate, z is the old vertical coordinate, H(x,y, t) is the total water depth, and ((x, y, t) is the surface elevation (measured from the mean water). The advantage of using the rgrid is that the same vertical resolution can he obtained in both the shallow water region and the deep water region. However, it has the disadvantage of creating larger numerical diffusion than the zgrid does, because the agrid enables the communication between gird points in the shallow and deep regions, which is not necessarily physically correct. The zgrid was used in the threedimensional models of Leenderste and Liu (1975), Lang et al. (1989), and Casulli and Cheng (1992). The agrid has been used in the 3D models of Blumberg and Meller (1987), Sheng and Butler (1982), and Sheng et al. (1990 and 1991). In order to significantly reduce the numerical error due to the use of agrid, Sheng et al. (1989c) modified Sheng's rgrid model to solve the salinity equation on the zplane (the horizontal gradient terms were calculated by first finding the salinity values at the zplane using linear interpolation from the salinity values at a points and then solving the difference equation on the zplane) in the 3D simulation of circulation in Chesapeake Bay. They were able to significantly reduce the numerical error of the agrid model and successfully simulate the stratification destratification cycle of salinity in the mid bay. Lee and Sheng (1993) and Choi and Sheng (1992) used the curvilinear grid 3D model (CH3D) originally developed by Sheng to simulate the circulations in Lake Okeechobee and James River, respectively. The present study uses a Cartisiangrid 3D model (EHSM3D) originally developed by Sheng to simulate the winddriven circulation in Lake Okeechobee. 49 3.2 Wave Boundary Layer 3.2.1 Pure Wave Motion Waves in lakes and estuaries are primarily generated by wind. High frequency oscillatory currents are usually associated with wind waves. Slowly varying currents are the result of tides, density gradients, or wind. Because of its high frequency oscillatory behavior, the wave boundary layer is much thinner than the current boundary layer. The thickness of a fully developed wave boundary layer is on the order of VA.T/(2r), where T is the wave period and A, is the vertical eddy viscosity near the bottom. If A. is 1 cm2/s and T is 2 seconds, the thickness is less than 2 cm. Consequently, the bottom shear stress induced by wind wave motion is much larger than that induced by slowly varying current. Waves can significantly affect the resuspension of bottom sediment in shallow lakes and estuaries (Sheng et al., 1989b). This is not only because a wave is capable of inducing a higher bottom shear stress than that by a comparable current, but also because the bed can undergo fluidization process under wave action. Various studies (e.g., Maa, 1986) indicated that waves are primarily responsible for the fluidization of bottom sediments. The erosional strength of the bed is weakened by the dynamic loading of waves which transmits normal and shear stresses downward to buildup excess pore pressure in the bed and thereby results in the rupture of interparticle bonds. As shown in Figure 3.1, the critical shear stress for erosion under wave action is much smaller than that without wave action. Because of the significant effect of waves on sediment dynamics in lakes, estuaries, and ocean, many researchers have investigated wave boundary layer dynamics. For example, a series of laboratory experiments of oscillatory flows were conducted by Jonsson (1966) and Jonsson and Carlsen (1976). By assuming that the velocity distribution is logarithmic, Jonsson and Carlsen (1976) calculated the bottom shear stress near a rough wall by integrating the momentum equation from the bottom to SoWithout Waves (Parchurs,1984) .With Waves (Maa,1986) w 0 z 77< S0.2 / Wave I Effect  z 0 "' ^  o . ' 0 0 5 10 15 PERIOD OF CONSOLIDATION (days) Figure 3.1: Critical shear stresses for erosion with and without wave action (from Maa, 1986). the top of the wave boundary layer. The friction factor, f,, was then calculated from the of the quadratic shear stress law: ^ = /PfU (3.2) where p is water density and Uoom is the maximum bottom orbital velocity computed from the linear wave theory. For an oscillatory motion with 2 m/sec amplitude and 8.39 sec period, they reported the logarithmic layer thickness to be approximately 6 cm, and the maximum bottom shear stress to be on the order of 400 dynes/cm2. Based on the experiments, they derived a semiempirical formulation for f,. Kamphuis (1975) also did extensive laboratory measurements of the maximum bottom shear stress induced by waves. He produced a diagram of the friction factor for different values of bottom roughness and Reynolds' number. Kajiura (1968) solved the vertically onedimensional timedependent wave bound ary layer equation using a threelayer eddy viscosity formulation. He assumed a time invariant eddy viscosity in studying the turbulent boundary layer dynamics: constant 51 eddy viscosity in the small inner layer (z < k./2, where k, is the Nikuradse roughness height), linearly increasing eddy viscosity in an overlap layer (k./2 < z < b, where 6 = U/(20), U,,(= V4( p) is the maximum shear velocity, and w is the angu lar frequency), and another constant eddy viscosity for the region outside b. From this assumed profile of eddy viscosity and by using the equation of motion, Kajiura obtained a theoretical solution and developed a theoreticalnumerical method to cal culate the friction factor in equation (3.2). Brevik (1981) simplified Kajiura's eddy viscosity concept by avoiding the inner layer to reduce the calculation in Kajiura's model. His results were very similar to Kajiura's. One common weakness in the Kajiura and Brevik models is the assumption of timeindependent eddy viscosity. Studies by Horikawa and Watanabe (1968) showed that eddy viscosity is a function of time and can vary significantly during the wave period. Bakker (1974) adopted a mixing length model similar to that of Prandtl to take into account the time dependence of eddy viscosity. He combined the mixing length model with the equation of motion to obtain an equation for frictional velocity (/'/p). However, his model did not show much improvement over Kajiura's or Brevik's model, because his mixing length assumption is questionable for oscillatory flow. To represent the temporal and spatial dependence of eddy viscosity in the wave boundary layer more precisely, Sheng (1982) used a Reynolds stress model to calculate the turbulence and the eddy viscosity. His results for currents and shear stresses compared well with the measurements of Jonsson and Carlsen. More recently, Sheng and Villaret (1989) used a turbulent kinetic energy (TKE) model to calculate the turbulence and the eddy viscosity. They also achieved good agreement with Jonsson and Carlsen's experiment. 3.2.2 WaveCurrent Interactions Wind generated waves in lakes and estuaries typically have a wave period ranging from 2 to 10 seconds. Such short period waves can feel the influence of the bottom 52 in approximately 10 30 meters of water, where waveinduced nearbottom orbital velocities are of the same order of magnitude as those due to winddriven and/or tidedriven currents in lakes and estuaries. For example, in Lake Okeechobee, the maximum winddriven current is on the order of 25 cm/sec, while the nearbottom orbital current can be equally strong. However, as discussed in the previous sec tion, because of the high frequency oscillatory nature of waveinduced current, the wave boundary layer is much thinner and waveinduced bottom shear stress is much higher than those associated with a slowly varying current of comparable magnitude. Therefore, waves are more capable of causing sediment resuspension when a current of comparable magnitude may be too weak to initiate sediment motion. Although waves can induce strong bottom shear stresses, their ability to transport sediments is very weak. Because the wave boundary layer is very thin, resuspended sediments near the bed may not be easily transported upward. Horizontally, because the oscillatory wave motion produces little net flow over a wave cycle, the net transport of sediments is only of second order importance. However, the combined presence of a current and a wave can cause significant sediment transport in both the vertical and horizontal directions. The upward transport of resuspended sediment is mainly due to the turbulent mixing processes in the current boundary layer which is generally quite thick (sometimes the current boundary layer extends throughout the entire water column). In order to consider the combined effects of waves and currents on sediment transport processes, it is necessary to first study the combined flow motion of waves and currents and quantify their interactions. Because the bottom shear stress and the eddy viscosity inside the wave boundary layer are associated with both the wave and the current, the combination of waves and current is nonlinear. The bottom shear stress for the combination of a wave and a current is not the linear addition of that due to the pure wave action and that due to pure current action. It is known that a 53 slowly varying current will experience more resistance if a wave is added to the flow (Grant and Madsen, 1979). In addition, current inside the wave boundary layer will change direction if the current is not in the same or opposite direction as that of the wave propagation (Davies et al., 1988). This is because the mean bottom shear stress over the wave period should always be in balance with the force driving the motion of the current outside the wave boundary layer. The governing equations for the combined wave and current motion are the ver tically onedimension equations of motion: Du 1pp (9 OT pOX + 5 p (3.3) Dv_ 1 09p+a (3.4) Tt _p0 OyY where u and v include the wave orbital velocities u, and v, as well as the slowly varying current velocities u, and v,, the pressure p includes p, and pc, and r. and ry are turbulent shear stresses in x and ydirections. The turbulent shear stresses T.. and ry are related to the vertical shear: (rT.) =A,, (u,v) (3.5) where A. is the vertical turbulent eddy viscosity. Grant and Madsen (1979) assumed the turbulent eddy viscosities (A,,) inside and outside the wave boundary layer as linear functions of the distance from the bed. By defining the characteristic shear velocities of respective current and wave boundary layers using a combined wavecurrent friction factor, the turbulent eddy viscosities were assumed to be A,, = ru:z z < <,. (3.6) A = Ku*z z > 6, (3.7) where K is the von Karman constant, 6, is the wave boundary layer thickness, and u:, and u* are the characteristic shear velocities of wave and current boundary layers, respectively: u(=f= V2)"ub1 (3.8) u Tb = = (f, ) ,1/2ub (3.9) where ub is the maximum bottom orbital velocity computed from linear wave theory, f is the combined wavecurrent friction factor, i, is the timeaveraged shear stress over the wave cycle, Ti,,ra is the maximum boundary shear stress due to the combined wave and current, and V2 and a are functions of current velocity, wave orbital velocity, and the angle between the current and the wave. In Grant and Madsen's model, it was assumed that the thickness of the logarith mic layer is constant and the eddy viscosity is timeinvariant. As already pointed out by Sheng (1984), the thickness of the logarithmic layer and the eddy viscosity change significantly with time during a wave cycle. Because of the timeinvariant eddy viscosity, Grant and Madsen's model can not predict the phase shift in the wave boundary layer correctly (Sheng, 1982). Davies et al. (1988) used a turbulent closure scheme to estimate the eddy viscosity shown in Equation (3.5). The eddy viscosity calculated is timedependent, but may not be accurate for the surface layer of the water column because they assumed a linear increase of mixing length with the distance above the bed, which is not correct (Nezu and Rodi, 1986). Although some model results were shown in their paper, no comparison between model results and measured data was made. Wavecurrent models that compare well with laboratory data were developed by Eidsvik (1990) and by Sleath (1991). Eidsvik used Prandtl's turbulent closure model for the region far above the wave boundary layer (&b) and a oneequation turbulent closure model for the region significantly below 6,. The effect of the oscillatory motion on the current is considered solely by a scalar valued bottom roughness (zoo), which is on the order of the wave boundary layer thickness, zoo is estimated by matching the current calculated from above and below b, at a height (a4b,) above the bed, where 55 the coefficient a4 should be determined by fitting the model results with measured data. Eidsvik's model is relatively simple, but requires very fine vertical resolution (Az < H) and a very small time step (At < ) because of the explicit treatment of the turbulent closure. The coefficient a4 may depend on the characteristics of the wave and the current, as well as the angle between them. In conclusion, numerous models for wavecurrent interaction have been developed in the last two decades. Most simple models contain simplifying assumptions and hence some limitations. For example, Grant and Madsen's model is only applicable for weak current situations. Onishi et al. (1993) pointed out that for strong current situations, Grant and Madsen's model results of bottom drag coefficient and bottom shear stress under wavecurrent interaction are smaller than those under pure current action. Eidsvik's model requires a lot of measured data to determine model param eters and is thus limited in its ability for prediction. On the other hand, numerical models with a sophisticated turbulent closure model are more robust because of their ability to accurately estimate the eddy viscosity. This study uses the TKE model developed by Sheng and Villart (1989) to calculate the the bottom shear stresses induced by the combined wave and current action (see Section 4.3 for detail). 3.3 Turbulent Mixing Turbulent mixing is a very important physical process which can affect the flow pattern and the horizontal and vertical distributions of sediments and nutrients. The eddy viscosity mentioned in the last section represents the vertical mixing produced by the turbulent motion. As can be seen, the basic differences among different wave boundary layer or wavecurrent interaction models are the different methods for pa rameterization of the eddy viscosity. The concept of eddy viscosity and eddy diffusivity has been widely used in mod eling the turbulent mixing of momentum, heat, and mass (salinity and other species). The basic assumption of the eddy viscosity/diffusivity model is that the turbulent 56 flux of momentum u'u'j (correlations of velocity fluctuations) is the product of mean velocity gradients and "eddy viscosities," and the turbulent mass flux is the products of mean mass concentration gradients and "eddy diffusivities": A ('u= + au, ii (3.10) =B, (3.11) where i,j can be 1,2, and 3, u, and u, are the mean velocity components, u' and u' are the fluctuations of u, and uj, respectively, represents the mean concentration of suspended sediment or nutrients, 0' is the fluctuation part of 0, ij is the Kronecker delta function, q2 = Uu = twice of turbulence energy, and At and BA are turbulent eddy viscosity and diffusivity, respectively. The horizontal turbulent eddy coefficients represent the mixing due to small scale subgrid scale variations of velocity and concentration which can not be resolved by the large grid size. Therefore, the values of horizontal eddy viscosities/diffusivities in a 3D model are dependent on the grid size used for a particular simulation. Generally, constant values can be used in a threedimensional simulation (e.g., Sheng et al., 1991; Lang et al., 1989). For vertical turbulent mixing, Lehfeldt and Bloss (1988) used a mixinglength based algebraic model: A,, = p,,A 2 m~o = ooV z)P + (LTV) 2 (3.12) A, = f(Ri)A,,, (3.13) B, = g(Ri)A, (3.14) where A,, is the eddy viscosity/diffusivity in nonstratified flow situation, A,, is the mixing length and is assumed to be a linear function of z, increasing with distance above the bottom or below the free surface with its peak value at middepth not ex ceeding a certain fraction of the local depth (in the presence of strong waves, turbulent 57 mixing in the upper layers may be significantly enhanced. In such a case, the length scale A. throughout the upper layers may be assumed to be equal to the maximum value at middepth) and f(Ri) and g(Ri) are stability functions dependent on the Richardson number Ri defined as: Ri 9 p/z (3.15) Po (9u/9z)2 + (v/z)2 (3.15) The stability functions f(Ri) and g(Ri) are traditionally assumed to have the following forms: f(Ri) = (1 +aRi)" (3.16) g(Ri) = (1 + a2Ri)R (3.17) As discussed by Sheng (1983), there exist great discrepancies among the many existing empirical forms of the stability functions. To determine the coefficients in these stability functions, one needs sufficient data to achieve the "best fit" between modeled and measured results. Therefore, specific formulae and coefficients will vary with the problem, the environment, and the available data. Munk and Anderson (1948) used the following formula in a study on the marine thermocline: f(Ri) = (1 + lORi)1/2 (3.18) g(Ri) = (1+33.3Ri)3/2 (3.19) Others (Kent and Pritchard, 1959; Blumberg, 1975; Bowden and Hamilton, 1975) used f(Ri) and g(Ri) in the form of Equations (3.16) and (3.17) but with coefficients significantly different from those in Equations (3.18) and (3.19). The secondorder closure models, on the other hand, resolve the dynamics of tur bulence by including the differential transport equations for the turbulence variables; i.e., the secondorder correlations (u'. u', j' and W). In the most complicated 58 case, i.e., the Reynolds stress model, differential transport equations of the following forms are solved (Lewellen, 1977; Sheng, 1982; Sheng and Villaret, 1989) u'f +, u ,OU9, uj 9Ou u *4 ; +=  +, + +gu (3.20) at Oxk xk o (ap, " 2enlk1uu 2aijflu'iu[ + v,57 ( qAP Xk OXI (^uj 2t< 2avuu A 3 5j' 'i 12A A2 Aujo u090 a 6U.1 u ,v' 3i t 7" OW at uuj0 4u+g, (3.21) O jXj = + ( axj ) , o9 au4 ~ u~k + vo2j ( a ') A A Aqe 2au77 X~u~2 + V% +,+,, a,01 ... ao o9 (901O01 +t OUj 2' + v ; O u (3.22) 2bsq 2as A A T where x, are coordinate axes, t is time, (u,, uj, uk) are the mean velocity components, (u', u', u') are the fluctuating velocity components, 4 and 0' are the mean and fluc tuating values of density or temperature or salinity, eik is the alternating tensor, i1 is earth rotation, Sj is Kronecker delta, q = (uIu)1/2 is the total rootmean square fluctuating velocity, and A is the turbulence macroscale which is a measure of the average turbulent eddy size. A total of five model coefficients: a, A, b, v., and s are contained in the above equations. By matching model results to data from critical laboratory experiments where only one or two of the turbulent transport processes are important, these model constants were found to be (Lewellen, 1977) a = 2.5, A = 0.75, b = 0.125, v, = 0.3 and s = 1.8. An additional equation for A is needed to close the system. 59 Despite the complexity of the Reynolds stress models, such models do contain more physics and hence the model constants remain unchanged for all model ap plications. The above Reynolds stress model was used to produce successful sim ulations of wave boundary layer (Sheng 1982; Sheng 1984; Sheng 1986; Sheng and Villaret, 1989), wavecurrent boundary layer (Sheng 1984), and vegetated boundary layer (Sheng 1982). Because of the complexity of Reynolds stress models, simplified secondorder clo sure models have been developed. These include the turbulent kinetic energy (TKE) closure model (Sheng and Villaret, 1989) and the equilibrium closure model. Both of these models contain the assumption that the turbulence time scale is much less than the mean flow time scale. The TKE closure model resolves the turbulent dy namics by including the differential transport equation for q2 (= u'u' + v'7 +w'w ), while the other secondorder correlation equations are simplified to algebraic equa tions. The TKE closure model is similar to the Level3 model of Mellor and Yamada (1982) and the k e model of Rodi (1980), although both models favor the solution of turbulence dissipation (e) instead of length scale (A). The equilibrium closure model (Sheng, 1985; Sheng and Chiu, 1986; Sheng et al., 1989d) is slightly more simplified than the TKE closure model, since only algebraic equations are used to resolve the secondorder correlations. The equilibrium closure model is similar to Mellor and Yamada's Level2 model. Instead of solving the differential transport equations for A which consists completely of modeled terms, simplified secondorder closure models often solve a set of integral constraints for A (Sheng 1985; Sheng and Chiu 1986). The equilibrium closure model computes the vertical turbulence based on the assumption that the local equilibrium condition is valid when the time scale of mean flow is much larger than that for turbulence and when turbulence varies little over the turbulence macroscale. The equilibrium closure model is significantly simpler than the complete secondorder closure model (Reynolds stress model) and has been found 60 to give very good results for mean flow, salinity and temperature with little or no tuning on model coefficients. In the equilibrium model, a set of algebraic equations are solved for the second order correlation quantities to obtain the stability functions f(Ri) and g(Ri) in terms of the mean flow variables. These equations, when written in dimensional and tensor forms, are (Sheng et al., 1989d) '' 0 U 'j gu""  uiP (3.23) ox xk T Po gPo 2eiktflkU~u' Crtk~u'kU; ^ I ) q A u ,3) 1i2A 0 T'P__ p __ D 7u f p" (3.24) 2eijkju7 0.75q! 0 P + 0.45qp (3.25) 49j A where the subscripts i,j, k can be 1, 2, or 3. Hence Equation (3.23) represents six equations and Equation (3.24) represents three equations in general. A total of five model coefficients are contained in the above equations. These coefficients were deter mined from comparing model results with data from critical laboratory experiments where only one or two of the turbulent transport processes are dominant. These co efficients remained "invariant" in application of the equilibrium closure model, the kinetic energy closure model (Sheng and Villaret, 1989) and the Reynolds stress model (Sheng, 1982; Sheng, 1984; Sheng and Villaret, 1989). As shown in Sheng et al. (1989d), q2 can be determined from the following dimensional equations when mean flow variables are known: 3A2bsQ4 + A[(bs + 3b + 7b2s)Ri Abs(1 2b)]Q2 (3.26) + b(s+3+4bs) Ri2+(bsA)(1 2b)Ri=0 where the model constants A, b, and s are "invariant" and have the values 0.75, 0.125, and 1.8, respectively, and Q = q (3.27) A Z F _9p Ri = z (3.28) (au (Ov 1' ()2 + (")2 The total rootmeansquare turbulent velocity q can then be obtained from the above equations after the mean flow variables are determined at each time step: q = QA au 2 + (3.29) Once q2 is computed, the vertical eddy coefficients can be computed from A + tiw'w' A = A A q (3.30) bs w' B (bsw)A q2 A q (3.31) where w = Ri/(AQ2) (3.32) wi = w/(1 w/bs) (3.33) w'w' = 1 2 (3.34) 3(1 2,b) q(.4 In the complete secondorder closure model, a differential transport equation for the turbulent macroscale A is derived. However, the A equation contains four model 62 coefficients that must be calibrated with experimental data. For ease of application, the turbulent macroscale A is assumed to satisfy a number of integral constraints. First of all, A is assumed to be a linear function of vertical distance immediately above the bottom or below the free surface. In addition, the turbulent macroscale A must satisfy the following relationships (Sheng et al., 1989d): dA < 0.65 (3.35) A < C, H (3.36) A < C H (3.37) A < C2'6,2 (3.38) A < (3.39) where C1 is a number usually within the range of 0.1 to 0.25; H is the total depth; Hp is the depth of the pynocline; C2, which is between 0.1 and 0.25, is the fractional cutoff limitation of turbulent macroscale based on 6q,2, the spread of turbulence determined from the turbulent kinetic energy (q2) profile; and N is the BruntVWis.la frequency defined as: N ( agP1/2 (3.40) \. p9z~ Although the simplified secondorder closure model presented above is strictly valid when the turbulence time scale (A/q) is much less than the mean flow time scale and when turbulence does not change rapidly over A, it has been found to be quite successful in simulating vertical flow structures in estuarine and coastal waters. Recent simulations of windinduced mixing of salinity in Chesapeake Bay (Sheng et 63 al., 1989c; Johnson et a!., 1989) and winddriven circulation and sediment transport in Lake Okeechobee (Sheng et al., 1989b) gave good results with little or no tuning of model coefficients. Compared to the equilibrium closure model, the kinetic energy closure model does not assume that the turbulence created at any second is in equilibrium with the turbulence disappearance due to the advection, diffusion and dissipation. It contains an extra equation for the turbulent kinetic energy, q2, which can be obtained directly (let i = j) from Equation (3.20): Nq2 tq2 2 t ui ui' .q2 q3 Suk = 2uukjO + 2g + 0.3 a(qA ) 3 (3.41) 9t x 9Xk xk O9xk 4A The vertically onedimensional version of this equation reads aq2 , ui u7 9 (qA q3 27 = 2ukk + 2g + 0.3 A) 4A (3.42) 5T dxk ifo Ix 9X3 4 A or 9q2 9u 9v 7?' Q 9Q2 q3 S= 2u'w' 2v'w' + 2g +0.3 (qA q) '(3.43) Wt 27 __2 9z_ p 9z 8z' 4A ' The secondorder correlation terms of fluctuating velocities and density, u'w7, y'w', and w'i' can be calculated from Equations (3.23), (3.24), and (3.25), or the solutions of the equilibrium closure model discussed above. They have the following forms: uw 1 = 8F+336 q u (3.44) 4A A Oz r= 81+336%P q36v (3.45) 4A A =^z (3.45) 4A A c9z 108F + 144T q3,9p (3.46) 4A A Oz where r = ()2 (3.47) g = 9_ (3.48) ozp 02 P A = 81()4 804 ()2 q2 + 928( )2q4 (3.49) A A 8z A az Substitute the secondorder correlation terms in Equation (3.43) with Equations (3.44) (3.45), the q2 equation can be solved with a finite difference method. Once q2 is solved from Equation (3.43), the vertical eddy viscosity (A.) and diffusivity (B.) can be calculated by Equations (3.30) and (3.31), where Q is calculated by Equation (3.27) and the turbulent macroscale A is also determined by the constraints mentioned above. Sheng and Villaret (1989) used the TKE closure model to simulate a wave bound ary layer, a thermally stratified boundary layer and a sedimentladen boundary layer. Mellor and Yamada (1982) generally favor the use of their Level3 model or Level 2.5 model. The k e models used by Rodi (1980) and others have also been quite successful. On the other hand, the equilibrium closure model has also been found to be very successful for simulating the vertical mixing in estuaries (Sheng et al., 1989c; Sheng et al., 1989d), despite the slightly further simplification. In this study, the TKE closure model is chosen to calculate the vertical eddy coefficients A. and B,, while the horizontal eddy coefficients AH and BH are assumed to be constant (both spatially and temporally). 3.4 Sediment Dynamics 3.4.1 Settling and Flocculation Settling is a process by which sediment particles fall through the water column because of density differences between sediment particles and water. This process is closely related to the flocculation process in the water column. In the flocculation process, sediment particles are flocculated together due to the turbulent shearing 65 of water, differential settling of sediment particles, Brownian motion, and electro chemical/biochemical attraction. The flocculation process results in the formation of flocs, which can generally fall through the water column faster than the primary sediment particles. The simplest parameterization of settling is the use of a constant settling velocity. For a single spherical sediment particle in the laminar flow, the settling velocity is described the Stoke's formula: I g(p. p.)d. w, 1 (3.50) w,18 r/ where p. and p, are densities of water and the sediment particle and d, is the diameter of the particle. Stoke's formula is a analytical solution for a very idealized case. In practice, sedi ment particles are not necessary spherical. Furthermore, instead of a single sediment particle, many particles are settling to the bottom at the same time. The fact that many particles are present in the water column not only suggests the influence of the settling behavior of one particle on the settling of other particles around it, but also the aggregation and flocculation of fine sediment particles (d < 60/gm) because of the significant electrostatic forces (cohesion) and the collision among them. The end product of the flocculation, flocs, can fall faster to the bottom than the primary sedi ment particles due to the relatively large size. However, if the sediment concentration is too higher (_> 10 g/l), collisions among sediment particles will actually retard the particle settling. This is the socalled "hindered settling". The cohesion of sediment particles depends on the type of clay minerals, cation exchange capacity (CEC), sodium adsorption ratio (SAR), and the salinity and pH of water. Gibbs (1983) studied the flocculation rates of clay minerals, and found that flocculation begins at low salinities of 0.05 to 1 ppt. Hunt (1982) performed labora tory experiments on flocculation of kaolinite and illite induced by Brownian motion and laminar shear, and obtained an equilibrium particle size distribution within a 66 few minutes. However, the equilibrium condition may not be valid for sediments in turbulent estuaries and lakes, where the flocculation time scale is typically on the order of a few hours. Most of flocculation models developed so far followed the basic theory proposed by Smoluchowski in 1918. Smoluchowski determined the rate of collisions among two particle groups due to Brownian motion (f3b), differential settling (/3d), and turbulent shearing (pt) to be as 2BT (r, + rj)2nnnj (3.51) 3pt rirj d = r(ri + rj)wi wjlnin, (3.52) 4G , S (ri + rj)3 nnj (3.53) 3( wherein i is the molecular viscosity, B is the Boltzman constant, T is the absolute temperature, r, and r, are radii of particles of size i and j, respectively, n, and n, are number densities of size i and j, respectively, wi and wj are settling velocities of size i and j, respectively, and G is the microshear which the particles experience. G is assumed uniform with a constant velocity gradient du/dz in Equation (3.53). In general, it is impossible to represent G by du/dz because there are turbulent velocity gradient superimposed upon the mean velocity gradient. Instead, a socalled 'absolute velocity gradient' defined as follows should be used for G: G= (3.54) wherein e is the total energy dissipation and v is the kinematic viscosity of water. Farley and Morel (1986) solved the particle size distribution equation following the basic work of Smoluchowski (1918). However, their studies treat the entire par ticle size distribution, and are too cumbersome for practical applications in a multi dimensional model. In addition, explicit accounting of the effect of turbulent motion 67 on flocculation was not included. Tsai et al. (1987) demonstrated the effect of fluid shear on flocculation. However, only shear stress, but not turbulence level or shearing rate, was measured in their study. Using Smoluchowski's equations for collision, Burban et al. (1989) developed a flocculation model with an additional breakup term induced by turbulent shearing. However, they could not fit their model results to their measured equilibrium par ticle size distribution in a laboratory device without modifying the breakup term and adjusting several model coefficients. Thus, the laboratory data of Burban et al. were used primarily for calibration of their flocculationbreakup model. Due to the lack of comprehensive data, the detailed dynamics of flocculation are still not well understood particularly in field conditions. Moreover, since their model is primarily a boxmodel, work is needed to extend it to multidimensional models for simulat ing field situations. In another study, Farley and Morel (1986) performed laboratory experiments on flocculationinduced mass removal in high sediment concentration en vironments (c > 0lg/i). They produced an empirical equation for the removal rate which depends on differential settling, fluid shear and Brownian motion: dt = Bc2.3 Bt c19 Bb C13 (3.55) where c is the volumeaveraged sediment concentration in the device, Bd, Bt and Bb are contributions of differential settling, fluid shear and Brownian motion on the collision frequency of particles, respectively. The empirical coefficients of 2.3,1.9 and 1.3 may require sitespecific parameter tuning. To apply their results to multidimensional modeling, one can derive an expression for the settling velocity in the following: W = A (B 23 + Btcls + Bbc13) (3.56) The above equation for w,, however, is only valid for high concentration situations where c is on the order of lOg/l or higher. Moreover, it does not contain any effect 68 of particle breakup on settling. This formula requires sufficient data to allow tuning of the model coefficients. Krone (1962) considered the flocculation at a constant collision rate and derived an analytical equation relating the settling velocity with the sediment concentration: w. = ac'" (3.57) where a is an empirical constant and m = 4/3. Such a power law of settling velocity has been supported by numerous field and laboratory measurements, but with slightly different exponent values (Table 3.1). For example, Thorn (1981) obtained the following settling velocity for mud from Severn Estuary, England: w = 0.513c1"29 c < 2g/l (3.58) For c > 2 g/l, he found that hindered settling dominated, and the settling velocity decreased with the increasing sediment concentration: w, = 2.6(1 0.008c)'465 (3.59) Sheng and Lick (1979) measured the settling velocity for sediments collected in Lake Erie in laboratory and used the measured settling velocity in a threedimensional simulation of sediment transport in Lake Erie. Gibbs (1985) measured the settling velocity of freshly formed flocs of illite, kaolinite and smectite in a laboratory and found that the settling velocity is related to the floc diameter through a nonStokesian relation which indicated that the floc density decreased and the drag increased as the flocs increased in size. Toorman and Berlamont (1993) measured the settling rate in a settling column, while considering the combined problem of settling and consolidation. He found two distinct settling modes in all his experiments: one corresponds to the hindered settling and the other corresponds to the consolidation, van Leussen et al. (1993) measured the settling velocity in the Ems Estuary in The Netherlands. Their Table 3.1: Measured exponents m in Krone's Equation. Exponent rnm Water body Reference San Francisco Bay (lab measurement) Chao Phya Estuary Thames (spring tide) Thames (neap tide) Thames (lab measurement) Severn Estuary Thames Elbe Estuary (limnetic zone) Elbe Estuary (brackish zone) Elbe Estuary (limnetic zone) Elbe Estuary (brackish zone) Ems Estuary (The Netherlands) Krone (1962) NEDECO (1965) Owen (1970) Owen (1970) Owen (1970) Thorn (1981) Burt (1986) Puls and Kuehl (1986) Puls and Kuehl (1986) Puls et al. (1988) Puls et al. (1988) 1.33 1.21.3 1.1 2.2 1.0 1.29 1.37 2.92 2.55 2.07 1.97 2.5 3.6 (1993) van Leussen et al. 70 measurements showed that there exists no unique relationship between the settling velocity, w,, and the suspended sediment concentration, c, over the entire estuary, although for a specific location there seems to be an unique relation between them. They used Equation (3.57) to fit their data with an exponent ranging from 2.5 at the mouth of the estuary to 3.6 at the upper stream of the estuary due to the different tidal, salinity, and flow conditions. They suggested that a locationdependent relation between the settling velocity and the suspended sediment concentration should be used in studying fine sediment transport in estuaries. Hwang and Mehta (1989) measured the settling velocity for sediments collected in Lake Okeechobee. They empirically parameterized the effects of flocculation and hindered settling on the settling velocity by a formulation proposed by Wolanski et al. (1989): Sac' (360) W (c2 + b2)" (3.60) where c is the sediment concentration and a,b,m, and n are empirical constants determined from laboratory experiments. The shape of w, as a function of c is shown in Figure 3.2 (Hwang and Mehta, 1989). As the sediment concentration increases, flocculation occurs as the result of increased collision among particles, thus leading to the increase of the settling velocity w,. As c exceeds 3000 mg/t, w, starts to decrease due to hindered settling, i.e., interference on settling due to the presence of other particles. It should be pointed out that such w, formulations expressed as functions of the sediment concentration c do allow the empirical parameterization of flocculation due to Brownian motion and differential settling. However, because of the use of the settling tube in the measurements, the real settling velocity in the field may be different from that obtained in the settling tube. Three factors may cause the deviation of measured settling velocity from the real one: (1) the lower turbulence level in the tube, (2) the small dimensions of the tube, and (3) the breakup of flocs during 0 1 10 CONCENTRATION, C(gL'1) Figure 3.2: Settling velocity as a function of suspended sediment concentration, with flocculation occurring at high concentrations (after Hwang and Mehta, 1989). 72 the sampling process in the field. Therefore, the measured empirical relationship between settling velocity and sediment concentration may require some adjustments if used in the numerical model. 3.4.2 Deposition Deposition is the process that governs the removal of sediment particles from the water column. The speed with which particles are deposited onto the bottom from the water column is the deposition velocity. Once the sediments are brought near the bottom by settling, deposition determines their removal from the water column. Deposition in estuaries and lakes are determined by the settling velocity, shear strength, and concentration of depositing aggregates as they encounter the higher velocity gradients near the bed (Krone, 1993). Many past studies (Lavelle et al., 1984; Bedford et al., 1987) assumed that the deposition velocity is the same as the settling velocity, which can be a constant or determined by Stoke's settling formula (Equation (3.50)). Krone (1962) observed that the occurrence of the final deposition of sediments is related to the flow condition just above the bed. Under calm flow conditions, most of the sediment particles in the bottom layer will deposit onto the bed. While the flow becomes stronger and stronger, however, more and more sediment particles will refuse to settle onto the bed. He introduced a probability function for deposition as followed: P=0 ) (3.61) where p is the probability of deposition, Tm is the bottom shear stress, and Ted is the critical shear stress for the deposition. Equation (3.61) shows that when rb = 0, the probability of deposition is 1, and when Tb Trd, there is no deposition. The probability of deposition of sediment particles in the bottom layer increases linearly as the bottom shear stress decreases. The deposition velocity, which is defined as the deposition rate D divided by the 73 sediment concentration near the bed cl, is then wp, or v= (1 + 1 ) (3.62) 2 Ted Ted The advantage of Krone's model lies in its simplicity. Laboratory experiments us ing sediments collected from the harbor at Redwood City, on South San Francisco Bay showed that results of Krone's model agree well with the experimental data (Krone, 1993). However, Krone's model lacks rigorous theoretical foundation. The empirical constant, Ted, incorporates many processes (such as the strength of aggregates, the flow condition near the bed etc.) associated with the deposition. Because of the empirical nature of the model, there is no guideline as to how red can be estimated. Based on rigorous theoretical consideration of the deposition process, Sheng (1986) considered the flow over a surface covered with several layers with the possible pres ence of vegetation, by expressing the deposition velocity as 1 ( = vd (Vdh)' + (Vd.)' + (Vd)' (3.63) where (w'c')o represents the deposition flux at the bed, cl is the concentration at the first grid point above the bottom, and Vdh, vd,, and vd, represent the deposition velocity within the hydrodynamic layer (logarithmic layer), the viscous sublayer, and the canopy layer, respectively Vdh = (3.64) 751n z' Vd8 \V ] Ut ud ().3 (D)O7u 3.5 + 0.l' [ exp(O.OSqT,/r)] + g Ut V 1 + LAI vdc = Vd, 1 + DT /D1 (3.66) 74 where z, is the distance of the first grid point above the bottom, DB = 2.176 x 109 /rp is the Brownian diffusion coefficient, v is the fluid kinematic viscosity, u. = v /f is the friction velocity, rb is the bottom stress, ul = (u./k)tn(z6/zo) is the fluid velocity outside the laminar sublayer, with zs = lOz/u. as the height of the viscous sublayer. 7, = w,/g is the relaxation time of sediment particles, q = (32)'/'u. is the total turbulent fluctuating velocity, LAI is the leaf area index (total wetted leaf area per unit bottom area), and Dp/DI is the ratio of the profile drag to the skin friction drag in the vegetation area. In the absence of a vegetation canopy, vd, = 0. The derivation of Sheng's deposition model is based on fundamental consideration of fluid dynamics principles and the fact that the inverse of vd is a "resistance" which is additive through the various layers. The advantage of Sheng's model is that vd can be explicitly computed if the turbulent and sediment parameters are known. The disadvantage, however, is that it is difficult to obtain the field data required for the validation of the model. It can be concluded that Sheng's deposition model is more general than Krone's one, because it considers several dominant processes and does not require the adhoc tuning of empirical constants. There is no Td in Sheng's model. In other words, the critical shear stress for deposition in Sheng's model is infinite. Krone's deposition model is simple and widely used, but requires the rather arbitrary choice of the em pirical constant of critical shear stress for deposition (Ted), which has little theoretical foundation. 3.4.3 Erosion Erosion is the process by which sediment is suspended from the bottom into the water column due to a large bottom shear stress induced by currents and waves. The erosion or resuspension of bottom sediment depends on the bed structure and the characteristics of the flow above the bed. The softer the bed, the more easily the bottom can be eroded. The stronger the flow above the bed, the higher the erosion 75 rate becomes. The mechanism of erosion for fine sediments is different from that of coarse sediments. For noncohesive sediments, erosion takes place when the lift force acting on a grain is larger than the downward force, while for cohesive sediments, the binding forces among flocs should be overcome before erosion occurs. Most of the existing erosion models can be grouped into two catalogues: (1) critical current model and (2) critical shear stress model. Markofsky et al. (1986) and Wolanski et al. (1988) calculated the erosion rate based on the depthaveraged current as follows: E = E(U I+ 1) (3.67) where E, is the erosion rate constant, U is the depthaveraged velocity, and uc, is the critical velocity for erosion. Because the upward motion of sediment particles is directly related to forces acting on the particles rather than the depthaveraged velocity, the critical current model as shown in Equation (3.67) is not very reasonable, although it is always possible to adjust the u, value to fit the data. On the other hand, Partheniades (1962) proposed a critical shear stress model for the erosion rate as follows: E = E.o( rb (.8 E= ( 1+ 1)" (3.68) where E. is the erosion rate constant, rn is the bottom stress, and Tr is the critical stress for erosion. Mehta et al. (1982) pointed out that the above critical shear stress model is appropriate for a dense bed in which the bed properties are vertically uniform. For a partially consolidated bed, both the density and the bed strength increase with the depth. The erosion rate for this kind of soft bed should be calculated using the following equation (Parchure and Mehta, 1985): E = Ef exp[a(b r)1/]'1 (3.69) wherein Ef is the floc erosion rate and r, is a measure of bed resistance to the shear 76 stress. re increases with the depth and has the unit of shear stress. When Tr = Te, the erosion rate is the floc erosion rate. Parchure and Mehta (1985) determined the coefficient E{ and a in Equation (3.69) by measuring the time history of the suspended sediment concentration and the known bed density variation in laboratory. They found that a = 18.4mN1/2 and Ef = 5 x 10'O cm2min1 for kaolinite in tape water, and a = 17.2mN1/2 and Ef = 1.4 x 10'g crn2minI for kaolinite in salt water. The above models for the erosion rate were developed for a steady unidirectional current. Under an oscillatory wave action, Maa (1986) found that the erosion rate can be calculated by a similar expression as that for the dense bed (Equation (3.68)), E = E. + ' 1 1 (3.70) 2* 1e TV where E, is the erosion rate constant under the wave action and rTb.. is the maximum bottom shears stress induced by the oscillatory flow in a wave cycle. Maa (1986) determined the erosion rate constant E,. and the critical shear stress Ter for Cedar Key mud. He obtained E, = 1.2 x 103g cmn2min1, which is about one order of magnitude higher than the erosion rate constant (E. = 9.7 x 10sg cm2min1) obtained by Villarent and Paulic (1986) under steady current over mud collected at the same place. Maa (1986) also determined the critical shear stresses for erosion under wave action for kaolinite beds. Critical shear stresses were found to be smaller than those determined by Parchure and Mehta (1985) under steady current conditions. For example, a r7 of 0.25 Nm 2 under steady current action can become only 0.03 Nm2 due to a wave action. In considering the relative merits of the different models for computing the erosion rate, the study done by Lavelle et al. (1984) should be mentioned. As shown in Figure 3.3, a summary by Lavelle et al. (1984) of erosion versus bottom stress for several locations indicated that a power law relationship exists between E and rTb, but no universal constants exists for the relationship. Thus, no a priori judgement can 77 be made as to which formula for the erosion rate is best. The best erosion rate for a particular application must therefore be selected on an ad hoc basis. Field validation of the erosion model must be conducted for each new sediment transport study. It should be clear that most of the erosion models are descriptive models which were developed with very simple concepts and calibrated with limited laboratory data and numerous ad hoc empirical assumptions. Although most of the erosion models did appear to reproduce the particular laboratory data with which the models were calibrated, applications of most models to field conditions have not been proven to be successful. Only a few of the models have been applied to the field successfully. In a recent study (Sheng et al., 1989a), field data and a watercolumn sediment mixing model were used to aid the determination of the erosion rate, which was found to be an order of magnitude higher than the laboratory determined erosion rate using the same sediment. This fact and other uncertainties associated with laboratory erosion studies strongly suggest the need to develop erosion rate from comprehensive field data, rather than laboratory data which represent very different hydrodynamic and sedimentary conditions. Even in laboratory experiments, there is a need to obtain more comprehensive data, particularly within the bottom 30cm of the water column and the top 10cm of the sediment column. In high concentration (c > l000mg/I) and low flow situations, a sharp density gradient can develop near the bottom to cause the formation of a lutocline, which is similar to the thermocline or pynocline in the water column. Although the simulation of lutocline has been attempted (Ross, 1988; Wolanski et al., 1988), the paucity of data near the bottom has limited the modeling to merely qualitative calibration. 3.4.4 Fluidization and Consolidation The erosion process is related to the sediment processes in the bed, including the fluidization process and the consolidation process. The fluidization is a process by which the bed is fluidized due to the excess pore stress induced by the flow movement 4 2 1 4 3 ' .6 Y 10 9 6 V V 7 V 0 .1 0 1 log1, T (dynes/cm2) Figure 3.3: Erosion rate versus bottom shear stress (after Lavelle et al., 1984; the numbers are for different sediments and are explained in Lavelle et al., 1984). in the water column (e.g., wave action). The opposite process of fluidization is the consolidation process, in which the bed is consolidated due to the overburden or selfweight. The fluidization process is one of the interactions between waves and the mud (other interactions between waves and mud include the wave energy dissipation due to the mud movement and the internal wave dynamics along watermud interface). It is a rheological response of fine sediments to the wave action. The basic mechanism of waveinduced fluidization is that the bed is shaken and water is pumped downward into the bed under highly oscillatory flows, so that internal pore pressure can be built up and the bed becomes fluidlike mud. The fluidization process has not been well studied due to the complexity of the rheological property of different sediments. Mallard and Dalrymple (1977) treated the mud as a pure elastic solid in studying the wavemud interaction. Voigt model (a linear viscoelastic model) was used by Hsiao and Shemdin (1980) and Maa and Mehta (1986) for muddy seabeds subjected to wave action. Mei and Liu (1987) used a Bingham model to simulate the response of muddy seabeds under long waves. Their model assumed that mud motion begins when the imposed shear stress is larger than the yield shear stress and stops when it is less than the yield stress. Their model results showed that mud motion is intermittent. But unfortunately, intermittent motions of mud were not observed in either the field and or laboratory (Chou et al., 1993). Maa (1986) studied the erosion process under wave action and obtained a rate expression for the erosion of fluidized mud. To study the waveinduced fluidization of a natural estuarine mud, Ross (1988) measured the pore and total pressures of Tampa Bay mud subjected to progressive, nonbreaking waves. By defining a small effective stress, he tracked the fluid mudbed interface as fluidization proceeds. His results showed that the the fluid mud thickness increased rapidly in the first 15 minutes and 80 continued to increase afterwards with a relative slow but constant rate. During the fluidization, the upper limit of the fluid mud, the socalled lutocline (Kirby and Parker, 1984), was found to increase only in the first 15 minutes and stayed almost stable later on. Ross' experiment showed the importance of waves during the fluidization process. However, since he was not be able to control the horizontal advective transport of the fluid mud, he could not obtain the steadystate fluid mudbed interface, which should be a function of the water depth and the wave parameters (wave height, wave period, or wave spectrum of wave group) for a given mud. Feng et al. (1992) also measured the porewater pressure under oscillatory loading in laboratory. The fluidization process measured by Feng et al. (1992) was subjected to a single progressive wave. They obtained some expressions for the fluidization rate for their experiment runs. However, no quantitative relations between fluidization rate and wave parameters or the final fluid mud thickness and the wave parameters were drawn. Tzang et al. (1992) recorded resonant amplification of pore pressure oscillation during their laboratory experiments of fluidization subjected to wave action. Mea surements by Clukey et al. (1983) showed that total or partial fluidization occurred in a soil trench below shallow water waves. Owing to the complexity of fluidization, most of the previous studies were limited in the laboratory experiment. Field measurements or model studies were barely done. Sakai et al. (1992) were among the few groups who measured the waveinduced transient porewater pressure in the field. Chou (1989) developed a fluidization model by assuming that the bed acts like an elastic solid at a low strain, a viscous fluid at a very high strain, and a viscoelastic material at an intermediate strain. His model predicted that the sediment fluidization depth increased with wave height, but did not show the effect of wave period on fluidization. 81 Like fluidization, the consolidation process is also a very complex one and has not been well studied. Alexis et at. (1992) did some in situ consolidation experiments for estuary mud. Some field and laboratory observations have been made by Toorman and Berlamont (1993) for dredged mud. Toorman found a hindered settling mode and a consolidation mode in all his experiments. Most of the consolidation models were developed from the point of view of soil mechanics, i.e., solve the void ratio described by Gibson's law. Alexis et al. (1992) proposed a new governing equation for the void ratio which takes into account the compressibility of fluid and grains, the history of material (erosion and deposition), and non pure Darcinian flow. Toorman and Berlamont (1992) developed a soil me chanics model and a hindered settling model to simulate the consolidation process he observed in experiments. His results showed that the soil mechanics model was more or less suitable for long term prediction and the hindered settling model only gives good result in the initial stage of consolidation. He found that the soil mechanics model requires a lot of experimental data for calibration and the hindered settling model requires longterm laboratory experiment data for optimal simulations. 3.4.5 Models for Multidimensional Sediment Transport in Lakes and Estuaries Although controlled laboratory experiments can be used to provide insight into certain sediment processes, sediment transport models developed from laboratory studies may not be applicable to field conditions due to scaling problems, both in terms of flow and sediment. Field data, on the other hand, are very difficult and ex pensive to obtain. Thus, mathematical models of fine sediment transport, in conjunc tion with field data, are necessary to provide quantitative understanding of sediment transport. Since the transport of fine sediment particles involves many complicated processes, there exists no fully validated and universal sediment transport model at the present time (Sheng, 1989). The difficulty of modeling 3D sediment transport in lakes and 82 estuaries is additionally compounded by the complicated flow conditions, the interac tions between sediment particles and flow field, the inhomogeneity of sediments over the entire bottom of the water body, the complicated geometry and bathymetry of the water body, and the boundary conditions at watersediment interface. In addition, 3D simulations of sediment transport require a lot of data for model calibration and validation, and data collection is very expensive. Due to the above problems, past comprehensive models for large scale sediment transport are rather scarce. Sheng and Lick (1979) developed the framework of a comprehensive sediment transport model, and validated it with data from the west ern basin of Lake Erie. Transport processes considered in their model were advection, turbulent diffusion, erosion, deposition, and gravitational settling. To quantify the transport process, they applied a 3D hydrodynamic model (Sheng et al., 1978) to compute the 3D lake circulation. To quantify the bottom stress acting on the lake bottom, they developed a wave model to compute the wave field and orbital veloc ity over the entire western basin of Lake Erie. Erosion and deposition rates were determined from laboratory experiments using actual lake sediments. Settling ve locity of sediments was also measured in the laboratory. Even though some of the individual process models contained some empiricism, Sheng and Lick were able to successfully simulate the large scale transport of fine sediments in the western basin of Lake Erie during a twoday episodic event. Ariathurai and Krone (1976) developed a 2D vertically integrated sediment transport model and applied it to study the dis persion of sediments in the Savannah River estuary. Despite the empiricism of their model, reasonable success was reported. Unfortunately, little data exist to validate their simulations. In addition, they considered deposition and erosion as a combined process and assumed that either net erosion or net deposition exists, depending on the bottom stress. Thus, for a certain range of bottom stress, their model predicts a zero net flux. Hayter and Pakala (1989) developed a 3D model to simulated the 83 contaminant transport in the Sampit River. They used the sediment removal rate obtained by Farley and Morel to calculate the settling velocity. This settling velocity may not be suitable for low sediment concentration in the water column (Sheng et al., 1991). No detailed comparison between model results and measured data was given by them. Lang et al. (1989) applied a 3D model developed by Markofsky et al. (1986) to simulate suspended sediment dynamics in Weser estuary in Germany. The model calculates eddy coefficients using a mixinglengthbased algebraic model developed by Lehfeldt and Bloss (1988). Equation (3.67) was used to calculated the erosion rate. Because of the explicit treatment in the hydrodynamics part of the model, the time step used was very small (about 15 seconds). Onishi et al. (1993) used a 3D model to simulate sediment transport processes during a storm event in New Bedford Harbor. They used Grant and Madsen's model to calculate bottom shear stresses due to wavecurrent interaction only during storm days and did not consider wave effects on sediment dynamics before and after the storm. The reason why they calculated bottom shear stresses in this way is because MrantMadsen's model failed to produce reasonable drag coefficients and bottom shear stresses before and after the storm. Obviously, this argument is questionable. 