UFLICOELTR/124
MUD BOTTOM EVOLUTION AT OPEN COASTS
by
HUGO N. RODRIGUEZ
Dissertation
2000
MUD BOTTOM EVOLUTION AT OPEN COASTS
By
HUGO N. RODRIGUEZ
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
2000
ACKNOWLEDGMENTS
First, I gratefully acknowledge my advisor, Dr. A.J. Mehta, who helped and guided
me for the past four years, giving me prompt answers to all my frequent inquiries.
I am also grateful to Dr. R.G. Dean, Dr. R.J. Thieke, Dr. Y.P. Sheng and Dr. U.H.
Kurzweg for their participation as supervisory committee members, and for the knowledge
given through their interesting lectures and discussions.
I am thankful to the staff of the Coastal and Oceanographic Engineering Program of
the Department of Civil and Coastal Engineering, especially the laboratory staff, for
assistance during the experimental phase of this research, to Helen Twedell for her
professional help in literature research, and to Becky Hudson for her help in dealing with the
university paperwork.
Thanks go to my fellow students in the program for their friendship, support, and
critiques and discussions during the past four years, especially to Kerry Anne Donohue, Peter
Seidle, Joel Melanson, Roberto Liotta, Justin Davis, Al Browder, and Jamie MacMahan.
Finally, I would like to express my heartfelt gratitude to my mother for her unfailing
love and for giving me all her support even during the hardest of times.
TABLE OF CONTENTS
page
ACKNOWLEDGMENTS ............................................ ii
LIST OF FIGURES ................................................... vi
LIST OF TABLES .................................................... xiii
LIST OF SYMBOLS .................................................. xv
ABSTRACT ........................................................ xxiii
CHAPTERS
1 INTRODUCTION ................................................ 1
1.1 Problem Statement ........................................... 1
1.2 Problem Approach .......................................... 4
1.3 Outline of Chapters ................ ..........................5
2 NEARSHORE HYDRODYNAMICS ................................ 7
2.1 Nearshore Wave Field ....................................... 7
2.1.1 TwoLayer Wave Propagation Theories ...................... 8
2.1.2 McPherson's TwoLayer Model ........................... 15
2.1.3 Wave Field Transformations ............................ 23
2.1.3.1 Wave ShoalingRefraction ........................ 23
2.1.3.2 Wave Breaking ................................ 26
2.2 W ave Modeling Scheme .................................. 29
2.3 Wave Field Simulations ..................................... 33
2.4 WaveInduced Current ...................................... 38
2.4.1 LonguetHiggins Approach ............................. 39
2.5 Tidal Current ............................................. 47
3 NEARSHORE SEDIMENT TRANSPORT ........................... 53
3.1 Sediment Transport Formulation .............................. 53
3.2 Sediment Mobilization ...................................... 56
3.2.1 Surf Zone ................ .......................... 56
3.2.2 Offshore Region ...................................... 63
3.3 Alongshore Sediment Transport ............................. 66
3.3.1 Surf Zone .......................................... 69
3.3.2 Offshore Region ...................................... 70
3.4 Depth of Closure ................ ......................... 70
3.5 CrossShore Sediment Transport .............................. 80
3.6 Sediment Transport M odel .................................. 83
4 CROSSSHORE SUSPENDED SEDIMENT DISTRIBUTION ............ 90
4.1 Introduction .............................................. 90
4.2 Laboratory Study .......................................... 90
4.2.1 Laboratory SetUp ................ .................. 90
4.2.2 Experimental Procedure and Conditions ................... 94
4.2.3 Data Analysis ....................................... 96
4.2.3.1 W ave Data ............... ................. 96
4.2.3.2 Erosion Rate ................................ 100
4.2.3.3 Suspended Sediment Concentration Data .......... 101
4.2.3.4 Settling Velocity ............................. 102
4.2.4 Test Results ....................................... 106
4.2.5 Alongshore Velocity ................................. 109
4.3 Concentration Profiles near Lian Island ......................... 114
5 FIELD EXAMPLES OF MUDDY COAST EVOLUTION ............... 122
5.1 Introduction ..............................................122
5.2 Lian Island Beach, China ................................. 123
5.2.1 Area Description .................................. 123
5.2.2 Hydrodynamic Characteristics .......................... 127
5.2.3 Sediment Properties ............................... 128
5.2.4 1983 M ud Transport ................................. 129
5.3 Mahin Beach, Nigeria ..................................... 130
5.3.1 Area Description .................................. 130
5.3.2 Hydrodynamic Characteristics .......................... 132
5.3.3 Sediment Properties ............... ................ 134
5.3.4 Recent Erosion History ............................... 135
5.4 Coast of Louisiana .. .....................................135
5.4.1 Area Description ....................................135
5.4.2 Hydrodynamic Characteristics ......................... 137
5.4.3 Sediment Properties ................................ 138
5.4.4 Beach Dynamics .................................... 139
5.5 Lake Ontario, Canada ................................... 139
5.5.1 Area Description .................................. 139
5.5.2 Hydrodynamic and Sedimentary Characteristics ............ 141
5.5.3 Beach Dynamics .................................... 141
6 MUD BOTTOM EVOLUTION SIMULATIONS ....................... 143
6.1 Introduction ............................................. 143
6.2 Test Scenarios ......................................... 144
6.2.1 Test Scenario 1 ................ .................... 144
6.2.2 Test Scenario 2 ................ .................... 149
6.2.3 Test Scenario 3 ..................................... 151
6.2.4 Test Scenario 4 ................ .................... 156
6.2.5 Test Scenario 5 .. ..................................160
6.3 Model Strengths and Limitations ............................. 166
7 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS FOR
FURTHER WORK ........................................... 169
7.1 Summary ............................................... 169
7.2 Conclusions ............................................. 170
7.3 Recommendations for Further Work ........................... 172
APPENDICES
A SURF ZONE EROSION RATE .................................. 173
A.1 Erosion Rate Expression ................................... 173
A.2 Laboratory Data .......................................... 174
A.3 Field Data ............................................... 177
A.4 Comparison of Data Sets ................ .................. 184
B LABORATORY AND FIELD SEDIMENT CONCENTRATION
DATA AND SIMULATIONS ................................. 187
C MUD BOTTOM EVOLUTION MODEL CODE ....................... 201
REFERENCES ...................................................... 220
BIOGRAPHICAL SKETCH ........................................... 231
LIST OF FIGURES
Figure page
2.1 Voigt viscoelastic model ................ ......................... 11
2.2 Maxwell viscoelastic model .................. ...................... 11
2.3 Rheological model used by Jiang (1993) .............................. 14
2.4 Definition sketch of a twolayer wave propagation system. ................ 15
2.5 Dimensionless wave damping coefficient as a function of dimensionless
frequency for different mud dynamic viscosities. ...................... 20
2.6 Dimensionless wave damping coefficient as a function of dimensionless
frequency for different mud elasticities. ........................... 20
2.7 Dimensionless wave number as a function of dimensionless frequency. ...... 21
2.8 Zero of a function using the secant method. ............................ 22
2.9 Shoreline related horizontal coordinate system. ......................... 24
2.10 Profile representation of the discretized grid. ......................... 31
2.11 Planform representation of the discretized grid. .......................... 31
2.12 Schematization of modeled domain ................................. 32
2.13 Wave rays over a rigid bottom (solid lines) and over a muddy bottom (dashed
lines).........................................35
2.14 Bathymetry at a sinusoidal embayment................................ 35
2.15 Wave rays over a rigid bottom (solid lines) and a muddy bottom (dashed
lines) at a sinusoidal embayment. ................................. 36
2.16 Energy flux ratio variation with deep water wave height and mud depth, h2 ... 37
2.17 Energy flux ratio variation with deep water wave height and bottom slope..... 38
2.18 Schematic sketch for an obliquely incident wave and associated velocity
components. ...................................................40
2.19 Theoretical forms of the nondimensional alongshore velocity Ox as a
function of y andr .............................................. 46
2.20 Tidal current with a lateral boundary layer. .............. .............. 48
2.21 Shear stress ratio R as a function of breaking wave height and tidal current
for awaveperiod T= 13 s. ................................... 51
2.22 Shear stress ratio R as a function of breaking wave height and tidal current
for a wave period T = 7 s. ....................................... 52
3.1 Sediment budget for a differential water element (control volume) .......... 54
3.2 Definitions related to the sedimentactive nearshore zone ................. 57
3.3 Schematic diagram for a breaking wave and associated velocities and impact
pressure ............. .................................... 58
3.4 a) Bed and fluid mud layers; b) Vertical forces on a particle of mass ni (after
Li and Mehta, 1997). ........................................... 73
3.5 Springdashpotmass mechanical analog of a mud bed. .................. 75
3.6 Definition of fluid mud thickness (after Mehta and Li, 1997). .............. 78
3.7 Comparison between water depth at initiation of liquefaction obtained using
equation (3.49) and the Fluidiza model. ............................ 79
3.8 Planform representation of the discretized grid ....................... 83
3.9 Profile representation of the onshore boundary condition .................. 87
3.10 Profile representation of the offshore boundary condition ................. 89
3.11 Plan view of the modeled domain with the lateral boundary conditions ....... 89
4.1 Plan view of the experimental basin .................................. 92
4.2 Elevation view of the experimental basin .............................. 92
4.3 W ave measurement locations ....................................... 97
4.4 Erosion rate function for AK sediment ............................... 101
1
4.5 Settling velocity, ws, as a function of suspended sediment concentration, C.
Solid line is the bestfit of data points. ws, is the free settling velocity
(after Jiang, 1999) ................ ........................... 105
4.6a Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No.8 ..................... 108
4.6b Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline, inside the surf zone, in test No.8 .... 108
4.7 Comparison of computed and measured alongshore current within the surf
zone inthe basin .............................................. 113
4.8 Measured versus calculated sediment discharge in the basin surf zone ....... 114
4.9 Comparison between computed crossshore (depthaveraged) suspended
sediment concentration distribution and data off Lian Island for a
measured wave height H5 =2m .................................. 118
4.10 Calculated crossshore distributions of suspended sediment transport for two
wave conditions offLian Island ............... ................. .121
4.11 Depthmean suspended sediment concentration as a function of crossshore
distance for wave height H5 = 2.0 m off Lian Island, China ............. 121
5.1 Global map showing locations of the muddy coasts of Jiangsu Province,
China, Ondo State, Nigeria, Louisiana, and Ontario, Canada ............ 122
5.2 Lian Island location in China ................ .................... 123
5.3 Mud discharge at Lian Island from dredging of navigation channel for
Lianyun Port (adapted from Yu et al., 1987) ....................... 124
5.4 Isopleths (thickness) of mud deposit at the end of discharge (September 5th,
1983) (after Yu et al., 1987) ..................................... 125
5.5 Isopleths of mud deposition two months after discharge (November 4th, 1983)
(after Yu et al., 1987) .......................................... 126
5.6 Bottom profile evolution in the NE direction from dredged material discharge
point at Lian Island (after Yu et al., 1987) ......................... 126
5.7 Bathymetric changes on the northern beach of Lian Island in 1984 and 1985
(after Yuetal., 1991) .......................................... 127
5.8 Mahin mud shore location in Nigeria ............................... 131
5.9 Shoreline recession in the vicinity of Awoye from 1972 to 1991 (adapted
from Eedy et al., 1994) ......................................... 133
5.10 Louisiana chenier plain and coastal sediments (after Wells, 1983) .......... 136
5.11 Surficial sediments, bathymetry and profile measurement station locations in
the study area (after Kemp, 1986) .............................. 137
5.12 Profile (below MSL) time series at Station 7 from February 1981 to May
1983 (after Kemp, 1986) ....................................... 138
5.13 Location of the till shores at Grimsby, Lake Ontario (after DavidsonArnott,
1986) ................................................... 140
5.14 Longterm glacial till profile change at Grimsby, Lake Ontario (after Coakley
et al., 1988) ..................................... ........... 142
6.1 Test scenario: straight shoreline forced by shorenormal waves ............ 144
6.2 Simulation of a laboratory AK mud profile (test A) considering sediment
conservation in the modeled domain .............................. 146
6.3 Simulation of a laboratory AK mud profile (test A) considering an offshore
sediment transport ............................................ 148
6.4 Simulation of a laboratory AK mud profile evolution (test B) considering an
offshore sediment transport ...................................... 149
6.5 Test simulation for Lake Ontario, Canada. The intermediate profiles are for
the years 1981 and 1982 ........................................ 151
6.6 Test scenario: shorenormal waves incident on a dredged channel. Awoye
shore erosion ................. ............................. 152
6.7 Simulated contour recession due to the sediment "sucking effect" of a
dredged channel .............................................. 155
6.8 Shoreline recession in the vicinity of Awoye in Nigeria from 1972 to 1991
(adapted from Eedy et al., 1994) ................................. 155
6.9 Test scenario: effect of seasonally changing waves on profile evolution. Mud
beach dynamics along coast of Louisiana ........................... 156
6.10 Profile evolution caused by seasonally varying wave climate .............. 158
6.11 Profile evolution caused by seasonally varying wave climate. More active
shoreward region ........................................... 159
6.12 Measured profiles at Station 7 on the Gulf coast of Louisiana (after Kemp,
1986) ...................................................... 159
6.13 Test scenario: convex shoreline forced by oblique waves. Lian Island mud
disposal problem ............................................. 160
6.14 Bottom contour evolution of a sinusoidal deposit forced by normal waves ... 162
6.15 Bottom contour evolution of a sinusoidal deposit forced by oblique waves ... 163
6.16 Contour comparison after two months for deepwater wave directions 9 = 0
and 45 ...................................................... 164
6.17 Simulated change in bottom profile along the centerline of the modeled
domain ..................................................165
6.18 Changes in bottom profiles in the NE direction from the dredged material
discharge point at Lian Island (after Yu et al., 1987) ................... 165
6.19 Contour comparison after two months with and without an alongshore (weak)
current (deep water wave direction = 0) ......................... 166
A.1 Definition sketch for surf zone mean erosion rate determination ........... 174
A.2 Erosion rate inside the surf zone as a function of breaking wave height for
laboratory data using AK sediment ............................... 176
A.3 Erosion rate inside the surf zone as a function of breaking wave height for
laboratory data using undisturbed samples from Lake Erie .............. 176
A.4 Erosion rate inside the surf zone as a function of breaking wave height for
data from the chenier plain of Louisiana ............................ 178
A.5 Erosion rate inside the surf zone as a function of breaking wave height for
data from Grimsby, Ontario, Canada ............................. 178
A.6 Erosion rate inside the surf zone as a function of breaking wave height for the
northern shore of Lake Erie using data of Kamphuis (1986) ............. 182
A.7 Erosion rate inside the surf zone as a function of breaking wave height for the
northern shore of Lake Erie using data of Gelinas and Quigley (1973) .... 183
B.1 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 1 ..................... 189
B.2 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 2. ..................... 189
B.3 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 3. ................... 190
B.4 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 4 .................... 190
B.5 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 5. .................... 191
B.6 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 6. .................... 191
B.7 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 7. .................... 192
B.8 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 8. ..................... 192
B.9 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 9 ..................... 193
B.10 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 10 .................... 193
B. 11 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 11 .................... 194
B.12 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 12. ................... 194
B.13 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 13 .................... 195
B.14 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 14 .................... 195
B.15 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 15 .................... 196
B.16 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 16 .................... 196
B.17 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 17 .................... 197
B.18 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 18 .................... 197
B.19 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 19 .................... 198
B.20 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 20 .................... 198
B.21 Measured and computed depthaveraged concentration as a function of
offshore distance from the shoreline in test No. 21 .................... 199
B.22 Comparison between computed crossshore (depthaveraged) suspended
sediment concentration distribution and data off Lian Island for a wave
height H= 2.0 m ............................................. 199
B.23 Comparison between computed crossshore (depthaveraged) suspended
sediment concentration distribution and data off Lian Island for a wave
height H= 1.5 m ............................................. 200
B.24 Comparison between computed crossshore (depthaveraged) suspended
sediment concentration distribution and data off Lian Island for a wave
height ofH5 = 0.5m ........................................ 200
LIST OF TABLES
Table pag
4.1 Water depth in cm at the five wave measuring locations .................. 97
4.2 Wave period, height and surf zone width. ........................... .99
4.3 Values of the damping coefficient k ................................. 100
4.4 Erosion rate constants for AK sediment .............................. 101
4.5 Suspended sediment concentration data for test No 8. .................. 102
4.6 Parameters in equations (3.18) and (3.31) for laboratory experiments ....... 106
4.7 Alongshore velocity data and calculation of c ......................... 111
4.8 Measured suspended sediment concentration offLian Island for H5 = 2 m. .. 116
4.9 Parameters for equations (3.18) and (3.31) for simulating Lian Island data.... 119
6.1 Parameters for test scenario 1 ...................................... 147
6.2 Parameters for test scenario 2. ..................................... 150
6.3 Parameters for test scenario 3. ..................................... 153
6.4 Parameters for test scenario 4. ..................................... 157
6.5 Parameters for test scenario 5. ..................................... 161
A. 1 Breaking wave height and erosion rate data from laboratory tests with AK
sediment, and Lake Erie samples ................................ .177
A.2 Breaking wave height and erosion rate data from chenier plain in Louisiana,
and Grimsby, Ontario .......................................... 179
A.3 Value of the constants obtained in the regression analysis of the data sets .... 179
A.4 18961975 average breaking wave height and erosion rate data for the
northern shore of Lake Erie ...................................... 183
A.5 Breaking wave height and erosion rate data for Lake Erie (Gelinas and
Quigley, 1973) ...................... ......... .... ..........184
B.1 Concentration data from test Nos. 1 through 11 ........................ 186
B.2 Concentration data from test Nos.12 through 21 ........................ 187
B.3 Parameters for equations (3.18) and (3.31) for simulating laboratory data .... 187
B.4 Concentration data offLian Island for H, = 2.0, 1.5 and 0.5 m ............. 188
B.5 Parameters for equations (3.18) and (3.31) for simulating Lian Island data ... 188
LIST OF SYMBOLS
a1 = Initial surface oscillation amplitude
a2 = Initial interfacial oscillation amplitude
aij = Weighting coefficient in equation (3.61)
am = Amplitude of the wave (horizontal) displacement at the bottom
A = Constant in the nondimensional alongshore velocity equation (2.47);
Volume per unit beach width eroded within surf zone
by = Weighting coefficient in equation (3.61)
BI = Constant in the nondimensional alongshore velocity equation (2.47)
B2 = Constant in the nondimensional alongshore velocity equation (2.47)
Bq = Sediment discharge per unit width in equation (3.57)
c = Wave celerity
co = Deep water wave celerity
ce = Sedimentspecific constant in the erosion rate equation (3.22) outside the
surf zone
cf = Bed friction coefficient
Cg = Wave group celerity
cg = Wave group celerity vector
ci = Weighting coefficient in equation (3.61)
C = Suspended sediment concentration
C = Depthaveraged suspended sediment concentration
C = Nondimensional, depthaveraged suspended sediment concentration
Cb = Sediment concentration at the bottom
Ck = Complex constant introduce in equation (3.14)
C, = Complex constant introduce in equation (3.14)
di = Weighting coefficient in equation (3.61)
D = Littoral zone width considered in equation (2.50)
D, = Mean rate of energy dissipation per unit water volume
D1j = Equilibrium value of the mean rate of energy dissipation per unit water
volume
E = Wave energy per unit bottom area
Ei. = Strain tensor
E /' = Deviatoric component of the strain tensor
fc = Current friction factor
few = Wave plus current friction factor
f, = Wave friction factor
g = Acceleration due to gravity
g/ = Reduced gravity
g = Gravity acceleration vector
G = Shear modulus of elasticity
G2 = Shear modulus of elasticity of fluid mud
hb = Water depth at breaking
h1 = Water depth
h2 = Mud layer thickness
hic = Depth of closure
hi = Water depth at initiation of mud liquefaction
H = Wave height
Hb = Wave height at breaking
Hbc = Critical value of the breaking wave height below which there is no
measurable erosion
Hb,max = Highest breaking wave within a certain period of time
HC = Wave height at closure
He,t = Nonbreaking significant wave height exceeded 12 hours in t years
Ho = Deep water wave height
Hsaxt = Maximum significant deep water wave height during time t
i, I = Discretization index in the xdirection
ind = Independent coefficient in equation (3.61)
j, J = Discrectization index in the ydirection
k = Complex wave number (k,+i ki)
ki = Wave damping coefficient
ki = Nondimensional wave damping coefficient
ki = Profileaveraged wave damping coefficient
kout = Proportionality constant in equation (6.1)
kq = Nondimensional proportionality constant in equation (3.54)
kr = Wave number
r, = Wave vector
kx = Wave vector component in the xdirection
ky = Wave vector component in the ydirection
k, = Nikuradse roughness parameter
K = Diffusion coefficient
K = Depthaveraged diffusion coefficient
KQ = Transport rate parameter in equation (3.51)
Kq = Activity factor in equation (3.53)
K = Proportionality constant relating shoreline recession rate to wave power
Q = Complex constant introduce in equation (3.14)
Lb = Wave length at breaking
Lc = Force on a floe due to cohesion
Li = Upward force on a floc due to wave motion
Limax = Maximum upward force on a floe due to wave motion
m = Crossshore bottom slope
m = Particle mass
n = Manning's bed resistance coefficient
xvii
N = Constant in the eddy viscocity equation
N = Nondimensional transport number
N2 = Nondimensional transport number
N3 = Nondimensional transport number
p, = Pressure in the water layer
Pb = Wave power at breaking
pi = Impact pressure at breaking
pic = Critical value of the impact pressure at breaking
4(x,y) = Suspended sediment mass transport rate
qut = Offshore sediment flux in equation (6.1)
q = Unit mass sediment transport rate in the xdirection
qx = Nondimensional unit mass sediment transport in the xdirection
qx = Sediment discharge due to a weak current in equation (3.60)
qxt = Water flux per unit length in equation (2.51)
qy = Unit mass sediment transport rate in the ydirection
Qx = Water flux between the shoreline and the offshore boundary
p = Waveinduced cyclic pressure
pO = Amplitude of the waveinduced cyclic pressure
P2 = Pressure in the mud layer
R = Ratio of wave plus current shear stress to wave shear stress;
shoreline recession rate
Ri = Richardson number
s = Coordinate axis in the direction of wave propagation
S.j = Mean momentum flux
Si = Radiation stress
S j = Integrated Reynolds' stress
Sk = Complex constant in equation (3.14)
S, = Complex constant in equation (3.14)
t = Time
xviii
T = Wave period
Tet = Wave period associated with He,t
Tsm,,t = Wave period associated with Hsm,,t
Tij = Stress tensor
T = Deviatoric component of the stress tensor
S = Total velocity vector at the bottom
ui = Fluid particle velocity vector in the water layer
u, = Depthaveraged velocity in the water layer
Ulb = Depthaveraged velocity in the water layer at breaking
i2 = Fluid particle velocity vector in the mud layer
ii, = Water particle velocity vector at the breaking wave crest
ui = Incident water particle impact velocity vector at breaking
Us = Sediment velocity vector
Ui = Reflected water particle impact velocity vector at breaking
u, = Maximum orbital velocity at the bottom
uw = Alongshore wave velocity
ux = Velocity component in the xdirection
fix = Oscillatory component of ux
ux = Turbulent component of Ux
Uy = Velocity component in the ydirection
ly = Oscillatory component of Uy
uy = Turbulent component of uy
Ut = Tidal velocity
U, = Depthaveraged maximum orbital velocity
Ux = Depthaveraged velocity in the xdirection
Uxb = Depthaveraged velocity at breaking in the xdirection
Uy = Depthaveraged velocity in the ydirection
V = Nondimensional alongshore velocity
W1 = Velocity component in the water layer in the zdirection
w2 = Velocity component in the fluid mud layer in the zdirection
w, = Settling velocity
Wq = Equilibrium value of the crossshore distance in equation (3.53)
x = Horizontal axis in the alongshore direction
, = Space vector in the [x,y] coordinate system
y = Horizontal axis in the crossshore direction
Yb = Value ofy at breaking
y = Nondimensional horizontal coordinate in the crossshore direction
y = Value of y at closure
Yf = Offshore extent of fluid mud layer
Yf = Nondimensional offshore extent of fluid mud layer
z = Vertical coordinate
ze = Fluid mud thickness
a = Constant in the depthaveraged longshore current equation (2.43)
ac = Proportionality constant between wave crest velocity and wave celerity
a, = Proportionality constant between impact pressure and kinetic energy
a, = Proportionality constant between impact pressure and breaking wave
height
a, = Proportionality constant in the alongshore transport rate equation (3.37)
a, = Nondimensional wave diffusion constant
ag = Characteristic modifier of the gravity effect
y = Shear strain
y, = Constant in the nondimensional alongshore velocity equation (2.49)
Y2 = Constant in the nondimensional alongshore velocity equation (2.49)
r = Constant in the nondimensional alongshore velocity equation (2.49)
6 = Sitespecific constant in the surf zone erosion rate equation (3.11)
6 = Thickness of the wave boundary layer in water
Ah = Discrete step in water depth
At = Time step
Ax = Discrete step in the x direction
Ay = Discrete step in the y direction
E = Arbitrarily small value
Eb = Mass erosion rate inside the surf zone
Ebo = Proportionality constant (mass erosion rate when Hb = 2 Hb)
ED = Rate of energy dissipation per unit area
ED* = Target value of ED
EM = Erosion rate constant
Eo = Mass erosion rate outside the surf zone
max = Amplitude of vertical acceleration of a particle
rI = Free surface oscillatoryy) displacement
rI = Free surface setup
112 = Interface oscillatoryy) displacement
0 = Angle between the wave vector and the positive xaxis
9 = Angle between the wave vector and the negative yaxis
K = Wave breaking index
9 = Dynamic viscosity
1i = Dynamic viscosity of water
P2 = Dynamic viscosity of fluid mud
PR = Dynamic eddy viscosity
v = Kinematic viscosity of water
v2 = Kinematic viscosity of fluid mud
ve = Apparent complex kinematic viscosity of fluid mud
p, = Water density
p2 = Fluid mud density
PD = Sediment dry density
a = Angular wave frequency
T = Shear stress; weighting coefficient
1b = Mean bottom stress
Tb = Wave bed shear stress
tS = Mean surface stress
T = Bed shear strength
,cw = Wave plus current shear stress
, = Wave shear stress
= Velocity potential in the water layer
= Phase angle between the surface and interface profiles
Q = Scalar phase function
xxii
j
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor in Philosophy
MUD BOTTOM EVOLUTION AT OPEN COASTS
By
Hugo N. Rodriguez
May 2000
Chairperson: Ashish J. Mehta
Major Department: Civil and Coastal Engineering
Nearshore evolution of muddy coasts by waves and weakcurrent forcing has been
examined. A mud transport model, which takes into account relevant physical processes in
the open coast environment, is developed. The model simulates the timeevolution of the
bathymetry due to sea forcing.
The model couples a sediment transport submodel with a hydrodynamic submodel.
A finite difference scheme over a staggered spatial grid is used. The sediment transport sub
model solves for the bottom change in each grid cell for each timestep due to changes in
crossshore and alongshore sediment fluxes. The hydrodynamic submodel, which drives
sediment transport, includes wave damping by fluid mud, wave shoaling and refraction, as
well as an imposed weak current in the alongshore direction. A depth of closure relevant to
mud profile response to waves is introduced. Deep water waves and alongshore current are
considered as input for each time step.
xxiii
The model is used to explain mud profile changes due to a variety of erosional and
accretional conditions at muddy coasts. Several field examples of such conditions and the
resulting bottom changes are presented. It is shown that accurate simulation of these changes
is sensitive to the hydrodynamic and sedimentary boundary conditions. Also, whereas in
analogous models for sandy profile evolution conservation of sediment mass within the
modeled domain is commonly assumed, such an assumption does not always apply in cases
involving fine sediments. This is because of the low settling velocities and long travel
distances associated with suspended fine sediment. For these reasons, synoptic and
synchronous data on hydrodynamic forcing and sediment fluxes are critical for simulating
profile evolution at muddy coasts.
xxiv
CHAPTER 1
INTRODUCTION
1.1 Problem Statement
The coastal zone has been recognized as a valuable resource, and many human
activities naturally depend on a comprehensive understanding of its behavior. Most of the
world's human population resides near the sea shore, where land is often highly developed
and in need of protection from marine elements. Many marine related activities such as
recreation, fishing and transportation mainly occur in nearshore waters. In order to sustain
such activities, several types of engineering works have been built in the coastal zone. The
impact these activities and structures have on the shoreline is the subject of major
engineering concern and has largely contributed to the development of coastal engineering.
Much progress has been made in the understanding of the nearshore region; most
work, however, has been focused on noncohesive, especially sandy, coasts. In recent years
interest in muddy coasts has increased, and is becoming a subject of great concern. Coastal
wetlands, primarily composed of mud, are diminishing due to forcing by waves, tides and
sea level rise (Kemp, 1986; Asangwe, 1993; Lee and Mehta, 1994). Port and harbor
developments near or on muddy coasts need to address the problem of sedimentation and
convenient disposal of fine grained dredged material. Shortage of sand sources for beach
2
nourishment has raised interest in the use of fine sediment at offshore locations as a mean
of shoreline protection by the high wave attenuation associated with mud beds (Lee, 1995;
Lee and Mehta, 1997).
A quantitative capability to model the hydrodynamic and sedimentary processes in
the nearshore muddy environment is essential to properly address the above issues. Physical
and mathematical modeling are possible approaches, and due to substantial advances in
numerical techniques and computer capability in recent years, numerical modeling is an
appropriate choice. Also, physical modeling is limited by the significant scale effects
associated with fine, especially cohesive, sediments. Consequently, the aim of this study is
to develop a predictive tool for the nearshore evolution of muddy coasts by means of
numerical simulations of the relevant hydrodynamic and sedimentary processes in the muddy
coastal environment.
Necessary and appropriate simplifications of the complex hydrodynamicsediment
interactions and the bathymetry of the coastal region are required to formulate the problem
in mathematical terms. A characteristic feature of muddy coasts is the generation of a fluid
mud layer due to liquefaction of the upper stratum by freesurface waves. This highly viscous
fluid layer in turn affects the wave field by dissipating part of its energy (Li and Mehta,
1997). The fluid mud layer can also migrate or be resuspended due to the combined action
of wave and tide induced currents and turbulence, sometimes causing the material to be
transported out of bounds of the littoral zone (Lee and Mehta, 1997; Tarigan et al., 1997).
This "unbalanced" sediment budget of the nearshore region characteristic of muddy coasts
further complicates computations of transport and bottom evolution.
3
The model developed in this study is meant to simulate mud transport and associated
bathymetric changes that reflect phenomena found in the prototype environment. Different
regions of the world where muddy coast evolution has been studied and data reported are
used as illustrative cases.
The first case is the sea coast of Lian Island, off Lianyun Harbor, in the Chinese
province of Jiangsu. Due to harbor development in the region, a large amount of dredging
was done there, and the dredged mud was disposed in coastal waters of the island through
pumping. A monitoring campaign was established in order to measure the "mesoscale"
evolution of the bottom affected by waves and tides (Yu et al., 1987).
Another case is one of shoreline recession due to wave action on the muddy coast of
Ondo State, Nigeria. A 20 km strip of muddy coast, centered in the village of Awoye, has
experienced considerable erosion in recent years (Eedy et al., 1994).
An interesting case is the Louisiana coast of the Gulf of Mexico, where mud stream
from the Mississippi River is responsible for a long shoreline composed of fine grained
sediment. From this site wave, current and sediment data along with seasonal bathymetric
changes are available (Kemp, 1986).
The last case considered is shoreline recession at Grimsby, Lake Ontario. At this
lakeshore site, an overconsolidated till bluff could not be reconstructed after the highly
cohesive bottom sediment was eroded and transported away from the profile. From this site,
limited wave data and profiles changes are available.
1.2 Problem Approach
A numerical model is developed for the simulation of bathymetric evolution of
muddy coast forced by waves and a weak tidal current. The model is based on a two
dimensional, horizontal, finite difference scheme using a staggered spatial grid. The
independent coordinates are the alongshore distance and the water depth, while the cross
shore distance is allowed to change. This grid arrangement makes it feasible to directly
follow bathymetric changes of each contour. The alongshore distance is divided into equal
segments, and the bathymetry is represented by specifying the depth of contours.
The model consists of two submodels: the first computes the wavecurrent field and
the second uses a combination of analytical and empirical equations to compute sediment
transport and bottom profile changes. In the first submodel, waves are introduced by
computing the wave field inside the modeled domain based on wave climate specified in
deep water. Wave shoaling, refraction, and energy dissipation by fluid mud are taken into
account in the conservation of wave number and energy equations. These two equations are
discretized and solved to obtain the wave field inside the domain. Tideinduced alongshore
current is considered as input at each time step.
In the second submodel, the effect of the hydrodynamic forcing on the bathymetry
is used to determine profile changes at the new time step by means of the continuity of
sediment mass.
5
1.3 Outline of Chapters
This study is presented in the following order. Nearshore hydrodynamics is described
in Chapter 2, including a review of twolayer wave propagation theories, followed by a
description of the theory used in the present work. Wave transformations in the muddy
bottom environment, including shoalingrefraction and breaking along with relevant wave
modeling, are included. This chapter also provides descriptions of relevant waveinduced and
tideinduced currents in the nearshore region.
In Chapter 3, nearshore sediment transport is described. Sediment mobilization
expressions for the surf zone and the offshore regions are developed along with an
expression for alongshore sediment transport. Crossshore sediment transport is formulated,
and a depth of closure for the muddy coast is introduced along with a sediment transport
model for the nearshore region.
Chapter 4 presents an analysis of crossshore suspended sediment transport
distribution developed in Chapter 3. Laboratory experimental procedure used to validate the
proposed distribution is described and the experimental data analyzed. Field data from the
Lian Island coast are also used to corroborate the crossshore distribution of suspended
sediment. Finally, experiments carried out to test the alongshore, waveinduced, current
formulation are described and results discussed.
Available field data on the hydrodynamics, sediment characteristics, and nearshore
bottom evolution of the muddy coasts of Lian Island, China, Mahin Beach, Nigeria, southern
Louisiana, and Lake Ontario, Canada are reported in Chapter 5.
6
Chapter 6 gives the results of a series of simulations carried out using the model
developed in Chapter 3. These simulations are critiqued in the light of field observations at
the sites described in Chapter 5. In Chapter 7, conclusions from this study and
recommendations for future research are presented.
Appendix A presents laboratory and field data analyses to verify the
phenomenological relationship between breaking wave height and erosion rate introduced
in Chapter 3. Appendix B presents the complete experimental data analysis of the crossshore
suspended sediment transport distribution developed in Chapter 4. The model code (in
Fortran) is given in Appendix C .
CHAPTER 2
NEARSHORE HYDRODYNAMICS
In order to simulate open muddy coast evolution, modeling of nearshore
hydrodynamics must be thoroughly examined. This includes the wave field, waveinduced
currents and tidal currents. These four topics are explored in what follows.
2.1 Nearshore Wave Field
Wave action is introduced by computing the wave field within the physical domain
modeled based on wave climate specified at the open boundaries of the domain. Wave
transformations which in general must be taken into account include refraction, shoaling,
diffraction, damping and breaking. A characteristic feature of muddy coasts is the generation
of a nearbottom fluid mud layer due to liquefaction of the upperstratum of the bottom by
surface waves. This soft mud layer in turn affects the wave field by dissipating part of its
energy. Hence, an appropriate wave theory that takes into account the coupling between
water and the soft mud layer must be used in contrast to the more common onelayer (water)
theory for a rigid bottom.
2.1.1. TwoLayer Wave Propagation Theories
Classical wave theories assume an inviscid fluid over a rigid, nonporous bottom. In
reality, the wave field interacts with the bottom resulting in wave height attenuation due to
bottom friction, percolation or viscous damping within the sediment bed, as well as a change
in the wave length and celerity, which in turn causes refraction. To account for these effects
accurately, different theories have been applied to simulate the propagation of waves over
a deformable bottom. These theories vary in the type of bottom considered. The main bottom
types include: elastic bottom (e.g., Mallard and Dalrymple, 1977; Dawson, 1978), porous
bottom (e.g., Liu, 1973), poroelastic bottom (e.g., Yamamoto and Takahashi, 1985), and a
wide range of viscous, viscoelastic and viscoplastic bottoms (Mehta et al., 1994).
Concerning the above formulations one notes that the various mathematical models
for wave propagation that include a layer of mud underlying the water layer mostly differ
in the rheology of the mud layer, and some of them in the way they treat the water layer, or
in the wave theory used. The behavior of mud under a dynamic load is characteristically very
complex. Mud rheologic properties are complex as well. Nevertheless, as noted, several
simple theological models have been introduced to represent bottom mud behavior, namely
viscous fluid, viscoplastic and viscoelastic. Other formulations including elastic and
poroelastic behavior are less commonly used in simulating muddy coast dynamics. The
elastic model does not predict wave damping, an important feature of clayey and silty bottom
environments. In the presence of such sediments, since they can be considered to be largely
1
9
impermeable to water flow, most of the energy dissipation is due to internal friction, rather
than percolation. Hence the applicability of the poroelastic model is not apparent.
Gade (1958) was among the first investigators to introduce a twolayer fluid model
for wave propagation. He considered the bottom layer to be a viscous fluid with a greater
density than water, and the water layer to be inviscid. He made a further restriction limiting
the calculation to shallow water linear waves. Dalrymple and Liu (1978) extended Gade's
work considering a linear wave (without the shallow water restriction) propagating over two
viscous layers. They called this model the "complete model", and also introduced a boundary
layer approach that is analytically simple and has explicit solutions.
Mei and Liu (1987) and Liu and Mei (1989) used a viscoplastic material, also called
a Bingham plastic. A viscoplastic material is characterized by a yield strength, meaning that
it can sustain low shear stresses, and when the shear stress reaches the yield value it flows
as a viscous fluid. In a strict sense, however, mud suspensions do not behave as viscoplastics,
as they do not have a well defined yield strength and, indeed, tend to move under any applied
shear stress. As Parker and Kirby (1982) have pointed out, yield strengths obtained by many
researchers are really extrapolations to zero shear strain rate based on measurements at higher
strain rates. Another drawback is that the Bingham plastic constitutive equation is nonlinear,
which is difficult to solve analytically except for a few simple cases which do not include the
wave propagation problem. Thus, for example, both models presented by Mei and Liu (1987;
1989) are limited to shallow water and comparatively thin mud layers.
Viscoelastic models are more appropriate because they account for energy loss due
to mud motion and also include the elastic properties, which cannot be neglected at sediment
10
concentrations found in mud beds (Golden et al., 1982). Viscoelastic materials show
properties common to both an elastic (Hookean) solid and a viscous (Newtonian) fluid. Thus
the constitutive equation for a simply described viscoelastic material can be obtained by
linearly adding the two effects. For the elastic behavior the stressstrain relationship is
T = Gy (2.1)
where is the shear stress, y the strain and G the shear modulus of elasticity. Equation (2.1)
can also be represented by a mechanical model of a linear spring with a rigidity constant G.
A Newtonian fluid can be mechanically represented as a linear dashpot with a
damping coefficient [t, and having the constitutive equation
r = Ut (2.2)
where the dot represent a time derivative and t is the dynamic viscosity.
The mechanical behavior of most materials can be represented by finite combinations
of the above two basic units. The two simpler representative models are the Kelvin or Voigt
model representing a solid, and the Maxwell model representing a fluid. Based on the
analogy of a springdashpot system, the Voigt model results from a spring and a dashpot in
parallel (Figure 2.1). For this model the strain in both elements (i.e., spring and dashpot) is
the same, and the total applied stress is the sum of the stresses in each element. Accordingly,
the resulting constitutive equation is
T = Gy + p' (2.3)
G I
Figure 2.1 Voigt viscoelastic model.
G ... 
Figure 2.2 Maxwell viscoelastic model.
12
The Maxwell model is one in which the spring and dashpot are connected in series
(Figure 2.2). In this case, the total strain is the sum of the strain in each element resulting in
the following constitutive equation
.i + t = t (2.4)
MacPherson (1980) developed a twolayer wave attenuation model considering mud
to be a Voigt continuum and water an inviscid fluid. He reduced the linearized equations for
the viscoelastic (Voigt) medium to the NavierStokes equations by assuming a sinusoidal
disturbance and introducing a complex viscosity of the form
.G
Ve = V + (2.5)
po
where p is the density of the Voigt material, v = tP/p is the corresponding kinematic
viscosity and a is the angular frequency of the sinusoidal perturbation. Thus, the study of
waves over a Voigt medium was reduced to that of waves over a viscous fluid. For a zero
value of G the Voigt element reduces to a viscous fluid, for an infinite G it becomes a rigid
body, and for a zero [t it transforms into a purely elastic material.
PiedraCueva (1993) extended MacPherson's model by introducing a boundary layer
at the watermud interface. MacPherson's model does not predict a shear stress at the
interface because the viscosity of water is neglected. By introducing the water boundary layer
at the interface, PiedraCueva provided shear stress continuity.
13
Following MacPherson's approach, Maa (1986) derived a complex viscosity for the
Maxwell element
v = 1 + ii (2.6)
1 + (jo/G)2 G
The apparent viscosity, Ve, in this case in frequencydependent. It is not defined for G = 0
and [t = 0, and for an infinite G behaves as a viscous fluid. In order to determine which
model better describes the stressstrain relation for muds, Maa (1986) carried out
experiments to observe mud behavior under a constant wavemean shear stress. These tests
provided evidence that in general muds do not behave as Maxwell elements, and that the
Voigt model gives a better description of mud rheology.
Jiang (1993) introduced a modification to simple mud rheology using the Burgers
viscoelastic model instead of a Voigt of Maxwell element. The constitutive equation for the
Burgers model is
J2 Gi G2 I2 G1
S+ = + (2.7)
GI + G2 G1 + G2 G, + G2
which corresponds to the mechanical model shown in Figure 2.3. This model can be reduced
to the Voigt model when G, approaches infinity. Jiang also extended the wave propagation
solution to second order using the perturbation method.
For the present study, MacPherson's model with a Voigt rheology will be used. Given
the scope of this work the modification introduced by PiedraCueva (1993) is not considered
14
essential, because the main dissipation occurs in the mud layer even after introduction of a
boundary layer in the water layer. Indeed, the main modification introduced, namely
continuity of shear stress at the mudwater interface, is not highly relevant to a 2D model
as the one implemented in this work. Modeling mud as a Voigt medium is justified following
Maa's (1986) work, which showed that it is a reasonable approximation to real mud
rheology. It also makes it feasible to treat the mud layer as a viscous fluid as a specific
theological subcase. The main modification introduced by Jiang (1993), namely second
order extension, is also not essential within the scope of the present work. Jiang's first order
solution with G1 oo is the same as the one obtained by MacPherson (1980) using the Voigt
model.
Figure 2.3 Rheological model used by Jiang (1993).
2.12 MacPherson's TwoLayer Model
Following MacPherson (1980), the governing equations and boundary conditions for
a wave train propagating over a twolayer system and their solution are introduced in this
section. A small amplitude wave is considered to propagate over a water layer above a mud
layer, with respective finite depths h1 and h2. Cartesian coordinates (s,z) are considered with
coordinate s at the undisturbed water surface and positive in the direction of wave
propagation, and z positive upwards (Figure 2.4) The water layer is assumed to be inviscid
with a density pi, and the mud layer is considered to be viscoelastic with a complex
viscosity ve given by equation (2.5) and density p2>Pl. The free surface wave profile is
defined as TI, and the interface wave profile as r12. Also, uil and ~, are the fluid particle
velocity vectors and p, and p2 are the corresponding pressures in the water and mud layers,
respectively.
z
(s)
Water
h, Layer
Mud
h2 Layer
Figure 2.4 Definition sketch of a twolayer wave propagation system.
16
The problem reduces to solving the linearized NavierStokes equations and continuity
in the mud layer, and the Laplace equation for the velocity potential in the water layer.
The equation of motion in the mud layer is:
a 1
2 Vp2 + VeV2ig2 (2.8)
at P2
where g is the gravity acceleration vector.
Continuity in the mud layer is:
Voi = 0 (2.9)
Laplace equation in the water layer is:
V2I = 0 (2.10)
where (, is the velocity potential in water, and is related to the particle velocity by:
U1 = V4)l.
The relevant linearized boundary conditions are as follows:
a) The combined dynamic and kinematic boundary condition at the free surface, which yields
the conventional freesurface boundary condition for irrotational waves at z = 0:
a21, 4),
+ g 0 on z = 0 (2.11a)
8t2 dz
17
b) The kinematic boundary condition at the interface, which requires continuity of vertical
velocity at z = h:
8112
a = w2 = w1 on z = h, (2.11b)
at
c) The dynamic boundary condition at the interface, which requires that both normal and
shear stresses be equal at the interface at z = h :
Sa2 W a41
P2 a + P2g2 + 2P2Ve aZ Pi at + Plgr12
at 8z at
on z = hi (2.11c)
P2 Ve a' + a W) 0
8z as
d) The noflow condition at the impermeable bottom and noslip condition at
z = (h,+h2):
i,2 = 0 on z = (h, +h2) (2.11d)
An interfacial wave of the form
12(s) = a2ekiscos(krs ot) (2.12)
which represents a progressive wave of amplitude a2e kis, is introduced. The term a2 is real
and represents the amplitude at s = 0, ki is the damping coefficient, kr is the wave number
18
and o is the frequency, all three common to both layers. Solving equations (2.8) through
(2.10) along with boundary condition (2.11), the following dispersion equation is obtained
(MacPherson, 1980):
p(04 g 2k2)tanhkh
+ p2gk +
gktanhkh, o2
+(2k2Vi(1)2[ (2k2i/ve)(QCkC kSkS,)2k2 +
(2k2 io/v)( SkC, kCkS,)
Sk 3V3(2k2 io/ve) 2k(kCkC, SkS,) 0
4P2k3e 2k([ 0kC
2k(QSkCAkCkS,)
(2.13)
where k = kr +iki is the complex wave number and
(2.14a)
k2 iY
Ve
Ck = cosh kh2;
Sk = sinh kh2;
CV = cosh Vh2
Se = sinh Vh2
(2.14b)
Expressing the free surface profile in the form
iq1(s) = ale kicos(krs (t+ +W)
(2.15)
the ratio of the interfacial amplitude to the free surface amplitude as well as the phase angle
i; between the free surface and interface profiles, respectively, are
coshkhi (gk/o2)sinhkhl
a1 (2.16)
j = arg[coshkhl (gk/o2)sinhkhl]
This model provides physically reasonable solutions for wave attenuation and the free
surface displacement to interface displacement ratio, and has the advantage of conforming
to a relatively simple dispersion equation (2.13). Figure 2.5 shows the variation of the
dimensionless wave damping coefficient, kih1, as a function of the dimensionless wave
frequency, o/g/h, for two different mud viscosities and zero elasticity. The values used
for these simulations are: hi =0.16m, h2 =0.18m, p1 =0.001 Pa.s, p = 1,000kg/m3 and
p2 = 1,200 kg/m 3. It is seen that wave attenuation decreases at low and high frequencies, with
a peak value at the resonance frequency. The wavemean rate of energy dissipation is given
approximately by (Dean and Dalrymple, 1991)
ED = 2 dz (2.17)
Sf az
Hence, maximum dissipation occurs at the maximum velocity gradient, when the wave
induced boundary layer within the mud nearly equals the mud thickness. At high frequencies,
the viscous boundary layer and viscous dissipation are reduced. At low frequencies, shallow
water effect dominates, and the waveinduced velocity is relatively uniform in depth, thus
reducing the velocity gradient (Jiang and Mehta, 1996).
0.4 0.6 0.8 1
a/Vg/
1.2 1.4 1.6
Figure 2.5 Dimensionless wave damping coefficient as a function of
dimensionless frequency for different mud dynamic viscosities.
1
10
S2
" 10
0.4 0.6 0.8 1
o/A/g/h.
1.2 1.4 1.6
Figure 2.6 Dimensionless wave damping coefficient as a function of
dimensionless frequency for different mud elasticities.
2
10
4
10
0.
2
103 2
0.2
SG, = 0 Pa
"G; = 3000 Pa
IL2 10 Pa.s
g2 I Pa.s
21
Figure 2.6 shows the variation of the dimensionless wave damping coefficient, kihl,
as a function of the dimensionless wave frequency, o/Jg/h for a constant mud viscosity
(,2 = 10 Pa.s) and two different elasticities. Figure 2.7 shows the dimensionless wave
number as a function of frequency for the following values of the constants: h, = h2 = 1 m,
u, = 0.001 Pa.s, p, = 1,000 kg/m3, p2 = 1,200 kg/m3, I2 = 10 Pa.s and G2 = 0. The
wave numbers for a single inviscid layer of thicknesses h, and h, +h2 (and linear wave
theory) are also shown. It is seen that the wave number for the viscoelastic model is in
between the two, but even for a relatively high viscosity of L2 = 10 Pa.s is closer to the
wave number for a single layer with depth equal to the depth of the water layer plus mud.
This behavior shows that the fluid mud layer has only a minor effect on the wave number,
and that the main influence of the mud layer is on wave damping (Wehausen and Laitone,
1960).
0.55
5 k, using viscoelastic bottom
0.5 
S *: k, using linear wave theory with h = h 
0.45 : k, using linear wave theory with h = h, + h2 "
0.4 /
0.35 
0.3
0.25
0.2
0.15 
0.1 
0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Figure 2.7 Dimensionless wave number as a function of dimensionless
frequency.
22
The wave number k = kr + i ki is computed by means of the dispersion equation
(2.13). As this is a complex transcendental equation, it has no analytic solution and must be
solved numerically. The technique used here to find the complex zeroes of a function is the
complex secant method (Hamming, 1973), a generalization of the secant method in the real
domain used to compute the zeroes of a real function. Given two points x1 and x2 close to
a zero of a function, the point of intersection of the secant, 12, and the xaxis, x3, is closer
to the zero than the two previous points (Figure 2.8). The method is used iteratively until
f(xi) is less than a certain arbitrarily chosen small value e.
y=f(x)
X3 X1
Figure 2.8 Zero of a function using the secant method.
The zero zi of the dispersion equation (2.13), a complex function of the form f(z)=0,
was obtained by means of
f(zi1)
zi = zi (i Zi2) (2.18)
f(zi1) f(zi2)
which was used iteratively with initial values z1 and z2 equal to the wave number obtained
from the linear wave theory for h = h1 and h = hi + h2, respectively, and e = 1010.
2.1.3 Wave Field Transformations
In addition to wave damping due to energy dissipation in the viscoelastic layer, the
wave field undergoes transformations in the nearshore region as a consequence of changes
in water and mud layer depths, which require consideration. Descriptions of these
transformations in a twolayer system including refraction, shoaling and breaking follow.
2.1.3.1 Wave ShoalingRefraction
Shoaling and refraction are two important transformations that waves undergo while
propagating in shallow water. Due to the dispersive characteristic of water waves,
represented in this case by equation (2.13), they propagate with different celebrities in
different water and mud depths. As the wave train arrives in shallower water it slows down,
and due to conservation of energy (specifically wave energy flux), the wave height increases
or shoals. In addition, there is a change in the wave direction due to the different velocities
of propagation of the same wave front over different water and mud depths. This
phenomenon of refraction causes a change in wave height as well due to energy conservation
along the direction of wave propagation. Thus, shoaling and refraction are both interrelated
j
24
through energy conservation. In the case of a twolayer system, damping due to viscous
dissipation must be taken into account as well in the energy balance.
The water surface profile described by equation (2.15) corresponds to unidirectional
wave propagation. In the actual nearshore situation the wave front will propagate over the
two horizontal dimensions. If we locate the coordinate system (as previously defined) in such
a way that the xaxis is along the coastline positive to the right when facing water and the y
axis positive offshore (Figure 2.9), the surface wave amplitude can be represented as
S= a e ki X2' cos (2.19)
where 0 is a scalar phase function such that the wave front crest occurs at 0 = 2nx for
integral values ofn (Dean and Dalrymple, 1991).
y
0
e
x
Shoreline
Figure 2.9 Shoreline related horizontal coordinate
system.
The wave vector is defined as
kr = VO (2.20)
in the direction of wave propagation, with its modulus equal to the wave number, i.e.,
 kr I = kr. Its components in each direction are kx and k If we define 0 as the angle
between the positive xaxis and the wave vector, positive counterclockwise, the phase
function can be expressed as (Dean and Dalrymple, 1991)
S= kxx + kyy ot = krcos6x + krsin0y ot = krX at (2.21)
The more commonly used angle, 0, is the angle between the negative yaxis and the
wave vector, positive counterclockwise. With the wave vector kr defined by equation (2.20)
and the coordinate vector K = [x,y], by making use of the property that the curl of a
gradient is identical to zero, the following equation is obtained
Vxkr = 0 (2.22)
or, by substituting the wave vector by its components, the resultant expression for the steady
state conservation of wave equation is obtained as
S(krsine) a(krcos) = 0 (2.23)
ax ay
26
The other governing equation, which together with equation (2.23) determines the
nearshore wave field, is the conservation of energy
V(EEg) = ED (2.24)
where E is the wave energy per unit horizontal area, Cg the group velocity vector and ED the
rate of energy dissipation per unit horizontal area. The rate of energy dissipation, ED, is given
by (Lee, 1995)
k,
ED = Eo (2.25)
kr
obtained by substituting E = pgH2/8 and the exponential decay law for wave height,
H(y) = He ky, into equation (2.24).
2.1.3.2 Wave Breaking
At numerous muddy coasts and under normal (i.e., nonstorm) conditions, due to the
high degree of dissipation caused by the fluid mud layer, the wave field is largely damped
out without much breaking (Azam, 1998). However, especially under storm conditions,
significant breaking can occur in spite of damping. This breaking generates a sizeable surf
zone and coastal currents which intensify sediment transport.
The wave breaking problem and the conditions at breaking have been given
substantial attention by investigators for rigid (or sandy) shorelines. In a muddy environment
breaking occurs once shoaling overcomes energy dissipation by damping. As the water depth
decreases the wave height increases, the wave celerity decreases, and as a result of inertia the
27
wave crest speed becomes faster than the trough speed and rolls down onto the wave front
or jets into the trough.
Different wave breaking criteria have been proposed for rigid bottoms based on
different wave theories. One of the simplest and early breaking criterion was proposed by
McCowan (1894), which states that the limiting wave height is a fraction of the local water
depth hb
Hb = Khb (2.26)
This "similarity" condition is applicable to mild slopes, and the value of the constant K
typically varies in the range 0.6 0.9, with the most widely used value being 0.78. A slightly
different and more general criterion was proposed by Miche (1951), in which he relates the
wave steepness to water depth at breaking
27ith
Hb = 0.142Lbtanh b (2.27)
Lb
where Lb is the local wave length. For the shallow water condition, Miche's equation (2.27)
reduces to McCowan's form (2.26) with the value of K equal to 0.89.
Other, more complex, expressions introduce the beach slope or relate specifically to
the wave theory used (Laitone, 1962; Lenau, 1966; Goda, 1975). As noted, all of these wave
breaking criteria have been developed for rigid bottoms and singlelayer wave theories. No
comprehensive study is believed to exist, at present, that specifically addresses the wave
breaking problem relevant to a twolayer inviscidviscoelastic model. For the sake of
28
simplicity, and until further studies are made, the McCowan criterion will be adopted here.
This is based on the observation that fluid mud density is typically no more than
approximately 20% larger than salt water density. The value of the constant K will be chosen
as 0.90 instead of the more common value of 0.78, in order to take into account the
"increased depth' effect of fluid mud.
Breaking begins when the shorewardpropagating wave reaches the breaking
condition as described. Inside the surf zone, the wave height decreases until the wave reaches
the shoreline where it effectively vanishes. Various models of wave decay inside the surf
zone have been proposed. The simplest one ("similarity model") assumes that the ratio of
wave height to water depth inside the surf zone is the same as the one at breaking (Horikawa,
1998).
Field evidence shows that the surf zone decay seldom follows such a simple rule,
even though it is widely used; hence other models have been proposed that take into account
energy dissipation in the surf zone in a more rational manner. For example, Dally et al.
(1985) proposed a model in which the rate of energy dissipation per unit horizontal area is
proportional (with a proportionality constant Ke) to the difference between the local energy
flux, EC and a stable value of the flux, (ECg)s that is,
aEc K
Ec= [Ecg (Ecg)2 (2.28)
ax h 1
Azam (1998) conducted laboratory experiments on wave breaking over mud bottoms
using sediment from mudflats in Malaysia. He concluded that the model proposed by Dally
29
et al. (1985) can be adopted provided suitable values of its constants are selected. However,
Azam's work is limited to waves propagating normally on a uniform beach, thus reducing
the problem to one of transport in the twodimensional vertical plane. This condition does
not address the breaking of a wave propagating over a nonuniform bathymetry, or of waves
propagating at a nonnormal angle with respect to the shoreline. Here, for the sake of
simplicity of treatment, the similarity model (2.26) is adopted, namely one in which the ratio
of wave height to water depth is assumed to remain constant (and equal to 0.90).
2.2 Wave Modeling Scheme
The wave modeling scheme selected is a finite difference scheme which solves the
conservation of waves and conservation of energy equations (equation 2.23 and 2.24) over
the grid shown in Figures 2.10 and 2.11.
This grid has been used by Perlin and Dean (1983) for an "nline" model of wave
propagation in water over a sandy bottom. It is compatible with a variable y with fixed x and
depth h (in contrast to the more common fixed grid system in x and y with a variable depth
h). I and J are the discretized variables in x and y, respectively, and the dotted lines in Figure
2.11 define the boundaries between discrete intervals along the xaxis. Given that y is a
dependent variable, two adjacent cells with the same discrete variable J are displaced with
respect to each other because they have different y and Ay, and therefore the equations in the
x and y directions are decoupled. To overcome this decoupling a weighting coefficient, t,
30
is used to attenuate the high frequency components. Accordingly, the difference form of the
conservation of wave equation (2.23) is
n = cos 1 [t(krcosO)i_lj+. + (1 2r)(krcos6)ij1 +
ij =kr ij
(2.29)
+ (krcosO)i+i+, y ((krsinO)il (ksin0),il)]}1
+ :(krCS 0)ijl 2Ax
This equation is solved iteratively for 60n", and after each iteration the wave angle
values, 0n", are averaged with those obtained in the previous iteration, On, in order to
improve convergence. The iteration is carried out until the average change in the angles is
less than 0.02, a small value.
For conservation of energy, the discretization scheme uses the same weighting
coefficient T. Expressing energy as a function of wave height from the linear wave theory as
0n", the difference form of the conservation of energy equation (2.24) is
Hnl 1 [1(H2cgsin0)ilj+1 + (1 2r)(H2cgSinO)ijl +
S (Cg sin0)i
+ (H2CSin0i+J + ((H 2CCOS0)ilj (H 2CCOS0)ilj) (2.30)
+ Ay ED ij] 1/2
Equation (2.30) is also solved by an iterative technique and the obtained wave
heights, H n+, are averaged with the previously obtained values, H in order to improve
convergence. The iteration is ended when the average change in the height is less than 0.01
m.
h, (I,J) Ah
y (I, J mu d
Sh
Figure 2.10 Profile representation of the discretized grid.
yII
(J+1)" contour,
....... .................
: Q,(I,J+1) I I
h(I,) J" contdur
Q.(Ij) / (I+)
Figure 2.11 Planform representation of the discretized grid.
32
The input parameters for the wave model are the water and mud properties, p p2, P2
and G2, the deep water waves characteristics Ho, 00 and T, the alongshore grid spacing Ax,
the maximum alongshore index value, Iax, the contours step Ah and the total number of
contours Jmax. For the modeled domain (Figure 2.12), the offshore distance y(I,J) and the
mud depth h2(I,J) at every grid point are also given as inputs.
yA Offshore boundary
Jmax contour
0 0
I I I 1
I I II
I=1 / Shoreline I = I x
Shoreward boundary
Figure 2.12 Schematization of modeled domain.
The offshore boundary condition, the condition at the lateral boundary of wave
approach, i.e., at I = 1 if 37/2 < 00 < 2i or at I =Imax if r < 00 < 37/2, as well as an initial
condition for the entire domain are needed to solve equations (2.29) and (2.30). A rigid
bottom outside the domain is assumed, and straight and parallel contours there. With these
assumptions, the alongshore variations (i.e., in the xdirection) in equation (2.23) are nil
33
outside of the modeled domain and the wave angle at the boundaries can be obtained from
Snell's law
sin sin90
si (2.31)
c co
where c and co are the wave celerity at the boundary of the modeled domain and in deep
water, respectively.
The wave height at the offshore boundary and the lateral boundary of wave approach,
Hbound, is obtained from linear wave shoaling and refraction from deep water, i.e.,
Co cos90
H bound = H (2.32)
S2Cg cosO
The initial values of Oi(I,J) and Hi(I,J) are also obtained using equations (2.31) and
(2.32). The submodel output includes values of wave height and angle of propagation inside
the domain, H(I,J) and O(I,J), respectively.
2.3 Wave Field Simulations
Some results are presented here to demonstrate the applicability of the wave model.
This analysis focuses on the effect of mud on wave refraction, and compares the results with
wave refraction over a sandy (rigid) bottom.
34
Figure 2.13 shows a comparison between wave rays over a rigid bottom and over a
1m thick mud at the bottom. In both cases the beach slope m = 1/160, the deep water wave
angle 00 = 2400, height Ho = 2 m and period T = 8 s. The fluid mud layer has a density
P2 = 1,200 kg/m3, dynamic viscosity t2 = 1 Pa.s and elastic modulus G2 = 0. The
domain is a 3 km stretch of a straightandparallel contoured nearshore region modeled from
the shoreline (h1 = 0) down to the 5 m bathymetric contour 800 m from the shore. As
expected from the results shown in Figure 2.7, there is less refraction at a coast with fluid
mud layer than over a rigid bottom for the same water depth (h1). Figure 2.13 in conjunction
with Figure 2.7 implies that the fluid mud layer decreases the wave number, due to an
effective increase in "water" depth. The wave rays in Figure 2.13, representing the trajectory
of the wave front, are used for illustration instead of model output, i.e., wave height and
angle at the grid points, for clarity of comparison.
Another comparison is made using a sinusoidal nearshore bathymetry as shown in
Figure 2.14. This bathymetry is of the form
y = 40(20 J)cos 2 + 800 + 80 J (2.33)
3,000)
where x and y are in meters and J varies from 1 at the shoreline to 20 at the offshoremost
contour at h, = 5 m. The other contours have depths (in meters) given by h1 (J) = [(J1)/4];
J=2,... 19.
S600
a 500
% 400
300 1 1 1
i ii i iii i i .i ii i i i i
200. .
100
0
0 500 1000 1500 2000 2500 3000
Alongshore distance, x(m)
Figure 2.13 Wave rays over a rigid bottom (solid lines) and over a muddy
bottom (dashed lines).
2500
2000
1500
1000
500
0 500 1000 1500 2000
2500
3000
Alongshore distance, x(m)
Figure 2.14 Bathymetry at a sinusoidal embayment.
*.:$ * Ns ' , .

...... ........... .
zz .//
N .
N 
N "Shoreline
____ I  I 
The results obtained for the sinusoidal enbayment with deep water wave and mud
characteristics identical to the previous simulation are shown in Figure 2.15. The same
behavior as discussed for Figure 2.13 can be observed. Again the fluid mud layer effect on
the wave refraction angle is that of an effective increase in depth over a rigid bottom.
A series of simulations are next carried out to analyze the effect of mud layer on the
wave field. A straight and parallel bathymetry is used from a water depth hi = 0.1 m to 5 m.
The deep water wave characteristics are, 00 = 3000, Ho = 2m and T = 8 s. Also,
p2 = 1,200kg/m 3, L2 = 0.1 Pa.s and G2 = 0 are selected.
2500...
E 2000 
0 '/oo *///////////// i i i
1 500
a 1000 / / / i
/ t i i i
U 500 I
.I
0
0 500 1000 1500 2000 2500 3000
Alongshore distance, x(m)
Figure 2.15 Wave rays over a rigid bottom (solid lines) and a muddy bottom
(dashed lines) at a sinusoidal embayment.
37
Figure 2.16 illustrates the effect of the fluid mud layer on wave damping. As
expected, the thicker the fluid mud layer the more the damping or, as plotted, lower the
energy flux at breaking, (EC )b, with respect to the same energy flux in deep water, (EC )o.
It is also observed that the damping effect is greater at a lower incident deep water wave
height and tends to a constant value for larger waves, irrespective of the thickness of the fluid
mud layer.
1 1.5 2 2.5 3 3.5 4
Ho(m)
Figure 2.16 Energy flux ratio variation with deep water wave height and mud
depth, h2.
Figure 2.17 presents the effect of bottom slope on wave damping. Milder slopes tend
to have a larger effect on damping because the bottom effect is exerted over longer distances
for the same water depth at the offshore boundary of the modeled domain.
0.7 : Slope m = 1:80 : Slope m = 1:320
0.7
0.6
i" 0.5
0
0.4
S0.3
0.2
0.1
0
0 1 2 3 4 5 6
Ho (m)
Figure 2.17 Energy flux ratio variation with deep water wave height and
bottom slope.
2.4 WaveInduced Current
Currents of different origin and character are present in coastal waters, including
waveinduced currents, tidal currents, windinduced currents and shallow ocean currents.
Shepard and Inman (1950) made the distinction between nearshore currents, which are
directly associated with the action of waves, and coastal currents comprising the effects of
ocean circulation, tides and wind. Here, wave and tideinduced nearshore currents will be
further explored.
39
Waveinduced currents are concentrated within the surf zone and near the breaker
line. The associated nearshore current system comprises of alongshore and crossshore
currents. For the alongshore (or longshore) current several formulas have been proposed
(Inman and Quinn, 1952; Shadrin, 1961; Eagleson, 1965; Sato and Tanaka, 1966).
One of the earlier approaches to simulate the alongshore current was made by Putnam
et al. (1949) using the momentum conservation principle. Assuming a uniform and steady
current along the shoreline inside the surf zone, the mean alongshore current Ux was
computed by equating the alongshore wave momentum flux at the breaker line to water mass
transport associated with Ux. Following this analysis, LonguetHiggins developed an
approach that is described next.
2.4.1 LonguetHiggins Approach
LonguetHiggins (1953) and LonguetHiggins and Stewart (1964) introduced the
concept of radiation stress, defined as the excess momentum flux induced by wave motion.
In subsequent papers, LonguetHiggins (1970a, 1970b, 1972) applied the radiation stress
notion to analyze the nearshore current system.
The equations of conservation of momentum are considered, but instead of using
their differential forms, they are used after integration over the total instantaneous depth
hi + il1 (where r,1 is the water free surface displacement), and averaging in time over a
wave period.
The mean waveinduced velocities in the horizontal directions, x and y (Figure 3.18),
can be defined as (Thornton and Guza, 1989)
I I
Ux = 1 fuxdz;
hI + \h
1 h
U h+ = fu dz
1 1l 1h
where Tr is the mean deviation of l, from the still water level. The corresponding velocity
components are
ux = U + i + Ux
Uy = U + il + uy
y y
in which superscript denotes the oscillatory component and the prime denotes the turbulent
component.
U
x
Uy
Shoreline
Figure 2.18 Schematic sketch for an obliquely incident wave and
associated velocity components.
(2.34)
(2.35)
41
Integrating the momentum equation, the components of the mean horizontal
momentum per unit area in the x and y directions, respectively, are obtained as
au au aux as,
p1(hi +1) I +U x +U +
a t ax a yy ax
(2.36a)
Sxy 11i s b
+ Plg(hl+1) + xTx
yy ax
a +U +U au + xy
p (h + t ax Y a y ax +
(2.36b)
asyy a"1 bs
+ = pg(hl +i~) +7y'T
ay ay
In equations (2.36a, b), ts and zb are the wavemean surface and bottom stresses,
respectively, and the first term on the right hand side is the horizontal force per unit area due
to the mean water level slope. The mean momentum flux S.i (with ij = xx, yy or xy) is
defined as
Sij = ij + S' (2.37)
The contribution due to turbulent motion, Si, is the integrated Reynolds stress and can be
parameterized using the eddy viscosity. The wave motion contributions are the radiation
stresses, which are
Sx Pi1 '2dz + f dz
hi hi
1pg2
2p, gh2
Syy = Pif dz + fpdz pgh12
hI h,
'ni
Sxy = Pl f ixiydz
h,
(2.38a)
(2.3 8b)
(2.38c)
where p, is the dynamic pressure. Using the linear wave theory, the following expressions
are obtained (Dean and Dalrymple, 1991)
x E 2C Ccos2
xx 2 c
Y E 2 c sin2
" 2 c
+ 2c
C
1)]
+ 2 c
c
(2.39a)
(2.39b)
(2.39c)
c
S, = E sinO cosO
C
Equations (2.39) show that the wave momentum fluxes (radiation stresses) are
proportional to the wave energy. Therefore any change in the energy, such as by dissipation
due to breaking, will cause the radiation stresses to change. These changes must be balanced
by external forces so that the momentum equations (2.36) are satisfied. That is, inside the
surf zone forces are induced that in turn account for water mass transport in this zone.
In the crossshore direction, as there is a boundary at the shore, the balancing force
is the pressure gradient associated with the wave setup rlT, and, therefore, Ux is equal to
zero. In the case of an open and unbounded (in the alongshore direction) shore, there is no
adverse pressure gradient in the alongshore direction capable of balancing the radiation
stress. As a result, a steady, wavemean current occurs, which in turn generates a bottom
shear stress. This stress balances the radiation stress.
For a long and straight shoreline that can be considered to be infinite in extent, all
derivatives with respect to x must be zero, so that the balance of forces in the alongshore x
direction, when the coupling term is neglected, is
xy (2.40)
ay
Given the turbulent nature of flow in the surf zone, the wavemean bottom shear
stress (due to turbulent flow) is usually expressed by the quadratic law
=b = 1CfJUJU
(2.41)
where cf is the bed friction coefficient and ii is the total velocity at the bottom. For the case
of small angles of wave incidence, LonguetHiggins (1970) simplified and linearized the
expression for the bed shear stress to
b P Cf
S= pc UmUx (2.42)
47r
where um is the maximum bottom orbital velocity. Solving for Ux, the following expression
is obtained
5rr am sine0
U gh1 (2.43)
8 c, co
where a is a constant of order 0.4, m is the bottom slope, hi is the water depth, 90 is the
angle between the wave ray and the shorenormal in deep water, and co the wave celerity.
In the case of a bottom of uniform slope (m = constant), the water depth can be
expressed as h1 = my, thus obtaining an alongshore velocity that is linear with y. This
velocity is zero at the shoreline, has a maximum value at the breaker line, and drops to zero
again outside the surf zone. The maximum velocity at the breaker line, Uxb, is
Uxb = am  g sineb (2.44)
8 Cf
45
The introduction of Reynolds' stresses in equation (2.40) reduces Uxb in (2.44) by
a factor of 0.2 according to field data analyzed by LonguetHiggins (1970a and 1970b).
Adopting and eddy viscosity of the form
le = Npy ghl (2.45)
where N is a constant, using the linear wave theory and the following nondimensional
variables
Y. UU
9 U ,x U (2.46)
yb Uxb
where subscript b indicates the value at breaking, the nondimensional alongshore velocity
becomes
Ux = A + B1'1, 0<,<1
(2.47)
x = B2.Y2, 1 <
where A, B1, B2, y and y2 are constants depending on the parameter r
p mN
N' (2.48)
2 ac,
and are given by
1
A=
(1 5r/2)'
Y2 1
B, 2 A;
Y Y2
3 9 1
Y= 2 42 1 +
'2 4 16 F '
Y1 1
B, A
Y1 Y2
The case described by equation (2.43), in which Reynolds stresses have been
neglected, corresponds to r = 0, and increasing F implies increasing turbulent dissipation by
lateral friction (Figure 2.19). For increasing turbulent dissipation the maximum velocity
decreases and moves shoreward, and outside the surf zone the velocity distribution extends
further offshore.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
y
Figure 2.19 Theoretical forms of the nondimensional alongshore velocity
ix as a function of 9 and r.
(2.49)
2.5 Tidal Current
Tidal currents often flow along the shoreline at open coasts and are strongly related
to the local morphology. Thus, they are easily influenced by the sea bottom and the shoreline
(Horikawa, 1988). Given the noslip flow condition at the shoreline for the alongshore depth
averaged tidal current, for simplicity of treatment Jiang and Mehta (1998) assumed a
parabolic form of the tidal velocityUt as
Ut(y) = 2( ) y )2 (D) (2.50)
where y is the crossshore coordinate system and Ut(D) is a reference tidal velocity at a
distance D from the shoreline (Figure 2.20). Equation (2.50) satisfies the noslip condition
at the shoreline [Ut(0) = 0] and zero stress [i.e., dUt/dy = 0] at D.
In the present study, in order to satisfy flow continuity in the modeled domain, the
water flux between the shoreline and the offshore boundary, Qxt, at each crossshore section
is assumed constant. Similar to equation (2.50), the crossshore distribution of the water flux,
qx(y), is obtained as
qxt(y) = )2]qxt(D) (2.50a)
48
The integral of equation (2.51) from y=0 to y=D is equal to Qxt. Note that the
offshore extent of the lateral boundary layer D is assumed to be larger than the surf zone
width yb.
Figure 2.20 Tidal current with a lateral boundary
layer.
Only cases of weak tidal currents will be considered in this study. A weak current is
defined here as one which does not influence the resuspension process, but is capable of
advecting sediment suspended by waves. In order to evaluate the maximum tidal velocity that
can be considered to be "weak", the following analysis of the resuspension capacity of waves
and of waves plus current is carried out.
The ratio of wave plus current shear stress, Tw, to wave shear stress, w,, is used as
the relevant parameter
R = w (2.51)
w
The wave shear stress is calculated by introducing the wave friction factor f,
(Jonsson, 1966)
1 2
S= fwPlm (2.52)
w here um is the amplitude of the wave horizontal velocity at the bottom (given by the linear
wave theory). The value of the wave friction factor, fw, used is obtained from
1 + log 1 0.08 + log (2.53)
4 4w ks
where am is the amplitude of the horizontal waveorbital displacement at the bottom (given
by the linear wave theory) and ks is the Nikuradse bottom roughness parameter.
The wave plus current shear stress, tCw, is computed from
w 1 fwPi(Ut + u,)2 (2.54)
2
where a compounded wave plus current friction factor, fw, is introduced (Jonsson, 1966;
Bruun et al., 1978)
f Ut+ fwu,
f (2.55)
Ut + uw
with fc the current friction coefficient obtained from Manning's n by (Vanoni, 1975)
2
fC = 2g (2.56)
hi1/3
The tidal current, Ut, is assumed to be parallel to the shoreline and the velocity, u,, in
equations (2.54) and (2.55) is the alongshore wave velocity component
u, = Ux + umcos (2.57)
where U, is the alongshore velocity defined in Section 2.5.1.
A value of Manning's n = 0.020, an angle at breaking of 0 = x/4, and a bottom slope
of 1/1,000 are selected. Wave heights at breaking are considered to range from 0 to 1.5 m.
Figures 2.21 and 2.22 respectively show plots of R for different values of the
breaking wave height and tidal current for two different wave periods that cover the range
of typical prototype values. From this comparison it can be deduced that a tidal current of 0.2
~ 0.25 m/s inside the surf zone does not increase R to more than 1.2 within the
characteristically common ranges of wave height and period encountered at open coasts. A
value of R less than 1.2 will be considered to be acceptable for computational purposes
within the scope of this study. This is because the wave direction (r/4 radians with respect
51
to the shoreline) represents a relatively severe situation, and also because a 20%
enhancement in the shear stress can be accounted for merely by slightly increasing the wave
friction coefficient.
1.4 R= 1 R1.1 R= 1.2
1.2
R= 1.5
1
S0.8.
.E 0.6
0.4
R= 2
0.2
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Tidal current, U, (m/s)
Figure 2.21 Shear Stress ratio R as a function of breaking wave height and
tidal current for wave period T = 13 s.
To reiterate, the above analysis is meant to include a tidal current which can advect
sediment eroded by waves, but cannot influence the erosion process. At higher tidal current
speeds the erosion process becomes more complex because currents and waves differ in the
way in which they erode a finegrained bed (Mehta, 1988).
R =1.5
S0.8
.H 0.6
0.4
0.2 R2
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Tidal current, U, (m/s)
Figure 2.22 Shear stress ratio R as a function of breaking wave height and
tidal current for wave period T = 7s.
CHAPTER 3
NEARSHORE SEDIMENT TRANSPORT
3.1 Sediment Transport Formulation
The timeevolution of muddy coasts can be evaluated from the spatial distribution of
sediment transport and associated bottom changes. At these coasts, sediment transport may
be conveniently divided into two main components by direction, namely crossshore and
alongshore. Crossshore transport rate is dependent mainly on wave orbital motion and
streaming, while alongshore transport can be considered to be primarily due to waveinduced
alongshore current, and also due to tidal current.
Bathymetric contours change their position as a result of changes in the crossshore
and alongshore sediment transport. For a sediment budget, the basic equation required is the
conservation of sediment mass, i.e., sediment continuity. Considering an elemental volume
of sides dx and dy, and depth hl(x,y,t) (Figure 3.1), the mass of sediment per unit length into
the differential element in time increment dt in the x direction is qx(x,y)dydt, and in the y
direction q (x,y)dxdt, where qx and qx are the components of unit mass sediment transport
rate in the x and y directions, respectively, in units of mass/length/time. The corresponding
components out of the element in the x and y directions are qx(x + dx,y)dydt and
qy(x,y + dy)dxdt respectively, obtained by means of Taylor expansion for small dx and dy.
L
Thus, the net mass of sediment transported out of the element in time dt in the x direction
is
[qx(x + dx,y) qx(x,y)] dydt dxdydt
ox
(3.1)
qx(x+dx)
Figure 3.1 Sediment budget for a differential water element (control
volume).
In an analogous way, the net mass of sediment leaving the control volume in the y
direction over time dt is
[qy(x,y + dy) qy(x,y)] dxdt y dydxdt (3.2)
By
55
The net loss of sediment within the control volume over time dt is equal to the
decrease in the bottom level times the bottom area times the dry density of the sediment, PD,
i.e.,
dpDh1
[PDhi(x,y,t + dt) pDhl(x,y,t)] dxdy = a dxdydt (3.3)
at
Equating the net loss of sediment mass in the control volume to the net mass
transporting out of the element, the conservation of sediment equation is obtained as
ah, qx aq
PD x y (3.4)
at ax ay
As the model considered in the present work simulates changes in contour position,
y, instead of changes in depth, h1, it is convenient to transform equation (3.4) in order to
yield h, as the independent variable and y as the dependent one. After introducing the local
crossshore bottom slope m = ahl/ay, the crossshore position of the contour line change
is obtained from
By_ 1 aqx 1 aqy
y 1 (3.5)
at mpD ax pD ahl
In order to simulate contour changes by means of equation (3.5), the alongshore and
crossshore components of sediment transport rate, qx and qy, are required. To calculate
56
these components, modes of sediment mobilization must be evaluated, as well as the relation
between current velocity and sediment transport velocity.
3.2 Sediment Mobilization
According to the analysis given in Chapter 2, wave action is the only bottom
sediment mobilizing factor considered in this study. The effect of waves in muddy
environments is complex. Waves can liquefy the upper stratum of the muddy bottom, which
in turn can interact with the wave field itself, thus causing a measurable damping of the wave
height by viscoelastic dissipation in the liquefied (fluid) mud layer. During storms, despite
this large dissipation, significant breaking may also occur close to the shore, thus exposing
the nearshore bottom to wave attack. In the surf zone the development of a stable fluid mud
layer is hindered because there this layer tends to be mobile or entrained, and consolidated
sediment beneath is exposed to breaking wave attack. In the offshore region, seaward of the
surf zone where a stable fluid mud layer can be present, entrainment and deposition tend to
occur at the watermud interface. Figure 3.2 is a schematic description of the surf zone and
offshore zone relevant to what follows.
3.2.1 Surf Zone
Inside the surf zone, breaking waves cause an impact pressure over the bottom
sediment. This impact leads to bottom scour and sediment resuspension which can be related
to the wave characteristics at breaking (Yamanishi et al., 1998). Accordingly, one can
57
introduce the breaking criterion considered in Section 2.2.3.2, Hb = K hb, and consider that
the water particle velocity at the breaking wave crest Ui (Figure 3.3) is a fraction of the wave
celerity, c, at breaking
JI = a c (3.6)
where ac is a proportionality constant which depends on, among other factors, the breaker
type, i.e., plunging, spilling, or surging, hence on bottom slope, m. Next, assuming the
shallow water condition at breaking the wave celerity depends only on the water depth at
breaking (c = Jg^) and substituting in equation (3.6)
lcl = ca gh (3.7)
Figure 3.2 Definitions related to the sedimentactive nearshore zone.
58
Furthermore, the incident and reflected fluid particle impact velocities, iii and Uir
can be considered to be proportional to Uic, with the proportionality constant dependent on
bottom slope.
Figure 3.3 Schematic diagram for a breaking wave and associated velocities
and impact pressure.
From energy considerations at the impact location, it can be stated that the pressure
energy at impact is a fraction of the incident kinetic energy
2 2 2
Pi ui uc Pi uc
__ oc_ oc = = a.
Pig 2g 2g pig 2g
(3.8)
and taking into account equation (3.7) the impact pressure at wave breaking is obtained as
Pi = %p I~gHb
(3.9)
where ap is a proportionality constant which, in general, depends on the bottom slope m.
59
Following the shear stressbased model of, among others, Kandiah (1974), the
erosion rate at the bottom inside the surf zone, Eb, is considered to be proportional to a
power of the excess impact wave pressure
Eb (Pi Pic)a (3.10)
where pi is the critical value of the impact pressure below which no measurable erosion
occurs. Taking into account equation (3.9), a simple heuristic erosion rate expression inside
the surf zone is proposed
Eb = H bc (3.11)
where Eb is the mass erosion rate inside the surf zone, EbO is a proportionality constant
(equal to the value of Eb when Hb = 2Hbc), Hb is the breaking wave height, Hbc is the
critical value of the breaking wave height below which there is no measurable erosion, and
6 is a sitespecific constant. Note that the erosion rate is a function of the wave height at
breaking rather than the local wave height, and thus represents a mean erosion rate over the
surf zone. Equation (3.11) is similar to the stressbased equation for the rate of bottom
surface erosion (Lee and Mehta, 1994)
E = EM 1 (3.12)
I S
60
where EM is a rate constant, Tbis the bed shear stress, and Tsthe bed shear strength. For
consolidated bottoms with uniform properties over depth, under laboratory conditions, and
under a steady or quasisteady (e.g., tidal) current, 6 = 1 in equation (3.12) (Mehta, 1988).
For bottoms that exhibit an increase in density and shear strength with depth, 6 reduces to
0.5 (Parchure and Mehta, 1985; Lee and Mehta, 1994). Under oscillatory flows induced by
windgenerated waves, 6 ranges from 0.95 to 1.82 (Mehta, 1996). Both 8 and EM are
dependent on sediment properties, and thus are specific to each location. Appendix A
presents an analysis of surf zone erosion based on equation (3.11) using several field and
laboratory data sets.
The suspended sediment concentration distribution along the vertical zdirection can
be evaluated considering steady state conditions from the vertical settlingdiffusion equation
(Mehta and Li, 1997)
aC
K + wsC = 0 (3.13)
Sz
where K is the upward diffusion coefficient for sediment mass, ws is the sediment settling
velocity, and C is the suspended sediment concentration. The selection of equation (3.13) for
the present study implies that the suspended sediment concentration at any position along the
offshore profile is assumed to be due to a local balance of the vertical fluxes of sediment. In
other words, timedependent effects as well as advection are ignored.
The diffusion coefficient, K, for wave motion was derived by, among others, Hwang
and Wang (1982) under the assumption that K is proportional to the vertical excursion of the
61
water particle motion (length scale), and also proportional to the vertical velocity of the wave
motion (velocity scale)
sinh2krz
K = aH 2o (3.14)
2 sinh2krhi
where a, is a nondimensional wave diffusion constant. Assuming the shallow water
approximation to be a reasonable choice in the present case and averaging in the z direction,
the depthmean value of K is obtained as
w H2o
K = (3.15)
6
Choosing a constant value of ws, equation (3.13) can be solved for C, and averaging
in the water column the mean concentration value obtained is
SCbK l e,hK) (3.16)
hi ws
where Cb is the concentration at the bottom, which can be obtained from the erosion flux
boundary condition at the bottom. At z = 0, assuming sedimentary equilibrium, the net
resuspension flux is zero, i.e., w Cb = Eb or, from equation (3.11)
SbO
Cb = Hb Hb (3.17)
ws Hbc
Introducing the similarity model H = Khi and the bottom slope m, so that the water
depth is h1 = my, the final expression of the steady, depthaverage sediment concentration
inside the surf zone as a function of offshore distance becomes
SEOb6w aK2m H H Hb 6w 1
ECb= wK2m b b y exp  0
6w2 Hbc ) awK2m y
where yb is the width of the surf zone.
Equation (3.18) can be nondimensionalized using the following variables
:= C Y
E0b Hb Hbc Yb (3.19)
Ws Hbc
Then, introducing the following dimensionless transport number
awOK2myb
N1 = (3.20)
6wS
the nondimensional expression of the depthaverage sediment concentration becomes
C = N 1 e 0<.<1 (3.21)
63
3.2.2 Offshore Region
Outside the surf zone the wave height no longer varies linearly with the water depth,
and resuspension of fluidlike mud occurs by entrainment at the waterfluid mud interface.
The entrainment mechanism can be simulated by essentially the same approach as that used
in the study of thermally stratified and salinity stratified flows by considering that
entrainment is controlled by shear flow instability at the interface (Scarlatos and Mehta,
1990). Entrainment is usually quantified as a function of the Richardson number, Ri, which
represents the ratio of buoyancy to inertia forces, and, accordingly, the rate of entrainment
can be expressed as (Mehta and Srinivas, 1993)
Ce P2UI
E = e 2u (3.22)
Rin
where e is the erosion rate outside the surf zone, ce and n are sediment specific constants, u1
is the mean velocity in the water layer, and Ri is defined as
Ri = (P2P)g 6,
R = 2 (3.23)
Pi ut
where 6w is the thickness of the wave boundary layer in water. Cb in the offshore region can
be obtained from the boundary condition at the bottom as before
Ce P2UI
Cb = (3.24)
WsRi
j
64
The local wave height outside the surf zone, H, is obtained by shoaling and damping
of the deep water wave height Ho
H = Ho e ki(Yfy) CO (3.25)
\ 2c8
where Yf is the offshore extent of the fluid mud layer. Assuming shallow water for
simplicity, H is obtained from
H ek(Yfy) 1/4
H = kf) g (3.26)
At the seaward edge of the surf zone, the sediment concentrations obtained inside and outside
this zone must be equal
H0 b Hbc Pl U2n+1 (3.27)
EOb Hb Ce P2 (P2 )g Ulb
where ulb is the mean velocity in the water layer at breaking. Assuming shallow water
e 0b Hb Hbc (P2P) g w 22n+1
P2 Hbc Pl (gK Hb)(2n+l)/2
Substituting equations (3.23) and (3.28) in equation (3.24)
C E 0 Hb Hbc H2n+1
b W Hbc (Khl Hb)(2n+l)/2
Substituting equations (3.15) and (3.29) in equation (3.16), the offshore depthaverage
sediment concentration becomes
C = 60b Wb Hb H 2n+3 (6w,/ao)hl/H2 (3.30)
6w2 (K H,)(2n+l)/2 Hbc h)(2n+3)/2
where H is given by equation (3.26), and h1 = my. Thus, the final expression of C in the
offshore region as a function ofy is obtained as
SEOb w(Hog 1/4)2n+3 Hb Hbc e (2n+3)ki(Yf)
12w2 (2 Hb(2n)2 Hbc (my)(6n+9)/4
s b(3.31)
yb
1 exp[(12wm3/2/a,gl/2H2) y3/2 /e 2k(Yfy)]}
Introducing the same nondimensional variables defined in equation (3.19), the
following dimensionless transport numbers can be defined
a,(H g 1/4)2n+3 1
N2 w (3.32)
12ws (2oKHb)(2n+l)/2 (m yb)(6n+9)/4
12 w, (m Yb)32
N3 = (3.33)
aw g1/2 HO
With these definitions, the following nondimensional equation for the depthaveraged
sediment concentration outside the surf zone is obtained
S= (2n+3)kN(Yf)1 N3 32/e 2ii(' ) 1
2 (6n +9)/4
where ki = kiYb and Yf = Yf/ Yb. The deep water wave height, H0, and the wave height at
breaking, Hb, are related through equation (3.26) by assuming the breaking criteria
Hb = Khb. Hence
5/4 k,(YfHb,/mK) HO
Hb e (Kg)/4 (3.35)
3.3 Alongshore Sediment Transport
The mechanism of sediment movement along the open coast is complex, and
governing equations based on the fundamental physics of individual sediment particles have
not been fully established. Consequently, formulas for the alongshore sediment transport rate
have been based on macroscopic approaches. Owing to difficulties in measuring alongshore
67
sediment movement, early studies focused on the total alongshore transport rather than its
crossshore distribution, and were conducted for noncohesive, mainly sandy, sediments
(Savage, 1962; Kraus et al., 1982). Kraus (1987) has made one of the few attempts to
measure the crossshore and vertical distributions of alongshore transport rate of non
cohesive sediment (sand) using streamer type traps.
Many total alongshore transport expressions simply state that the total alongshore
sediment transport rate is proportional to the alongshore wave energy flux evaluated at the
line of wave breaking (Inman and Bagnold, 1963). Several studies have dealt with the
determination of the proportionality constant, and the relationship of this constant with the
governing physical parameters (Bodge, 1989).
Another concept used considers alongshore transport as the product of sediment
volume mobilized by some mechanism, and a sediment velocity which transports material
alongshore. This socalled energetics based approach" allows for a description of the cross
shore distribution of the alongshore transport, and has a more rational basis. In any event,
different approaches diverge in the mobilizing mechanism used and the degree of interaction
between the mobilizing and transporting forces (Bodge, 1989).
As previously stated, efforts to describe alongshore sediment transport have focused
on noncohesive, sandy sediments. Given this restriction, using the mobilizationtransport
approach an expression for muddy alongshore transport will be developed here.
The suspended sediment mass transport rate can be formally expressed as (Horikawa,
1988)
T TI(x,y)
(x,y) fdt f C(x,y,z,t)ii(x,y,z,t)dz (3.36)
0 hl(x,y)
in which T is the wave period, r1, is the water surface elevation, C is the concentration of
sediment in motion (in units of mass of dry sediment per unit volume of suspension), and Us
is the sediment velocity vector. From the momentum equation for an individual sediment
particle it is possible to relate the sediment velocity vector of a particle to the corresponding
fluid velocity (Hinze, 1975). From the velocity value for a particle the analytical functions
C(x,y,z,t) and Uis(x,y,z,t) needed in equation (3.36) can be found by summing the movement
of all particles.
Even though the above approach is a rational general technique, at present it is not
possible to evaluate individual particle velocities with an acceptable degree of accuracy.
Moreover, this particlebyparticle approach neglects particleparticle interactions (Horikawa,
1988). In the case of noncohesive sediments these interactions do not always affect
significantly the final result; however, for fine sediments they constitute one of the main
factors affecting C(x,y,z,t) and Uis(x,y,z,t). Fine, cohesive sediments do not move as
individual grains but coagulate and form flocs of different sizes composed of many
individual grains. The velocities of these flocs, i f(x,y,z,t), must be considered in equation
(3.36). These flocs continuously interact with each other under characteristically turbulent
flow conditions, changing sizes and making it difficult to track their movements. Efforts
have been made recently to track the movement of flocs (McAnally, 1999); however, at
69
present no practical results are available in order to apply equation (3.36) to cohesive
sediments.
As an alternative to a complete formulation dealing with equation (3.36), the
alongshore rate of transport, qx (x,y), can be evaluated in a semiempirical and approximate
way from a recognition of global patterns of movement. Accordingly, taking into account the
energetic based approach and equation (3.36), qx can be expressed as
qx = aq CUxh (3.37)
in which aq is a proportionality constant, C is the depthaveraged mass sediment
concentration and Ux is the depthaveraged alongshore velocity current. The proportionality
constant aq is introduced in order to empirically account for different influential factors
including the relationship between the sediment alongshore velocity and the fluid velocity,
flocculation, and the use of depthaveraged values as opposed to vertical distributions.
3.3.1 Surf Zone
Considering waves to be the sole sediment mobilizing factor, the depthaveraged
concentration C can be obtained from equation (3.18). Also, in the absence of a tidal current,
the waveinduced alongshore current, which transports sediment, can be evaluated using the
LonguetHiggins approach summarized in Section 2.5.1. Using the nondimensional
expressions given by equations (3.21) and (2.47), and considering that h1 = my, the non
dimensional expression of the alongshore sediment mass transport rate inside the surf zone
is
qx=aqNlm2 (1 e N'y (A +B1 ) 0<<1 (3.38)
where the alongshore mass sediment transport rate has been nondimensionalized as
=qx
ob HbHb (3.39)
Uxb bb
Ws H bc
3.3.2 Offshore Region
Using the same approach as in the previous section, the nondimensional expression
for the alongshore sediment transport rate outside the surf zone can be found by means of
equations (2.47) and (3.34) as
e (2n+3)ii(f9) N3 3/2 2 f)(3.40)
qx = qN2m B2 e e N3<
N(6n+5)/4y2
in which the sediment transport rate is given by equation (3.39).
3.4 Depth of Closure
The depth of closure, hi, (Figure 3.2), is the offshore depth at which no measurable
change in bottom elevation occurs (Birkemeier, 1985). For a sandy bottom, this depth is
71
related to the seaward limit of the alongshore transport. Therefore it is slightly larger than the
maximum breaker depth because alongshore transport characteristically decreases rapidly
outside the breaker line. The depth of closure is not an absolute seaward boundary, and
bottom changes tend to occur farther offshore (Hallermeier, 1978), although they are
typically much less important there than those landward of this depth. Thus, the concept of
a depth of closure is essentially empirical, and is based on the observation that with
increasing depth the vertical change of bottom decreases and, at a certain point, becomes
insignificant. It has been used in sediment budget estimates, and to provide an offshore
boundary in numerical modeling of coastal sediment transport (Kraus and Harikai, 1983).
The depth of closure is best determined by direct measurement of profile change at
the site. In most cases, however, profile data are not available, or are scarce, hence the depth
of closure is estimated from observationbased formulas.
Based only on crossshore conditions (i.e., neglecting the alongshore gradients) and
field data, Hallermeier (1978) proposed an analytical expression to estimate the depth of
closure from extreme wave conditions. In a generalized, timedependent form (Nicholls et
al., 1996) this expression is
2
hi, = 2.28He,t 68.5 (3.41)
gTe,
in which hict is the predicted depth of closure measured from mean low water level over t
years, He,t is the nonbreaking significant wave height which is exceeded 12 hours in t years,
and Tet is the associated wave period. Because the second term on the right hand side of
72
equation (3.41) has a relatively small contribution, the depth of closure is approximately
twice the significant wave height, leading Bruun and Schwartz (1985) to propose a simpler
rule that the limiting depth for active movement is 2Hb,max, where Hb,max is the height of the
highest breaking waves within a certain period of time. Down to this depth 90% of the total
profile movement can be accounted for over that time period, with the remainder extending
down to approximately 3.5 Hb,max
Hallermeier's equation has been applied with reasonable results (Kraus and Harikai,
1983), and its longterm validity has been evaluated (Nicholls et al., 1998). These results,
however, are applicable to sandy bottoms only. The analysis of Hallermeier (1978) is based
on the assumption that the only governing factor is a sediment entrainment parameter, which
is the ratio of the wave energy per unit sediment grain volume to the energy needed to raise
an immersed grain. Thus, this parameter relates the hydrodynamics to the settling of a grain.
Mud profile response to wave forcing is characteristically more complex, and the
depth of closure cannot be related only to particle or floc settling. Also, unlike sandy profiles,
no systematic studies addressing this issue appear to have been made. Nevertheless, an
analogous "operational" definition of closure depth can be introduced as the depth seaward
of which the influence of waves on the bottom can be neglected. Given the observation that
the "active profile" is often underlain by fluid mud (Lee and Mehta, 1997), this depth can
be taken as the depth where the bottom liquefaction potential of waves practically vanishes.
Fluid mud can be considered to be generated when the maximum upward force due
to wave motion, Limax, equals the sum of the downward forces, i.e., the submerged weight
of the bed floc, mpg/, and cohesion, Le (Figure 3.4)( Li and Mehta, 1997).
73
k Fluidized Layer
Bed
(a) (b)
Figure 3.4 a) Bed and fluid mud layers; b) Vertical forces on a particle
of mass mp (after Li and Mehta, 1997).
Although the bed matrix can withstand a certain degree of compression, it cannot
resist significant stretching forces. Thus, when at any elevation the maximum upward force
equals or just exceeds the sum of downward forces, the bed is critically stretched and fluid
mud is generated above that elevation. Note that g = g [p2 l] / p is the reduced gravity
and mp is the particle mass. Under these conditions, the criterion for fluid mud generation
can be formulated as (Li, 1996):
Limax = mpg +Le (3.42)
Then, by considering the resistive effect to be a modifier of submerged weight, equation
(3.41) can be reformulated as:
max = ag / (3.43)
where max is the amplitude of vertical acceleration of a particle, and a is a characteristic
modifier of the gravity effect. Equation (3.43), which is based on forces on a single particle
at a given elevation, can be considered to reflect the corresponding condition causing the
separation of the entire fluid mud layer from the bed at the same elevation.
For a known bed density, p2 and a, the fluid mud thickness, zc, can be calculated
from equation (3.43) provided max is determined by considering mud response to wave
loading. For a given wave, this response depends on the theological behavior of mud. In that
regard, the shear viscoelasticity of a mud bed can be represented by the (shear) Voigt model
(Maa, 1986), according to
T = 2GEi + 2E (3.44)
where Ti. = deviatoric component of stress the tensor Ti, Ei = deviatoric component of the
strain tensor Ej, subscripts i and j denote directions and ij denotes a second order tensor, G
= elastic modulus, [t = viscosity and the dot denotes derivative with respect to time.
Following the work of Li and Mehta (1997) a numerical model, Fluidiza, which
makes use of the theological properties of the bed, is used to determine Cmaxe Considering
the bed as an impermeable continuum represented as an equivalent mechanical spring
dashpotmass analog (Figure 3.5), the application of a constant normal stress, Tn, to the bed
surface is considered first, obtaining the displacement at the bed surface
0
m
Sk
Figure 3.5 Springdashpotmass mechanical analog of a mud bed.
2(0,t) = 1 ec) (3.45)
Equation (3.45) represents the extensional response of a Voigt solid to loading by Tn,
with a characteristic retardation time constant c/k (Barnes et al., 1989). The response at any
bed depth, z' measured from the bed surface, can be obtained as
2( ,t) = C2(0,t)( 1 (3.46)
h2^
