Effects of hydrodynamics and sediment transport processes on nutrient dynamics in shallow lakes and estuaries

MISSING IMAGE

Material Information

Title:
Effects of hydrodynamics and sediment transport processes on nutrient dynamics in shallow lakes and estuaries
Physical Description:
Book
Creator:
Chen, Xinjian, 1962-
Publication Date:

Record Information

Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 030582149
oclc - 31920333
System ID:
AA00022412:00001

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 3-D 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 parents-in-law 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 Short-Term Survey ....................... 139

7 MODEL APPLICATIONS TO LAKE OKEECHOBEE .......... 143
7.1 Modeling of the 1993 Storm Event .................... 143
7.2 1-D Simulations of 1989 Spring Survey .................. 150
7.3 3-D 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 3-D MODEL SYSTEM .............. 248

C NUMERICAL SOLUTION ALGORITHMS IN HYDRODYNAMIC AND
SEDIMENT TRANSPORT MODELS ...................... 250
C.1 3-D Hydrodynamics Model ........................ 250
C.2 1-D Hydrodynamics Model ........................ 252
C.3 3-D 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 1-D Simulation of Lake Okeechobee Phosphorus Dynamics ...... 307
H.2 3-D Simulation of Lake Okeechobee Phosphorus Dynamics ...... 308
H.3 1-D 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 two-film 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 u-velocities, v-velocities, 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 root-mean-square 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 three-day
event in the 1989 Spring Survey in Lake Okeechobee . . ... .. 165

7.14 SSC profiles in the lutocline simulation at Station C during a 3-day
event in Week 2 of the 1989 Spring Survey in Lake Okeechobee. 166

7.15 A Cartesian grid system for the 3-D 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 wave-induced bottom shear stress from the 1-D TKE
m odel . . . . . . . . . . . . . . . . .. 172

7.18 Comparison of model (1-D TKE) results of bottom shear stresses
induced by a single wave to those induced by group waves. . 174

7.19 Simulated u-velocities, v-velodcities, and suspended sediment con-
centrations at Station C of Lake Okeechobee form May 21 to June
16, 1989 . . . . . . . . . . . . ...... .. 176

7.20 Simulated u-velocities, v-velocities, and suspended sediment con-
centrations at Station D of Lake Okeechobee form May 21 to June
16, 1989 . . . . . . . . . . . . . . . .. 177

7.21 Simulated u-velocities, v-velocities, 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 3-D 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 3-D
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 u-velocity, v-velocity, 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 u-velocity, v-velocity, and sediment concentration
at Station A during the storm event of 1993 . . . . ... .. 220

9.6 Model results of u-velocity, v-velocity, 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 3-D model system of hydrodynamics, sediment
transport, and nutrient dynamics in lakes and estuaries . . 248

B.2 Flow chart of the 3-D 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 u-velocities, v-velocities, and suspended sediment con-
centrations at Station A from May 21 to June 16, 1989 . . 261

D.7 Measured u-velocities, v-velocities, and suspended sediment con-
centrations at Station B from May 21 to June 16, 1989 . . 262

D.8 Measured u-velocities, v-velocities, and suspended sediment con-
centrations at Station C from May 21 to June 16, 1989 . . 263

D.9 Measured u-velocities, v-velocities, and suspended sediment con-
centrations at Station D from May 21 to June 16, 1989 . . 264

D.10 Measured u-velocities, v-velocities, and suspended sediment con-
centrations at Station E from May 21 to June 16, 1989 . . 265

D. 1 Measured u-velocities, v-velocities, 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 u-velocities, v-velocities, 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 u-velocites and v-velocities 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 u-velocity, v-velocity, 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 u-velocity, v-velocity, 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 u-velocity, v-velocity, 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 1-D 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 3-D nutrient model . .197

H.I Results of some sensitivity tests for Lake Okeechobee phosphorus
dynamics . . . . . . . . . . . . . . ... .. 308

H.2 Model coefficients used in the 3-D simulations of Lake Okeechobee
phosphorus dynamics . . . . . . . . . . .... .. 309








H.3 Model coefficients used in the 1-D 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 1-D and

3-D 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 1-D and 3-D 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 wind-induced 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 real-time simulation

(10-15 min.) is generally much smaller than the time scale of desorption-adsorption

reactions of nutrients, equilibrium models for desorption-adsorption 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 60-90

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 real-time

nutrient variations caused by a short-term physical process such as an episodic storm

event because these models are unable to resolve such a short-term event.

Most of the previous water quality models (e.g., Chen and Orlob, 1975; Ambrose

et al., 1991) have been used to simulate long-term 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 water-sediment 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 short-term nutrient dynamics in shallow lakes and

