Erosion of soft muds by waves


Material Information

Erosion of soft muds by waves
Physical Description:
Maa, P.-Y., 1948-
Publication Date:

Record Information

Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
oclc - 15300522
System ID:

Table of Contents
    Title Page
        Page i
        Page ii
        Page iii
    Table of Contents
        Page iv
        Page v
        Page vi
    List of Tables
        Page vii
    List of Figures
        Page viii
        Page ix
        Page x
        Page xi
        Page xii
        Page xiii
        Page xiv
    List of symbols
        Page xv
        Page xvi
        Page xvii
        Page xviii
        Page xix
        Page xx
        Page xxi
    Chapter 1. Introduction
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
    Chapter 2. Background study
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
    Chapter 3. Numerical one-dimensional Bingham fluid model
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
    Chapter 4. Two-dimensional multi-layered hydrodynamic model
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
    Chapter 5. Viscoelastic property tests
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
    Chapter 6. Erosion tests: Facility and procedure
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
    Chapter 7. Erosion tests: Results and discussion
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
        Page 173
        Page 174
        Page 175
        Page 176
        Page 177
        Page 178
        Page 179
        Page 180
        Page 181
        Page 182
        Page 183
        Page 184
        Page 185
        Page 186
        Page 187
        Page 188
        Page 189
    Chapter 8. Summary, conclusions and recommendations
        Page 190
        Page 191
        Page 192
        Page 193
        Page 194
        Page 195
        Page 196
        Page 197
        Page 198
        Page 199
    Appendix A. Boundary conditions for multi-layered model
        Page 200
        Page 201
        Page 202
        Page 203
        Page 204
        Page 205
        Page 206
        Page 207
    Appendix B. Comparison between prediction and measurement: Dynamic pressure and horizontal velocity
        Page 208
        Page 209
        Page 210
        Page 211
        Page 212
        Page 213
        Page 214
        Page 215
        Page 216
        Page 217
        Page 218
        Page 219
        Page 220
        Page 221
    Appendix C. Miscellaneous considerations for the wave
        Page 222
        Page 223
        Page 224
        Page 225
        Page 226
        Page 227
        Page 228
        Page 229
        Page 230
        Page 231
        Page 232
        Page 233
        Page 234
        Page 235
        Page 236
        Page 237
    Appendix D. Data on erosion tests
        Page 238
        Page 239
        Page 240
        Page 241
        Page 242
        Page 243
        Page 244
        Page 245
        Page 246
        Page 247
        Page 248
        Page 249
        Page 250
        Page 251
        Page 252
        Page 253
        Page 254
        Page 255
        Page 256
        Page 257
        Page 258
        Page 259
        Page 260
        Page 261
        Page 262
        Page 263
        Page 264
        Page 265
        Page 266
        Page 267
        Page 268
        Page 269
        Page 270
        Page 271
        Page 272
        Page 273
        Page 274
        Page 275
        Page 276
    Biographical sketch
        Page 277
        Page 278
        Page 279
        Page 280
Full Text








The author would like to express his sincerest appreciation to his

research advisor and supervisory committee chairman, Dr. A.J. Mehta,

Associate Professor of Civil Engineering and of Coastal and Oceano-

graphic Engineering, for his continuous guidance and encouragement

throughout this research. Appreciation is also extended for the valu-

able advice and suggestions of the supervisory committee cochairman,

Dr. B. A. Christensen, Professor of Civil Engineering, as well as the

guidance received from the other committee members: Dr. R. G. Dean,

Graduate Research Professor of Coastal and Oceanographic Engineering,

Dr. B. A. Benedict, Professor of Civil Engineering, and Dr. A. K.

Varma, Professor of Mathematics.

Sincere thanks also go to Dr. R. A. Dalrymple, Dr. J. T. Kirby,

Dr. M. C. McVay, Dr. L. E. Malvern and Mr. M. Ross for their sugges-

tions and help in this study.

Special thanks go to Mr. Vernon Sparkman, Mr. C. Broward, and

other staff of the Coastal Engineering Laboratory for their assistance

with the experiments performed during this research. In addition, the

author thanks Ms. L. Lehmann and Ms. H. Twedell of the Coastal Engine-

ering Archives for their assistance.

Gratitude is due to the U.S. Army Corps of Engineers Waterway

Experiment Station, Vicksburg, Mississippi, for their financial support

of this research under Contract DACW39-84-C-O013.

Finally, the author thanks his wife, Tai-Fang, for her love,

encouragement and patience, and his parents for their love and support.



ACKNOWLEDGEMENTS ............................................. ii

LIST OF TABLES ............................................... vii

LIST OF FIGURES .............................................. viii

LIST OF SYMBOLS .............................................. xv

ABSTRACT ..................................................... xx


1 INTRODUCTION ......................................... 1

1.1 Significance of Wave Erosion of Soft Mud ........ 1
1.2 Factors Characterizing the Wave Erosion Process. 3
1.3 Objectives and Scope of the Present Study ....... 7
1.4 Outline of Presentation ......................... 8

2 BACKGROUND STUDY ...................................... 10

2.1 Introduction .................................... 10
2.2 Wave Erosion/Resuspension ....................... 10
2.3 Constitutive Models for Soft Mud ................ 17
2.3.1 Nature of Clayey Soil .................... 18
2.3.2 Elastic Model ............................ 21
2.3.3 Poro-elastic Model with Coulomb Friction. 21
2.3.4 Viscous Fluid Model ...................... 22
2.3.5 Viscoplastic Model ....................... 26
2.3.6 Viscoelastic Model ....................... 32
2.4 Bed Shear Stress ................................ 36
2.5 Wave Diffusion Coefficient ...................... 39


3.1 Introduction .................................... 44
3.2 Problem Formulation ............................. 44
3.3 Results ......................................... 50
3.4 Conclusion ...................................... 55


4.1 Introduction .................................... 57
4.2 Formulation ..................................... 57
4.3 Solution Technique .............................. 64
4.4 Input Data ...................................... 66

4.5 Model Results ................................... 69
4.5.1 Velocity ................................. 70
4.5.2 Pressure ................................. 73
4.5.3 Shear Stress ............................. 75
4.5.4 Water Wave Decay ......................... 77
4.5.5 Interfacial Wave Amplitude ............... 78
4.6 Resonance Characteristics ....................... 80

5 VISCOELASTIC PROPERTY TESTS .......................... 87

5.1 Introduction .................................... 87
5.2 Relaxation Test ................................. 87
5.2.1 Apparatus ................................ 88
5.2.2 Procedure ................................ 92
5.2.3 Results .................................. 93
5.3 Tests for Material Constants of Voigt Element... 97
5.4 Tests for Rheological Properties of Bingham
Fluid Model ..................................... 102


6.1 Introduction .................................... 105
6.2 Wave Flume ...................................... 105
6.2.1 Wave Reflection .......................... 107
6.2.2 Wave Horizontal Velocity ................. 108
6.2.3 Wave Decay for Rigid Bed ................. 108
6.3 Sediment ........................................ 109
6.3.1 Kaolinite ................................ 109
6.3.2 Cedar Key Mud ............................ 110
6.4 Eroding Fluid ................................... 113
6.5 Instrumentation ................................. 116
6.5.1 Wave Gauges .............................. 116
6.5.2 Pressure Transducers ..................... 116
6.5.3 Light Meter .............................. 118
6.5.4 Data Acquisition System .................. 120
6.5.5 Suspended Sediment Sampler ............... 122
6.5.6 Bed Sampler .............................. 124
6.5.7 Current Meter ............................ 124
6.6 Procedure ....................................... 125
6.6.1 Process of the Digital Data .............. 126
6.6.2 Suspension Samples ....................... 128
6.6.3 Wave-averaged Mud Surface Elevation ...... 130


7.1 Introductory Note ............................... 131
7.2 Results ......................................... 131
7.2.1 Bed Density .............................. 131
7.2.2 Pressure Response ........................ 137
7.2.3 Instantaneous Sediment Concentration ..... 142
7.2.4 Time-averaged Sediment Concentration ..... 146
7.2.5 Wave-averaged Mud Surface Elevation ...... 151

7.3 Wave Erosion and Entrainment .................... 156
7.3.1 Erosion .................................. 156
7.3.2 Erosion Rate ............................. 165
7.3.3 Entrainment .............................. 174
7.4 Dye Diffusion Tests ............................. 182
7.5 Concluding Comments ............................. 185


8.1 Summary and Conclusions ......................... 190
8.2 Recommendations for Future Study ............... 197




EROSION EXPERIMENT ....................................... 222

D DATA ON EROSION TESTS .................................... 238

REFERENCES ................................................... 270

BIOGRAPHICAL SKETCH .......................................... 277



2-1 Diffusion Coefficients in Wave Flow ...................... 42

4-1 Input Thickness for Each Layer .......................... 68

4-2 Interfacial Shear Stresses and Mud Wave Amplitudes ...... 78

5-1 Vane Dimensions ......................................... 91

5-2 Coefficients of the Correlation Equations ............... 101

6-1 Sediment Properties ..................................... 113

6-2 Composition of Clay Fraction ............................ 113

6-3 Chemical Composition of the Eroding Fluid ............... 115

6-4 Eroding Fluid Properties ................................ 115

6-5 Wave Conditions for Erosion Tests ....................... 129

7-1 Average Erosion and Entrainment Rates ................... 163

7-2 Coefficients for the Erosion Function ................... 171

7-3 Measured Diffusion Coefficient in the Upper Layer ....... 179

7-4 Constants for Steady Flow Erosion Function .............. 187



1-1 Schematic Depiction of Response of Wave-Mud
Interaction ............................................. 5

2-1 A Plot of the Suspended Solids and Bed Shear Stress
for a Wave Resuspension Test ............................ 12

2-2 Measured Suspended Sediment Concentration Profiles
from a Wave Resuspension Test ........................... 14

2-3 Shear Stress-Strain Loop for Clay under Cyclic Load ..... 19

214 Energy Dissipation Ratio X versus Shear Strain .......... 19

2-5 Shear Stress-Strain Loops and Effects of Cyclic Load .... 20

2-6 A Plot of the Horizontal Velocity Amplitude Profiles
by Using the Two-Layered Viscous Fluid Model ............ 25

2-7 Relationship between Yield Strength and
Sediment Concentration .................................. 28

2-8 Relationship between Apparent Viscosity and
Sediment Concentration .................................. 29

2-9 System Diagram for Simple Viscoelastic Models and
Response under Constant Loading ......................... 33

2-10 Viscoelastic Behavior of Clay ........................... 33

2-11 Wave Friction Factor .................................... 38

3-1 System Diagram for the 1-D Numerical Bingham
Fluid Model ............................................. 45

3-2 Horizontal Velocity Profiles for Two Steady Flows.
(a) Two-Layerd Viscous Fluid; (b) Bingham Fluid with
Exponentially Increasing K and v ........................ 51

3-3 Comparison of the Wave Velocity Amplitude Profiles from
1-D model and Dalrymple and Liu's Model with P, -1000
kg/m3, p2 1160 kg/m3, v, -lxl06 and v2-2x10- m2/s... 52

3-4 Wave Velocity Profiles for the Water Bingham Fluid
System with T 3 sec. (a) H 2.5 cm and kim-
0.0011 m-1; (b) H 5 cm and kim:O.016 m-1 .............. 54

3-5 Shear Stress versus Shear Strain Rate for Kaolinite
Bed with Two-Day Consolidation Period ................... 56

3-6 Wave Velocity Profiles for Water-Bingham Fluid
System with Measured Yield Strength and Viscosity ........ 56

4-1 Schematic Figure for the Multi-Layered Model ............ 58

4-2 Layout of the Coefficient Matrix for the Multi-Layered
Model ................................................... 65

4-3 Comparison of Computed Bed Shear Stress (model and
Kamphuis) at the Interface to Determine the Eddy
Viscosity of Water ...................................... 68

4-4 Comparison of the Model Prediction and Measurement
for Run 1-2. (a) Velocity; (b) Pressure ................. 71

4-5 Comparison of the Model Prediction and Measurement
for Run 5-2. (a) Velocity; (b) Pressure ................. 72

4-6 An Example of the Model-Predicted Shear Stress Profile.. 76

4-7 Comparison of the Measured and Predicted Wave Decay
Coefficient ............................................. 79

4-8 Non-Dimensional Interfacial Wave Amplitude versus
Complex Reynolds Number ................................. 82

4-9 Resonance Phenomenon of Water-Mud System.
(a) with d2 = 0.05 m; (b) with d2 = 0.14 m .............. 84

4-10 Frequency Response of the Water-Mud System. Mud is
Assumed as (a) Viscoelastic Material; (b) Viscous Fluid. 85

5-1 Apparatus for Relaxation Tests and Viscoelastic
Constant Measurement. (a) System Diagram; (b) Top View.. 88

5-2 Structure of Brookfield Viscometer ...................... 89

5-3 Miniature Vanes and Sample Container .................... 90

5-4 A Plot of the Residual Torque and Angular Displacement
versus Elapsed Time for Run 1, Kaolinite with One-Week
Consolidation Period .................................... 94

5-5 A Plot of the Residual Torque and Angular Displacement
versus Elapsed Time for Run 4, Cedar Key Mud with
One-Week Consolidation Period ........................... 95

5-6 A Plot of the Residual Torque and Angular Displacement
versus Elapsed Time for Run 5, Cedar Key Mud with
Two-Day Consolidation Period ............................ 96

5-7 A Plot of Shear Stress against Shear Strain for
Determination of Viscoelastic Constant .................. 99

5'8 Relationship between Viscoelastic Constants and Mud
Dry Density. (a) Viscosity; (b) Shear Modulus ........... 100

5-9 Shear Stress versus Shear Strain Rate obtained from
Kaolinite Bed, Run 1, 6.2 cm below Mud Surface .......... 104

6-1 Wave Flume for the Erosion Experiments .................. 106

6-2 Water-Mud Interface Elevations during Deposition
and Consolidation of the Kaolinite Slurry ............... 111

6-3 Dispersed Grain Size Distribution for Sediment
in Runs 4, 5, and 6 ..................................... 114

6-4 Pressure Transducer, Druck Model PDCR 135/A/F ........... 117

