Mud bottom evolution at open coasts

UFLICOEL-TR/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 Two-Layer Wave Propagation Theories ...................... 8
2.1.2 McPherson's Two-Layer Model ........................... 15
2.1.3 Wave Field Transformations ............................ 23
2.1.3.1 Wave Shoaling-Refraction ........................ 23
2.1.3.2 Wave Breaking ................................ 26
2.2 W ave Modeling Scheme .................................. 29
2.3 Wave Field Simulations ..................................... 33
2.4 Wave-Induced Current ...................................... 38
2.4.1 Longuet-Higgins 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 Cross-Shore Sediment Transport .............................. 80
3.6 Sediment Transport M odel .................................. 83

4 CROSS-SHORE SUSPENDED SEDIMENT DISTRIBUTION ............ 90

4.1 Introduction .............................................. 90
4.2 Laboratory Study .......................................... 90
4.2.1 Laboratory Set-Up ................ .................. 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 two-layer 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 non-dimensional 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 sediment-active 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 Spring-dashpot-mass 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 best-fit of data points. ws, is the free settling velocity
(after Jiang, 1999) ................ ........................... 105

4.6a Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No.8 ..................... 108

4.6b Measured and computed depth-averaged 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 cross-shore (depth-averaged) suspended
sediment concentration distribution and data off Lian Island for a
measured wave height H5 =2m .................................. 118

4.10 Calculated cross-shore distributions of suspended sediment transport for two
wave conditions offLian Island ............... ................. .121

4.11 Depth-mean suspended sediment concentration as a function of cross-shore
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 Davidson-Arnott,
1986) ................................................... 140

5.14 Long-term glacial till profile change at Grimsby, Lake Ontario (after Coakley
et al., 1988) ..................................... ........... 142

6.1 Test scenario: straight shoreline forced by shore-normal 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: shore-normal 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 deep-water 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 depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 1 ..................... 189

B.2 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 2. ..................... 189

B.3 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 3. ................... 190

B.4 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 4 .................... 190

B.5 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 5. .................... 191

B.6 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 6. .................... 191

B.7 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 7. .................... 192

B.8 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 8. ..................... 192

B.9 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 9 ..................... 193

B.10 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 10 .................... 193

B. 11 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 11 .................... 194

B.12 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 12. ................... 194

B.13 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 13 .................... 195

B.14 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 14 .................... 195

B.15 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 15 .................... 196

B.16 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 16 .................... 196

B.17 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 17 .................... 197

B.18 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 18 .................... 197

B.19 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 19 .................... 198

B.20 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 20 .................... 198

B.21 Measured and computed depth-averaged concentration as a function of
offshore distance from the shoreline in test No. 21 .................... 199

B.22 Comparison between computed cross-shore (depth-averaged) suspended
sediment concentration distribution and data off Lian Island for a wave
height H= 2.0 m ............................................. 199

B.23 Comparison between computed cross-shore (depth-averaged) suspended
sediment concentration distribution and data off Lian Island for a wave
height H= 1.5 m ............................................. 200

B.24 Comparison between computed cross-shore (depth-averaged) 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 1896-1975 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 non-dimensional 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 non-dimensional alongshore velocity equation (2.47)
B2 = Constant in the non-dimensional alongshore velocity equation (2.47)
Bq = Sediment discharge per unit width in equation (3.57)
c = Wave celerity
co = Deep water wave celerity
ce = Sediment-specific 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 = Depth-averaged suspended sediment concentration
C = Non-dimensional, depth-averaged 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 = Non-breaking 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 x-direction
ind = Independent coefficient in equation (3.61)
j, J = Discrectization index in the y-direction
k = Complex wave number (k,+i ki)
ki = Wave damping coefficient
ki = Non-dimensional wave damping coefficient
ki = Profile-averaged wave damping coefficient
kout = Proportionality constant in equation (6.1)
kq = Non-dimensional proportionality constant in equation (3.54)
kr = Wave number
r, = Wave vector
kx = Wave vector component in the x-direction
ky = Wave vector component in the y-direction
k, = Nikuradse roughness parameter
K = Diffusion coefficient
K = Depth-averaged 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 = Cross-shore bottom slope
m = Particle mass
n = Manning's bed resistance coefficient

xvii

N = Constant in the eddy viscocity equation
N = Non-dimensional transport number
N2 = Non-dimensional transport number
N3 = Non-dimensional 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 x-direction

qx = Non-dimensional unit mass sediment transport in the x-direction
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 y-direction

Qx = Water flux between the shoreline and the offshore boundary
p = Wave-induced cyclic pressure

pO = Amplitude of the wave-induced 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, = Depth-averaged velocity in the water layer
Ulb = Depth-averaged 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 x-direction
fix = Oscillatory component of ux

ux = Turbulent component of Ux
Uy = Velocity component in the y-direction
ly = Oscillatory component of Uy

uy = Turbulent component of uy
Ut = Tidal velocity
U, = Depth-averaged maximum orbital velocity
Ux = Depth-averaged velocity in the x-direction