estuaries because of the following reasons:

1. Wave boundary layer dynamics and wave-current 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 3-D 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 wind-induced 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 near-bed 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 set-up 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, set-up in shallow lakes is enhanced (Harleman and Shanahan,

1983). Because of the set-up, 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, wave-induced 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 wave-current 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 well-aerated because of the shallowness and the turbulent mixing induced

by wind (Somlyody, 1983). Because the depth-intergal 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 real-time 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 short-term and

long-term 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 3-D model and a 1-D model for hydrodynamics and sed-

iment transport. This chapter contains (1) assumptions used in the 3-D and 1-D

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 3-D model and a 1-D 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
four-week survey in the spring of 1989 and a three-day 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: 1-D simulations of sediment and phosphorus dynamics at L002 (a

station) during a three-day storm event in February 1993, 1-D simulations of sediment

dynamics at several stations during the 1989 Spring Survey, and 3-D simulations

of hydrodynamics, sediment transport, and phosphorus dynamics of the four-week

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 1-D model of hydrodynamics, sediment

transport, and nutrient dynamics to Tampa Bay. These applications include (1) a
one-day storm event observed in February 1992, (2) a one-day storm event observed

in March 1993, and (3) a one-week 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 9--z.:- 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, blue-green, 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) sorption-desorption 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 organic-bound phosphorus (Jensen and

Andersen, 1992).

The mineralization rate of DOP is usually modeled by a first-order 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). NaOH-Extractable P, commonly assumed to correspond to

aluminum- and iron-bound phosphorus, is the sediment form most readily available






16
to algae (Bostrom et al., 1988). HCl-Extractable P, corresponding to calcium-bound

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- Oc-3 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 right-hand 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 Sorption-Desorption Reactions

Sorption-desorption 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.
Sorption-desorption 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

log-log 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 sorption-desorption 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 sorption-desorption reactions are accomplished instantaneously.

The frequently used kinetics models for sorption-desorption 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 intraparticle-diffusion model to
simulate sorption-desorption 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 first-order 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 first-order 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 ammonium-nitrogen,

2. nitrification of ammonium-nitrogen in the presence of oxygen; in this process,

ammonium is first converted to nitrite and then converted to nitrate,

3. volatilization of ammonia-nitrogen,

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 first-order 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) zero-order equation,
(2) first-order 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 zero-order kinetics and

the first-order kinetics based on the substrate concentration of ammonium. For N,

much lower than H,, it approaches the first-order equation; and for N, much higher
than H,2, it approaches the zero-order equation. Since the H, value is typically

smaller than N, in natural aquatic systems, zero-order 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, 76-100% of nitrogen flux

across the sediment-water interface is N2, but only 5-70% 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 zeroth-order rate equation, a

first-order rate equation, or the standard Michaelis-Menten equation. Dawson and






26
Murphy (1972) found that the denitrification is of zeroth-order
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 first-order at low nitrate

nitrogen concentration. Reddy et al. (1978) also showed that under carbon-limiting
conditions the denitrification followed the first-order 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 half-saturation 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 free-living bacteria or blue-green algae.

It is an energy-intensive 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

sea-water compared to freshwater may inhibit the ability of nitrogen-fixers 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 so-called two-film model (J0rgensen,
1983c), which is often used in chemical engineering. The basic idea of a two-film

model assumes that a stagnant liquid film of thickness dL and a stagnant gas film of
thickness dG exist at the water-air 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


Air-Water N dG
Interface -.- -------- N; -..-- -- ------

Water N4




Figure 2.2: The two-film 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 half-saturation 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 half-saturation 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 blue-green

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

particle-size 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 half-saturation algal level for raptors and R-., is the maximum

snatching rate of raptors. uR_., is a function of temperature.

2.4.3 Algal Non-Predatory 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 20-50 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 non-predatory mortality of algae can be modeled by a first-order equation

(Gargas, 1976):








OCa
--- = a ,C, (2.48)

where K., is the mortality rate of algae and is temperature-dependent.

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) non-algal 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

blue-greens 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 blue-greens sank as detritus. The cell structure of the blue-green 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. Blue-green 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 high-density

storage products that act as ballast and affect cell buoyancy (Spencer and King, 1989).

Some blue-green algae that use one or more of these mechanisms are Anabaena sp.,