6&5 System Diagram of the PhotoSensing Light Meter ......... 119

6-6 Suspended Sediment Sampler .............................. 123

6-7 Wave Loading for the Erosion Test, Run 3 ................ 127

6-8 A Plot of Wave Height Decay in the Wave Flume ........... 129

7-1 Measured Mud Density Variation due to Wave Action,
Run 1 ................................................... 132

7f 2 Volumetric Swelling of Mud Bed Caused by Wave Action.
(a) Run 1; (b) Run 4 .................................... 134

7-3 Depth-Averaged Bed Density Variation Caused by
Wave Action ............................................. 135

7-4 Dimensionless Mud Dry Density Profiles. (a) Kaolinite;
(b) Cedar Key Mud ....................................... 136

7-5 Instantaneous Pressure Response in the Kaolinite Bed,
Run 1. (a) Wave Profile; (b) Pressure ................... 138

7,6 Average Pressure Response in the Kaolinite Bed, Run 1.
(a) Dynamic Pressure Fluctuation; (b) Wave-Averaged
Pressure ................................................ 140

7-7 Variation of Dimensionless Apparent Bulk Density Caused
by Wave Action, Run 1, Normalized by the Initial Value..141

7-8 Instantaneous Sediment Concentration Response in
the Water Column. (a) Run 3; (b) Run 5 ..................143

7-9 Wave-Averaged Response from Light Meter during
Wave Erosion Process .................................... 144

7-10 Average Sediment Concentration Fluctuation Detected
by the Light Meter during Wave Erosion Process .......... 145

7-11 Sample Profiles for the Suspended Sediment
Concentration for Run 1. (a) at STA. B; (b) at STA. D... 147

7-12 Sample Profiles for the Suspended Sediment
Concentration for Run 5. (a) at STA. B; (b) at STA. D ...148

7-13 Longitudinal Variation of the Depth-Averaged
Sediment Concentration, Run 4 ........................... 150

7-14 Mud Surface Profiles at Selected Times.
(a) Run 1 ............................................... 152

7-14 (continued). (b) Run 2; (c) Run 3 ....................... 153

7-15 Response of the Spatially-Averaged Mud Surface
Elevation. (a) Kaolinite Bed; (b) Cedar Key Mud Bed ..... 155

7-16 Wave Erosion/Entrainment Behavior for Kaolinite.
(a) Run 1 ............................................... 157

7-16 (continued). (b) Run 2 .................................. 158

7-16 (continued). (c) Run 3 .................................. 159

7-17 Wave Erosion/Entrainment Behavior for Cedar Key Mud.
(a) Run 4; (b) Run 5 .................................... 160

7-17 (continued). (b) Run 6 .................................. 161

7-18 Shear Strength Profiles for Resistance of Wave Erosion.
(a) Kaolinite Bed; (b) Cedar Key Mud Bed ................ 167

7-19 Influence of the Period of Consolidation on Erosion
Resistance. (a) Kaolinite; (b) Cedar Key Mud ............ 168

7-20 Erosion Rate versus Dimensionless Excess Shear Stress.
(a) Kaolinite; (b) Cedar Key Mud ........................ 170

7-21 Dimensionless Erosion Rate versus Dimensionless
Excess Shear Stress ..................................... 172


7-22 Erosion Rate versus Bed Shear Stress, Run 5 ............. 173

7-23 Suspended Sediment Concentration Profiles at
Steady State. (a) Run 1; (b) Run 2 ...................... 175

7-23 (continued). (c) Run 3; (d) Run 4 ....................... 176

7-23 (continued). (e) Run 5; (f) Run 6 ....................... 177

7-24 Settling Velocity versus Sediment Concentration for
Kaolinite ............................................... 179

7-25 Spatially-Averaged Variation of (a) the Distance
between the Lowest Tap and Mud Surface;
(b) Sediment Concentration at the Lowest Tap ............ 181

7-26 Diffusion of a Line-Source Dye in a Wave Flow with
T = 1.2s, H 7 cm. (a) t = 0.5s; (b) t 5s;
(c) t 14 s; (d) t 20 s .............................. 183

B-I Comparison between Observed and Predicted Horizontal
Velocity Amplitudes in the Bed .......................... 210

B,2 Comparison of Model Prediction and Measurement
for Run 1-1. (a) Velocity; (b) Pressure ................. 211

B-3 Comparison of Model Prediction and Measurement
for Run 2-1. (a) Velocity; (b) Pressure ................. 212

B4 Comparison of Model Prediction and Measurement
for Run 2-2. (a) Velocity; (b) Pressure ................. 213

BA5 Comparison of Model Prediction and Measurement
for Run 3-1. (a) Velocity; (b) Pressure ................. 214

B-6 Comparison of Model Prediction and Measurement
for Run 3-2. (a) Velocity; (b) Pressure ................. 215

B-7 Comparison of Model Prediction and Measurement
for Run 4-1. (a) Velocity; (b) Pressure ................. 216

B-8 Comparison of Model Prediction and Measurement
for Run 4-2. (a) Velocity; (b) Pressure ................. 217

B,9 Comparison of Model Prediction and Measurement
for Run 5-1. (a) Velocity; (b) Pressure ................. 218

B-10 Comparison of Model Prediction and Measurement
for Run 6-1. (a) Velocity; (b) Pressure ................. 219

B,11 Comparison of Model Prediction and Measurement
for Run 6-2. (a) Velocity; (b) Pressure ................. 220


B-12 Comparison of Model Prediction and Measurement
for Run 63. (a) Velocity; (b) Pressure ................. 221

C-i Normalized Frequency Response of the Current Meter......225

C-2 Circuit Diagram for the LED Driver of Light Meter ......227

C-3 Light Meter Probe ....................................... 228

C-4 Circuit Diagram for the Pre-Amplifier of Light Meter .... 230

C-5 Circuit Diagram for the Main Processor of Light Meter.. 231

C16 Light Meter Calibration Device .......................... 232

C-7 Light Meter Calibration Curves for Kaolinite ............ 233

C-8 Light Meter Calibration Curves for Cedar Key Mud ........ 234

C-9 Response Function of the Low-Pass Filter ................ 236

C.I10 Comparison of Data Record with/without the
Low-Pass Filter ......................................... 237

D-1 Wave Loading for the Wave Erosion Test, Run 1 ........... 253

D-2 Wave Loading for the Wave Erosion Test.
(a) Run 2; (b) Run 4 .................................... 254

D-3 Wave Loading for the Wave Erosion Test.
(a) Run 5; (b) Run 6 .................................... 255

D-4 Wave-Averaged Pressure Response in the Mud Bed.
(a) Run 2; (b) Run 3 .................................... 256

D-5 Wave-Averaged Pressure Response in the Mud Bed.
(a) Run 4; (b) Run 5 .................................... 257

D-6 Wave-Averaged Pressure Response in the Mud Bed, Run 6.. .258

D-7 Average Dynamic Pressure Fluctuation in the Mud Bed,
Run 2 ................................................... 258

D-8 Average Dynamic Pressure Fluctuation in the Mud Bed.
(a) Run 3; (b) Run 4 .................................... 259

D-9 Average Dynamic Pressure Fluctuation in the Mud Bed.
(a) Run 5; (b) Run 6 .................................... 260

D-10 Average Apparent Bulk Density Variation.
(a) Run 2; (b) Run 3 .................................... 261

D-11 Average Apparent Bulk Density Variation.
(a) Run 4; (b) Run 5 .................................... 262

D-12 Average Apparent Bulk Density Variation, Run 6 .......... 263

D-13 Variation of Kaolinite Bed Density Profile Due to
Wave Action, Run 2 ...................................... 264

D-14 Variation of Kaolinite Bed Density Profile Due to
Wave Action, Run 3 ...................................... 265

D-15 Variation of Cedar Key Mud Bed Density Profile
Due to Wave Action, Run 4 ............................... 266

D-16 Variation of Cedar Key Mud Bed Density Profile
Due to Wave Action, Run 5 ............................... 267

D-17 Variation of Cedar Key Mud Bed Density Profile
Due to Wave Action, Run 6 ............................... 268

D-18 Wave-Averaged Response from Light Meter, Run 6.
(a) Sediment Concentration; (b) Sediment Concentration
Fluctuation ............................................. 269


a: Wave amplitude in the x wave propogation direction.

a : Wave amplitude at x = 0.
ab: Semi-excursion distance of a water particle just outside 6.

A/D: Analog to digital.

boi: The i-th interfacial mud wave amplitude at x 0, a complex no.

C: Suspended sediment concentration (g/l).

c: Complex coefficient matrix c(i,j) or wave celerity.

cd: Damping coefficient in the disk-shaft system.

cg: Wave group velocity.

: Average suspended sediment concentration at the 3 lowest taps.

d i: Thickness of i-th layer.

d : Total thickness of mud bed.

D: Rate of deformation tensor, D I (- j )"
S2 axj Dxi
D': Deviator part of the rate of deformation tensor, Dj= D -- D
i ii 3 r
Dm: Apparent Diffusion coefficient (molecular diffusion + turbulent
diffusion + dispersion).

D : Diameter of the miniature vane.
D90 Grain size through which 90% of the total soil mass is finer
than this size.

E': Deviator part of the small strain tensor, El = Eij - E
1 3a ij 3 rr
E : Small strain tensor, Eij= + (- j -){.

j: Rate of strain.

j': Rate of shear strain.

En: Wave energy.

Ed : Time-averaged energy dissipation per unit surface area.





















Wave friction factor.

Acceleration due to gravity.


Shear Modulus of elasticity.

A constant coefficient.

Water depth.

Wave height.

Height of the miniature vane.

Free index.

Free index or V --T.

Grid number at the water surface.

Complex wave number, k kr + j kim

Imaginary part of the complex wave number, also the wave
decay coefficient, defined as a(x)= a exp(-k imAx).

The real part of the complex wave number, k r= 27/L.

Yield strength of Bingham Fluid.

Wave length.

Constant coefficient, or free index.

Sediment mass entrained into water column measured by water
samples withdrawn from water column.

Sediment mass eroded by waves.

Constant coefficient.

Number of layers in the hydrodynamic model presented in Cha

Dynamic pressure.

Correction term for the static pressure due to density diff

Total pressure;

Dynamic pressure amplitude.

Radius, or free index.

pter 4.


R: Complex Reynolds number, R Rr + j Rim.

Rw: Wave Reynolds number, R = 2 / o-
w 0w 0V
RMS: Root mean square.

S : Undrained shear strength.

SWL: Still water level.

t: Time.

T: Wave period or stress tensor.

T': Deviator part of the stress tensor, Tij= Tij Trr
Tdc: Period of Consolidation of a sediment slurry.

j': Rate of shear stress.

ui Velocity in the xi direction; horizontal velocity for i-th layer.

Ub: Horizontal velocity amplitude just outside the boundary layer.

u1: Horizontal velocity amplitude.

w Vertical velocity.

wi: Vertical velocity amplitude.

wl: First derivative of wi with respect to z.
w": Second derivative of wi with respect to z.
w : Settling velocity of a sediment particle or floc.

x: Longitudinal direction, or xI.

xi: Cartesian coordinate using index notation, i 1, 2, 3.

X: Column matrix containing Al, BI, .... CN, DN, and b.,N-1.

y: Laterial direction, or x2.

z: Vertical coordinate, or x3, starting from the water surface,
positive upward.

z': Vertical coordinate, starting from mud surface, positive downward.

z": Vertical coordinate, starting from bed, positive upward.

z : An elevation where a pressure transducer was mounted.

z : An elevation where a pressure transducer was mounted.

: Average distance from the 3 lowest taps to the bed surface.

At: Increment in the time domain.

Ax: Increment in the x direction.

Az: Increment in the z direction.

a: Constant coefficient.

a : Constant coefficient.

B: Constant coefficient.

a: Wave frequency, a = 2n/T.

6: Thickness of wave boundary layer.

6 L: Thickness of the viscous wave boundary layer.

p: Material density.

PB: Bulk density of mud.

PD: Dry density of mud.
P: Depth-averaged dry density of mud at a time.

PDo: Depth-averaged dry density of mud before wave loading.

Pe: Average apparent bulk density of mud between two sensors.

Peo: Pe value at the beginning of wave loading.

Ps: Density of sediment grain.

pwt Density of water.

p': Normalized apparent bulk density between two sensors, = p / p
e e
T: Shear stress.

;: Shear stress amplitude.

Tb x): Bed shear stress amplitude at x=x.

Tbol Bed shear stress amplitude at x-0.

: Spatially-averaged bed shear stress amplitude.


T Critical shear stress. The minimum bed shear stress that will
cause erosion of sediment bed.