Uxb = Depth-averaged velocity at breaking in the x-direction
Uy = Depth-averaged velocity in the y-direction
V = Non-dimensional alongshore velocity

W1 = Velocity component in the water layer in the z-direction

w2 = Velocity component in the fluid mud layer in the z-direction
w, = Settling velocity
Wq = Equilibrium value of the cross-shore 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 cross-shore direction

Yb = Value ofy at breaking
y = Non-dimensional horizontal coordinate in the cross-shore direction
y = Value of y at closure
Yf = Offshore extent of fluid mud layer
Yf = Non-dimensional offshore extent of fluid mud layer
z = Vertical coordinate
ze = Fluid mud thickness
a = Constant in the depth-averaged 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, = Non-dimensional wave diffusion constant
ag = Characteristic modifier of the gravity effect

y = Shear strain
y, = Constant in the non-dimensional alongshore velocity equation (2.49)
Y2 = Constant in the non-dimensional alongshore velocity equation (2.49)
r = Constant in the non-dimensional alongshore velocity equation (2.49)
6 = Site-specific 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 x-axis
9 = Angle between the wave vector and the negative y-axis
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 weak-current 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 time-evolution of the

bathymetry due to sea forcing.

The model couples a sediment transport sub-model with a hydrodynamic sub-model.

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 time-step due to changes in

cross-shore and alongshore sediment fluxes. The hydrodynamic sub-model, 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 non-cohesive, 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 hydrodynamic-sediment

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 free-surface 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 "meso-scale"

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 sub-models: the first computes the wave-current field and

the second uses a combination of analytical and empirical equations to compute sediment

transport and bottom profile changes. In the first sub-model, 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. Tide-induced alongshore

current is considered as input at each time step.

In the second sub-model, 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 two-layer wave propagation theories, followed by a

description of the theory used in the present work. Wave transformations in the muddy

bottom environment, including shoaling-refraction and breaking along with relevant wave

modeling, are included. This chapter also provides descriptions of relevant wave-induced and

tide-induced 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. Cross-shore 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 cross-shore 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 cross-shore distribution of suspended

sediment. Finally, experiments carried out to test the alongshore, wave-induced, 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 cross-shore

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, wave-induced

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 near-bottom fluid mud layer due to liquefaction of the upper-stratum 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 one-layer (water)

theory for a rigid bottom.

2.1.1. Two-Layer Wave Propagation Theories

Classical wave theories assume an inviscid fluid over a rigid, non-porous 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 two-layer 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 non-linear,

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 stress-strain 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 spring-dashpot 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 two-layer 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 Navier-Stokes 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.

Piedra-Cueva (1993) extended MacPherson's model by introducing a boundary layer

at the water-mud 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, Piedra-Cueva 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 frequency-dependent. 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 stress-strain relation for muds, Maa (1986) carried out