Aphanizomenon flos-aque, 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[Rl--T (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) = K20OT-20 (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 temperature-dependency 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 Beer-Lambert'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 non-point sources accounted

for only one-third of the nitrogen released from the sediments. Adding inputs from

the point and non-point 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 in-situ 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 wind-wave induced sediment resuspension in Lake Erie,

J = kp, P' w, (2.57)
p. p-Pw
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

one-dimensional model to simulate ammonia concentrations in an idealized estuary.

Two cases were studied: one with the real-time tidal current, and the other with

through-flow velocity due to freshwater inflow. He found that ammonia concentra-

tions predicted in the latter case (through-flow) differ by as much as 100 percent from

the real-time tidal averages even though the same biochemical model is used. This

indicates that the consideration of short-term hydrodynamics and sediment transport

processes is very important in a water quality model.

To study phosphorus dynamics in Lake Okeechobee with real-time hydrodynamic

and sediment transport processes, Dickinson et al. (1992) developed a 3-D phosphorus

model, which coupled with a 3-D model of hydrodynamics and sediment dynamics

developed by Sheng et al. (1991). The time step used in the 3-D phosphorus model

was the same as that in the 3-D 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, blue-green, 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 short-term event. Some models may fit field data well, but

are generally site-specific 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, wave-current 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 wind-driven. In a large lake, differential heating over

the lake can sometimes result in horizontal density-driven 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 well-mixed estuary, tidal range is generally large and the

river discharge is generally low. Therefore, a well-mixed 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 semi-diurnal 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 time-averaged tidal

currents are ebb-directed on one side of an estuary cross section and flood-directed
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 time-averaged 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.
Wind-driven 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. Wind-driven






46
circulation is particularly pronounced in shallow lakes and lagoons. Depending on the

wind and the water body, wind-driven circulation can be very important even when

gravitational and tidal circulations are important in an estuary. For example, in the
event of an oil-spill, 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, wind-driven circulation is very complicated. Sheng and Lick

(1979) developed a three-dimensional model for simulating wind-driven circulation in
Lake Erie. Sheng et. al (1991) developed Cartesian-gird and curvilinear grid three-

dimensional models to study wind-driven circulation in Lake Okeechobee. Their study

showed that wind-driven 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 wind-driven 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 three-dimensional free-surface

time-dependent circulation models.

Earlier free-surface three-dimensional models (e.g., Leenderste and Liu, 1975)
used explicit finite-difference methods to solve the three-dimensional 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 time-consuming. Leenderste and Liu (1975) used a time step of 10

seconds for modeling circulation in Chesapeake Bay and San Francisco Bay. The
three-dimensional 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 3-D free-surface time-

dependent models, Simons (1974) and Sheng et al. (1978) applied the time-splitting

technique to models of Lake Ontario and Lake Erie, respectively. The time-splitting

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 vertically-integrated

flow). Using the time-splitting 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 3-D 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 vertically-integrated 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 three-dimensional free-
surface models. One is the z-grid and the other is the a-grid. In z-grid 3-D models,

the vertical grid system is composed of fixed straight lines (layers), while in a-grid






48
3-D models, grids are vertically stretched through the following transformation:

0, = ((, Y,(3.1)
z-H(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 r-grid 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 z-grid does, because
the a-grid enables the communication between gird points in the shallow and deep

regions, which is not necessarily physically correct.

The z-grid was used in the three-dimensional models of Leenderste and Liu (1975),

Lang et al. (1989), and Casulli and Cheng (1992). The a-grid has been used in the
3-D 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 a-grid, Sheng et al. (1989c) modified Sheng's r-grid model to solve the
salinity equation on the z-plane (the horizontal gradient terms were calculated by first

finding the salinity values at the z-plane using linear interpolation from the salinity

values at a points and then solving the difference equation on the z-plane) in the 3-D

simulation of circulation in Chesapeake Bay. They were able to significantly reduce

the numerical error of the a-grid 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 3-D model (CH3D) originally developed by

Sheng to simulate the circulations in Lake Okeechobee and James River, respectively.
The present study uses a Cartisian-grid 3-D model (EHSM3D) originally developed

by Sheng to simulate the wind-driven 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 semi-empirical 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 one-dimensional time-dependent wave bound-

ary layer equation using a three-layer 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 theoretical-numerical 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
time-independent 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 Wave-Current 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 wave-induced near-bottom orbital

velocities are of the same order of magnitude as those due to wind-driven and/or

tide-driven currents in lakes and estuaries. For example, in Lake Okeechobee, the

maximum wind-driven current is on the order of 25 cm/sec, while the near-bottom

orbital current can be equally strong. However, as discussed in the previous sec-

tion, because of the high frequency oscillatory nature of wave-induced current, the

wave boundary layer is much thinner and wave-induced 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 one-dimension 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 y-directions. 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 wave-current 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 wave-current friction factor, i-, is the time-averaged 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 time-invariant. 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 time-invariant
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 time-dependent, 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.
Wave-current 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 one-equation 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 wave-current 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 wave-current 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 wave-current 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, i-i- (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 = U-u = 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
sub-grid 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 3-D model are dependent on the grid size used for a particular simulation. Generally,
constant values can be used in a three-dimensional simulation (e.g., Sheng et al., 1991;
Lang et al., 1989).
For vertical turbulent mixing, Lehfeldt and Bloss (1988) used a mixing-length-
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 non-stratified 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 mid-depth 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 mid-depth) 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 second-order 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 second-order 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

6-U.1 u ,v' 3i t -7" OW
at -uuj0- 4u-+g, (3.21)
O --jX-j = -+ ( axj )
, o9 au4
-~ u~k + vo-2-j ( a ')

A A
Aq-e 2au77

-X-~u~2 + V%
+,+,, a,-0--1 ... ao o9 (901O-01
+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 root-mean-
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), wave-current boundary layer (Sheng 1984), and vegetated boundary

layer (Sheng 1982).
Because of the complexity of Reynolds stress models, simplified second-order 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 second-order correlation equations are simplified to algebraic equa-
tions. The TKE closure model is similar to the Level-3 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

second-order correlations. The equilibrium closure model is similar to Mellor and
Yamada's Level-2 model. Instead of solving the differential transport equations for A
which consists completely of modeled terms, simplified second-order 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 second-order 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)
o-x 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+(bs-A)(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 root-mean-square 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 (bs-w)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 second-order 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 cut-off

limitation of turbulent macroscale based on 6q,2, the spread of turbulence determined
from the turbulent kinetic energy (q2) profile; and N is the Brunt-VWis.la frequency

defined as:
N ( agP1/2 (3.40)
\. p9z~
Although the simplified second-order 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 wind-induced mixing of salinity in Chesapeake Bay (Sheng et






63
al., 1989c; Johnson et a!., 1989) and wind-driven 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 = -2uuk-jO + 2g--- + 0.3 a-(qA ) -3 (3.41)
9t x 9Xk xk O9xk 4A
The vertically one-dimensional 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 second-order 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 second-order 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 sediment-laden boundary layer.
Mellor and Yamada (1982) generally favor the use of their Level-3 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 so-called "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 micro-shear 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 so-called '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 break-up 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 break-up term

and adjusting several model coefficients. Thus, the laboratory data of Burban et al.
were used primarily for calibration of their flocculation-breakup 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 box-model, work is needed to extend it to multi-dimensional models for simulat-

ing field situations. In another study, Farley and Morel (1986) performed laboratory
experiments on flocculation-induced 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 volume-averaged 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 site-specific parameter tuning. To apply their results to multi-dimensional

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)'4-65 (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 three-dimensional

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 non-Stokesian
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.2-1.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 location-dependent 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 10-9 /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 ad-hoc
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 non-cohesive 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 depth-averaged
current as follows:
E = E(U -I+ 1) (3.67)

where E, is the erosion rate constant, U is the depth-averaged 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 depth-averaged 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.4mN-1/2
and Ef = 5 x 10-'O cm-2min-1 for kaolinite in tape water, and a = 17.2mN-1/2 and

Ef = 1.4 x 10-'g crn-2min-I for kaolinite in salt water.

The above models for the erosion rate were developed for a steady uni-directional
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 10-3g cmn-2min-1, which is about
one order of magnitude higher than the erosion rate constant (E. = 9.7 x 10-sg
cm-2min-1) 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 Nm-2 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 water-column 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

self-weight.
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 water-mud interface). It

is a rheological response of fine sediments to the wave action. The basic mechanism

of wave-induced 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 fluid-like 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 wave-mud 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 wave-induced fluidization of

a natural estuarine mud, Ross (1988) measured the pore and total pressures of Tampa
Bay mud subjected to progressive, non-breaking waves. By defining a small effective

stress, he tracked the fluid mud-bed 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 so-called 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 steady-state fluid mud-bed 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 pore-water 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 wave-induced

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 long-term laboratory experiment data for optimal simulations.

3.4.5 Models for Multi-dimensional 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 3-D 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 water-sediment interface. In addition,

3-D 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 3-D hydrodynamic model (Sheng et al., 1978) to

compute the 3-D 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 two-day episodic event. Ariathurai and Krone (1976) developed

a 2-D 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 3-D 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 3-D model developed by Markofsky et

al. (1986) to simulate suspended sediment dynamics in Weser estuary in Germany.
The model calculates eddy coefficients using a mixing-length-based 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 3-D 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 wave-current 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 Mrant-Madsen's

model failed to produce reasonable drag coefficients and bottom shear stresses before
and after the storm. Obviously, this argument is questionable.