Te Excess shear stress, Te =T b Ts(z').

Ts: Erosion resistance or bed shear strength.

Tao: Erosion resistance at the bed surface.

Ts: Depth-averaged (top 5 mm) erosion resistance of mud bed.

n: Water surface displacement.

C: A constant eddy viscosity coefficient.

c Erosion rate constant for wave flow.
Cf Erosion rate constant for steady flow.

c n Rate of upward entrainment of sediment.

E Erosion rate.
E: Interfacial wave displacement, 0 b J(kx-at).

Ci: Displacement in the x direction.

U: Dynamic viscosity.

Pe: Apparent dynamic viscosity, Pe = PBe"

Vo: Constant coefficient.

pw: Dynamic viscosity of water.

Im: Dynamic viscosity of mud.

V: Kinematic viscosity.

v e: Apparent kinematic viscosity.

e: Angular shear strain.

de/dt or 6: Rate of angular shear strain.

IID,: Second invariant of the deviator part of D'.

: Diameter.

0: Angular speed.

A: Energy dissipation ratio.

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



P.A-Y. Maa

August 1986

Chairman: Dr. A. J. Mehta
Cochairman: Dr. B. A. Christensen
Major Department: Civil Engineering

The influence of water waves on soft muds can result in high

turbidity and significant wave attenuation. A multi-layered 2hD hydro-

dynamic model was developed and experimentally verified in order to

evaluate the bed shear stress at the water-mud interface, for inter-

preting results on interfacial sediment erosion. In the model, the mud

bed was discretized into layers of constant density, viscosity and

shear modulus of elasticity. The last two parameters characterize mud

as a linear viscoelastic material. Layering was essential for simu-

lating the depth-varying bed properties. For given non-breaking,

regular waves, the model gives profiles of velocities, pressure and

shear stress as well as wave attenuation coefficient.

Laboratory flume tests of sediment erosion under waves were con-

ducted using a commercial kaolinite and an estuarine mud. Sediment

concentration response in the water column and the variation of total

eroded sediment mass with elapsed time under different wave loading

were measured. Additional tests were conducted to evaluate the visco-

elastic constants of the selected sediments.


Two types of erosion process were identified. The first is

applicable to a partially consolidated bed, in which erosion resist-

ance, Ts(Z'), increases with depth, z', and the erosion rate, Er,

decreases as erosion proceeds. The proposed erosion rate function is

r Lb T, >? br sZ
E-- = Or [ I ] r>r (z')

E0 5 Z'

where Tb is the bed shear stress, co is a rate constant and a is an

empirical exponent. The average value of a for all test data was close

to one, which reveals a simple, linear relationship. The second type

of erosion process is for a newly deposited bed, in which the erosion

resistance is a small constant, and the erosion proceeds at a constant


Measured suspended sediment concentration profiles indicated

relatively low concentrations in the upper 80% of the water column and

the formation of a high density fluid mud layer near the bed. Low

upward diffusion, high longitudinal dispersion and advective transport

characterized the near-bed layer. Bed softening and fluid-mud-gene-

rating capacity of waves appear to be significant features of the wave

erosion process.


1.1 Significance of Wave Erosion of Soft Muds

Most estuaries and several reaches of shorelines have soft muddy

bottoms consisting primarily of silt and clay comprising a cohesive

sediment. With the rapid development of coastal and estuarine harbors

and increasing concerns from the environmental protection point of

view, the need to predict the erosion of soft mud has increased.

Basically, in these environments, tidal currents and water waves are

the two forces that cause erosion of sediments. Extensive studies of

the erosion behavior of cohesive sediment under tides, in quasi-steady

flow, have been conducted by various researchers. Expressions of the

erosion rate have been proposed (Partheniades, 1965; Krone, 1962; Mehta

et al., 1982; Parchure, 1984). Mathematical models for cohesive sedi-

ment transport by current are also available (Odd and Owen, 1972;

Ariathurai, 1974; Hayter, 1983). However, the influence of water waves

on the erosion of cohesive sediment has by and large remained unex-

plored. Only a few important laboratory studies have been made

(Alishahi and Krone, 1964; Anderson, 1972; Jackson, 1973; Thimakorn,

1980). However, data obtained have been scattered and incomplete. A

comprehensive and systematic study of wave erosion behavior is neces-

sary to elucidate the mechanism of the wave erosion process. Here the

physical processes which have been observed at muddy coasts are intro-

duced first.

The most obvious feature of the muddy coastal zone is the lack

of large waves due to the great damping ability of the soft mud beds

(Wells, 1983). This has been noticed since the early 1800s in South-

west India and northeastern South America. The damping ability comes

from the motion of the mud bed, a phenomenon observed in both labo-

ratory and field studies (Tubman and Suhayda, 1976; Schuckman and

Yamamoto, 1982). An immediate benefit of this high damping ability

is protection of shorelines. Fishermen in the Gulf of Mexico have

taken refuge during stormy weather in waters over the "mud hole," a

localized accumulation of soft mud bottom (Morgan et al., 1953).

High suspended sediment concentration (0.1 -2 g/l) is another

feature of soft muddy coasts even in calm weather (wave height 10

cm). In severe wave conditions, the entire soft mud can be stirred

up resulting in great turbidity.

The mechanism of wave transformation when propagating over soft

mud beds and the resulting erosion of mud are important to both coastal

and environmental engineers because waves and the associated sediment

transport are major factors for design of most coastal waterways and

harbors, e.g., the maintenance of navigation channels, and the control

of pollutants. However, studies on this subject are limited apparently

because the change of characteristics of surface water waves caused by

the dissipation of wave energy in the mud and resuspension of mud are

interlinked and complex. Added complexities arise from the influences

of physicochemical properties of the sediment and the fluid (Mehta,


1983), the effect of bed structure (Parchure, 1984), and the time

dependence of bed structure (Schuckman and Yamamoto, 1982) on the

overall problem.

1.2 Factors Characterizing the Wave Erosion Process

A summary of the parameters which influence the wave erosion

problem is presented in Fig. 1-1. Thick arrows show the dominant

processes, and thin arrows are the less important processes. Dashed

lines represent processes or parameters that are known to have influ-

ences on the dynamics of wave-mud interaction, but which were not

considered in detail in this study. The responses of this wave-mud

system, contained in the dashed rectangular box, subjected to dynamic

wave loading, are dependent on the characteristics of both wave and

mud. Pressure and shear stress at the water-mud interface are the two

forces that have linked the water and mud bed together. Two immediate

responses are wave damping and erosion, shown in the circles. The

contribution of pressure fluctuation at the interface is mainly res-

ponsible for the motion of the mud beds, and therefore, wave damping.

Shear stresses mainly contribute to the erosion process, a conclusion

from steady flow studies. Turbulence in the water flow influences the

bed shear stress and the sediment transport process as well. Erosion

is dependent on the properties related to waves, i.e., bed shear

stress, as well as erosion resistance of the mud material, i.e., shear

strength. Therefore, development of a hydrodynamic model that can

simulate the water wave decay and the interfacial shear stress (also

called bed shear stress, Tb) is necessary before an erosion rate


expression can be developed. With Tb known, the results from wave

erosion tests can be better interpreted.

Defining the properties of soft mud beds is in itself an inter-

esting subject. Such beds were assumed to be a viscous fluid by Gade

(1958) and Dalrymple and Liu (1978); an elastic material by Mallard and

Dalrymple (1977) and Dawson (1978); a poro-elastic material by Yamamoto

and Takahashi (1985); a viscoelastic material by MacPherson (1980) and

Hsiao and Shemdin (1980); a Bingham fluid by Migniot (1968), Krone

(1965), Pazwash and Robertson (1971), Gularte (1978), Vallejo (1979,

1980) and Englund and Wan (1984). However, what is important in this

case is to find a model that can reasonably simulate the two important

parameters: bed shear stress and mud motion. These two parameters are

directly related to wave erosion and wave attenuation. This is the

first step in studying the wave erosion process.

Initially a 1-D numerical Bingham fluid model was developed by the

author to simulate these two observed phenomena. However, the results

demonstrated that it could not appropriately simulate wave-mud interac-

tions. A linearized multi-layered model was then developed to predict

the wave decay and the bed shear stress. In this model, the mud bed

was arranged as layers of linear viscoelastic material to meet the

depth variation of mud properties. The water layer was characterized

by having a constant viscosity. Velocities, pressure, shear stress

profiles and wave attenuation coefficient were the products.

Following the selection of a model for the soft mud, an effort was

made to measure the mud properties by using the viscometer and other

auxiliary equipment.



r -- .. .__ _

Water: Turbulence Wave '
Tu uk cea q rre n4u
Newtonian Fluid -'- ..urr....
ft I 2 1 T rso Sedirnm 1 --

5 Soft Mud:

Non-Newtonian 'TmI
I Material I ime- Deperdent\
L Jl 1 Long -Term Bed
\ProPerty Change /

1. Pressure and shear stress at the water-mud interface contribute
to the motion and erosion of the soft mud.

2. Momentum transfer from waves to muds alters the kinematics of water

3. Turbulence structure of the water determines the suspended sediment

4. Bed shear stress and the mud resistance to erosion determine the
rate of erosion.

5. Mud motion and associated energy dissipation plus the minor energy
loss from the water viscosity cause the attenuation of water waves.

6. Eroded sediment mass serves as a sediment source in the sediment

7. Wave-induced second order current contributes to the sediment
transport process.

8. Successive oscillatory shear deformation in the mud bed alters
the mechanical properties of the mud, and therefore, the dynamic
response of the wave-mud system.

Legend: Box with dashed line represents the system, ovals show the
processes, circles show the results.

Fig. 1-1. Schematic Depiction of Response of Wave Mud System.

Because of the complexity of the erosion process itself, labora-

tory tests to observe wave erosion under clearly specified conditions

were absolutely necessary. Alishahi and Krone (1964), Jackson (1973),

and Thimakorn (1980) have all done some initial studies. However,

their data were not complete enough to give an expression for the

erosion rate function. Therefore, wave erosion tests were performed to

provide information on the actual erosion and entrainment behavior. In

the present study, quantities of sediment mass eroded were obtained and

incorporated with the bed shear stress calculated from the 2-D mathe-

matical model to develop an erosion rate function.

The diffusion process in the wave flow field was studied by ob-

serving the spreading of dye and the suspended sediment concentration

profiles. A conclusion concerning the possible profile of the diffusion

coefficient is proposed.

Marine sediment beds tend to change their physical properties

under wave action. Alishahi and Krone (1964) reported a sudden loos-

ening of the bed. Schuckman and Yamamoto (1982) reported the reduction

of shear strength of a bentonite bed. Both observations point to the

degradation and possible liquefaction of mud bed under wave loading.

This may result in a change of flow structure and a change in the

interfacial shear stress and energy dissipation as well. Therefore,

this long-term bed property change was also monitored in the wave

erosion tests.

In addition to the above, there are other complexities which need

to be briefly addressed: the effects of salinity differences between


the eroding fluid and the pore water; the construction of a mattress-

like layer by organic material; the effects of wave-induced current,

etc. Although these were not here studied, they may be considered in

future investigations.

1.3 Objectives and Scope of Study

The objectives of this study were as follows:

1. to develop a hydrodynamic model that can be used to simulate the two

important parameters of wave-mud interaction: mud bed motion and bed

shear stress,

2 to conduct laboratory experiments to obtain necessary data as input

for the model,

3. to conduct laboratory experiments to verify model results,

4. to conduct laboratory experiments on the wave erosion and entrain-

ment of soft mud beds, and

5. to propose an erosion rate function.

The scope of this study was defined as follows:

1. The hydrodynamic model developed was limited to the primary forcing

mechanism only, i.e., linear waves; secondary effects were excluded.

2. Water and mud were considered immiscible layers in the model.

3. The equations of motion were linearized by neglecting the convective

acceleration terms.

4. Verification of model results was limited to comparing measurements

of wave attenuation, velocities, and pressures to calculated values.

5. Wave erosion tests were conducted in a wave flume to be described in

Chapter 6. Mud bed thicknesses were limited to within 10 15 cm.

Waves were restricted to regular (monochromatic), progressive and

non-breaking types.

6. The pore fluid and the eroding fluid were the same. This was

achieved by equilibrating the sediment with the eroding fluid.

7. Vertical sediment concentration profiles were observed but no

attempt was made to develop a theoretical approach for simulating

the depth-distribution.

8. Three consolidation periods, 2 days, I week, and 2 weeks, for the

mud beds were selected, and limited to self-weight consolidation.

9. Erosion rate studies were limited to two flocculated sediments and

one fluid. Kaolinite and Cedar Key (Florida) mud were the two

sediments. Tap water without added salt was the eroding fluid.

1.4 Outline of Presentation

This study is presented in the following order. A brief review

of previous studies on wave erosion of soft muds is given in Chapter

2, followed by a discussion of problems in those studies, including

those related to constitutive models for the soft mud, bed shear

stress, and diffusion features in the wave flow. Chapter 3 describes

an effort to develop a 1-D numerical Bingham fluid model to simulate

the water-mud system and the conclusion that the 1-D model is not

necessarily a good approach. Chapter 4 presents a 2-D hydrodynamic

model which treats muds as a viscoelastic material. Bed shear stresses

needed for correlation with the erosion data from Chapter 6 are also

presented. Chapter 5 describes the experimental procedures used to

obtain the viscoelastic constants for the mud beds used in the wave


erosion tests. Experiments to measure the erosion of kaolinite and

Cedar Key mud are given in Chapter 6 and Chapter 7. Erosion rates are

also presented. Summary and conclusions of the study are given in

Chapter 8. The derivation of boundary conditions for the multi-layered

hydrodynamic model is presented in Appendix A. In Appendix B the

horizontal velocity and pressure predicted by the model are compared

with measurements. Calibration of the current meter in an oscillatory

flow, design of the light meter, and development of a numerical low-

pass filter are given in Appendix C. In Appendix D, data from the wave

erosion experiments are presented.


2.1 Introduction

Previous studies on the wave erosion, constitutive models for mud

and their behaviors, and wave diffusion coefficients will be reviewed.

These subjects are important in the understanding of the wave-mud system

and provide necessary information to simulate the dynamics of this


2.2 Wave Erosion/Resuspension

Water waves apply periodic shear stresses at the water-mud inter-

face. Whenever the shear stress exceeds the bonding force (cohesive

forces and self weight) between particles, erosion occurs. The eroded

material is transported and dispersed into the water column and high

turbidity is subsequently observed.

Investigations of wave erosion phenomena have been made primarily

on cohesionless material, i.e., sand beds. Only a few experiments deal

with cohesive sediment. In this section, erosion studies with cohesive

sediment under non-breaking wave conditions are discussed.

Alishahi and Krone (1964) conducted the earliest wind wave erosion

experiments for the San Francisco Bay mud. The test section, where the

mud was placed, was 1.2 m in length and located at the downstream side

of the wave flume (18 m in length). The average thickness of the mud

beds was about 1 cm. Water samples were taken at several locations

via a special device that gave only depth-averaged sediment concentra-

tion. Total suspended sediment mass was then obtained by summation of

the product of sediment concentration and water volume. Figure 2-1

shows one of their two test results. Their most significant result was

"the demonstration of the existence of critical wave induced shear for

suspension of the bed surface." No or negligible erosion occurred for

a bed shear stress less than the critical shear. See the response in

the first 4 hr in Fig. 2-1. Another interesting response reported was

"a sudden loosening of the bed and direct movement of sediment into

suspension." This phenomenon occurred in both runs and was responsible

for a significant increase of suspended sediment. See the response

after 4.5 hr in Fig. 2-1. However, possible reasons for this sudden

loosening were not given. Their Run 1 showed that an equilibrium state

had been reached during the erosion process. Unfortunately, Run 2 was

not sufficiently long for an equilibrium condition to have been reached,

as was noted by the investigators.

Anderson (1972) performed field experiments in The Great Bay

estuarine system in New Hampshire. He took records of tide, wind

speed, wind direction, wave height, and suspended sediment concentra-

tion. He found a linear relationship between suspended sediment con-

centration and water wave height at flood tide, while at ebb tide

these two factors could not be correlated. He also concluded that wave

force is most important for introducing sediment into the water column.

#. 0.48

g W 0 6 -. " '

CO 0.08-

0 2 4 6 8


Fig. 2-1. A Plot of the Suspended Solids and Bed Shear Stress for a Wave Resuspension Test
(After Alishahi and Krone, 1964).

Jackson (1973) conducted laboratory experiments to study the

effects of water waves on the resuspension of muds on a beach. Breaking

waves were involved because of the sloping beach. He found that waves

produced little erosion or resuspension when the bottom orbital velocity

was less than 10 12 cm/s. For orbital velocities greater than 20

cm/s, particle-by-particle erosion behavior shifted to mass erosion.

Thimakorn (1980) conducted non-breaking wave resuspension tests of

mud taken from the Sakohon River mouth, an estuary in central Thailand.

His wave flume had dimensions of 45 cm width, 60 cm depth, and 45 m

length. The grain size distribution and clay composition were not

given but the investigator noted that all sediments were washed through

a U.S. standard sieve no. 200. The mud was rinsed with fresh water

and the experiments were conducted in fresh water. Mud slurry was

uniformly distributed in the entire flume, and was allowed to settle for

two weeks. The final bed thickness was 2.5 cm. He took water samples

at nine elevations at only one location. Figure 2-2 shows examples of

the concentration profile. The data show a strong vertical mixing

mechanism illustrated by the uniform concentration profile for the top

25 cm and gradient over the bottom 5 cm. He integrated these profiles

and obtained the time variation of total suspended sediment concentra-

tion. From these data, he concluded that the wave erosion process

did reach a steady state after 3 or 4 hours. However, the possibility

of a strong sediment concentration gradient in the longitudinal direc-

tion was not considered. This gradient, for example, could be induced

by the longitudinal variation of bed shear stress. Therefore, his

statement regarding the steady state condition may not be fully



U. 30
z- t a~in

0 00

00 500100

t 00min


0) 500- 1ID


h= 29, 4 cm
H = 5.4 cm
T= 1.7 s
Rw =4900

Fig. 2-2. Measured Suspended Sediment Concentration Profiles
from a Wave Resuspension Test (after Thimakorn, 1980).


justified in case of a strong longitudinal gradient of sediment concen-


Thimakorn (1984) also presented a mathematical model to simulate

the suspended sediment concentration profile. This was composed of (1)

a shear model for predicting the sediment mass eroded by wave shear

stress (which had a similar function to that of the erosion rate) and

(2) a diffusion model for distributing the eroded sediment into the

water column. The diffusion model, based on a constant diffusivity,

should have had the well-known log-linear distribution curve for the

sediment concentration profile at steady state. However, Thimakorn

did not demonstrate this characteristic but introduced a depth-varying

diffusivity profile and asserted that the results were good compared to

measured data.

Some relevant conclusions based on the aforementioned studies

are as follows:

1. The critical shear stress concept for erosion by steady flow

(Partheniades, 1965; Mehta et al., 1982) can also be applied for

wave erosion. In other words, wave-induced bed shear stress must

exceed a critical value to cause noticeable erosion.

2. The mud beds used in these studies were quite thin. They were on

the order of 1 to 2.5 cm. Therefore, the interaction between the

mud bed and water waves was limited. Water wave attenuation was not

mentioned in these studies, possibly because it was not noticeable

for such a small mud thickness. However, this is not the case in

the prototype, where significant wave decay has been reported

(Tubman and Suhayda, 1976; Wells, 1983). Therefore, a thicker mud

bed should be employed for wave erosion/resuspension tests.

3. Wave erosion rates were not mentioned even though this information

is most important for estimation of sediment transport.

4. Data on the vertical structure of suspended sediment profiles are

rare. One practical reason is the huge amount of sediment samples

that must be gathered for this purpose. An easy and reliable

measurement technique for rapidly recording cohesive sediment con-

centration is not available yet. It is believed that light attenu-

ation caused by individual sediment particles or flocs can be

utilized for this purpose, e.g., the Iowa Sediment Concentration

Measuring System for suspended sand (Glover et al., 1969). However,

unresolved problems still remain in measuring silt and clay suspen-


5. It remains unclear whether the wave erosion process reaches a

steady state or not. Suppose a steady state can be reached, then

it must reflect a downward increase of erosion resistance in the

bed, provided no deposition is occurring. A similar behavior has

been observed for steady flow erosion processes (Mehta et al.,

1982). However, more data are required to support the premise

of a steady state for wave erosion.

6. The sudden movement of bed material may be caused by the cyclic

mobility or liquefaction of the bed by water waves. The former,

for clays subjected to oscillatory shear loads, has been demon-

strated by Schuckman and Yamamoto (1982) for bentonite beds. The

latter, for silt or fine sand under cyclic loads, also has been

verified by Turcotte et al. (1984). Pamukcu et al. (1983) claimed

that these two terminologies (cyclic mobility and liquefaction)

actually describe the same phenomenon but for a different type of

soil. The important results from this "loosening" are the changing

flow structure (for both water and mud) as well as the erosion


7. In the previous studies noted, the bed layers were assumed to be

rigid. Bed shear stresses were evaluated from linear wave theory

and the wave friction factor concept. Such an assumption is valid

for a hard soil or compacted sediment of low water content. In

general, particularly for soft muds, motion within the bed can be

significant, and the shear stress between the eroding fluid and the

bed must be evaluated by considering the dynamics of mud motion.

However, the question of which bed properties should be selected has

not been resolved.

2.3 Constitutive Models for Soft Muds

Elastic, poro-elastic, viscous fluid, viscoplastic, and visco-

elastic behavior are the five behavioral models which have been proposed

for simulating mud properties. They are discussed in this section. It

may be unrealistic to say that any simple model, such as those described

below, is good enough to represent the rather complex behavior of mud.

However, in general, the goal is not to obtain an absolutely correct

answer but to find a relatively reasonable predictive model.

Before discussing any model, it may be helpful to examine the

nature of a clayey soil under a dynamic shear load. Although the


characteristics of the clay bed discussed are for quite low water

contents, it does give some insight into soft bed behavior.

2.3.1 Nature of Clayey Soil

Hardin and Drnevich (1972) measured clayey soil response under

cyclic load. They followed the definition of energy loss ratio given by

Kolsky (1963) as the ratio of the energy dissipation in taking a spec-

imen through a stress cycle and the elastic energy stored in the spec-

imen when the strain is a maximum. See Fig. 2-3. For a stress-strain

cycle not too far from the linear range, the energy dissipation is

small (see the small hysteresis loop). However, when the strain is

far off the linear range, the stress-strain path follows the large

hysteresis loop and results in high energy dissipation. Kovacs et al.

(1971) indicated that the ratio of energy dissipation, X, increases

at large shear strains. See Fig. 2-4. This implies that, for a given

shear load, soft muds have high energy dissipation rates due to the

large shear strain. The more dense muds, however, have smaller strains

and consequently small energy dissipation as well. This also implies

that energy dissipation is velocity-dependent (strickly speaking,

shear displacement-dependent) since high velocities usually mean large

shear displacements.

Thiers and Seed (1968) demonstrated that the shear strength of

clays can be reduced when subjected to a sufficient number of cyclic

shear loads. See Fig. 2-5. This phenomenon was also observed by

Schuckman and Yamamoto (1982) for marine sediment from measurement of

the undrained shear strength, Su, of a bentonite bed before and after

Strain 7

/-Hysteresis Loop

Fig. 2-3. Shear Stress-Strain Loop for Clay under
Cyclic Loading.


0.04 .

0.1 0.4 1.0 40 10.0

Fig. 2-4. Energy Dissipation Ratio A versus Shear
Strain (after Kovacs et al., 1971).

Fig. 2-5. Shear Stress-Strain Loops and Effects of Cyclic
Loading (after Thiers and Seed, 1968).


wave loading. Su (N/m2) was empirically related to the shear modulus

G (N/m2) according to G = a Sum, with a = 0.66 and m= 1.75, applicable

to G values on the order of 102 to 104 N/m2. A maximum reduction of Su

of about 55% was reported. The reduction was fully recovered 9 days

after stopping the wave action. The reduction of shear modulus is

typically due to dynamic softening experienced by the sediment when

subjected to cyclic shear stresses.

2.3.2 Elastic Model

For an isotropic and homogeneous linear elastic body, the consti-

tutive equation for an incompressible material can be represented by

T'= 20 E' (2-1)


E. = (- -) (2-1a)
E. 2 3x i ax

where T and E are the stress and strain tensor, respectively, the prime

indicates the deviator parts of those two tensors, i,j, are free in-

dices, and ri is the displacement in the xi direction. Mallard and

Dalrymple (1977) used this model to study the stress and displacement

of a soil bed. Dawson (1978) also used this model but included the

inertia term to improve upon Mallard and Dalrymple's results.

A serious drawback of this model is that no damping effects can

be obtained because no viscous or friction terms are involved. There-

fore, this model can not justifiably represent the interaction of the

wave-mud system, where energy dissipation is important.


2.3.3 Poro-elastic Material with Coulomb Friction

Yamamoto and Takahashi (1985) suggested that the Internal loss of

energy is due to Coulomb friction and is independent of the velocity of

soil movement. Because of the nature of the energy loss mechanism they

chose, a direct measure of the energy dissipation of mud was selected.

This direct measure does not require any assumption concerning the

nature of internal friction. However, the results are dependent on the

strain amplitude as well as the previous stress experience. This

approach complicates the modeling of soil behavior because more pa-

rameters are involved.

Yamamoto (1982) incorporated the inviscid fluid flow assumption

for the water waves (therefore, the momentum is transferred to the mud

by normal pressure only). He then solved the equations of motion for a

poro-elastic soil skeletal frame and Darcy's equation for the pore water

pressure. Energy dissipation was evaluated by using a complex shear

modulus, in which the imaginary part stands for the friction damping

coefficient. Because of their inviscid fluid assumption for water, the

bed shear stress (an important parameter for erosion) was ignored. The

applicability of "poro-elastic material" to muds (which have low per-

meability and can be considered as a continuum) is not apparent.

2.3.4 Viscous Fluid Model

In this model the sediment bed is treated as a fluid which has a

greater viscosity and higher density than water. The constitutive

equation is

T' = 2p D' (2-2)

I. i + ._u. (2-2a)
mij i x

where D' is the deviator part of the rate of deformation tensor D which

is close to E for small strains, and ui is the velocity in the xi


Gade(1958) may have been the first to adopt this model to simulate

energy dissipation in a two-layered system. He assumed the overlying

water to be an inviscid fluid and the driving force to be restricted to

long waves only. In addition to theoretical derivation, he also con-

ducted experiments to verify his results. Unfortunately, his experi-

ments did not provide any information on real mud behavior because

he used a sugar solution instead of mud. Since a general model based

on the same viscous fluid assumption was developed later by Dalrymple

and Liu (1978), Gade's work is not further discussed.

Dalrymple and Liu (1978) published their."complete model," in

which they used two layers of (laminar) viscous fluid to represent the

overlying water and the soft mud bed. They then solved the linearized

Navier-Stokes equations for wave propagation. The boundary conditions

were specified as follows. First, the no-slip condition was satisfied

at the rigid bottom. Second, at the mud-water interface the horizontal

velocity, vertical velocity, shear stress, and normal stress were

matched. Third, the zero shear stress and zero normal stress conditions

were specified at the free surface. Finally, the linearized kinematic

boundary condition at the free surface and the water-mud interface were

invoked. This resulted in 10 equations and 10 unknowns. After some


manipulations, the 10 equations were reduced to a single one, and a

numerical method was utilized to solve it. Velocities, pressure,

and the interracial wave amplitude were typical outputs. Although the

shear stress profile was not shown it can be worked out easily from the

velocity gradients.

Based on their model, Fig. 2-6 shows the horizontal velocity

amplitude profiles for four mud viscosities, v2. All other parameters

are the same: wave height, H = 0.026 m, period, T = 1.75 sec, water

density, pl, and mud density, P2, equal to 1000 and 1260 kg/m3, respec-

tively, and water viscosity v, = 10-6 m2/sec. The interfacial shear

stress, Tb, and the decay coefficient, kim, are also displayed. It is

observed that because of the viscous fluid assumption, the interfacial

shear stress turns out to be small. Figure 2-6 indicates that the

horizontal velocity in the bottom layer is of the same order as that in

the upper water layer provided the kinematic viscosity of mud is less

than 0.01 m2/sec (1000 times higher than that for water). The dif-

ference in the velocity phase for these two layers is found to be about

45 degrees. This means that the velocity gradient at the water-mud

interface is small, which explains why the interface shear stress was


Also notice in Fig. 2-6, that for small v2, e.g., less than 10-3

m2/sec, the viscous wave boundary layer 22a is limited to the mud

layer. For high mud viscosity, the wave boundary layer thickness

extends up to the interface and a second wave boundary layer in the

water layer is developed. This second wave boundary is characterized

by the velocity overshoot above the interface. This overshooting is due

0.0001 0.017 0.018 i I
0.001 0.057 0.019 1 f j
0.01 0.235 0.099
0.05 0.148 0.179

H /


0- 41 6 81
2' 1'B0t


Fig. 2-6. A Plot of the Horizontal Velocity Amplitude Profiles
by Using the Two-Layered Viscous Fluid Model.


to the influence of the viscosity within the second boundary layer.

This phenomenon has been explained by Dean and Dalrymple (1984) and

observed by Sleath (1970) for a smooth bed, and Jonsson and Carlsen

(1976) for a rough bed.

The velocity in the lower layer is significantly reduced for high

viscosity because the entire mud layer is inside the viscous wave bound-

ary layer. However, the fact that measured viscosity data (see Chapter

5) indicate typically low viscosities, e.g., 10-3 to 10-4 m2/sec, for

the soft mud prohibits the use of high viscosity in the model. However,

model results from a high viscosity are closer to what has been observed

(Migniot, 1968; Nagai el al., 1982; Schuckman and Yamamoto, 1982).

This hints at another velocity restraining mechanism, which can be in-

corporated, for instance, into the viscosity term to create a high

apparent viscosity. This can then be used to simulate the mud motion.

Either the Bingham model's yield strength or the viscoelastic model's

elasticity could be used. These two models will be discussed next.

Despite the drawback in the selection of a viscous fluid model

to simulate wave-mud interaction, Dalrymple and Liu's model is indeed

the one with the most complete approach. Efforts to generalize their

approach as well as to incorporate other restraining mechanisms have

been made and are presented in Chapter 4.

2.3.5 Viscoplastic model

Viscoplastic material, also called a Bingham plastic, is charac-

terized by a special feature recognized as the yield strength. Basi-

cally, a viscoplastic material can sustain a shear stress even when it


is at rest, but when the shear stress intensity reaches a critical

value the material flows, with viscous stresses proportional to the

excess stress intensity. The constitutive equation (Malvern, 1969) is

+ K(+ )' (2-3)

D' 2 i D! (2-3a)
I' ij lj

where IID, is the second invariant of the deviator part of the rate

of deformation tensor, K and 4 are the yield strength and apparent

viscosity. Notice that the extra term, i.e., the second invariant in

the constitutive equation, results in a nonlinear equation of motion

which is difficult to solve analytically except in a few simple steady-

flow cases. Unfortunately, the wave propagation problem is not in-

cluded. Although many researchers have studied the theological proper-

ties of mud by using this model, none has been able to solve the equa-

tions of motion of a Bingham material under an oscillatory force.

Krone (1965) used a fairly extensive experimental procedure to

measure the rheological properties, i.e., yield strength and apparent

viscosity, of cohesive sediment. The maximum sediment concentration in

his experiments was about 110 g/l. The period of consolidation was not

mentioned. However, it can be inferred that he did not allow the

sediment to settle down because he remolded the mud quite frequently

to keep the sediment in suspension. He found that the yield strength

for all the samples in his experiments was proportional to the 2.5th

power of sediment concentration. See Fig. 2-7. The apparent visco-

sity for mud, pm, was 1 to 10 times higher than those for water, 1w.

See Fig. 2-8.

O San Francisco Bay Mud
() Wilmington Mud
B runswick Mud
10 Koolinite Suspension

10 10


16 _

Fig. 2-7. Relationship between Yield Strength and Sediment
Concentration (Data from Krone, 1965; Gularte, 1978;
Engelund and Wan, 1983).


@ Son Francisco Bay Mud
(Z) Wilmington Mud
() Brunswick Mud
@ Grundite
(5) Koolinite Suspension




'0 200 400 600 800 1000
Fig. 2-8. Relationship between Apparent Viscosity and Sediment Concentration.
(Data from Krone, 1965; Gularte, 1978; Engelund and Wan, 1983).





Migniot (1968) also made a series of laboratory experiments to

study the physical properties of cohesive sediments and their behavior

under hydrodynamic forces. He observed that muds have an "initial rigi-

dity" and that it is proportional to the 5th power of sediment concen-

tration, ranging from 120 to 800 g/l. The initial rigidity was defined

as the yield strength extrapolated from a curve (not a straight line)

which fitted the measured data on shear stress versus rate of shear

strain. The apparent viscosities were small for low sediment concen-

tration but increased in significance when the sediment concentration

became greater than 200 g/l. For example, m 20 pw when the sediment

concentration was around 200 g/l. However, pm 150 Vw when the concen-

tration went up to 500 g/l.

Gularte (1978) used a vane viscometer to measure the viscosity

and yield strength of grundite, which is an equal mixture of illite

and silt. The sediment concentration ranged from 600 g/l to 1000 g/l.

He found that the Bingham yield strength was proportional to the 10th

power of sediment concentration. See Fig. 2-7. The apparent viscosity,

pm, was 1000 times higher than pw. See Fig. 2-8.

Englund and Wan (1984) studied the instability of the hyper-

concentrated flow. They found that the Bingham model is quite satis-

factory for the description of the instability of the fluid surface

elevation. They found the yield strength to be proportional to the 3rd

power of sediment volume concentration. The correlation between Um and

sediment concentration is also shown in Fig. 2-8.

Parker and Kirby (1982) pointed out that, although mud suspensions

have been characterized as Newtonian fluid at low solid concentration,

to Bingham plastic at high concentration, the true behavior should

follow a pseudoplastic model, and that the yield strengths obtained by

many researchers are really extrapolations (to zero shear strain rate)

based on measurements at higher shear strain rates. They also presented

data from Williams and James (1978) to support their view point.

To summarize the afore-described studies on the viscoplastic

model, some conclusions may be stated as follows:

1. The yield strengths and apparent viscosities observed vary signi-

ficantly. No unique, correlative relationships have so far been


2. Mud suspensions at high sediment concentration do not follow

Newtonian fluid behavior. Although a Bingham model is commonly

suggested, a pseudoplastic model may be more appropriate.

3. In previous studies, attention does not seem to have been paid to the

effects of self-weight consolidation. This may be important, how-

ever, because muds possess a shear strength and a behavior which

approaches that of an elastic solid. Field-measured acceleration

data (in two horizontal and one vertical direction) of mud do not

seem to show the quiescence-moving-quiescence behavior predicted

by the Bingham fluid model (Tubman and Suhayda, 1976). Instead, a

continuously varied acceleration, which implies an elastic behavior,

is observed.

4. The Bingham model predicts no motion in the mud whenever the shear

stress is less than the yield strength. Therefore, no energy loss

can be predicted under this condition. However, unless the equation


of motion for a Bingham fluid can be solved, it would be difficult

to say whether there is a movement or not.

5. To further examine the suitability of the Bingham model for the

wave-mud system, solving the velocity field for both water and mud,

even one-dimensionally, is important. The velocity and associated

energy dissipation data can provide information for the justifica-

tion of using this model. Efforts to evaluate the one-dimensional

water-mud system have been made and the results are discussed in

Chapter 3.

2.3.6 Viscoelastic Model

Muds can have elastic properties when the sediment concentration

exceeds a specific value. Golden et al. (1982) reported that kaolinite

suspensions can have a shear strength when the solid fraction exceeds

0.04 (sediment concentration over 107 g/l) and it increases markedly

with higher sediment concentration. Mud beds usually have a sediment

concentration greater than 200 g/l. Therefore, elasticity can not be

neglected. Furthermore, in order to account for the energy loss due to

mud motion, a dissipation function must be included. Therefore, visco-

elastic models may be used.

Viscoelastic models account for the energy loss in an elastic

solid in an indirect manner. It is assumed that the restoring forces

are proportional to the amplitude of vibration, while the dissipative

forces are proportional to the velocity.

Two simple linear models, namely those based on the VoIgt and

Maxwell elements, are available. These are based on the analogy of a



a) Voigt Element

2G 2.

b) Maxwell Eement