experiments to observe mud behavior under a constant wave-mean 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 Piedra-Cueva (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 mud-water interface, is not highly relevant to a 2-D 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 sub-case. 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 Two-Layer Model

Following MacPherson (1980), the governing equations and boundary conditions for

a wave train propagating over a two-layer 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 two-layer wave propagation system.

16

The problem reduces to solving the linearized Navier-Stokes 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:

V-oi = 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 free-surface 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 no-flow condition at the impermeable bottom and no-slip condition at

z = -(h,+h2):

i,2 = 0 on z = -(h, +h2) (2.11d)

An interfacial wave of the form

12(s) = a2e-kiscos(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

+(2k2V-i(1)2[ (2k2-i/ve)(QCkC -kSkS,)-2k2 +
(2k2 io/v)( SkC, -kCkS,)

Sk 3V3(2k2 -io/ve) -2k(kCkC, SkS,) 0
-4P2k3e 2k([ 0-kC
2k(QSkCA-kCkS,)

(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 wave-mean 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 wave-induced 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

S-2
" 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, 1-2, and the x-axis, 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(zi-1)
zi = zi- (i Zi-2) (2.18)
f(zi1) f(zi-2)

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 = 10-10.

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 two-layer system including refraction, shoaling and breaking follow.

2.1.3.1 Wave Shoaling-Refraction

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 inter-related

j

24

through energy conservation. In the case of a two-layer 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 x-axis 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 x-axis and the wave vector, positive counterclockwise, the phase

function can be expressed as (Dean and Dalrymple, 1991)

S= kxx + kyy ot = krcos6x + krsin0y ot = kr-X at (2.21)

The more commonly used angle, 0, is the angle between the negative y-axis 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., non-storm) 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 single-layer wave theories. No

comprehensive study is believed to exist, at present, that specifically addresses the wave

breaking problem relevant to a two-layer inviscid-viscoelastic 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 shoreward-propagating 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 two-dimensional vertical plane. This condition does

not address the breaking of a wave propagating over a non-uniform bathymetry, or of waves

propagating at a non-normal 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 "n-line" 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 x-axis. 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 de-coupled. To overcome this de-coupling 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)i-lj+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 x-direction) 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 sub-model 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 straight-and-parallel 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 offshore-most

contour at h, = 5 m. The other contours have depths (in meters) given by h1 (J) = [(J-1)/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 Wave-Induced Current

Currents of different origin and character are present in coastal waters, including

wave-induced currents, tidal currents, wind-induced 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 tide-induced nearshore currents will be

further explored.

39

Wave-induced currents are concentrated within the surf zone and near the breaker

line. The associated nearshore current system comprises of alongshore and cross-shore

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, Longuet-Higgins developed an

approach that is described next.

2.4.1 Longuet-Higgins Approach

Longuet-Higgins (1953) and Longuet-Higgins and Stewart (1964) introduced the

concept of radiation stress, defined as the excess momentum flux induced by wave motion.

In subsequent papers, Longuet-Higgins (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 wave-induced 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) + x-Tx
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 wave-mean 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
2-p, 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 cross-shore 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, wave-mean 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 wave-mean 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, Longuet-Higgins (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 shore-normal 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 Longuet-Higgins (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 non-dimensional

variables

Y. UU
9 U ,x U (2.46)
yb Uxb

where subscript b indicates the value at breaking, the non-dimensional 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 non-dimensional 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 no-slip 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 cross-shore 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 no-slip 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 cross-shore section

is assumed constant. Similar to equation (2.50), the cross-shore 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 wave-orbital 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 fine-grained 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 time-evolution 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 cross-shore and

alongshore. Cross-shore transport rate is dependent mainly on wave orbital motion and

streaming, while alongshore transport can be considered to be primarily due to wave-induced

alongshore current, and also due to tidal current.

Bathymetric contours change their position as a result of changes in the cross-shore

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

cross-shore bottom slope m = ahl/ay, the cross-shore 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

cross-shore 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 water-mud 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 sediment-active 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 stress-based 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 site-specific 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 stress-based 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 quasi-steady (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

wind-generated 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 z-direction can

be evaluated considering steady state conditions from the vertical settling-diffusion 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, time-dependent 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 non-dimensional wave diffusion constant. Assuming the shallow water

approximation to be a reasonable choice in the present case and averaging in the z direction,

the depth-mean 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, depth-average 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 non-dimensionalized 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 non-dimensional expression of the depth-average 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 fluid-like mud occurs by entrainment at the water-fluid 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 = (P2-P)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(Yf-y) 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 e-k(Yf-y) 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 depth-average

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(Yf-y)]}

Introducing the same non-dimensional 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 non-dimensional equation for the depth-averaged

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,(Yf-Hb,/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

cross-shore distribution, and were conducted for non-cohesive, mainly sandy, sediments

(Savage, 1962; Kraus et al., 1982). Kraus (1987) has made one of the few attempts to

measure the cross-shore 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 so-called 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 non-cohesive, sandy sediments. Given this restriction, using the mobilization-transport

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 particle-by-particle approach neglects particle-particle interactions (Horikawa,

1988). In the case of non-cohesive 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 semi-empirical 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 depth-averaged mass sediment

concentration and Ux is the depth-averaged 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 depth-averaged values as opposed to vertical distributions.

3.3.1 Surf Zone

Considering waves to be the sole sediment mobilizing factor, the depth-averaged

concentration C can be obtained from equation (3.18). Also, in the absence of a tidal current,

the wave-induced alongshore current, which transports sediment, can be evaluated using the

Longuet-Higgins approach summarized in Section 2.5.1. Using the non-dimensional

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 non-dimensionalized as

=qx

ob Hb-Hb (3.39)
Uxb bb
Ws H bc

3.3.2 Offshore Region

Using the same approach as in the previous section, the non-dimensional 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(f-9) -N3 3/2 -2 f-)(3.40)
qx = qN2m B2 e e N3< N(6n+5)/4-y2

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 observation-based formulas.

Based only on cross-shore 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, time-dependent 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 non-breaking 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 long-term 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-

dashpot-mass 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 Spring-dashpot-mass 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^

MILLISECOND	CLASS.METHOD	MESSAGE
0	sobekcm_page_globals.constructor
0	sobekcm_page_globals.constructor	Application State validated or built
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.constructor	Navigation Object created from URI query string
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.display_item	Retrieving item or group information
0	sobekcm_page_globals.get_entire_collection_hierarchy	Retrieving hierarchy information
0	sobekcm_assistant.get_entire_collection_hierarchy
0	cached_data_manager.retrieve_item_aggregation
0	cached_data_manager.retrieve_item_aggregation	Found item aggregation on local cache
0	item_aggregation_builder.get_item_aggregation	Found 'all' item aggregation in cache
0	system.web.ui.page.page_load (ufdc.page_load)
0	sobekcm_page_globals.constructor.on_page_load
0	html_echo_mainwriter.add_style_references	Adding style references to HTML
0	html_echo_mainwriter.add_text_to_page	Reading the text from the file and echoing back to the output stream
44	html_echo_mainwriter.add_text_to_page	Finished reading and writing the file