Fig. 2-9. System Diagram for Simple Viscoelastic Models and
Response under Constant Loading.


e (rad,/,ec)
o .206
2 o .0265
a .0122
F- .00135
0I I I
0 .02 .O1 ,rJ .08
Fig. 2-10. Viscoelastic Behavior of Clay (after Stevenson, 1973).
The rate of angular displacement is denoted by .


spring-dashpot system. They are plotted in Fig. 2-9 and the constitu-

tive equations are presented as

Voigt element: T' = 20 E + 2P E' (2-4a)

Maxwell element: T' + p/G T' = 2p E' (2'4b)

where T and E are the stress and strain tensor, respectively, the prime

indicates the deviator part of those two tenors, the dot indicates the

derivative with respect to time, and G and p are constants. Figure 2-9

shows the shear strain as function of time under a constant shear load.

For the Voigt element, the response of elastic shear strain is delayed,

but it finally reaches the value that the spring alone would reach.

However, the Maxwell element will keep creeping because of the nature

of the dashpot. The Maxwell element is often used to describe the

relaxation of shear stresses under a constant strain. In the Voigt

element, the spring and dashpot are subject to the same displacement

but different stresses. It is easily shifted between a fluid and an

elastic solid and the behavior is frequency-independent. The Maxwell

element, however, behaves as an elastic solid at a high frequency load

and as a fluid at a low frequency load.

Hsiao and Shemdin (1980) studied the wave decay problem by assuming

that soft muds have Volgt properties. They solved the linearized two

dimensional equation of motion for the Voigt element given by Kolsky

(1963) and assumed the overlying water to be an inviscid fluid. Kolsky


2 2 a 2 .
p ax 3 G + (
P t - 9t x x 2 ( at (2-5


where i is the displacement in the xi direction, t is the time, p is

the pressure, and p is the material density. Hsiao and Shemdin (1980)

gave results showing a rapid wave energy dissipation.

Macpherson (1980) also studied the wave attenuation problem by

assuming that muds have the properties of a Voigt element. He had a

different treatment of Eq. 2-5 and arrived at a linearized Navier-

Stokes equation but with a complex viscosity. His apparent viscosity

for the mud is given in Eq. 2-6a. Here he represented the shear modulus

as the imaginary part of the apparent viscosity term and kept the

kinematic viscosity as the real part. Although his model could not

predict a bed shear stress because he neglected the viscosity of the

water layer, his approach greatly simplified Eq. 2-5.

Voigt element: ve V + jG(2-6a)
(1 0
Maxwell element: ve ( + j P- (2-6b)

Following Macpherson's approach for an oscillatory load, the

apparent viscosity for a Maxwell element is also given in Eq. 2-6b.

However, unlike the Voigt element, the apparent viscosity is frequency-

dependent. For an infinite G, the Maxwell element behaves as a viscous

fluid but the Voigt element becomes a rigid body. For zero V, the

Voigt element is a perfect elastic material while the Maxwell element

is not defined. Also the Maxwell element is not defined for zero G;

however, the Voigt element is a viscous fluid at this condition.


Questions concerning which element simulates mud behavior more

closely are not discussed by Macpherson (1980) or by Hsiao and Shemdin

(1980) and still remain unclear. Stevenson's (1973) data, see Fig.

2-10, showed that the Voigt element was applicable for his soil sam-

ples. However, his samples had a high shear strength and were not like

a soft mud. Therefore, experiments to observe the response of mud

under constant shear loads were conducted in the present study. The

objective was to find whether the Voigt model or the Maxwell model was

better in modeling the stress-strain relation, and also to find the

material constants G and p. Details are given in Chapter 5.

2.4 Bed Shear Stress

The bed shear stress is well known as the most important factor in

the erosion process in steady flow as well as in an oscillatory wave

flow. Unlike steady flow, for a smooth bed (ratio of the bed roughness

r to the semi-excursion distance just outside the boundary layer ab,

r/ab < 0.01), the wave bed shear stress has a period which is basically

equal to the wave period. For a rough bed, r/ab > 0.01, higher order

harmonics may be induced by the rhythmic formation and release of vor-

tices (Lofquist, 1980). Most papers dealing with bed shear stress in

the oscillatory flow assume that the bed is rigid and that the water is

free of sediment. This is a reasonable assumption for a sandy bed or a

clayey bed with low water content since there is no appreciable motion

in the bed. However, for a cohesive sediment bed with high water

content, this assumption is not always applicable because the bed may

also move (Tubman and Suhayda, 1976; Nagai et al., 1982). This would

mean that the bed shear stress would be different because the velocity

gradient is altered at the interface. Therefore, to accurately estimate

the bed shear stress one must consider the dynamic response of mud as

well. In this section, only the wave bed shear stress for a rigid bed

is briefly summarized.

Jonsson (1966) proposed the following equation to evaluate the

maximum bed shear stress (which is one of the most popular equations

still used):

1 2
= 2 p fw Ub (2-7)

where ub=abu is the maximum particle velocity just outside the bottom

wave boundary layer which can be predicted by classical wave theory, and

fw is the wave friction factor. This is a simple form and the only

unknown is fw.

The wave friction factor is a function of the surface roughness

and wave Reynolds number Rw, defined by Jonsson (1966):

Rw= u2/ va (2-8)

Here v is the kinematic viscosity of the water. For the case when

Rw < 1.25 x 104, the viscous wave boundary condition, Jonsson (1966)


fw= 21 Aww (2-9)
w w

At higher wave Reynolds numbers, the bottom roughness is involved.

This is the transition zone between viscous and full turbulent flow.

Kamphuis (1975) presented a chart (Fig. 2-11) to evaluate fw. For

.- Rough Turtbient 0"

- 80

r --- 40---
N. ---80
Z -Ierpoloted 20
o ld----Extrpolated 400 i
r Lowen Limit of Rough %be
** l lt I I ,itt*l g *pttl| llllplI a i flit
.21 I0 1 1 0,1 1 1 111 1 1 0111 C


Fig. 2-11. Wave Friction Factor (after Kamphuis, 1975). The
grain size, which is larger than 90% of the total
mass of soil sample, is denoted as D90.


higher Rw, only the roughness is important. Jonsson (1966-, 1975, 1980)

also gave a formula to evaluate the wave friction factor under this


2.5 Wave Diffusion Coefficient

Fisher (1968) pointed out that the concept of the apparent diffu-

sion coefficient Dm (turbulent diffusion and dispersion) has been

widely employed to predict the transport of pollutants. This concept

has been used for sediment transport as well. However, the problem of

evaluation of the diffusion coefficient for a non-breaking wave environ-

ment has not been fully clarified yet. Previous studies relative to

this coefficient, although having no consistent results, are briefly

summarized here.

Kennedy and Locher (1972) gave a good summary of the early

studies. Key points are given here. The pioneering experimental work

done by Shinohara et al. (1958) showed that a constant diffusivity

is good for sand with a mean diameter of 0.2 mm. However, specifying

only one constant diffusivity for pulverized coal with a mean diameter

of 0.3 mm was not sufficient. This was because of the discontinuity in

the concentration profile near the bed. Hattori (1969) also employed

the constant diffusion concept and showed good agreement with his

experiments for the top 80% of the water column.

Das (1971) concluded, from his oscillating bed experiment, that

the diffusion coefficient varies linearly with distance above the bed.

Hom-ma and Horikawa (1963) gave Eq. 2-10 for the diffusion coef-

ficient under a wave field.

D 1 Ho sinh3kz"
m 2B sinh kb cosh 2 kz (2-10)

where is a constant related to the ripple geometry, H is wave height,

h is water depth, k is wave number, z" is elevation above bed, and c is

wave celerity. This formula uses the velocity profile from the small

amplitude wave theory and, therefore, the diffusion coefficient in-

creases slowly upward since the velocity profile does not vary sharply

with elevation. However, their sediment concentration data (for the

upper 70% of water layer) showed quite good linearity on a log-linear

plot which implies a constant Dm.

Horikawa and Watanabe (1970) measured the turbulent velocity

fluctuations, u' and w', rather than deduce the variation of diffusion

coefficient from the vertical distribution of the mean sediment con-

centration. They concluded that the root mean square (RMS) value of w'

is almost independent of z, whereas the RMS value of u', which is

greatest near the bed, decreases first just above the bed, and then

becomes nearly constant for the most part.

Jonsson and Carlsen (1976) measured wave velocity profiles and

then calculated the profile of eddy viscosity, c. Their flow had a high

Rw and a rough bed. Their results showed that the maximum E is located

within the turbulent wave boundary layer. Outside the boundary layer,

c decreases drastically and approaches zero for the upper water column.

Here c = 0 reflects the fact that turbulence is not important in the

calculation of velocity for wave flows. However, turbulence is indeed

important in sediment transport. For cohesive sediment, momentum

transfer coefficients are almost the same as mass transfer coefficient,

diffusion coefficient, (Jobson and Sayre, 1970). Therefore, it can be

inferred that the turbulent diffusion coefficient should be also the

largest near a rough bed.

Wells et al. (1978) presented Eq. 2-11 as "the most justifiable

estimate for suspended-sediment concentration" along the Louisiana

coast as a result of wave resuspension:

C(z)= 5.31 exp(-1.8 z") (2-11)

where C is the sediment concentration in kg/m3 and applies to a water

depth h < 5 m. Their result implies a constant diffusion coefficient

since there is a log-linear relationship between C and the vertical

coordinate, z".

Hwang and Wang (1982) summarized the models available for the

turbulent diffusion coefficient outside the wave boundary and suggested

that a modified version of their model, Eq. 2-12, would be the most

plausible one:

D =oH2a sinh2kz" (2-12)
m 2 sinh2kh

where a is a constant.

Thimakorn (1984) also gave an equation for the diffusion coef-

ficient profile, which was close to that given by Hwang and Wang (1982),

based on the concept of turbulent energy.

The above information is summarized in Table 2P1. The wave

Reynolds number, Rw, which indicates the flow regime near the mud sur-

face, is also included. Notice that some of the laboratory studies


have had a wide range of Rw which indicates that both the yiscous and

turbulent wave boundary layers were included. Two types of diffusion

coefficient profiles have been introduced. Data from Shinohara et al.

(1958), Hattori (1969), Horikawa and Watanbe (1970), and Wells et al.

(1978) suggest a constant coefficient for the top 80% of the water

layer, the remainder apparently having another constant value. However,

to the contrary, the others (Hom-ma and Horikawa, 1963; Das, 1971;

Hwang and Wang, 1982; Thimakorn, 1984) suggest a minimum diffusivity at

the bottom and increasing upward.

Table 2-1. Diffusion Coefficients in Wave Flow

Researcher Sediment Rw Suggested Dm

Shinohara et al. Sand, < 6700 Constant

Hattori Sand < 24000 Constant

Hom-ma &
Horikawa Sand < 26400 Eq. 2-10

Horikawa & > 12000 Constant

Wells et al. Mud > 12000 Constant

Hwang & Wang Sand Eq. 2-12

Thimakorn Mud < 10000 similar to
Eq. 2-12

It is interesting to note that none of these studies has mentioned

a possible longitudinal dispersion phenomenon. In actual shear flow,

fluid layers near the bottom tend to exhibit a sharp velocity gradient.

This results in a much greater longitudinal dispersion. Taylor (1954)


gave the example of steady uniform pipe flow and showed that the disper-

sion coefficient can be about 200 times larger than the mean value of

the turbulent diffusion coefficient. However, Awaya (1969) pointed out

that in an oscillatory flow the dispersion coefficient can not be that

high, because of the limited time-scale associated with this type of

flow, tidal flow in his case. It can still be inferred, however, for

wave-induced flow, that there must be a strong longitudinal mixing

layer above the bed surface because of sharp velocity gradients in the

bottom boundary layer. It is further noted that the bed surface is the

main source of turbulence (Nielsen, 1984). This may explain what was

observed by Shinohara et al. and others mentioned. To further under-

stand the wave diffusion behavior, dye diffusion tests were conducted

in the present study. Suspended sediment concentration profiles were

also collected. Details are discussed in Chapter 6 and Chapter 7.


3.1 Introduction

In this chapter a numerical solution of the velocity field and

energy dissipation in a vertical one-dimensional water-mud system is

presented. The mud is assumed to have Bingham plastic properties with

constant density but variable viscosity and yield strength. The over-

lying water is assumed to have constant density and viscosity. The

objectives are to examine the velocity profile and the associated

energy dissipation in this system, and to justify this 1-D model for

simulating the water-mud system.

3.2 Problem Formulation

Figure 3-1 shows the system diagram. A variable grid system, with

grid size Azj, was chosen to give better resolution at the bottom, the

interface, and the free surface. The grid number starts from I at the

rigid bottom, and ends at the free surface with number J. The total

depth is h. The fluid properties, e.g., density p, kinematic viscosity

v, and yield strength K, as well as the horizontal velocity u, are

specified at the center of each grid element. Shear stresses, T, are

specified at the node points. The one-dimensional linearized equation

of motion is


._Azj I
kJi T_ -

Fig. 3-1. System Diagram for the 1-D Numerical Bingham
Fluid Model.

du dI 1 dT (3-I)


+ K du (3-1a)

where dn/dx is the slope of the water surface profile, t is time, g is

gravity, and z is the vertical direction. The boundary conditions are

as follows:

1. no-slip at the rigid bottom, u=0 at z=O.

2. velocity is maximum at the free surface, du/dz=0 at z=h, or,

3. The initial condition says that the fluid is at rest at the

beginning, u=0 at t=O.

To ensure that the flow velocity is maximum at the free surface, a

fictitious layer (from J to J+1) is added at the top of the water

column. Equation 3"1, with the boundary conditions, can be solved by

the double sweep method (Abbott, 1980). The finite difference formula-

tion is

n+1 n n+1 n+1 n n
uj Uuj = dn )n+ /2+ 1 (P+ lA 1= ) + .L ( J+ I J)
At dx 2pj Az. 2p. z (32)

k n 1 un+1
n+1 j1(Az++ Az+1) 2 uj+1 i (3-2a)
2 luj+1 -u I CAZ + Azj+I)

n+1 n+1
k. (Azj+ Az _) 2 (u. un+I
In- =[P + n j-1 u j-1 (3-2b)
2 jun ju11 (AZj + AZjI)

Notice that the velocities at time n are also used to evaluate

the velocity gradient in the denominator for the shear stress term at

time n+1. The superscript n+1/2 implies the center of each time step.

Equation 3-2 can be rearranged in terms of un+1 and rewritten as

n+1 un+1 n+1
Au .B uj. = D (3-3)
J1+ .j j-1

K (Azi + Az +I)
Au-+1 _un 1 I (Az ) 7Z. (3-3a)
2 + AZ i Z

K 1(Azj + Azj I)
B = 1 + ( j + I + ) j A t
K. ( Azj

+(.j 3 Azj1) At
2 n n p (Az + Az 1 ) Azj (3-3b)
2 1 uj un_1P1 (zi+A J1A

K.(Azj + Azj1 ) At
2- (uj n _u n P (Azi + Az j) AZ (3-3c)
2lu n I jAj z

D= -gAt (dnfn-1/ +

Kj+ (Azj + Azj+1 At Un
(Pj+I 2 Jun1 n pj(Azj + Az ) Az. j+1
j j j j

Kj (AZj ) Az
+U 1 -uj Aj+
+Kj(Az + AZjI At U
+ 31-(1+1 + 2 Jun+ unI P .(Az. + Azj~ ) AZ

(Ij + nuKj(Azj + AzA]
2 u n -u11 PJ(AZ + AZ1) A1 j

Kj (Azj +z
2 n U (3-3d)
2 lu' -ujI pj(Azj + Azj I) Azj j- 3-d

Now introduce the auxiliary equation, Eq. 3-4, and substitute into

Eq. 3-3 to obtain Eq. 3-5

n+1 n+1
u+1 =E i+ F (3-4)
uj+I = j 3 3

n+1 -C n+1 (D A F.)
j A E.+ -B uj-1 + (A E j+ B)

=Ej_ ujI + IF (3-5)
i- J-1 J-i

Comparing Eq. 3-4 and Eq. 3-5 reveals that Ej and Fj, from j=J-1 to 1,

can be obtained because Ej=1 and Fj=O are specified by the zero velocity

gradient condition. Therefore, the first sweep yields all the Ej and

Fj. The second sweep, from Eq. 3-4 and the no-slip boundary condition,

u1=O, solves for the velocity for all the grid elements.

Because there is a velocity gradient term in the denominator

(Eq. 3-1a) for evaluating the shear stress of the Bingham fluid, a

problem of calculating A, B, C, and D was encountered because of the

zero velocity gradient at the onset. However, letting the first sweep

stop at the uppermost element, S, where the yield strength is zero,

and moving the no-slip condition to the element one cell below element

S, allows calculation of the velocity profile above this element. By

comparing the shear stress acting at the bottom of the moving fluid and

the yield strength one layer below, a decision can be made as to whether

the next layer will move or not and therefore to stop the first sweep


for the next time step. If this is done, the first sweep can go down

step-by-step to a cell where the shear stress and shear strength are

equal, or to the rigid bottom provided the driving force is large


Wave energy components are evaluated separately as kinetic energy

and potential energy. The time-averaged total energy per unit surface

area is

E= -4 pgdz dt + 1 T hpu2dz dt


a a exp( 2k.X) (3-6)
0 m
where a is a constant coefficient, kim is the wave decay coefficient,

and ao is the wave amplitude at x = 0. The time-averaged energy dis-

sipation per unit surface area is calculated as follows:

I T h du zd
Ed 1 0 hd dz

1 n n n n
= I I (T j1 u + -/2 u I/2)]At (3-7)
n-1 j=1 J-

where T and u are given by the model and N is the total number of time

steps for a complete wave cycle. The energy dissipation can also be

written as

Sc DE = 2 cga kima exp(2k X)
d ax g 0in im

= 2 cgkimEn (3-8)
where Cg is the wave group velocity and can be replaced by vg- because

of the long wave assumption. Therefore, the decay coefficient can be

obtained as

k im E / 2V9 En) (3-9)

3.3 Results

The numerical scheme was tested for the steady-state solutions

first, simply because analytical solutions are available. The water

surface gradient was arbitrarily selected to be 10-5. The time incre-

ment was also arbitrarily chosen as 50 seconds since the implicit

scheme is unconditionally stable. Figure 3-2a shows the velocity

profile for a two-layered viscous fluid at selected times. Squares

show the analytical solution. Figure 3-2b shows the velocity profile

for an assumed fluid which has exponentially distributed yield stress

and viscosity. The results are satisfactory in both cases.

The pressure forcing term (dn/dx) in Eq. 3-1 was then changed to

a cosine function to simulate a periodic wave load. The time increment

was also reduced to 0.05 second because the wave period was relatively

short and also to minimize possible error introduced by the step-by-step

downward computational procedure. The scheme was compared with the

results from Dalrymple and Liu's two-layered viscous fluid model and

showed reasonable agreement. See Fig. 3-3. The observed discrepancy

was due to the long wave assumption in the numerical model. In this

example, the long wave assumption could not account for the viscosity

variation at the interface. This may be interpreted as follows. Long

waves are characterized by a constant dynamic pressure in the entire

water depth (i.e., pressure is hydrostatic) and by the absence of

t=166 min t- 7 mi
66 mi

n t=332 min t=75 mi
t=500 min tal50 mi

15 t=666 min

LwA 10

a: Top Bottom K=0 I N/m2
S p(kg/m3) 1000 1160 ixlO' m2/s
ci pl60 kg/n3
(a) v(m2/s) x10-6 2xO-4 (b)
0 s0 oo SO 200 250 0 20 40 so 80 100


Fig. 3-2. Horizontal Velocity Profiles for Two Steady Flows. (a) Two;Layerd Viscous
Fluid; (b) Bingham Fluid with Exponentially Increasing K and U.

W 2

'15 iInterface Elevation_

W 10
() I
z I

05 I H T
7 Legend Model (cmXsec
I-D 2.5 3
-DaL Z53
-- D L 5 40
5 10 15 20 25

Fig. 3-3. Comparison of Wave Velocity Amplitude Profiles from
1-D Model and Dalrymple and Liu's Model with P=11000
kg/m3, P2=1160 kg/m3, \i=IxIO-6 and v2=2x10-4 rm2/s.


vertical velocity. For cases where the viscosity of the bottom layer

is not too large, the viscous boundary layer, of thickness 6L, is far

below the interface. Then the equation of motion, for the part outside

6L, can be simplified by neglecting the shear stress term, the last

term in Eq. 3-1. Although the shear stress in the bottom layer is

larger than in water, because v2 > v1, it is still small when compared

with the pressure and inertia terms. Therefore, a single equation can

be used from the water surface down to a particular elevation which is

below the interface but above the viscous boundary layer. This phe-

nomenon was demonstrated by another run of Dalrymple and Liu's model

(see Fig. 3-3) but with a long wave period, T = 40 sec.

A simulation of water-Bingham fluid system was then performed.

Input data were specified as follows. The water layer had a constant

density (1000 kg/m3) and viscosity (1 x 10-6 m2/sec). The bottom

Bingham layer also had a constant density, 1160 kg/m3, viscosity, 0.0002

m2/sec, and yield strength, 0.1 N/m2. Fig. 3-4a,b shows the velocity

profiles under two wave conditions.

It was expected that the low layer would not move for a small wave

height (Fig. 3-4a). This is because the bed shear stress for this wave

was less than the yield strength. For larger waves, the lower layer is

observed to move with a zero velocity gradient in much of the lower

layer, indicating a plastic sliding. Nevertheless, the results show

that the velocities in the mud bed were either zero or close to that in

the water layer.

Energy dissipation and total energy at each time step were evalu-

ated. Time-averaged values over one wave period were calculated next.

tsw LWL - -- -


0tS aat-
Q t6 a6a
c~20 1!


o Interfacei 7T


Cr 2

(a) (b)

-5 0 10 is -5 0 5 10 is 20


Fig. 3-14. Wave Velocity Profiles for the Water-Bingham Fluid System with T =3 see.
(a) H -2.5 cm and ki,= 0.0011 m'1; (b) H = 5 cm and kim=0.016 m'1.

The decay coefficient was then obtained based on Eq. 3-9 and is given

in Fig. 3-4.

Yield strength and dynamic viscosity data were obtained from

viscometer measurement (see Chapter 5 for detail) for a kaolinite bed

with a two-day consolidation period, which was the softest mud bed in

the wave erosion tests (described later in Chapter 6 and Chapter 7).

Results are given in Fig. 3-5. The yield strength and dynamic viscosity

varied from 2.5 to 4.4 N/m2 and 0.1 to 0.2 N-s/m2, respectively. The

corresponding kinematic viscosity was on the order of 1 x 10-4 m2/s.

Because of the relatively high yield strength, the model-simulated mud

bed showed no motion at all. See Fig. 3-6. The energy dissipation was

significantly low because there was no motion in the bottom layer.

3.4 Conclusion

It must be realized that muds can dissipate energy only if they

are moving. The 1-D numerical model presented here shows no motion in

the mud layer, therefore no energy dissipation, when the yield shear

strength is higher than the bed shear stress. However, such a condition

of no motion no dissipation was not observed in experiments described

later, and thus indicated that the 1-D Bingham model is not generally

applicable for the simulation of watermud interaction in erosion


Leptd Elev. Below Mud
Surface (cm)
o 1.5

w + 9.3
1 K=4.4N/n,2


2 I I I I I I
2 4 6

Fig. 3-5. Shear Stress versus Shear Strain Rate for Kaolinite
Bed with Two-Day Consolidation Period.

10 1



i .



-5 0 5 10 15 20

Fig. 3-6. Wave Velocity Profiles for Water-Bingham Fluid System
with Measured Yield Strength and Viscosity.



4.1 Introduction

A multi-layered hydrodynamic model to simulate the two important

parameters, wave damping coefficient and interracial shear stress, is

presented in this chapter. The water layer was assumed to be charac-

terized by a constant kinematic and eddy viscosity. A simple linear

viscoelastic behavior, the Voigt element, and the linearized equations

of motion were selected for modeling mud response under wave loading.

This selection was based on the discussion presented in Chapter 2,

Chapter 3 and the results from Chapter 5.

Depth variation of mud properties can be simulated by further

dividing the mud bed into layers. The equations of motion for the

water layer are also linearized to simplify the problem. The water

and mud layers are considered incompressible and immiscible such that

density, viscosity, and shear modulus are constant in each layer.

4.2 Formulation

Figure 4-1 is a definition sketch of the water-mud system. Water

waves propagate in the x-direction. The water layer has a thickness,

dl, and the underlying mud bed can be considered as composed of several

layers of linear viscoelastic material with thicknesses di. The free



Pi ,vi ,Gi diI

Fig. 4-1. Schematic Figure for the Multi-Layered Model.

surface displacement is n. Here i 2, 3, ..., n. The displacement at

each interface is i. The displacements n and i can be expressed as

n = ao exp[ J(kx-ot)] (4-1a)

i= boi exp[ j(kx-ot)] (4-1b)


where ao is the given water wave amplitude at x = 0, boi a e unknown

complex variables for interface wave amplitudes and phases at x = 0,

J /CT, o=2ir/T, T is the wave period, t is time, and k is a complex wave

number. The real part, kr=27/L, is related to the wave length, L. The

imaginary part, kim, represents the decay coefficient defined as

a(x) = aoexp(-kimx) (4-2)

This expression correctly represents the wave decay in a flume

(Schuckman and Yamamoto, 1982; Nagai et al., 1982).

The linearized equations of motion for the incompressible upper

fluid and incompressible, viscoelastic Voigt element layers under

oscillatory load are

u a t a 2 U 2u
__. = 1 + Vei( ) (4-3a)
at P 3x ei x2 z 2

t 2i 2i
Dw 1 1 Pi + e ( -i + -i )(4-3b)
at -P -z ei x2 3z2

where u and w are velocities in x and z direction, respectively. The

apparent viscosity, vei, actually contains two terms. For the water

layer, it contains the kinematic viscosity, vI, and a constant eddy

viscosity, e. The purpose of involving e in the model is for the bed

shear stress calculation. The approach presented later demonstrates

that this choice of e is not connected with the turbulent diffusion

coefficient, Dm. Each mud layer is characterized by a kinematic

viscosity, vi, and a shear modulus, Gi, as discussed in Chapter 2. The


subscript i1, 2 ..., N indicate the top, second, and subsequent

layers, respectively, pt is the total pressure defined as

t 0
Pi = Pi Pigz Pi


0, if i = 1
o ii iil (4-4)
Pi I(mP-)g E },d] if i > 1

where m and r are dummy indices. The continuity equation is

Sui w i
ax + = 0 (4-5)

The solutions of these variables are assumed separable according to

ui(x,z,t) = u1 (z) exp( j(kx-at))

wi(x,z,t) = wi(z) exp( j(kx-at)) (4-6)

Pi(xz,t) = pi(z) exp( j(kx-at))

Substititing u i wi into the continuity equation one obtains

u = J ~ (4-7)
i k i

where the prime indicates differentiation with respect to z. Intro-

ducing this equation for u, into the horizontal momentum equation

(4-3a) yields an expression for pi, i.e.,

'ei w 2
w'A ) (4-8)
;i k2 1 ii
2 2 -1
Ai k jvei

Substituting P into Eq. 4-3b yields


2 2 22
k X + k A wi 0 -(4-9)

The solutions for the water layer and the subsequent mud layers are
in the form:

w1(z)=A sinh Z, + B10osh Z+ C1exp(A1z) +D1exp(-A1(z+dl)) (4-10a)

wi(z)=A sinh Zi + Bicosh Zi + Cisinh Ni + Dicosh Ni (4-l0b)

with i=2,3 ,...,N, Zi-k(d1+d2+...+di+z), and Ni=Ai(d1+d2+...+di+z). For

only one mud layer, there are eight unknown complex coefficients (Al,

...,D2) for the vertical velocities and the two unknown complex varia-

bles, i.e., wave number, k, and interfacial wave amplitude bo. There-

fore, 10 boundary conditions are required for solving for the 10 un-

known coefficients. The 10 boundary conditions will be summarized here

to show how the coefficient matrix is constructed. Details are shown

in Appendix A.

At the free surface there is no external driving force, i.e.,
2w (2u 2w
normal stress, an + 2.e a, and shear stress, T = p L + L).

are zero. Because the actual free surface is unknown a prior, these

two boundary conditions are applied at the mean or still water surface.

Taylor expansion technique is applied as shown in Appendix A to obtain

an approximate solution. Only the basic harmonic term is considered.

All high harmonic terms are neglected. The zero normal stress condi-

tion at the free surface gives
M ( A cosh kd1+ B1sinh kd ) 2veIX1C = plgao (4-11)

The linearized kinematic free surface boundary condition yields

A 1sinh kd1+ B 1cosh kd1+ CI -joa0 (4-12)

Equations 4-11 and 4-12 are combined to give Eq. 4-13

c1,1 A,+ c1,2 B, + c1,3 C = 0 (4-13)

The coefficients cij are given in Appendix A. The zero shear stress

condition at the mean free surface gives

c2,1 A, + c2,2 B, ++c2,3 CI = 0 (4-14)

At the water-mud interface E, (a prior unknown elevation), the hori-

zontal and vertical velocities for the water and mud layer are matched.

See Eq. 4-15 and 4-16. The linearized kinematic boundary condition at

the interface gives Eq. 4-17. The matched normal stresses and shear

stresses conditions between the water layer and the mud layer give Eqs.

4-18 and 4119. For details, see Appendix A.

c Ac Dc A2+c_ .c Cc 0 (4-15)
c3,1A c 3,4 D1+ 3,5 2 3,6B c3,7 C2+ c3,8D 2

c 4,2Bi+ c 4,4 c 4,5A2+ 4,6B2 c 4,7C+ c 0 4,8D2 0 (4-16)

c5,2 B1 + c5,4 D1+ c5,9 b 10 0 (4-17)

c6,1A, + c6,4D+ c 6,5 A2 + c6,6 B2+ c6,7C2 +

c6,8 D2 + c6,9 b 10 0 (4-18)

c7,2B+ c7,4D+ c 7,5 A2+ c7,6 B2+ c7,7 C2+ c7,8D2 0 (4-19)

At the rigid bottom (z=-dl-d2), the no-slip condition (u-w=O0) must be

satisfied, i.e.,

08,5 A2 + c8,7 C2. 0 (4-20)
c9,6 B2 + c9,8 D2= 0 (4-21)

These nine equations, Eq. 4-13 through Eq. 4-21, can be written as


a homogeneous matrix equation c X = 0, where c is a 9 by 9 matrix

with non-zero terms specified from Eqs. 4-13 through 4-21. Here X is a

column matrix which contains coefficients A,, BI, C1, D1, A2, B2, C2,

D2, and bol. The unknown complex wave number k is involved in the

coefficient matrix c. This is an Elgenvalue problem since the system

has solutions only for special k values such that the matrix c is

singular, i.e., the determinant of matrix c is zero.

Adding one more mud layer in the bottom gives five more unknowns:

four coefficients (A3, B3, C3, D3) in the added vertical velocity equa-

tion, and one (bo2) for the unknown new interracial amplitude and

phase. Five more boundary conditions must also be specified similar to

Eq. 4-15 through 4-19. The newly added layer should be specified as

the lowest one, directly above the rigid bottom. Then, the added

equations can be listed as

c8,5 A2+ c8,7C 2+ c8,10A3 c 8,11 B3 + c8,12 C3+ c8,14D3' 0 (4-22)

c9,6 B 2+ c09,8D C9,10 A3+ c9,11B 3 c9,12 C + c9.3 D 3 0 (4-23)

c0,6 B2 + c10,8D 2 c10,14 b2o' 0 (4-24)

c11,5A2 + c11,7 C2+ c11,10 A3 + c11,11B 3 + c11,12C3 +

c11,13 D + c11,14 b2o 0 (4-25)

c12,6B2 + c12,8D 2+ c12,10A 3 + c12,11B 3+ c12,12C3 +

c 12,13D3 0 (4-26)

The no-slip condition, now at z=-dl-d2-d3, is replaced by

c13,10 A3 + c13,12 C3 0 (4-27)

c14,11 B3 + c14,13 D3' 0 (4-28)


The size of the new coefficient matrix is then expanded to 14 by

14, and five more coefficients (A3 D3, and b,2) are appended to the

column matrix X. The homogeneous matrix equation, however, remains in

a similar form, and the same technique may be utilized to solve it.

When adding one more layer, matching conditions at the new interface

and the no-slip condition at the rigid bottom are similar to Eqs.

4-22 through 4,28 because of the same equation for w. However, the

subscripts are increased by 5. The matrix c expands to 19 by 19.

Figure 4-2 shows the layout of matrix equation, c X = 0, for N layers.

The calculation of matrix c may be performed in an iterative manner

(a DO-loop in the computer program). This approach allows greater

flexibility to add layers in meeting the depth variation of mud proper-

ties. Any number of mud layers may be assigned provided this is ne-

cessary. Basically, the fewer the layers, the faster the rate of con-


4.3 Solution Technique

To find the complex wave number k for a zero determinant, the

Newton-Raphson method was used. The principle of this method for a

complex function is the same as that for a real function. The first

guess of k value used is the wave number obtained from Airy's (linear)

wave theory for a water depth dl. The second guess is arbitrarily set

to be 1% off from the first one. These two values of k are real

numbers, but then k shifts to a complex number because the coefficient

matrix is a complex function. Examining the trajectory of the complex

wave number and the value of the determinant suggests that the function

xI I I I
x x x Ilx iX I

x x xix xix i i
x x X x x -
x x xx xl
x x KXX
X X X X___ _

x Xl X I -L -L

i~~ ~ I [ I lEI x

le Loyer 2ndLayer 3,d Layer



Fig. 4-2. Layout of the Coefficient Matrix for the Multi-Layered Model.

_Vj 11
thLayer Nth Layer


has a tornado-like shape. To hit the target (exact k value for zero

determinant), high precision in the computer program is required.

Sometimes, an absolute zero determinant is difficult to obtain due to

the large matrix and limited accuracy of the computer. Therefore,

another convergence criterion was desirable. Comparing the results,

e.g., velocity and pressure, obtained from a finite determinant value

but with Ak < 10-14 and that from a zero determinant showed no account-

able difference. Therefore, Ak < 10-14 was selected as the new con-

vergence criterion. This criterion requires a reasonable number of

iterations, i.e., less than 20, to obtain the results.

The elements of the column matrix X, after solving the matrix

equation, do not represent the true answer but are proportional to the

true value. Equation 4-11 (or Eq. 4-12) is invoked to obtain the pro-

portionality constant and then to obtain the true Ai,Bi,..., and bi-1.

Pressure and velocity profiles are then determined easily from Eqs.

4-7, 4-8, and 4-10. Shear stress profile is obtained as

= (u + 9w)
Pe~z ax

SPe(W" + k2 w) exp[ j(kx ot)] (4-29)

4.4 Input data

The model was run through all the experimental cases (see Chapter

6 for detail). Wave amplitude, period, and the thickness, bulk densi-

ty, viscosity, and shear modulus for each layer are the input data

needed. The number of layers and thicknesses di were selected to meet

the variation of measured density profile. For the water layer, the


density was selected to be 1000 kg/m3. The kinematic viscosity, vI,

was 1 x 10-6 m2/sec. The eddy viscosity, e, is dependent on the wave

Reynolds number and is discussed later. Details for measuring Gi and

vi for the selected muds are given in Chapter 5. In this study, G and

p for the six experimental runs were correlated with the dry density.

The correlative equations are presented in Chapter 5. See Table 5-2.

Bed density information follows in Chapter 7, see Eq. 7-1 through

7-3. Wave conditions for the six experimental runs are presented in

Chapter 6 (Table 6-4). The thickness of each layer for the six experi-

mental runs are listed in Table 4-1. Each run had two wave loadings

except Run 6 which had three. See Chapter 6 for detail.

An assumption concerning the eddy viscosity was that E (water

layer) is only a function of the water-wave parameters and is indepen-

dent of mud bed rigidity. In other words, the eddy viscosity in the

water column for a rigid bed was assumed to be the same as for a soft

movable bed, given the same wave conditions. This assumption allows

the determination and use of eddy viscosity in the model. Details are

described next.

Given a high shear modulus for the mud layer, e.g., G > 10000

N/m2, the model may be used to simulate the wave bed shear stress for a

rigid bed. By comparing the model-predicted bed shear stress and that

calculated from Jonsson's friction factor concept (using the friction

factor diagram presented by Kamphuis, 1975: see Fig. 2-11 in Chapter

2), a value of the apparent eddy viscosity, Vel = "I + E, can be esti-

mated. From this, c can be obtained given vI. Figure 4-3 shows an

example with d, 3 m, ao = 0.25 m, a = 0.79 rad/sec, Pi = 1035 kg/m3,

Table 4-1. Input Thickness for Each Layer (cm)

Run dI d2 d3 d4 d

1-1 21.65 2.5 2.5 4.0 5.0
2 24.15 1.5 3.0 3.0 4.0
2-1 19.15 1.5 3.0 3.0 4.0
2 19.65 1.0 3.0 3.0 4.0
3-1 16.15 2.5 3.0 4.0 5.0
2 18.15 1.5 2.0 4.0 5.0
4-1 26.35 1.3 2.0 3.0 3.0
2 28.65 1.0 1.0 2.0 3.0
5-1 19.65 2.0 3.0 5.0 6.0
2 21.05 1.6 2.0 5.0 6.0
6-1 24.65 1.0 2.0 3.0 5.0
2 25.15 0.5 2.0 3.0 5.0
3 25.15 0.5 2.0 3.0 5.0

Fig. 4-3. Comparison of Computed Bed Shear Stress (Model and
Kamphuis) at the Mud Surface to Determine the Eddy
Viscosity of Water.


vi = 1x10-6 m2/s, bed surface roughness r = 20 Om and Rw = 3.0 x 105.

Values of d2, P2, and P2 were immaterial because of the high G2 value.

With a fw from Kamphuis' diagram, Eq. 2-7 gives a maximum bed shear

stress of 0.6 N/m2 which is equivalent to the model having a total

viscosity of 2.4 x 10-6 m2/sec. Therefore, the corresponding eddy

viscosity would be 1.4 x 10-6 m2/sec and may be used to calculate the

interfacial shear stress for a soft movable bed under the same wave

load. This procedure illustrates a practical method for estimating the

water-mud interfacial shear stress for any wave flow condition. For

the six experimental runs, however, the value of c was set to be equal

to zero because all Rw < 104 (presented later in Table 4-2). This is

because Rw < 104 indicates a viscous sub-layer just above the mud

surface and consequently zero eddy viscosity for the purpose of calcu-

lating the bed shear stress.

In oscillatory flows, the pressure force is balanced mainly by the

inertia force. The shear force plays little role when outside the

bottom boundary layer, 6, even though turbulence is involved. This can

be understood by noting that a high value of Reynolds number, R=udl/vel

300 2000, implies that the shear stress term can be neglected

outside the boundary layer. This explains why the potential flow ap-

proach works quite well for water wave problems. For wave flows in-

volving viscosity, an example of the horizontal velocity profile given

by Dean and Dalrymple (1984, p. 264) indicates that the only influence

of increased vel is an increase in 6. The velocity profile above 6

does not change at all. Therefore, the value of c for the water layer

is primarily necessary for the determination of the bed shear stress.

4.5 Model Results

Model results in terms of velocities, pressure, and shear stress

as well as the corresponding phase angles for the six experimental

runs are summarized as follows. Details are presented in Appendix B.

Five layers were chosen in all runs.

4.5.1 Velocity

Two types of vertical water velocity profiles are obtained by

model simulation. In the first example (see Fig. 4-4a), the vertical

velocity decreased to a minimum a few centimeters above the mud surface

and then increased to a maximum at the mud surface. The vertical

velocity profile in the mud started from maximum at the mud surface and

decreased to zero at the rigid bottom. This profile tends to show that

there is a "false bottom" above the mud surface. The tendency to have

a false bottom is due to the vertical motion of mud being about 90

degrees out of phase with the surface wave. Mallard and Dalrymple

(1977) have shown mathematically that a false bottom can exist for a

perfectly elastic bed which has a vertical motion 180 degrees out of

phase with the surface wave. In the present example the viscosity

damps out the response and therefore only shows a tendency. In the

second type (see Fig. 4-5a for Run 5-1) the vertical velocity does not

have a phase lag between the water and mud layers.

Horizontal velocities were reduced significantly in the mud layers

and also shifted 90 to 120 degrees. This implies that the velocity

PHASE LAG (degree)

0 4 8 12


16 20 24 28

& U (cm/s)

PHASE LAG (degree)
-180 -120 -60 0 60

120 180






(9- measured

0.60 0.70 0.80 030 1.00 1.1,0

P Pig

Fig. 4-4. Comparison of the Model Prediction and Measurement for Run 1-2. (a) Velocity; (b) Pressure.

PHASE LAG (degree) PHASE LAG (degree)

(a) Run S-2 (b)
SUL gW __R_ __ _
3S : ^--


t '- "Ph a
L~J 25
4, : Phase
/ dt
- 20 ^U "+
,I/ -Phase -;'

I" i t d ,-4

d+I 12
(0- measured u) ds (e-mesured

0 5 to is 20 2S 30 3S 0.50 0.60 0.2o 0.80 0.0 1.00, 1.10

UELOCITY U & U (cm/9) 0/Plga
Fig. 45. Comparison of the Model Prediction and Measurement for Run 5-2. (a) Velocity; (b) Pressure.

gradient at the water-mud interface is larger than what meets the eye

in Figs. 4-4a and 4-5a. Notice that the boundary layer thickness for

the mud layers, Vve/2o 10 cm, was greater than the thickness of each

layer and resulted in the build-up of a horizontal velocity overshoot

above the water-mud interface. The measured horizontal velocities (see

Chapter 7 for details) in the water layer are also displayed. Measured

values in these two examples show a reasonable agreement with model-

simulated values.

Although the model gave information on the motion in the mud

layers, experimentally measuring the motion was difficult. Here the

horizontal mud displacements were measured by visual observation from

the side wall for experimental runs 5 and 6. The amplitudes of move-

ment of sediment particles near the side wall were converted to the

corresponding velocity amplitudes (by using the formula u = ao) and

compared with model results. The observed mud velocities, however,

were reduced, due to the side boundary effects; therefore, they were

smaller than predicted (by a factor of about 1:3.7, see Appendix B).

In justifying the measurements, a relationship for the lateral distri-

bution of the horizontal velocity, u(y), where y is the lateral

coordinate, must be worked out. Efforts to interpret the results by

considering the side wall influence are given in Appendix B. The

resulting value of the ratio of u in the center of the flume and u near

the side wall ranged from about 4 to 7, which is not far from the

average value of 3.7 obtained by a direct comparison of prediction and


4.5.2 Pressure

The pressure is always in phase both in the water and mud layers

since it is the driving force. This phase relationship was also veri-

fied from experiments, see Section 7.2.2 in Chapter 7. The amplitude

decreased in the water column, dropped at each interface between the

layers, and increased within the mud layers. Figures 4-4b and 4-5b

show two examples. Note that the measured pressure amplitudes (see

Chapter 6 for details) were normalized by p1ga. At each interface, the

normal stresses were matched from the two adjoining layers. Therefore,

the pressure, -P on 2p Dw/3z dropped between layers because of a

step-increase in p and almost the same 3w/az. Within a mud layer, the

shear stress term becomes important because of the high ve value. The

behavior in a vertical pressure profile can be depicted easily by

substituting Eq. 4-6 and 4-7 into the vertical momentum equation (Eq.

4-3b). The vertical pressure gradient can be obtained as shown in

Eq. 4-30. For low viscosity, e.g., water with v = 1 x 10-6 m2/s, the

first term of the right hand side dominates, and therefore yields a

positive pressure gradient. For large ve the second term is also

important and may result in a negative pressure gradient as shown in

Figs. 4-4b and 4-5b.

T p [(ja-k v ) w w" ] (4-30)

As already noted, the model uses layers of mud which have constant

properties to simulate continuously varying properties. Therefore, the

true pressure profile may be estimated by a smooth curve, e.g., the

broken line in Fig. 4-5b.


Run 5 shows a reasonable agreement; however, Run 3 indicated that

the measured pressure was larger (about 30%) than predicted-. For other

data, see Appendix B. Generally, the measured dynamic pressure was

higher than that predicted by the model for kaolinite beds. The same

was true for the Cedar Key mud bed, but only during the initial period

of the wave loading. One possible reason is that the pressure trans-

ducer measured the total dynamic pressure, which also included the

effect of time-varying mud structure and generation of pore water pres-

sure. The model, however, does not account for these effects. After

the completion of these processes, the measured dynamic pressure

principally reflected the wave-induced pressure.

4.5.3 Shear Stress

The shear stress amplitude profile in this system showed a linear

increase with depth characteristic of the mud bed over a rigid bottom,

and negligible shear stresses in the water column except near the

water-mud interface. See Fig. 4-6 for an example. The small total mud

thickness, less than half of the wave length, is the likely reason for

a linearly increasing shear stress profile. For an infinite mud thick-

ness, the shear stress would eventually reach a maximum and then de-

crease (Mallard and Dalrymple, 1977). Shear stresses in the water

layer are negligibly small because of the small velocity gradient

except near the water-mud interface. Notice the nearly 90 degree shift

of shear stress in the mud layers. This implies that the maximum shear

stress occurred between the wave crest and trough. It may further

indicate that the shear stress in the mud was induced mainly by the


PHASE LAG (degree)
4-8 10 -60 0 60 .120 180 240
Run 1-1


LU 2s

" 20
T .b=025N/r'

U- 4

Cr d3


0 20 40 so 80 t0 120 140


Fig. 4-6. An Example of the Model-Predicted Shear Stress
Amplitude Profile.


longitudinal pressure variation at the mud surface rather than the

water-mud interracial shear stress. The maximum longitudinal pressure

gradient occurs near the mean water elevation (90 degrees after wave

crest), while the interfacial shear stress, being proportional to the

velocity, is minimum. The figure shows oscillations of the shear

stress just above mud surface. This would not be expected in reality

and may be due to the calculation of phase lag by using small shear

stress resulting in a numerical oscillation. Finally, it should be

noted that the large shear stresses in the mud layer are due to the

large viscosity, ve, and finite shear strain.

The model-predicted bed shear stress amplitude at x = 0, Tbo, is

b e (w" + k2 w) (4-31)
bo el k 1 1

The values of Tbo for the six runs are summarized in Table 4-2. Bed

shear stress amplitudes based on the rigid mud bed assumption and Eq.

2-7 are also listed in the parentheses for comparison. It is observed

that the model predicted generally larger Tbo (up to 30%) than those

from a rigid bed. This is because of the opposite (100 120 degree

phase lag) motion between the mud and the water. However, for a very

soft mud, e.g., Run 5, the model predicted a smaller value. It should

also be noted that the longitudinal shear stress profile also exhibits

an exponential decay law for Rw < 104. The bed shear stress amplitude

at any location can therefore be evaluated as Tb(x) = Tbo exp(-kimx).

The spatially-averaged bed shear stress amplitude is

b 1 L Tbo
<'b> T W bX dx [1 -exp(-kimL)] (4 -32)
b L b im Li


where L = 8 m is the length of mud bed section and kim is the measured

wave decay coefficient (see Chapter 6 for details). The corresponding

Rw are also included in Table 4-2.

4.5.4 Water Wave Decay

A comparison of predicted and measured wave decay coefficient,

kim, is an easy way to check the validity of the model. A viscous

damping mechanism implies that the greater the mud velocity, the higher

the rate of energy dissipation. Therefore, a highly consolidated

(rigid) bed can not be expected to have a high energy dissipation.

Both the predicted and measured kim show this behavior. See Fig. 4-7

and also Table 6-5, Chapter 6. Except for Run 6-I, the predicted and

measured wave decay coefficient show satisfactory agreement.

Table 4-2. Interfacial Shear Stresses and

Mud Wave Amplitudes

Run T bo Rw bol

N/m2 N/m2 cm

1-I 0.247 (0.187) 0.170 3600 0.24
2 0.386 (0.304) 0.242 4040 0.23
2-1 0.166 (0.163) 0.146 2360 0.05
2 0.291 (0.290) 0.245 3070 0.09
3-I 0.210 (0.178) 0.124 2830 0.23
2 0.293 (0.258) 0.158 2720 0.27
4-i 0.196 (0.192) 0.178 2510 0.06
2 0.297 (0.294) 0.251 2660 0.05
5-I 0.272 (0.274) 0.191 5130 0.14
2 0.421 (0.425) 0.314 6220 0.15
6-1 0.226 (0.200) 0.149 3450 0.07
2 0.352 (0.306) 0.248 3590 0.13
3 0.431 (0.376) 0.323 3810 0.21


Digits Indicate Run Numbers

6 1
6 5 .

4 2




aO 0.2 03

Fig. 417. Comparison of the Measured and Predicted Wave
Decay Coefficient.



E 0.2




t ....