3rd DRAFT

EFDC Technical Memorandum Theoretical and Computational Aspects of Sediment and Contaminant Transport in the EFDC Model
Prepared for: US Environmental Protection Agency, Office of Science and Technology 401 M Street SW Washington, DC 20460

Prepared by: Tetra Tech, Inc. 10306 Eaton Place Suite 340 Fairfax, Virginia 22030

May 2002

Table of Contents
1. 2. 3. 4. 5. 6. 7. 8. 9. Introduction Summary of Hydrodynamic and Generic Transport Formulations Solution of the Sediment Transport Equation Hydrodynamic and Sediment Boundary Layers Sediment Bed Mass Conservation and Geomechanics Noncohesive Sediment Settling, Deposition and Resuspension Cohesive Sediment Settling, Deposition and Resuspension Sorptive Contaminant Transport References 3 4 9 11 14 26 34 47 53

1.

Introduction

This report summarizes theoretical and computational aspects of the sediment and sorptive contaminant transport formulations used in the EFDC model. Theoretical and computational aspects for the basic EFDC hydrodynamic and generic transport model components are presented in Hamrick (1992). Theoretical and computational aspects of the EFDC water quality-eutrophication model component are presented in Park et al. (1995). The paper by Hamrick and Wu (1997) also summarized computational aspects of the hydrodynamic, generic transport and water quality-eutrophication components of the EFDC model. The EFDC model has been extensively applied to estuaries (Fredricks and Hamrick, 1996; Shen and Kuo, 1999; Shen et al., 1999; Ji et al., 2001), lakes (Jin et al., 2000; 2002), reservoirs (Hamrick and Mills, 2000), rivers (Ji, et al., 2002), and wetlands (Moustafa and Hamrick, 2000). The model has also been used for a number of fundamental process studies (Hamrick, 1994; Kuo, et al., 1996; Yang, et al., 2000). This report is organized as follows. Chapter 2 summarizes the hydrodynamic and generic transport formulations used in EFDC. Chapter 3 summarizes the solution of the transport equation for suspended cohesive and noncohesive sediment. A discussion of near bed boundary layer processes relevant to sediment transport is presented in Chapter 4. Sediment bed mass conservation and methods for representation of the beds geomechanical properties are discussed in Chapter 5. Chapters 6 and 7 summarize noncohesive and cohesive sediment settling, deposition and resuspension process representations. The final chapter, Chapter 9, documents the EFDC model's sorptive contaminant transport and fate formulations.

2.

Summary of Hydrodynamic and Generic Transport Formulations

This section summarizes the hydrodynamic and transport equations used by the EFDC model. Reference is made to Hamrick (1992) and Hamrick and Wu (1997) for details of the computational procedure. This section does however describe modifications to the solution procedure when the model operates in a geomorphologic mode. The EFDC model's hydrodynamic component is based on the three-dimensional hydrostatic equations formulated in curvilinear-orthogonal horizontal coordinates and a sigma or stretched vertical coordinate. The momentum equations are:
t

mx my Hu my H
x

my Huu my

mx Hvu
* x b

mx my wu
z

f e mx my Hv Av H
1/ 2 z

patm

z xH u

mx my

(2.1)

my
x

mx

HAH xu

mx HAH my
y

mx my c p Dp u 2 v 2

mx my Hv mx H
y

my Huv mx

mx Hvv
* y b

mx my wv
z

f e mx my Hu mx m y Av H
1/ 2 z

p HAH

patm v

z v

(2.2)

my
x

mx

mx HAH my

mx my c p D p u 2 v 2

mx my fe
xz

mx my f
,
yz

u y mx v x my
1 z

(2.3) (2.4)

Av H

u ,v

where u and v are the horizontal velocity components in the dimensionless curvilinearorthogonal horizontal coordinates x and y, respectively. The scale factors of the horizontal coordinates are mx and my. The vertical velocity in the stretched vertical coordinate z is w. The physical vertical coordinates of the free surface and bottom bed are zs* and zb* respectively. The total water column depth is H, and is the free surface potential which is equal to gzs*. The effective Coriolis acceleration fe incorporates the curvature acceleration terms, with the Coriolis parameter, f, according to (2.3). The Q terms in (2.1) and (2.2) represents optional horizontal momentum diffusion terms. The vertical turbulent viscosity Av relates the shear stresses to the vertical shear of the horizontal velocity components by (4.4). The kinematic atmospheric pressure, referenced to water density, is patm, while the excess hydrostatic pressure in the water column is given by:

gHb

gH

1 o o

(2.5)

where and o are the actual and reference water densities and b is the buoyancy. The horizontal turbulent stress on the last lines of (2.1) and (2.2), with AH being the horizontal turbulent viscosity, are typically retained when the advective acceleration are represented by central differences. The last terms in (2.1) and (2.2) represent vegetation resistance where cp is a resistance coefficient and Dp is the dimensionless projected vegetation area normal to the flow per unit horizontal area. The three-dimensional continuity equation in the stretched vertical and curvilinearorthogonal horizontal coordinate system is:
t

mx my H

my Hu

mx Hv

mx my w

QH

0 QSS

QSW

(2.6)

with QH representing volume sources and sinks including rainfall, evaporation, and lateral inflows and outflows having negligible momentum fluxes. The terms QSS and QSW are the net volumetric fluxes of sediment and water between the bed and water column, defined as positive from the bed to the water column, when the model operates in a geomorphologic mode. The delta function, (0) indicates these fluxes enter the bottom layer of the water column. Integration of (2.6) over the depth gives
t

mx my H

my Hu

mx Hv

QH

QSS

QSW

(2.7)

In the geomorphologic mode, the water column continuity equation is coupled to a bulk volume conservation equation for the sediment bed.
t

mx my B

QGW

QSS

QSW

(2.8)

where B is the total thickness of the resolved sediment bed and QGW is the volumetric ambient groundwater inflow at the bottom of the sediment bed. The bed surface elevation is defined by
* B zbb

(2.9)

Where zbb* is the time invariant elevation at the bottom of the sediment bed. Using (2.9), equation (2.8) can be written as
t

mx my

QGW

QSS

QSW

(2.10)

Adding (2.7) and (2.10) gives

t

mx my

my Hu

mx Hv

QH

QGW

(2.11)

where the water surface elevation, , is defined by

* zs

(2.12)

The EFDC model solves the external mode continuity equation (2.11) using a two-step procedure. The first step corresponding to the standard implicit external mode hydrodynamic solution is
mx m y
*

mx m y
y

2
n 1

my Hu
y

n 1

2
n

m y Hu

(2.13)

mx Hv

mx Hv

n QH 1/ 2

where is the time step between n and n+1. The intermediate time level notation, n+1/2, denotes an average between the two time levels. The second step is taken after the bed volumetric continuity equation is updated to time level n+1 and is

mx my

n 1

mx my

n QG

1/ 2

(2.14)

Combining (2.13) and (2.14) gives the equivalent full step.

mx m y 2
y n 1

mx m y
n 1

2 2
y

my Hu
n

n 1

2
n QH 1/ 2

m y Hu
n QG 1/ 2

(2.15)

mx Hv

mx Hv

The water column depth is then updated by

Hn
1 n 1 n 1

(2.16)

prior to the next hydrodynamic time step. The EFDC model includes the ability to simulate drying and wetting of shallow areas. Drying and wetting is iteratively determined during the implicit solution of equation (2.13) after the time discrete depth average horizontal momentum equations have been inserted to form an elliptic equation for the water surface elevation. The solution procedure is as follows. A preliminary solution for the water surface elevation is determined by solving (2.13) with all horizontal grid interior horizontal cell faces open. The resulting cell center water depth in each cell is then compared to a small dry depth Hdry. If the depth is greater than the dry depth, the cell is defined as wet. If the depth is less than the dry depth and less than the depth at the previous time step, the cell is defined as wet and its four flow faces are blocked. If the depth is less than the dry depth, but greater than the depth at the previous time step, the direction of flow on each cell face is checked and faces having outflow are block. Following this checking and blocking,

(2.13) is solved again, followed by the same checking procedure. This iteration is repeated until wet or dry status of each cell does not change from that of the subsequent iteration. Typically two or three iterations are required. This implementation of drying and wetting is fully mass conservative and does not produce negative water column depths. The generic transport equation for a dissolved or suspended material having a mass per unit volume concentration C, is
t

mx my HC my
x

my HuC
x

mx HvC
y

mx m y wC
z

mx m y wscC Qc

(2.17)

mx

HK H

mx HK H my

K mx m y v H

where KV and KH are the vertical and horizontal turbulent diffusion coefficients, respectively, wsc is a positive settling velocity went C represents a suspended material, and Qc represents external sources and sinks and reactive internal sources and sinks. The solution of the momentum equations, (2.1) and (2.2) and the transport equation (2.17), requires the specification of the vertical turbulent viscosity, AV, and diffusivity, Kv. To provide the vertical turbulent viscosity and diffusivity, the second moment turbulence closure model developed by Mellor and Yamada (1982) and modified by Galperin et al., (1988) and Blumberg et al., (1988) is used. The MY model relates the vertical turbulent viscosity and diffusivity to the turbulent intensity, q, a turbulent length scale, l, and a turbulent intensity and length scaled based Richardson number, Rq, by:

Av
A

Ao ql

1 R1 1 Rq 1 R2 1 Rq 1 R3 1 Rq A1 1 3C1 B2 3 A2 R1
1

Ao

6 A1 B1

1 1/ B1 3
(2.18)

3 A2

6 A1 1 B1 1 3C1

3C1 B2 6 A1 6 A1 B1

R2 1 R3 1

9 A1 A2

3 A2 6 A1 B2

Kv
K

K o ql

1 1 R3 1 Rq A2 1 6 A1 B1

(2.19)

Ko

Rq

gH z b l 2 q2 H 2

(2.20)

where the so-called stability functions, A and K, account for reduced and enhanced vertical mixing or transport in stable and unstable vertically density stratified environments, respectively. Mellor and Yamada (1982) specify the constants A1, B1, C1, A2, and B2 as 0.92, 16.6, 0.08, 0.74, and 10.1, respectively. The turbulent intensity and the turbulent length scale are determined by the transport equations:
t

mx m y Hq 2
z

m y Huq 2 mx m y v Aq H
2 2 zq

mx Hvq 2 2 mx m y Hq 3 B1l

mx m y wq 2

(2.21)

2 mx m y

Av H

2 z

p p

c Dp u 2 v 2

3/ 2

gK v

Qq

mx my Hq 2l Aq H Av H
z

m y Huq 2l

mx Hvq 2l l Hz
p p 2

mx m y wq 2l
2

mx my

ql
2

Hq 3 mx my 1 E2 B1
z

E3

l H 1 z
3/ 2

(2.22)

mx my E1l

gK v

c Dp u 2 v 2

Ql

where (E1, E2, E3) = (1.8, 1.33, 0.25). The second term on the last line of each equation represents net turbulent energy production by vegetation drag where p is a production efficiency factor having a value less than one. The terms Qq and Ql may represent additional source-sink terms such as subgrid scale horizontal turbulent diffusion. The vertical diffusivity, Aq, is set to 0.2ql following Mellor and Yamada (1982). For stable stratification, Galperin et al., (1988) suggest limiting the length scale such that the square root of Rq is less than 0.52. When horizontal turbulent viscosity and diffusivity are included in the momentum and transport equations, they are determined independently using Smagorinsky's (1963) subgrid scale closure formulation.

Vertical boundary conditions for the solution of the momentum equations are based on the specification of the kinematic shear stresses, equation (2.4), at the bed and the free surface. At the free surface, the x and y components of the stress are specified by the water surface wind stress
xz

yz

sx

sy

2 cs U w Vw2 U w ,Vw

(2.23)

where Uw and Vw are the x and y components of the wind velocity at 10 meters above the water surface. The wind stress coefficient is given by:
cs 0.001
a w 2 0.8 0.065 U w Vw2

(2.24)

for the wind velocity components in meters per second, with a and w denoting air and water densities respectively. At the bed, the stress components are related to the near bed or bottom layer velocity components by the quadratic resistance formulation
xz

yz

bx

by

cb u12 v12 u1 , v1

(2.25)

where the 1 subscript denotes bottom layer values. Under the assumption that the near bottom velocity profile is logarithmic at any instant of time, the bottom stress coefficient is given by
2

cb

ln(

/ 2 zo )

(2.26)

where , is the von Karman constant, 1 is the dimensionless thickness of the bottom layer, and zo=zo*/H is the dimensionless roughness height. Vertical boundary conditions for the turbulent kinetic energy and length scale equations are:
q2 q2 B12 / 3 s B12 / 3 b : z 1 : z 0

(2.27) (2.28) (2.29)

: z 0,1

where the absolute values indicate the magnitude of the enclosed vector quantity. Equation (2.28) can become inappropriate under a number of conditions associated with either or both high near bottom sediment concentrations and high frequency surface wave activity. The quantification of sediment and wave effects on the bottom stress is discussed in Chapter 4.

3.

Solution of the Sediment Transport Equation

This section describes the solution of the transport equations for suspended sediment. The general procedure follows that for the salinity transport equation, which uses a high order upwind difference solution scheme for the advective terms, described in Hamrick (1992). Although the advection scheme is designed to minimize numerical diffusion, a small amount of horizontal diffusion remains inherent in the scheme. Due the small inherent numerical diffusion, the physical horizontal diffusion terms in (2.17) are omitted as to give:
t

mx my HS j
z

my HuS j
z

mx HvS j
z

mx my wS j Q
I sj

(3.1)

mx my wsj S j

K mx my V H

Sj

E sj

where Sj represents the concentration of the jth sediment class and the source-sink term has been split into an external part, which would include point and nonpoint source loads, and internal part which could include reactive decay of organic sediments or the exchange of mass between sediment classes if floc formation and destruction were simulated. Vertical boundary conditions for (3.1) are:
KV wsj S j J jo : z 0 zS j H KV wsj S j 0 : z 1 zS j H

(3.2)

where Jjo is the net water column-bed exchange flux defined as positive into the water column. The numerical solution of (3.1) utilizes a fractional step procedure. The first step advances the concentration due to advection and external sources and sinks having corresponding volume fluxes by

H n 1S * my Hu
n 1/ 2

H nS n

mx my mx Hv

E Qsj

n 1/ 2

(3.3)

mx my

Sn

n 1/ 2

Sn

mx my wn

1/ 2

Sn

where n and n+1 denote the old and new time levels and * denotes the intermediate fractional step results. The portion of the source and sink term, associated with volumetric sources and sinks is included in the advective step for consistency with the continuity constraint. This source-sink term, as well as the advective field (u,v,w), is defined as intermediate in time between the old and new time levels consistent with the temporal discretization of the continuity equation. Note that the sediment class subscripts

have been dropped for clarity. The advection step uses the anti-diffusive MPDATA scheme (Smolarkiewicz and Clark, 1986) with optional flux corrected transport (Smolarkiewicz and Grabowski, 1990). The second fractional step or settling step is given by

S **

S*

n 1

ws S **

(3.4)

which is solved by a fully implicit upwind difference scheme

** Skc * Skc z ** Sk * Sk 1

n 1

ws S ** ws S ** ws S **

(3.5)
kc

n 1

ws S ** S1**

k 1 z

n 1

: 2 k

kc 1

S1*
z

n 1

marching downward from the top layer. The implicit solution includes an optional antidiffusion correction across internal water column layer interfaces. The third fractional step accounts for water column-bed exchange by resuspension and deposition
S1*** S1**
z *** Lo J o

(3.6)

n 1

Where Lo is a flux limiter such that only the current top layer of the bed can be completely resuspended in single time step. The representation of the water column bed exchange by a distinct fractional step is equivalent to a splitting of the bottom boundary condition (3.2) such that the bed flux is imposed intermediate between settling and vertical diffusion. For resuspension and deposition of suspended noncohesive sediment, the bed flux is given by
*** Jo

ws

Seq

S1***

(3.7)

which will be further discussed in Chapters 4 and 6. Inserting (3.7) into (3.6) gives

Lo ws n 1 zH

S1***

S1**

Lo ws n 1 zH

(3.8)

Seq

For cohesive sediment resuspension, the bed flux is specified as a function of the bed stress and bed geomechanical properties. For cohesive sediment deposition, the bed flux is typically given by
*** Jo

Pd ws S1***

(3.9)

where Pd is a probability of deposition which will be further discussed in Chapter 7. Inserting (3.9) into (3.6) gives

Pd ws S1*** n 1 zH

S1**

(3.10)

The remaining step is an implicit vertical turbulent diffusion step corresponding

S
n 1

*** z

KV H2

n 1 z

(3.11)
S
n 1

with zero diffusive fluxes at the bed and water surface.

4.

Hydrodynamic and Sediment Boundary Layers

Both two-dimensional and three-dimensional applications of the EFDC model require parameterization of near bed boundary layer processes. In the absences of high frequency surface gravity waves and when sediment transport is not being simulated, this parameterization is made through the bottom friction coefficient, (2.26) and the bottom turbulence intensity boundary conditions (2.28). The presence of high frequency surface gravity waves and near bed gradients of suspended sediment requires additional parameterization since the sediment and wave boundary layers cannot be directly resolved by typical vertical grid resolution. Approximate parameterizations of hydrodynamic and sediment boundary layer appropriate for representing the bottom stress and the water column-bed exchange of sediment under conditions including ambient flow, high frequency surface waves and high near bed suspended sediment gradients can be derived form simplified forms of the momentum and sediment transport equations and the turbulent kinetic energy equation.

4.1

Boundary Layer Equations

First consider the horizontally homogeneous momentum equation written in vector form
t

p g

1 z

AV z u

(4.1)

The horizontal velocity, pressure and water surface elevation can be decomposed into components associated with the current or mean flow and the high frequency surface gravity wave motion
u uc u w p
c

pw
w

(4.2)

where the current pressure in excess of hydrostatic pressure has been set to zero. Assuming the current is steady with respect to the time scale of the wave motion and inserting (4.2) into (4.1) gives
t

uw

pw

1 z

AV

u w uc

(4.3)

On non-geophysical scales where the bottom current boundary layers does not exhibit Ekman effects, equation (4.3) can be vectorially split into components aligned with the wave and current directions

t w

pw g

Av H2 Av H2

z w

(4.4)

cos

z c

cos

t w

pw Av H2

Av H2 0

z w

(4.5)

z c

where c and w are the directions of the current and wave propagation, respectively, and for simplicity in notation uw and uc are the wave and current velocities in these two directions. Subtracting the wave period average of (4.4) from (4.4) gives an equation for the wave motion

t w

pw

Av H2 Av H
2

z w

Av H2
z c

z w

(4.6)

cos

Av

Averaging (4.5) over the wave period gives an equation for the mean current

Av
c c z

z c

cos

Av H2

z w

(4.7)

Wave-current boundary layer models formulated for use with numerical circulation models typically neglect variations in the vertical turbulent viscosity at the wave time scale (Styles and Glenn, 2000) allowing (4.6) and (4.7) to be reduced to
u pw g Av H2 u 0

(4.8)

t w

z w

Av H2

z c

(4.9)

Above the wave boundary layer, the wave velocity field is inviscid and (4.8) reduces to
t w

pw

(4.10)

which is subtracted from (4.8) to give the wave boundary layer equation
u Av H2 u u

(4.11)

t w

z w

t w

The boundary conditions for (4.11) are

Av H

z w

wb

(4.12)

or uw 0

As z goes to the roughness height zo, and

uw uw

(4.13)

as z becomes large. Integrating of (4.9) over the bottom hydrodynamic model layer and subtracting the results from (4.9) gives the current boundary layer equation

Av H

z c 1

c1

cb

(4.14)

where the c1 and cb subscripts denote the shear stresses at the top and bottom of the bottom grid layer. Integration of (4.14) gives
Av H u z
cb c1 cb 1

(4.15)

z c

Where 1 is the dimensionless thickness of the bottom grid layer. For small z near the bed, (4.15) is approximated as a constant stress layer

Av H
The boundary condition for (4.16) is

z c

(4.16)
cb

uc

(4.17)

as z goes to the roughness height zo. In the bottom hydrodynamic layer the integral condition

(4.17)

uc dz
0

1 1

is imposed where u1 is the current velocity in the bottom grid layer. The sediment boundary layer equation can be derived form the horizontally homogeneous approximation to the sediment transport equation (3.1).
HS ws S KV H S 0

(4.18)

Integrating (4.18) over the bottom hydrodynamic layer gives

HS1 J sb
1

J s1

(4.19)

Where S1 is the bottom layer sediment concentration and Jsb and Js1 are the sediment fluxes at the bed and the top of the bottom grid layer. Subtracting (4.19) from (4.18) gives
HS HS1 ws S KV H S J s1 J sb
1

(4.20)

Assuming that the temporal derivative in (4.20) is small and can be neglected, (4.20) is integrated to give
ws S KV H S J sb J s1 J sb z
1

(4.21)

For small z near the bed, (4.21) is approximated as a constant flux layer

ws S

KV H

J sb 1

z
1

(4.21x)

ws S

KV H

J sb

(4.22)

The bottom boundary condition for (4.22) is

S Sr

(4.23)

as z goes to the dimensionless sediment reference height zr, which can be the roughness height. In the bottom hydrodynamic layer the integral condition
1

(4.24)

Sdz
0

1S1

is imposed. The near bed wave, current and sediment boundary layer equations (4.11, 4.16, and 4.22) require specification of the near bed forms of the vertical turbulent viscosity and diffusion coefficients. Near the bed, the turbulent kinetic energy equation (2.21) can be approximated by its equilibrium form

q3 B1l

Av H2

2 z

Kv H

(4.25)

where the vegetation term has been dropped since the horizontal velocity components approach zero. Introducing the definitions of Av and Kv given by (2.18) and (2.19) and solving for the turbulent intensity gives
q2 B1 Ao 1 B1 K o
A K

Rq

l2 H2

zu

zv

(4.26)

Equation (4.25) can be also be written in terms of the shear stresses after multiplying by Av, inserting the definitions of Av and Kv given by (2.18) and (2.19), and using (2.20), to give
q
2

B1 Ao A

1/ 2

1 B1 K o

K Rq

1/ 2

2 xz

2 1/ 2 yz

(4.27)

When (4.27) is evaluated at the bed, the results

q
2 b

B1 Ao A

1/ 2

1 B1 K o

K Rq

1/ 2

2 bx

2 1/ 2 by

(4.28)

is equivalent to (2.28) under neutral conditions where Rq is equal to zero. High near bed sediment concentrations and associated vertical gradients can result in nonzero values of Rq immediately above the bed. The buoyancy gradient near the bed is primarily due to gradients in suspended sediment concentration with the effect of sediment on density given by

1
j

Sj
w sj

Sj
sj sj

(4.29)

where Sj is the mass concentration of sediment class j per unit volume of the watersediment mixture. The buoyancy is expressed in terms of the sediment concentration using
b
w w j sj w w

Sj
j sj j

Sj

(4.30)

which can be used to evaluate the buoyancy gradients. When high frequency surface waves are present, the velocity components in (4.25) and (4.26) and the shear stress components in (4.26) and (4.27) can be decomposed into

u uc cos v uc sin
xz yz cz cz

c c

uw cos uw sin
c c wz wz

w w

(4.31)

cos sin

cos sin

w w

(4.32)

where uc and uw are the current and wave velocities and c and w are the current and wave shear stress magnitudes, each aligned with the current and wave directions denoted by c and w. Using (4.32) and (4.32), the shear and bed stress terms can be written as
z

2 z

z c

z w

2cos

z c

z w

(4.33) (4.34)

2 xz

2 yz

2 cz

2 wz

2cos

cz wz

Assuming the wave velocity and shear stress to be periodic

uw uwm sin
wz w wm wzm

t t t

sin

(4.35)

sgn sin

the wave period averages of (4.31) and (4.32) are

z c

z w

z c

1 2

z wm

cos

wm

z c

z wm

(4.36)

2 cz

2 wz

2 cz

1 2

2 wzm

cos

(4.37)
c wm cb wzm

Analytical solutions of the wave, current and sediment boundary layer equations (4.11, 4.16, and 4.22), as exemplified most recently by Styles and Glenn (2000), typically assume tractable forms of the vertical turbulent viscosity and diffusivity inside the wavecurrent and the current boundary layers. The following sections discuss boundary layer parameterization for neutral and stratified boundary layers in absences and presences of waves.

4.2

Neutral Current and Sediment Boundary Layers

For neutral conditions, the turbulent intensity (4.27) and the vertical turbulent exchange coefficients (2.18) and (2.19) can be written as

q2 Avn
K vn

B12/ 3 Ao ql

2 xz

2 1/ 2 yz

(4.38) (4.39) (4.40)

2 xz

2 1/ 4 yz

l
l

K o ql

Ko Ao

2 xz

2 1/ 4 yz

For three-dimensional, multiple vertical layer applications equation (4.16) becomes

l H

z c

(4.41)
cb

Letting l/H = z, and using (4.17) gives the logarithmic profile

uc

cb

ln

z zo

(4.42)

Applying the integral condition (4.17) over the bottom layer gives
cb

cb u1 u1
2

(4.43)

cb

ln(

/ 2 zo )

For two-dimensional depth average applications (4.15) becomes

l H

z c

cb

1 z

(4.44)

For consistency with the subsequent solution of the sediment boundary layer equation, the length scale is chosen as

l H
With the solution of (4.44) becoming
uc
cb

z 1 z

(4.45)

ln

z zo

(4.46)
z zo

Applying the integral constraint (4.14) to (4.46) gives

cb

cb u1 u1
2

cb

(4.47)

ln(1/ 2 zo )

For three-dimensional multiple layer, applications, the sediment boundary layer equation (4.22) can be written as
z R S S J sb ws

(4.48)

where

Ao Ko

ws
cb

(4.49)

is the Rouse parameter. The solution of (4.48) is

S J sb ws C zR

(4.50)

For noncohesive sediment, the constant of integration is evaluated using

Seq : z

zeq

and

J sb

(4.51)

which sets the near bed sediment concentration to an equilibrium value, Seq, defined just above the bed under no net flux condition. Using (4.51), equation (4.50) becomes

zeq z

Seq

J sb ws

(4.52)

For nonequilibrium conditions, the net flux is given by evaluating (4.52) at the equilibrium level
J sb ws Seq Sne

(4.53)

where Sne is the actual concentration at the reference equilibrium level. Equation (4.53) indicates that when the near bed sediment concentration is less than the equilibrium value a net flux from the bed into the water column occurs. Likewise when the concentration exceeds equilibrium, a net flux to the bed occurs. For the relationship (4.53) to be useful in a numerical model, the bed flux must be expressed in terms of the model layer mean concentration. For a three-dimensional application, (4.53) and the constraint (4.24) give
J sb ws Seqe S1

(4.54)

where

Seqe

ln

zeq1 zeq1 1
1 R

Seq : R 1
(4.55)

zeq1 Seqe 1 R

1 Seq : R 1

zeq1 1

defines an effective bottom layer mean equilibrium concentration in terms of the near bed equilibrium concentration. The corresponding quantities in the numerical solution bottom boundary condition (3.7) are

Wr Sr Wd

ws Seqe ws

(4.56)

If the dimensionless equilibrium elevation, zeq exceeds the dimensionless layer thickness, (4.54) and (4.55) can be modified to
J sb ws Seqe S

(4.57)

Seqe

ln M zeq1 M zeq1 1 M zeq1

1 R

Seq : R 1
(4.58)

1 Seq : R 1

Seqe

1 R M zeq1 1

where the over bars in (4.57) and (4.58) implying an average of the first M grid layers above the bed. When multiple sediment size classes are simulated, the equilibrium concentration, Seq, in (4.55) and (4.58) are reduced from their uniform values by multiplying by the sediment class volume fractions at the bed surface. For cohesive sediment resuspension, the flux is presumed known, and the constant of integration in (4.48) is determined by the integral constraint with the resulting sediment concentration distribution being

J sb ws S J sb ws

1 R z
R 1 R 1 1

zr z
1 R r

S1 S1

J sb ws J sb ws

: R 1
(4.59)

zr
1 r

z ln

: R 1

For cohesive sediment deposition, the bed flux is given by

J sb Pd ws S r

(4.60)

where Pd is the probability of deposition. Evaluating (4.59) at the reference level, inserting into (4.60) and solving, gives the deposition flux in terms of the bottom layer concentration
1

J sb

1 Pd

Pd 1 R zrR
1 R 1

zr
R 1

1 R zrR
1 1 R 1

zr z1 r
R

z1 r
1

Pd ws S1 : R 1

(4.61)

J sb

1 Pd

Pd zr ln

zr
1 1 zr

zr
1 1 zr

zr ln

Pd ws S1 : R 1

The sediment concentration profile under depositional conditions is also give by (4.59) using the flux from (4.61). For depth average applications, the sediment boundary layer equation (4.21) can be written as

1 z l R H

S zS

J sb 1 z ws

(4.59)

A closed form solution is possible by choosing

l H
with (4.59) becoming
S R S z

z 1 z

(4.60)

J sb R 1 z ws z

(4.61)

The solution of (4.61) is

S 1 Rz 1 R J sb ws C zR

(4.62)

Evaluating the constant of integration using (4.51) gives

zeq z

Seq

Rz 1 R

J sb ws

(4.63)

For nonequilibrium conditions, the net flux is given by evaluating (4.63) at the equilibrium level

J sb

ws

1 R 1 R 1 zeq

(4.64)

Seq

S ne

where Sne is the actual concentration at the reference equilibrium level. Since zeq is on the order of the sediment grain diameter divided by the depth of the water column, (4.64) is essentially equivalent to (4.54). To obtain an expression for the bed flux in terms of the depth average sediment concentration, equation (4.63) is integrated over the depth to give

J sb

ws

21 R 2 R 1 zeq

(4.65)

Seqe

where

Seqe

ln zeq1 zeq1 1
R zeq 1 1

Seq : R 1

(4.66)
Seq : R 1

Seqe

1 R zeq1 1

When multiple sediment size classes are simulated, the equilibrium concentration, Seq, in (4.66) is reduced from its uniform value by multiplying by the sediment class volume fractions at the bed surface. The corresponding quantities in the numerical solution bottom boundary condition (3.7) are

Wr Sr

ws

21 R 2 R 1 zeq 21 R 2 R 1 zeq

Seqe
(4.67)

Wd

ws

For cohesive sediment resuspension, the flux is presumed known, and the constant of integration in (4.62) is determined by the integral constraint with the resulting sediment concentration distribution being
R 1 zr 1 21 R
1 1 1 1

zr 1 R
1 R

zr
1 R

1 R

1 zR

Rz 1 R

J sb ws

(4.68)

zr 1 R
1 R

zr
1 1

S1 zR
1

: R 1 zr 4 z 2 J sb ws

zr z
1

ln zr z

zr ln
1

1
1

zr

S1

: R 1

For cohesive sediment deposition, the bed flux is given by

J sb Pd ws S r

(4.69)

where Pd is the probability of deposition. The depositional flux can be determined by evaluating (4.68) at the reference level, inserting the results into (4.69), and solving for the flux. The sediment concentration profile under depositional conditions is also give by (4.68) using the depositional flux.

4.3

Stratified Current and Sediment Boundary Layers

Analytical solutions for stratified current and sediment boundary layers are difficult to obtain unless tractable expressions are assumed for the near bed distribution of the vertical turbulent viscosity and diffusion coefficients. An alternative is a numerical solution of the boundary layer equations using a sub-grid embedded in the bottom hydrodynamic grid layer. The distribution of the vertical turbulent viscosity and diffusion coefficients is presumed known form the sub-grid layer solution at the previous time step using (4.26) or (4.27) to determine the turbulent intensity. The sub-grid layer solution proceeds by writing equation (4.16) in finite difference form as
u
k 1 c

k c

cb

H Av
s

(4.70)

where k denotes the sub-grid layer and

1 s

zo ks

(4.71)

is the thickness of the sub-grid layers with ks being the number of sub-grid layers embedded in the bottom grid layer. The integral constraint (4.17) becomes
ks

uck
k 1

k s uc1

(4.72)

where uc1 is the current velocity in the bottom grid layer. Solving the recursion (4.70) and substituting into (4.72) gives
1 c

1 ks

ks

ks
1

H Av
s

k cb

(4.73)
uc1

The velocity profile in the bottom half of the near bed sub-grid layer is assumed to be logarithmic
1 uc cb

(4.74)
ln
s

Inserting (4.74) into (4.73) gives

1 u 1 ks
1 c

ks

ks k
1

H Av
s

(4.75)
1 uc

ln

uc1

which can be solved iteratively for the current velocity in the bottom sub-grid layer when the distribution of the turbulent viscosity at the sub-grid interfaces is known. The recursion (4.70) can then be solved for the velocity in the remaining sub-grid layers. The finite difference form of the sediment boundary layer equation (4.22) is

Kv sH

k s

Kv sH

(4.76)

k s

k 1

J Sb

where equals 1 for upwind settling and 0.5 for central difference settling. The constraint equation is
ks

Sk
k 1

k s S1

(4.77)

For noncohesive sediment transport, the sub-grid near bottom sediment concentration S1 is specified as a function of the bed stress and the bed composition. The sediment flux and primary bottom grid layer concentration, S1, must then be determined. This is accomplished by introducing a dimensionless sediment variable
k

w1 S k s J Sb

(4.78)

Into (4.76) to give

k k 1 k k

(4.79)

where
k

Kv w1 s H s Kv w1 s H s
k

wsk w1 s 1 wsk w1 s

(4.80)

Since S1 is known, the first equation becomes

(4.81)

and (4.78) now represents a closed system of ks-1 equations. The of solution of (4.79) can be written as
k

(4.82)

Which is the sum of the solutions of the two simpler linear systems
k

k
1

(4.83) (4.84)

1 2

k 1 : k

The solution of (4.83) can be written as

k k 1 k k 1 k 1 1

(4.85)

while (4.80) is solved numerically. The dimensionless form of the constraint (4.77) is

1 ks

ks k 1

(4.86)

and can be written as

1 1

(4.87)

where

1 ks
k k 1

ks k k 1

(4.88)

1: k 1
k 1 k 1

:k

2
(4.89)

1 ks

ks k 1

Reverting to the original variables gives the bed flux

J Sb

w1 s

S 1 S1

(4.90)

where S1 can be interpreted as the equilibrium sediment concentration for the bottom layer of the primary vertical grid. The flux relationship (4.90) is used to determine the sediment concentration, S1 in the bottom grid layer, using
S1 S1old
1

J sb

(4.91)

where is the time step for integration of the primary grid equations. The flux is then evaluated and used to determine the vertical sediment concentration distribution in the sub-grid layers using

Sk

S1

J Sb w1 s

: k

(4.92)

Which follows from (4.76), (4.82), and (4.85). The sediment concentration is used to determine the buoyancy distribution in the sub-grid layers. For cohesive sediment resuspension, the bed flux is known as a function of the bed stress and geomechanical bed properties. The sediment concentration in the bottom grid layer, S1, can be determined using (4.91). The ks-1 equations (4.76) supplemented by (4.77) then form a tri-diagonal system linear system, with a zero lower diagonal, supplemented by a full last row. The system is readily solved using the Sherman-Morrison formula (Press, et al., 1992) for the vertical distribution of sediment in the boundary layer subgrid. For cohesive sediment deposition, the bed flux can be represented by

J Sbd

Pd w1S 1 s

(4.93)

where Pd is the probability of deposition which depends on the bed stress and a critical depositional stress. Inserting (4.93) into (4.76) gives a system of ks-1 equations which must be supplement the equation formed by introducing (4.93) and (4.91) into (4.74) or
ks

Sk
k 1

ks

Pd w1S 1 s
b

ks S1old

(4. 94)

The resulting system of linear equations is of tri-diagonal form, with a zero a zero lower diagonal, supplemented by a full first column and a full last row. The system is readily solved using the Sherman-Morrison formula (Press, et al., 1992) for the vertical distribution of sediment in the boundary layer sub-grid.

4.4

Neutral Wave, Current, and Sediment Boundary Layers

Analytical solutions of the wave, current and sediment boundary layer equations (4.11, 4.16, and 4.22), as exemplified most recently by Styles and Glenn (2000), typically assume tractable forms of the vertical turbulent viscosity and diffusivity inside the wavecurrent and the current boundary layers. Closed form solutions, using special mathematical functions, are possible for the neutral case where the sediment concentrations are low enough to assume that the buoyancy is zero. An alternate approach is to extend the numerical sub-grid approach of the previous section to include a numerical solution for the wave boundary layer. with the resulting formulation being applicable to both neutral and sediment stratified conditions. The sub-grid formulation for the wave boundary layer, which is applicable to both neutral and sediment stratified conditions will be presented in the following section, while this section presents a semianalytical solution appropriate for neutral conditions. For both the semi-analytical and sub-grid solution of the wave, current and sediment boundary layers, the turbulent viscosity and diffusion coefficients are assumed to be time invariant with (2.18) and (2.19) written in terms of the root mean square turbulent intensity

Av

Ao

q2 l

(4.95) (4.96)

Kv

Ko

q2 l

Equations (4.26) and (4.36) used to determine the mean square turbulent intensity
u
2

q2

B1 Ao 1 B1 K o

A K

Rq

l2 H2

z c

1 2
wm

z wm

(4.97)
z uc z uwm

cos

Converting the shears in (4.97) to stresses using (4.95) gives

Av H
K 2 z c

2 2 A

Ao

B1 B1 A Ao

1 Av 2 H Av H
z c

2 z wm

(4.98)
u Av H
z wm

K o Rq

cos

wm

Which for neutral conditions reduces to

q2

B14 / 3

Av H 4 cos
c

2 z uc

1 Av 2 H Av H
z c

2 z u wm

(4.99)
u Av H
z wm

wm

The neutral version of the Styles and Glenn (2000) wave current boundary layer formulation defines two regions for the turbulent intensity
q
2

qwc

1/ 3 1

2 cb

1 2

2 wbm

1/ 4

cos

wm

cb wbm

:0

qwc qw

wc

(4.100)

q2

qc

1/ B1 3

2 1/ 4 cb

: z

qwc qw

wc

and three regions for the length scale

l l l
w

zH : zo H :
wc

z z qwc qc

wc

qwc qc
wc

wc

(4.101)

kzH : z

where wc is a characteristic thickness of the wave-current boundary layer relative to the water column depth. The resulting turbulent viscosity distribution is
Avn H Avn H Ao qwc z : zo
wc

z z qwc qc

wc

Ao qwc Avn H

wc

qwc qc
wc

(4.102)
wc

Ao qc kz : z

with corresponding distributions for the vertical turbulent diffusivity. Rather than solving the wave boundary layer velocity distribution using special mathematical functions and then approximating these functions by series expansions, the solution proceeds by introducing an approximate velocity distribution in the lower region
u Re U w1 ln z exp i t zo U w2 z zo exp i t zo

(4.103)
: zo z
wc

And an exact distribution in the constant viscosity central region

uw uw Re U w3 exp z
wc

(4.104)
wc

exp i t

wc

where Uw1, Uw2, and Uw3 are complex constants and

2

H wc Ao qwc

(4.105)

Since

is of order unity, the wave boundary layer scale is on the order of

wc

O O

Ao qwc H Av H2

(4.106)

2 wc

The solution the lower region is obtained by a Galerkin procedure. Substitution of (4.103) into (4.11) gives a residual error:
E i U w1 ln z zo U w2 z zo zo Av U w1 H2 z U w2 zo

(4.107)
i Uw

The Galerkin weighted residual errors are then set to zero

wc

ln
zo
wc

z Edz zo

(4.108)

0 0

zo

z zo Edz zo

Expanding (4.107) and integrating the vertical stress gradient by parts gives
wc

i ln 2
zo

z zo i

Av 1 dzU w1 H 2 z2
wc

wc

i ln
zo

z zo

z zo zo Av H2 U
z

Av 1 dzU w 2 H 2 zo z

(4.109)

ln
zo

z dzU w zo

ln

wc

zo

wc

wc

i
zo

z zo z ln zo zo i
zo

Av 1 dzU w1 H 2 zo z
wc

wc

i
zo

z zo zo Av H2
z

Av 1 dzU w 2 2 H 2 zo

(4.110)

z zo dzU w zo

zo

zo

Uw
wc

Or

a11U w1 a12U w2 a21U w1 a22U w2

Where
wc

b1 b2

(4.111)

a11
zo
wc

ln 2 z zo

z zo

i z zo zo
2

Ao qwc H i

1 dz z Ao qwc H z 2 zo Twi H zo Twi H 1 dz zo

a12

a21
zo
wc

ln

a22
zo
wc

z zo zo

Ao qwc H i ln
wc

(4.112)

b1
zo
wc

ln

z dzU w zo

zo
wc

b2
zo

z zo dzU w zo Twi Av H
z

i U
z

zo

wc

The solution is
wc

U w1 U w2

a11 a21

a12 a22

1 zo
wc

z ln dz zo z zo dz zo

(4.113)

Uw

a11 a21

a12 a22

i ln

wc

zo
wc

zo

Twi H

zo

zo

Or in symbolic form

U w1 U w2

A11U w A21U w

A12 A22

Twi h Twi h

(4.114)

where the complex stress amplitude Twi at the interface between the lower and central regions remains to be determined as does the constant, Uw3, in the central region solution. The two constants, Twi and Uw3, must be determined such that the velocity and its vertical gradient are continuous at the interface between the two boundary layer regions by the solution of
A12 ln Uw A12
wc wc

zo A11U w

A22 ln
wc

wc

zo

zo zo U w3
wc

Twi h

U w3
wc

A21U w A11
wc

zo

(4.115)

zo A21 Uw zo

A22 zo

Twi h

The solution provide the interface stress in terms of the inviscid wave velocity amplitude and in turn allows Uw1 and Uw2 to be expressed in terms of the inviscid wave velocity amplitude. The maximum wave bed stress can then be determined by
wbm

Ao qwc U w1 U w 2

(4.116)

Note that in the absences of currents

qwc
With (4.116) becoming
2 wbm

1/ B1 3 21/ 4

1/ 2 wbm

(4.117)

2 cbw

U w1 U w 2
2

cbw U w
2

(4.118)

U w1 U w 2 Uw
2

The solution for the current velocity is

uc
in the lower region,
uc
cb

cb

Ao qwc

ln

z zo

(4.119)

z
wc

(4.120)
ln
wc

Ao qwc

zo

in the central region, and

uc

cb

Ao qwc ze e
wc

qwc z ln qc ze
qc qwc

(4.121)

qwc zo e qc wc

in the upper region. To enforce the integral constraint, the current profile is integrated over the three regions to give
wc

cb

Ao qwc

ln
zo

z dz zo
wc

(4.122)
cb

Ao qwc
qwc qc

wc ln

wc

zo

wc

zo

(4.123)
z
wc

cb

Ao qwc
cb wc

ln qwc qc
m

wc

wc

zo

1 dz 1 2

Ao qwc

1 qwc 2 qc

qwc 1 ln qc

wc

zo

qwc Ao qwc k qc
cb cb

ln
wc

qwc qc

z dz ze m ze qwc ze qc
wc

(4.124)

qwc m Ao qwc k qc qwc Ao qwc k qc

cb 2 wc

ln

ln

With the general integral constraint being

m

m
cb

zo
k 1 wc

uck zo
wc

(4.125)

Ao qwc
cb wc

wc

ln

zo qwc qc ln

wc

Ao qwc

1 qwc 2 qc
cb

qwc 1 ln qc m ze qwc ze qc
wc

zo

1 2

Ao qwc
cb

qwc m qc qwc qc
2 wc

Ao qwc

ln

For the thickness of the wave-current layer exceeding the lower hydrodynamic grid layer. When the wave-current boundary remains inside the bottom layer of the hydrodynamic grid, (4.125) reduces to
3 qwc 2 qc
wc

wc

wc

Ao qcuc1

cb

qc qwc qwc qc

wc

(4.126)
zo qwc ze qc
wc

2
wc

zo

ln

wc

zo

ln

ze

ln

Introducing (4.100) into (4.126) gives

cb

cbc uc1 uc1

2

(4.127)
2 2

cbc ln ze 1
wc

zo 2

qwc qc

qc qwc

zo

an expression for the current stress and bottom current friction coefficient.

4.5

Stratified Wave, Current, and Sediment Boundary Layers

In this case, the turbulent intensity and vertical turbulent transfer coefficients become

q2 Avn
K vn

B12/ 3 Ao ql

2 xz

2 1/ 2 yz

(4.37) (4.38) (4.39)

2 xz

2 1/ 4 yz

l
l

K o ql

Ko Ao

2 xz

2 1/ 4 yz

The neutral version of the Styles and Glenn (2000) wave current boundary layer formulation defines two regions for the turbulent intensity
q qwc B
1/ 3 1 2 cb

1 2 q

2 wbm

1/ 2

cos

wm

cb wbm

: 0

qwc qw

(4.39)
wc

qc

1/ B1 3

2 1/ 2 cb

: z

qwc qw

wc

and three regions for the length scale

l l l
w

z : zo :
wc

z z qwc qc

wc

(4.39)
wc

qwc qc
wc

kz : z

where wc is a characteristic thickness of the wave-current boundary layer. The resulting turbulent viscosity distribution is

Avn Avn

Ao qwc z : zo
w

z z qwc qc

wc

(4.39)
wc

Ao qwc Avn

wc

qwc qc
wc

Ao qc kz : z

with a corresponding distributions for the vertical turbulent diffusivity. For stratified boundary layers, Styles and Glenn modify the neutral transfer coefficients using a MoninObokov length based stability function which leads to none closed form solutions of the wave, current and sediment boundary layer equations.

Kv
Kwc

Kwc wc wc

q l

Ko 1 R3 1 Rqwc
gH
z 2 qwc

(2.19)

Rqwc

wc

2 lwc H2

(2.20)

4 qwc

B1
Awc

1 B1

Kwc

Rqwc

2 cb

1 2

2 wbm

(4.26)

cos

wm

cb wbm

and inside the current boundary layer above the wave boundary layer

Av
Ac

Ac c c

ql

Ao 1 R1 1 Rqc 1 R2 1 Rqc 1 R3 1Rqc

Kv
Kc Kc c c

(2.18)

ql

Ko 1 R3 1 Rqc
gH qc2
z

(2.19)

Rqc

lc2 H2
1 2 cb

(2.20)

qc4

B1
Ac

(4.26)

1 B1

Kc Rqc

5.

Sediment Bed Mass Conservation, Armoring and Consolidation (final revision 05/21/2003)

The general conservation of mass for bed sediment has the form
t

SiB

i k , kt J SB

i k , kt J PA

i k , kt 1 J PA

(5.1)

where S is the mass concentration of per total volume of a bed layer k, B is the layer thickness, JSB is the net sediment mass flux, mass per unit area and unit time, positive from the bed to the water column, A is an armoring parameter (1 for armoring, 0 otherwise), and JPA is the parent to armoring layer flux when the top or surface layer of the bed, kt, acts to simulate armoring. The superscript i denotes the ith sediment size-type class. The sediment concentration can also be defined by

Fi 1

i s

(5.2) is the

where F is the sediment volume fraction, s is the sediment particle density, and void ratio. The sediment volume fraction is defined by

i k i

SiB
i s k

SiB
i s k

(5.3)

Assuming that the sediment particles are incompressible (5.1) can be alternately expressed by

FiB 1

k , kt
k

i J SB i s A

k , kt

i J PA i s A

k , kt 1

i J PA i s

(5.4)

Summing (5.4) over the sediment size classes gives

B
t

k , kt
k i

i J SB i s A

k , kt
i

i J PA i s A

k , kt 1
i

i J PA i s

(5.5)

The conservation of water volume in a bed layer is given by

B
t

qw:k
k

qw:k
i dep i

k , kt
i

i kt

max

i J SB i s

,0
i J PA i s

min

i J SB i s

,0
i J PA i s

(5.6)

k , kt
A

k , kt 1

i kt 1 i

max

,0
i

i kt

min

,0

Where without the i superscript is the bulk void ratio of the bed layer, and s with superscripts i denote sediment class void ratios required by the mixed material consolidation formulation to be subsequently discussed. Equations (5.5) and (5.6) combine to give
t i kt

Bk
i s

qw:k
i J SB

qw:k 1
i i dep

k , kt
i

max

,0

min

i J SB i s

,0

k , kt
i

i kt 1

max

i J PA i s

,0
i

i kt

min

i J PA i s

,0

(5.7)

k , kt 1
i

i kt 1

max

i J PA i s

,0
i

i kt

min

i J PA i s

,0

The solution procedure for the bed uses a fractional step approach. The first step involves deposition and resuspension while the second step involves pore water flow and consolidation. 5. 1 Deposition, Resuspension, and Armoring

The discrete deposition and resuspension step for the sediment class i mass conservation equation (5.1) is
SiB
* k

SiB

n k

i k , kt J SB

i k , kt J PA

i k , kt 1 J PA

(5.8)

Or
FiB 1
*

FiB 1

k , kt
k

i J SB i s A

k , kt

i J PA i s A

k , kt 1

i J PA i s

(5.9)

The corresponding discrete forms of (5.5), (5.6) and (5.7) are

B 1

B 1

k , kt
k i

i J SB i s i J PA i i s A

k , kt
i

i J PA i s

(5.10)

k , kt 1

B 1
A

B 1 k , kt

k , kt
k i

i kt

max

i J SB i s

,0

i dep

min

i J SB i s

,0
(5.11)

k , kt 1
i
* Bk

i kt 1

max

i PA i s

,0

i kt

min

i PA i s

,0

Bkn
i s

k , kt
i

i kt

max
i kt 1

i J SB

,0

i dep

min
i kt

i J SB i s

,0

k , kt
i

max

i J PA i s

,0

min
i kt

i J PA i s

,0

(5.12)

k , kt 1
i

i kt 1

max

i J PA i s

,0

min

i J PA i s

,0

When the armoring option is inactive, the deposition and resuspension step operates only on the top layer of the bed with (5.8) solved for the new top layer sediment mass per unit area
SiB
* kt

SiB

n kt

i J SB

(5.13)

using a known sediment depositional or resuspension flux. If the flux in (5.13) is positive, representing resuspension, it is limited over the time step by
i J SB i min J SBR , 1

SiB

n kt

(5.14)

where the subscript SBR represents the predicted resuspension flux. Following the solution of (5.13) for each sediment class, equations (5.12) is solved for the new top layer thickness and (5.10) is solved for the new top layer void ratio. When the noncohesive sediment armoring option is active, equation (5.8) is applied to the top two layers of the bed

SiB

* kt * kt 1

SiB

n kt n

i J SB

i J PA i J PA

(5.15a) (5.15b)

SiB

SiB

kt 1

with the flux limiter (5.14) being applied to (5.15a) for resuspension flux from the top layer. Two options exist for determining the parent to active layer flux. One option is to require that the total mass of sediment in the surface, active layer remains constant during the deposition-resuspension step. The total parent to active layer flux is then given by (5.10) as
i J PA i i s i i J SB i s

(5.16)

The class fluxes can then be assigned by

i J PA i s

i kt 1

max
i

i J SB i s

,0

F min
i

i kt

i J SB i s

,0

(5.17)

allowing (5.15) to be updated. Equation (5.10) and (5.12) are then solved for the new thicknesses and void ratios of the parent and active layer. Another option is to require that the thickness of the active layer to be time invariant during the deposition and resuspension step. Equation (5.12) reduces to

1
i

i kt 1

max
i kt

i J PA i s

,0
i i SB i s

1
i dep

i kt

min J

i J PA i s i SB i s

,0
(5.18)

1
i

max

,0

min

,0

The sediment class fluxes can be assigned by

i J PA i s i Fkt 1 i kt 1

1 1
i

max QSW , 0 max J

i SB i s

i Fkt

1 1

i kt

min QSW , 0
(5.19)

QSW

i kt

,0

i dep

min

i SB i s

,0

allowing solution of equations (5.15), (5.10), and (5.12).

5. 2

Consolidation of Homogeneous Sediment Beds

This section discusses options for representing consolidation of sediment beds containing either cohesive sediment or a mixture of noncohesive sediments defined by multiple size classes. Mixed cohesive and noncohesive bed consolidation is discussed in the subsequent section. For the second, consolidation half step, the sediment mass per unit area and the sediment volume per unit area remain constant, with (5.1) and (5.5) giving
SiB
n 1 k

SiB

* k

(5.20)
*

B 1

n 1

B 1

(5.21)
k

The second half step for the water volume conservation equation (5.6) is

B 1

n 1

B 1

qw:k
k

qw:k

(5.22)

Equations (5.21) and (5.22) can be combined to give

Bkn
1 * Bk

qw:k

qw:k

(5.23)

an equation for the layer thickness, and

n 1 k * k

1 B

n 1

qw:k
k

qw:k

(5.24)

an equation for the void ratio. The EFDC model includes four options for consolidation and pore water flow. The first option is a constant porosity bed, with (5.24) giving
qw:k qw:k qGW

(5.25)

which indicates that the pore water specific discharge is equal to a specified groundwater specific discharge at the bottom of the lowest bed layer. The second option is a simple consolidation model based on relaxation of the vertical void ratio profile to a specified profile given by
m o m

exp

t to

(5.26)

where c is a consolidation rate coefficient, and m is an ultimate minimum void ratio, which can be dependent on the vertical position in the bed. Evaluating (5.26) at two successive time levels gives
n m o m

exp
exp
c

n
n

to
to

(5.27) (5.28)

n 1 m o m

Taking their ratio gives

n 1 m * m

exp

(5.29)
c

or
n 1 k m * k m

exp

(5.30)

Using the new void ratio given by (5.30), the new bed layer thickness is updated by (5.21). The pore water specific discharges are then given by recursively solving (5.23)
qw:k qw:k
1

Bkn

* Bk

(5.31)

From k = 1, kt using
qw:1 qGW

(5.32)

The third option for consolidation and pore water flow is based on finite strain consolidation theory. Use of this option requires the bed layers to be composed of either cohesive or noncohesive sediment, such that a single set of constitutive relationships are used over the entire thickness of the bed. The specific discharges in (5.23) or (5,24) are determined using the Darcy equation
q K g w u

(5.33)

where K is the hydraulic conductivity and u is the excess pore pressure defined as the difference between the total pore pressure, ut, and the hydrostatic pressure, uh.
u ut uh

(5.34) and

The total pore pressure is defined as the difference between the total stress effective stress e.

ut

(5.35)

The total stress and hydrostatic pressure are given by

zb

pb

g
z

1 1
uh pb
w

1
w

Fi
i

i s

(5.36)
dz

zb

(5.37)

where pb is the water column pressure at the bed surface zb. Solving for the excess pore pressure using (5.34) through (5.37) gives
zb

1 1

(5.38)
s w

w z

dz

where
s i

Fi

i s

(5.39)

is the average sediment density. The specific discharge (5.33), can alternately be expressed in terms of the effective stress

q
or the void ratio

K g w

s z e w

K 1

(5.40)

K g w

e z

s w

K 1

(5.41)

where d /d c is a coefficient of compressibility. For consistency with the Lagrangian representation of sediment mass conservation, a new vertical coordinate , defined by

d dz

1 1

(5.42)

is introduced, with (5.40) and (5.41) becoming

qw

K 1 g

e w

s w

K 1

(5.43)

and

qw
Where is a length scale

K 1

s w

K 1

(5.44)

1 g
w

(5.45)

The consistency of (5.43) and (5.44) at bed layer interfaces also requires consideration. The finite difference form of (5.33) in the transformed coordinate, defined by (5.42), at an interface between bed layers can be written as

q
below the interface and

2 g
w

K 1
k

ui uk
k

(5.46)

q
above the interface, where

2 g
w

K 1
k 1

uk

ui

(5.47)

k 1

Bk

(5.48)

is the transformed coordinate thickness of the bed layer. Solving (5.47) for the interface excess pore pressure and inserting the results into (5.46) gives
q K 1
k 1/ 2

2 g
w

uk

uk
k

(5.49)

k 1

where
1
k 1 k

1 K
k 1 k 1/ 2

1 K
k k 1

(5.50)
K
k

defines the hydraulic conductivity at the bed layer interface between layers k and k+1. The discrete from of (5.38) in the transformed coordinate is

u g
w k 1

u g
w k e

s w

2
e

k 1

(5.51)

k 1

which after introduction into (5.49) gives

K 1
k 1/ 2

2
k 1 k s k 1/ 2 w

e k 1

(5.52)

K 1
where
s w

1
k 1/ 2

1
k 1/ 2

1
k 1 k 1 k

s w

1
k 1

s k w

1
k

(5.53)

In terms of the void ratio, (5.52) is

K 1
k 1/ 2

k 1/ 2 k 1 k k

K 1
k 1/ 2

s w

1
k 1/ 2

(5.54)

k 1

where

1
k 1/ 2

e,k 1 w k 1

e,k k

(5.55) For

The effective stress and hydraulic conductivity are functions of the void ratio. cohesive material
e eo

exp
eo

(5.56)

exp

K Ko

exp
K

(5.57)

are the simplest functional relationships consistent with observational data. Figures 5.1 through 5.4, based on data presented in Cargill (1985) and Palermo et al., (1998) confirm these choices. However, they show essentially two regions of behavior, below and above a void ratio of approximately 6. For noncohesive material the linear relationships
e eo e eo

(5.58)

K Ko
are appropriate.

1
K

(5.59)

Given the unique dependence of the specific discharge on the void ratio, the void ratio form of the consolidation step, (5.24) is selected for the solution, with the thickness of the bed layers then determined by (5.23). The specific discharges at the top and bottom of layer k, follow from (5.54) and are given by
qw:k qw:k 2
k 1 k k k k 1

K 1 K 1
k k k 1 k 1 k k

s w s w

1
k

K 1 K
k k

(5.60)

2
k

For the bottom layer of the bed,

qw:1

qgwi

(5.61)

where qgwi is a known specific discharge due to groundwater interaction. For the top layer of the bed, two alternate formulations are possible. The first formulation assumes that the void ratio at the water column-sediment bed interface is specifed by dep, with (5.60a) modified to

qw:kt

kt kt

K 1
dep kt k

s w

1
kt

K 1
kt

(5.62)

The second formulation assumes that the excess pore pressure, u, at the water columnsediment bed interface is zero with (5.46) giving

qw:kt

K 1
Kt

2
Kt

u g
w Kt

(5.63)

Using (5.38) the excess pore pressure at the midpoint of the top layer is

u g
w

s w

(5.64)

Which combines with (5.63) to give

qw:kt

K 1
Kt

2
Kt

K 1
Kt

s w

(5.65)

1
Kt

Equation (5.65) can be expressed in terms of the void ratio at the new time level n+1 by expanding the effective stress at time level n+1 in a Taylor series
n 1 e n e n e n 1 *

(5.66)

Substituting (5.61) into (5.65) gives

qw:kt K 1 1
Kt

2
Kt

* n 1

K 1
Kt * * Kt n e

s w

1
Kt

(5.67)

K
Kt

2
Kt

The numerical values of the various parameters in the expressions for the specific discharge indicate that an implicit solution of (5.24) is necessary. This is done in two stages with an intermediate void ratio, denoted by **, determined by substituting the internal specific discharges, written as

qw:k qw:k

2
k 1

* k k * k k 1

K 1 K 1

* ** k 1 k * ** k k k 1 w s k w s

1
k

K 1 K
k

k *

(5.68)

2
k

and one of the surface specific discharges corresponding to (5.62)

qw:kt

K 1

* dep kt k

* ** s w

1
kt

K 1

(5.69)

kt

kt

or (5.67)

qw:kt

K 1

2
*

* ** kt

K 1
e

* s kt * w

1
(5.70)

kt

K 1
into
** k * k

2
kt

kt

kt

1 2 B

**

(5.71)

qw:k
k

qw:k

and solving the resulting tri-diagonal system of equations. The specific discharges are then exactly calculated using (5.68) and (5.69) or (5.70). The new time level thickness of the layers is determined by (5.23) with the void ratios determined from (5.24). The linearized form of this scheme is unconditionally stable.

5. 3

Consolidation of Mixed Cohesive and Noncohesive Sediment Beds

This section presents a methodology for representing consolidation of sediment beds containing both cohesive and noncohesive sediments. The methodology allows for both cohesive and noncohesive sediment in any bed layer and is based on the following assumptions. First, it is assumed that during the consolidation step, a fraction of the bed pore water volume per unit horizontal area is associated with each sediment type or
B 1
wc wn

(5.72)

where the subscripts wc and wn denote water associated with cohesive and noncohesive sediment, respectively. Likewise the volume of sediment per unit horizontal area can be fractionally partitioned between cohesive and noncohesive

B 1
sc sn

(5.73)

Following the Lagrangian formulation of the previous section, the total volume of sediment and the fractional sediment volume in a bed layer remain constant during a consolidation step.
t

sc

sn

(5.74)

Fractional void ratios can also be defined

wc c sc

(5.75)

wn n sn

(5.76)

And using (5.72) and (5.73), the void ratio of the mixture is
sc c sc sn n sn

(5.77)

Is the sediment volume weighted average of the void ratios of the two sediment types. The second assumption is that during the consolidation time step, the fraction of water associated with noncohesive sediment remains constant, as does the fractional void ratio. This is equivalent to the assuming that the portion of the bed layer associated with noncohesive sediment is incompressible, and that the pore water associated the noncohesive sediment is specified by n. Consistent with the preceding assumptions, the thickness of the bed layer can be divided into cohesive and noncohesive fractions, Bc and Bn, respectively.
Bc Bn
wc wn sc sn

B B

1 1

c n

sc

B B

(5.78)

sn

The hydraulic conductivity of the layer can be expressed by

Bc Bc Kc

Bn Bn Kn

(5.79)

Which is equivalent to an infinite number of alternating infinitesimal cohesive and noncohesive sublayers of proportional thickness comprising the mixed bed layer. Equation (5.79) can be written as

K 1 f sc 1 Kc
c

1 f sn 1 Kn
n

(5.80)

where

f sc
sc

sc sn sn sc sn

(5.81)

f sn

Are the time invariant total cohesive and noncohesive sediment fractions in the bed layer. Likewise, (5.77) can be write as
f sc
c

f sn

(5.82)

The final assumption for the mixed material consolidation formulation is that changes in effective stress are due entirely to changes in the cohesive void ratio. Under this assumption, the specific discharge given by (5.54) can be written as
q K 1
k 1/ 2

2 K 1

k 1/ 2 k s k 1/ 2 w

f sc 1

c k 1

f sc

c k

(5.83)

k 1

k 1/ 2

with (5.55) becoming

1
k 1/ 2 e,k 1 w e,k

f sc

c k 1

f sc

(5.84)
c k

The other layer interface quantities in (5.83) remain defined by (5.50) and (5.53). When the depositional void ratio is specified for the surface layer, (5.62) is modified to

qw:kt

kt kt

K 1
c dep kt c k

s w

1
kt

K 1
kt

(5.85)

When the zero excess pore pressure boundary condition at the bed surface is used, (5.67) becomes
qw:kt K 1
Kt

2
Kt

f sc

n 1
C

K
Kt

s Kt * c Kt w

1
*

1
Kt

(5.86)

K 1
Kt

2
Kt

n e

f sc

Equation (5.71) for updating the void ratio is modified using (5.82) to give

f sc

** c k

f sc

* c k

1 2 B

**

(5.87)

qw:k
k

qw:k

Thus the mixed bed layer consolidation formulation essentially solves of the space and time evolution of fsc c with the continuum constitutive relationship for given by

1 f sc

(5.88)
w

The formulation has the desirable characteristic of reducing to the well established coheasive formulation in the absence of noncohesive material. The solution for fsc c proceeds by introducing (5.83) and (5.85) or (5.86) into (5.87) and solving the resulting tridiagonal sytem of equations. The new specific discharges are then directly calculated using (5.83) and (5.85) or (5.86) and used to update the layer thickness using (5.23)
Bkn
1 * Bk

qw:k

qw:k

(5.23)

Equation (5.21) is then used to solve for the void ratio

B 1

n 1

B 1

(5.21)
k

Followed by the solution of (5.82) for the cohesive void ratio

f sn f sc
n

(5.82)

1e+02 1e+01 1e+00 1e-01 1e-02 1e-03 1e-04 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Void Ratio

Figure 5.1. Specific Weight Normalized Effective Stress Versus Void Ratio.
1e+02 1e+01 1e+00 1e-01 1e-02 1e-03 1e-04 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Void Ratio

Figure 5.2. Compress Length Scale, g

1 w

/d

, Versus Void Ratio.

1e-01 1e-02 1e-03 1e-04 1e-05 1e-06 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Void Ratio

Figure 5.3. Hydraulic Conductivity Versus Void Ratio.

1e-02

1e-03

1e-04

1e-05

1e-06 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Void Ratio

Figure 5.4. Hydraulic Conductivity/(1 + Void Ratio) Versus Void Ratio

6.

Noncohesive Sediment Settling, Deposition and Resuspension

Noncohesive inorganic sediments settle as discrete particles, with hindered settling and multiphase interactions becoming important in regions of high sediment concentration near the bed. At low concentrations, the settling velocity for the jth noncohesive sediment class corresponds to the settling velocity of a discrete particle:

wsj

wsoj

(6.1)

Useful expressions for the discrete particle settling velocity which depends on the sediment density, effective grian diameter, and fluid kinematic viscosity, provide by van Rijn (1984b) are:

wsoj g 'd j

: d 100 m 18 10 2 1 0.01Rdj 1 : 100 m d j Rdj 1.1 : d j 1000 m

Rdj

(6.2)

1000 m

where

g'

sj w

(6.3)

is the reduced gravitational acceleration and

Rdj

d j g 'd j

(6.4)

is a the sediment grain densimetric Reynolds number. At higher concentrations and hindering settling conditions, the settling velocity is less than the discrete velocity and can be expressed in the form
I

wsj

1
i

Si
si

(6.5)
wsoj

where s is the sediment particle density with values of n ranging from 2 (Cao et al., 1996) to 4 (Van Rijn, 1984). The expression (6.2) is approximated to within 5 per cent by
I

wsj

1 n
i

Si
si

(6.6)

wsoj

for total sediment concentrations up to 200,000 mg/liter. For total sediment concentrations less than 25,000 mg/liter, neglect of the hindered settling correction results in less than a 5 per cent error in the settling velocity, which is well within the range of uncertainty in parameters used to estimate the discrete particle settling velocity. Noncohesive sediment is transported as bed load and suspended load. The initiation of both modes of transport begins with erosion or resuspension of sediment from the bed when the bed stress, b, exceeds a critical stress referred to as the Shield's stress, cs. The Shield's stress depends upon the density and diameter of the sediment particles and the kinematic viscosity of the fluid and can be expressed in empirical dimensionless relationships of the form:
csj csj 2 u*csj

(6.7)

g 'd j

g 'd j

f Rdj

Useful numerical expressions of the relationship (6.5), provided by van Rijn (1984b), are:
2 0.24 Rdj/ 3 2 0.14 Rdj/ 3 csj 0.64 1 2 : Rdj/ 3

: 4 : 10 : 20

2 Rdj/ 3 10

(6.8)

0.04 R

2/3 dj

0.1

2/3 dj

20

2 0.013 Rdj/ 3

0.29

2 Rdj/ 3 150

2 0.055 : Rdj/ 3 150

A number of approaches have been used to distinguish wheather a particular sediment size class is transported as bed load or suspended load under specific local flow conditions characherized by the bed stress or bed shear velocity:

u*

(6.9)

The approach proposed by van Rijn (1984a) is adopted in the EFDC model and is as follows. When the bed velocity is less than the critical shear velocity
u*csj
csj

g 'd j

csj

(6.10)

no erosion or resuspension takes place and there is no bed load transport. Sediment in suspension under this condition will deposit to the bed as will be subsequently discussed. When the bed shear velocity exceeds the critial shear velocity but remains less than the settling velocity,

u*csj

u*

wsoj

(6.11)

sediment will be eroded from the bed and transported as bed load. Sediment in suspension under this condition will deposit to the bed. When the bed shear velocity exceeds both the critical shear velocity and the settling velocity, bed load transport ceases and the eroded or resuspended sediment will be transported as suspended load. These various transport modes are further illustrated by reference to Figure 1, which shows dimensional forms of the settling velocity relationship (6.2) and the critical Shield's shear velocity (6.10), determined using (6.8) for sediment with a specific gravity of 2.65. For grain diameters less than approximately 1.3E-4 m (130 um) the settling velocity is less than the critical shear velocity and sediment resuspend from the bed when the bed shear velocity exceeds the critical shear velocity will be transported entirely as suspended load. For grain diameters greater than 1.3E-4 m, eroded sediment be transported by bed load in the region corresponding to (6.11) and then as suspended load when the bed shear velocity exceeds the settling velocity. In the EFDC model, the preceding set of rules are used to determine the mode of transport of multiple size classes of noncohesive sediment. Bed load transport is determined using a general bed load transport rate formula:
qB sd g d ,

(6.12)
cs

where qB is the bed load transport rate (mass per unit time per unit width) in the direction of the near bottom horizontal flow velocity vector. The function depends on the Shield's parameter
b

g 'd j

u*2 g 'd j

(6.13)

and the critical Shield's parameter defined by (6.7) and (6.8). A number of bed load transport formulas explicitly incorporate the settling velocity. However, since both the critical Shield's parameter and the settling velocity are unique functions of the sediment grain densimetric Reynolds number, the settling velocity can also be expressed as a function of the critical Shield's parameter with (6.12) remaining an appropriate representation. A number of bed load formulations developed for riverine prediction (Ackers and White, 1973; Laursen, 1958; Yang, 1973; Yang and Molinas, 1982) do not readily conform to

(1) and were not incorporated as options in the EFDC model. Two widely used bed load formulations which do conform to (6.12) are the Meyer-Peter and Muller (1948) and Bagnold (1956) formulas and their derivatives (Raudkivi, 1967; Neilson, 1992; Reid and Frostick, 1994) which have the general form
,

(6.14)
cs cs cs

where
cs

or

Rd

(6.15)

The Meyer-Peter and Muller formulations are typified by

3/ 2 cs

(6.16)

while Bagnold formulations are typified by (6.17)

cs cs

with Bagnold's original formula having equal to zero. The Meyer-Peter and Muller formulation has been extended to heterogeneous beds by Suzuki et al. (1998), while Bagnold's formula has been similarly extended by van Niekerk et al. (1992). The bed load formulation by van Rijn (1984a) having the form
2.1 cs

0.053 1/ 2.1 Rd 5 cs

(6.18)

has been incorporated into the CH3D-SED model and modified for heterogeneous beds by Spasojevic and Holly (1994). Equation (6.18) can be implemented in the EFDC model with an appropriately specified . A modified formulation of the Einstein bed load function (Einstein, 1950) which conforms to (6.12) and (6.14) has been presented by Rahmeyer (1999) and will be later incorporated into the EFDC model. The procedure for coupling bed load transport with the sediment bed in the EFDC model is as follows. First, the magnitude of the bed load mass flux per unit width is calculated according to (6.12) at horizontal model cell centers, denoted by the subscript c. The cell center flux is then transformed into cell center vector components using

qbcx qbcy

u u 2 v2 v u 2 v2

qbc qbc

(6.19)

where u and v are the cell center horizontal velocities near the bed. Cell face mass fluxes are determined by down wind projection of the cell center fluxes

qbfx qbfy

qbcx qbcy

upwind

(6.20)

upwind

where the subscript upwind denotes the cell center upwind of the x normal and y normal cell faces. The net removal or accumulation rate of sediment material from the deposited bed underlying a water cell is then given by:
mx m y J b m y qbfx
e

m y qbfx

mx qbfy

mx qbfy

(6.21)
s

where Jb is the net removal rate (gm/m2-sec) from the bed, mx and my are x and y dimensions of the cell, and the compass direction subscripts define the four cell faces. The implementation of (6.19) through (6.21) in the EFDC code includes logic to limit the out fluxes (6.20) over a time step, such that the time integrated mass flux from the bed does not exceed bed sediment available for erosion or resuspension. Under conditions when the bed shear velocity exceeds the settling velocity and critical Shield's shear velocity, noncohesive sediment will be resuspended and transported as suspended load. When the bed shear velocity falls below both settling velocity and the critical Shield's shear velocity, suspended sediment will deposit to the bed. A consistent formulation of these processes can be developed using the concept of a near bed equilibrium sediment concentration. Under steady, uniform flow and sediment loading conditions, an equilibrium distribution of sediment in the water column tends to be established, with the resuspension and deposition fluxes canceling each other. Using a number of simplifying assumptions, the equilibrium sediment concentration distribution in the water column can be expressed analytically in terms of the near bed reference or equilibrium concentration, the settling velocity and the vertical turbulent diffusivity. For unsteady or spatially varying flow conditions, the water column sediment concentration distribution varies in space and time in response to sediment load variations, changes in hydrodynamic transport, and associated nonzero fluxes across the water column-sediment bed interface. An increase or decrease in the bed stress and the intensity of vertical turbulent mixing will result in net erosion or deposition, respectively, at a particular location or time.

To illustrate how an appropriate suspended noncohesive sediment bed flux boundary condition can be established, consider the approximation to the sediment transport equation (3.1) for nearly uniform horizontal conditions
HS Kv H S ws S

(6.22)

Integrating (6.22) over the depth of the bottom hydrodynamic model layer gives
t

HS

J0

(6.23) .

where the over bar denotes the mean over the dimensionless layer thickness, Subtracting (6.23) from (6.22) gives
HS Kv H S ws S J0 J

(6.24)

Assuming that the rate of change of the deviation of the sediment concentration from the mean is small
t

HS

HS

(6.25)

allows (6.24) to be approximated by

Kv H

S ws S

J0

(6.26)

Integrating (6.26) once gives

Kv H

S ws S

J0

J0

(6.27)

Very near the bed, (6.27) can be approximated by

Kv H

S ws S

J0

(6.28)

Neglecting stratification effects and using the results of Chapter 4, the near bed diffusivity is approximately

Kv H

Ko q

l H

u* z

(6.29)

Introducing (6.29) into (6.28) gives

S R S z R Jo z ws

(6.30)

where
R ws u*

(6.31)

is the Rouse parameter. The solution of (6.30) is

S Jo ws C zR

(6.32)

The constant of integration is evaluated using

Seq : z

zeq

and

Jo

(6.33)

which sets the near bed sediment concentration to an equilibrium value, defined just above the bed under no net flux condition. Using (6.33), equation (6.32) becomes
S zeq z
R

Seq

Jo ws

(6.34)

For nonequilibrium conditions, the net flux is given by evaluating (6.34) at the equilibrium level
Jo ws Seq Sne

(6.35)

where Sne is the actual concentration at the reference equilibrium level. Equation (6.35) clearly indicates that when the near bed sediment concentration is less than the equilibrium value a net flux from the bed into the water column occurs. Likewise when the concentration exceeds equilibrium, a net flux to the bed occurs. For the relationship (6.35) to be useful in a numerical model, the bed flux must be expressed in terms of the model layer mean concentration. For a three-dimensional application, (6.34) can be integrated over the bottom model layer to give
Jo ws Seq S

(6.36)

where

Seq

ln

zeq1 zeq1 1
1 R

Seq : R 1
(6.37)

zeq1 Seq 1 R

1 Seq : R 1

zeq1 1

defines an equivalent layer mean equilibrium concentration in terms of the near bed equilibrium concentration. The corresponding quantities in the numerical solution bottom boundary condition (3.6) are

wr Sr Pd ws

ws Seq ws

(6.38)

If the dimensionless equilibrium elevation, zeq exceeds the dimensionless layer thickness, (6.19) can be modified to

Seq

ln M zeq1 M zeq1 1 M z
1 1 R eq

Seq : R 1
(6.39)

1 Seq : R 1

Seq

1 R M zeq1 1

where the over bars in (6.36) and (6.38) implying an average of the first M layers above the bed. For two-dimensional, depth averaged model application, a number of additional considerations are necessary. For depth average modeling, the equivalent of (6.27) is

Kv H

S ws S

Jo 1 z

(6.40)

Neglecting stratification effects and using the results of Chapter 4, the diffusivity is

Kv H

Ko q

l H

u* z 1 z

(6.41)

Introducing (6.41) into (6.40) gives

zS

R z 1 z

R 1 z z

Jo ws

(6.42)

A close form solution of (6.42) is possible for equal to zero. Although the resulting diffusivity is not as reasonable as the choice of equal to one, the resulting vertical distribution of sediment is much more sensitive to the near bed diffusivity distribution than the distribution in the upper portions of the water column. For equal to zero, the solution of (6.42) is
S 1 Rz 1 R Jo ws C zR

(6.43)

Evaluating the constant of integration using (6.43) gives

zeq z

Seq

Rz 1 R

Jo ws

(6.44)

For nonequilibrium conditions, the net flux is given by evaluating (6.44) at the equilibrium level

Jo

ws

1 R 1 R 1 zeq

(6.45)

Seq

Sne

where Sne is the actual concentration at the reference equilibrium level. Since zeq is on the order of the sediment grain diameter divided by the depth of the water column, (6.45) is essentially equivalent (6.35). To obtain an expression for the bed flux in terms of the depth average sediment concentration, (6.44) is integrated over the depth to give
21 R 2 R 1 zeq

(6.46)
Seq S

Jo

ws

where
ln zeq1 zeq1 1 z
R 1 eq

Seq

Seq : R 1

(6.47)
Seq : R 1

Seq

1 R zeq1 1

The corresponding quantities in the numerical solution bottom boundary condition (3.6) are

wr S r

ws

21 R 2 R 1 zeq 21 R 2 R 1 zeq

Seq
(6.48)

Pd ws

ws

When multiple sediment size classes are simulated, the equilibrium concentrations given by (6.37), (6.39), and (6.47) are adjusted by multiplying by their respective sediment volume fractions in the surface layer of the bed. The specification of the water column-bed flux of noncohesive sediment has been reduced to specification of the near bed equilibrium concentration and its corresponding reference distance above the bed. Garcia and Parker (1991) evaluated seven relationships, derived by combinations of analysis and experiment correlation, for determining the near bed equilibrium concentration as well as proposing a new relationship. All of the relationships essential specify the equilibrium concentration in terms of hydrodynamic and sediment physical parameters
Seq Seq d ,
s

, ws , u* ,

(6.49)

including the sediment particle diameter, the sediment and water densities, the sediment settling velocity, the bed shear velocity, and the kinematic molecular viscosity of water. Garcia and Parker concluded that the representations of Smith and McLean (1977) and Van Rijn (1984b) as well as their own proposed representation perform acceptably when tested against experimental and field observations. Smith and McLean's formula for the equilibrium concentration is

Seq
where
o

0.65 oT 1 oT

(6.50)

is a constant equal to 2.4E-3 and T is given by

b cs

cs

2 2 u* u*cs u*2cs

(6.51)

where b is the bed stress and cs is the critical Shields stress. The use of Smith and McLean's formulation requires that the critical Shields stress be specified for each sediment size class. Van Rijn's formula is
Seq 0.015 d 3/ 2 1/ 5 T Rd * zeq

(6.52)

where zeq* ( = Hzeq ) is the dimensional reference height and Rd is a sediment grain Reynolds number. When Van Rijn's formula is select for use in EFDC, the critical Shields stress in internally calculated using relationships from Van Rijn (1984b). Van Rijn suggested setting the dimensional reference height to three grain diameters. In the EFDC model, the user specifies the reference height as a multiple of the largest noncohesive sediment size class diameter. Garcia and Parker's general formula for multiple sediment size classes is

S jeq

A
s

Zj

5 5

(6.53)

1 3.33 A
u* 3/ 5 Rdj FH wsj
dj d 50
1/ 5

Zj

(6.54)

(6.55)

FH

1
o

(6.56)

where A is a constant equal to 1.3E-7, d50 is the median grain diameter based on all sediment classes, is a straining factor, FH is a hiding factor and is the standard deviation of the sedimentological phi scale of sediment size distribution. Garcia and Parker's formulation is unique in that it can account for armoring effects when multiple sediment classes are simulated. For simulation of a single noncohesive size class, the straining factor and the hiding factor are set to one. The EFDC model has the option to simulate armoring with Garcia and Parker's formulation. For armoring simulation, the current surface layer of the sediment bed is restricted to a thickness equal to the dimensional reference height.

7.

Cohesive Sediment Settling, Deposition and Resuspension

The settling of cohesive inorganic sediment and organic particulate material is an extremely complex process. Inherent in the process of gravitational settling is the process of flocculation, where individual cohesive sediment particles and particulate organic particles aggregate to form larger groupings or flocs having settling characteristics significantly different from those of the component particles (Burban et al., 1989,1990; Gibbs, 1985; Mehta et al., 1989). Floc formation is dependent upon the type and concentration of the suspended material, the ionic characteristics of the environment, and the fluid shear and turbulence intensity of the flow environment. Progress has been made in first principles mathematical modeling of floc formation or aggregation, and

disaggregation by intense flow shear (Lick and Lick, 1988; Tsai, et al., 1987). However, the computational intensity of such approaches precludes direct simulation of flocculation in operational cohesive sediment transport models for the immediate future. An alternative approach, which has met with reasonable success, is the parameterization of the settling velocity of flocs in terms of cohesive and organic material fundamental particle size, d; concentration, S; and flow characteristics such as vertical shear of the horizontal velocity, du/dz, shear stress, Avdu/dz, or turbulence intensity in the water column or near the sediment bed, q. This has allowed semi-empirical expressions having the functional form
Wse Wse d , S , du ,q dz

(7.1)

to be developed to represent the effective settling velocity. A widely used empirical expression, first incorporated into a numerical by Ariathurai and Krone (1976), relates the effective settling velocity to the sediment concentration:
ws wso S So
a

(7.2)

with the o superscript denoting reference values. Depending upon the reference concentration and the value of , this equation predicts either increasing or decreasing settling velocity as the sediment concentration increases. Equation (7.2) with user defined base settling velocity, concentration and exponent is an option in the EFDC model. Hwang and Metha (1989) proposed
ws aS n S 2 b2
m

(7.3)

based on observations of settling at six sites in Lake Okeechobee. This equation has a general parabolic shape with the settling velocity decreasing with decreasing concentration at low concentrations and decreasing with increasing concentration at high concentration. A least squares for the paramters, a, m, and n, in (7.3) was shown to agree well with observational data. Equation (7.3) does not hav a dependence on flow characteristics, but is based on data from an energetic field condition having both currents and high frequency surface waves. A generalized form of (7.3) can be selected as an option in the EFDC model. Ziegler and Nisbet, (1994, 1995) proposed a formulation to express the effective settling as a function of the floc diameter, df
ws ad b f

(7.4)

with the floc diameter given by:

1/ 2

df

(7.5)
2 xz

2 xz

where S is the sediment concentration,

is an experimentally determined constant and

xz and yz are the x and y components of the turbulent shear stresses at a given position in the water column. Other quantities in (7.4) have been experimentally determined to fit the relationships:
0.85

B1 S

2 xz

2 xz

(7.6)

0.8 0.5 log S

2 xz

2 xz

B2

(7.7)

where B1 and B2 are experimental constants. This formulation is also an option in the EFDC model. A final settling option in EFDC is based on that proposed by Shrestha and Orlob (1996). The formulation in EFDC has the form

ws

S exp

4.21 0.147G

(7.8)

0.11 0.039G
where

2 z

(7.9)

is the magnitude of the vertical shear of the horizontal velocity. It is noted that all of these formulations are based on specific dimensional units for input parameters and predicted settling velocities and that appropriate unit conversion are made internally in their implementation in the EFDC model. Water column-sediment bed exchange of cohesive sediments and organic solids is controlled by the near bed flow environment and the geomechanics of the deposited bed. Net deposition to the bed occurs as the flow-induced bed surface stress decreases. The most widely used expression for the depositional flux is:

(7.10)
ws S d J od 0 :
b cd cd cd b

wsTd S d

cd

where b is the stress exerted by the flow on the bed, cd is a critical stress for deposition which depends on sediment material and floc physiochemical properties (Mehta et al., 1989) and Sd is the near bed depositing sediment concentration. The critical deposition stress is generally determined from laboratory or in situ field observations and values ranging form 0.06 to 1.1 N/m**2 have been reported in the literature. Given this wide range of reported values, in the absence of site specific data the depositional stress and is generally treated as a calibration parameter. The depositional stress is an input parameter in the EFDC model. Since the near bed depositing sediment concentration in (7.10) is not directly calculated, the procedures of Chapter 5 can be applied to relate the the near bed depositional concentration to the bottom layer or depth averge concentration. Using (6.14) the near bed concentration during times of deposition can be determined in terms of the bottom layer concentration for three-dimensional model applications. Inserting (7.10) into (6.14) and evaluating the constant at a near bed depositional level gives

Td

1 Td

R zd Sd zR

(7.11)

Integrating (7.11) over the bottom layer gives

ln zd 1 zd 1 1 zeq1 Sd Td 1 R
1 R 1

Sd

Td

1 Td 1 1 Td

S :R 1
1

(7.12)

zd 1 1

S :R 1

The corresponding quantities in the numerical solution bottom boundary condition (3.6) are

Pd ws

Td

ln

zd 1 zd 1 1

1 Td 1 1 Td

ws : R 1
1

(7.13)
ws : R 1

z Pd ws Td 1 R

1 1 R eq

zd 1 1

For depth averaged model application, (7.10) is combined with (6.25) and the constant of integration is evaluated at a near bed depositional level to give
S 1 Rz 1 R Td Sd 1 1 Rzd 1 R Td Sd
R zd zR

(7.14)

Integrating (7.14) over the depth gives

Sd

2 R 1 zd 21 R 2 R 1 zd 21 R Td

Td

ln zd 1 zd 1 1
R zd 1

1 1

1 R 1 zd 1 R 1 1 Rzd 1 R

Td

S :R 1
(7.15)
1

Sd

1 R zd 1 1

Td

S :R 1

The corresponding quantities in the numerical solution bottom boundary condition (3.6) are

Pd ws

2 R 1 zd 21 R 2 R 1 zd 21 R Td

Td

ln zd 1 zd 1 1
R zd 1

1 1

1 R 1 zd 1 R 1 1 Rzd 1 R

Td

ws : R 1
1

(7.16)

Pd ws

1 R zd 1 1

Td

ws : R 1

It is noted that the assumptions used to arrive at the relationships, (7.12) and (7.15) are more teneous for cohesive sediment than the similar relationships for noncohesive sediment. The settling velocity for cohesive sediment is highly concentration dependent and the use of a constant settling velocity to arrive at (7.12) and (7.15) is questionable. The specification of an appropriate reference level for cohesive sediment is difficult. One possibility is to relate the reference level to the floc diameter using (7.5). An alternative is to set the reference level to a laminar sublayer thickness

zd

S Hu*

(7.17)

where (S) is a sediment concentration dependent kinematic viscosity and the water depth is include to nondimensionlize the reference level. A number of investigators, including Mehta and Jiang (1990) have presented experimental results indicating that at high sediment concentrations, cohesive sediment-water mixtures behave as high viscosity fluids. Mehta and Jain's results indicate that a sediment concentration of 10,000 mg/L results in a viscosity ten time that of pure water and that the viscosity increases logrithmically with increasing mixture density. Use of the relationships (7.12) and (7.16) is optional in the EFDC model. When they are used, the reference height is set using (7.17) with the viscosity determined using Mehta and Jain's experimental relationship between viscosity and sediment concentration. To more fully address the deposition prediction problem, a nested sediment, current and wave boundary layer model based on the near bed closure presented in Chapter 4 is under development. Cohesive bed erosion occurs in two distinct modes, mass erosion and surface erosion. Mass erosion occurs rapidly when the bed stress exerted by the flow exceeds the depth varying shear strength, s, of the bed at a depth, Hme, below the bed surface. Surface erosion occurs gradually when the flow-exerted bed stress is less than the bed shear strength near the surface but greater than a critical erosion or resuspension stress, ce, which is dependent on the shear strength and density of the bed. A typical scenario under conditions of accelerating flow and increasing bed stress would involve first the occurrence of gradual surface erosion, followed by a rapid interval of mass erosion, followed by another interval of surface erosion. Alternately, if the bed is well consolidated with a sufficiently high shear strength profile, only gradual surface erosion would occur. Transport into the water column by mass or bulk erosion can be expressed in the form
r Jo

wr Sr

mme

(7.18)

Tme

where Jo is the erosion flux, the product wrSr represents the numerical boundary condition (3.6), mme is the dry sediment mass per unit area of the bed having a shear strength, s, less than the flow-induced bed stress, b, and Tme is a somewhat arbitrary time scale for the bulk mass transfer. The time scale can be taken as the numerical model integration time step (Shrestha and Orlob, 1996). Observations by Hwang and Mehta (1989) have indicated that the maximum rate of mass erosion is on the order of 0.6 gm/sm**2 which provides an means of estimating the transfer time scale in (4.10). The shear strenght of the cohesive sediment bed is generally agreed to be a linear function of the bed bulk density (Metha et al., 1982; Villaret and Paulic, 1986; Hwang and Mehta, 1989)
s

as

bs

(7.19)

For the shear strength in N/m**2 and the bulk density in gm/cm**3, Hwang and Mehta (1989) give as and bs values of 9.808 and -9.934 for bulk density greater than 1.065 gm/cm**3. The EFDC model currently implements Hwang and Mehta's relationship, but can be readily modified to incorporated other functional relationships. Surface erosion is generally represented by relationships of the form
J or wr S r dme dt
b ce ce

ce

(7.20)

or

J or

wr Sr

dme exp dt

b ce

ce

ce

(7.21)

where dme/dt is the surface erosion rate per unit surface area of the bed and ce is the critical stress for surface erosion or resuspension. The critical erosion rate and stress and the parameters , , and are generally determined from laboratory or in situ field experimental observations. Equation (7.20) is more appropriate for consolidated beds, while (7.21) is appropriate for soft partially consolidated beds. The base erosion rate and the critical stress for erosion depend upon the type of sediment, the bed water content, total salt content, ionic species in the water, pH and temperature (Mehta, et al., 1989) and can be measured in laboratory and sea bed flumes. The critical erosion stress is related to but generally less than the shear strength of the bed, which in turn depends upon the sediment type and the state of consolidation of the bed. Experimentally determined relationships between the critical surface erosion stress and the dry density of the bed of the form
ce

d s

(7.22)

have been presented (Mehta, et al., 1989). Hwang and Mehta (1989) proposed the relationship
ce

b b l

(7.23)

between the critical surface erosion stress and the bed bulk density with a, b, c, and l equal to 0.883, 0.2, 0.05, and 1.065, respectively for the stress in N/m**2 and the bulk density in gm/cm**3. Considering the relationship between dry and bulk density
b d s s w w

(7.24)

equations (7.22) and (7.23) are consistent. The EFDC model allow for a user defined constant critial stress for surface erosion or the use of (7.23). Alternate predictive expression can be readily incorporated into the model. Surface erosion rates ranging from 0.005 to 0.1 gm/s-m**2 have been reported in the literature, and it is generally accepted that the surface erosion rate decreases with increasing bulk density. Based on experimental observations, Hwang and Mehta (1989) proposed the relationship

log10

dme dt

0.23exp

0.198 1.0023 b

(7.25)

for the erosion rate in mg/hr-cm**2 and the bulk density in gm/cm**3. The EFDC model allow for a user defined constant surface erosion rate or predicts the rate using (7.25). Alternate predictive expression can be readily incorporated into the model. The use of bulk density functions to predict bed strength and erosion rates in turn requires the prediction of time and depth in bed variations in bulk density which is related to the water and sediment density and the bed void ratio by

1
b

(7.26)
s

Selection of the bulk density dependent formulations in the EFDC model requires implmentation of a bed consolidation simulation to predict the bed void ratio as discussed in the following chapter.

7.

Sediment Bed Geomechanical Processes

This chapter describes the representation of the sediment bed in the EFDC model. To make the information presented self contained, the derivation of mass balance equations and comparison with formulations used in other models is also presented. Consider a sediment bed represented by discrete layers of thickness Bk, which may be time varying. The conservation of sediment and water mass per unit horizontal area in layer k is given by:
s t

Bk
k

(7.1)

J s:k

J s:k

k , kb J sb
(7.2)
k max J sb , 0 b min J sb , 0

w k t

Bk
k

J w:k

J w:k

k , kb

w s

where is the void ratio, s and w are the sediment and water density and Js and Jw are the sediment and water mass fluxes with k- and k+ defining the bottom and top boundaries, respectively of layer k. The mass fluxes are defined as positive in the vertical direction and exclude fluxes associated with sediment depostion and erosion. The last term in equation (7.1) represents erosion and deposition of sediment at the top of the upper most bed layer, k=kb, where

k , kb
Consitent with this partitioning of flux,
J s:k

1: k kb 0 : k kb

(7.3)

0:k

kb

(7.4)

The last term in (7.2) represents the corresponding entrainment of bed water into the water column during sediment erosion and entrainment of water column water into the bed during deposition. The water flux, Jw:k+, at the top of the upper most layer, kb, is not necessarily zero, since it can include ambient seepage and pore water explusion due to bed consolidation. Assuming sediment and water to be incompressible, (7.1) and (7.2) can be written as:

Bk
t

1
k s

J s:k

J s:k

k , kb

J sb
s

(7.5)

k t

Bk
k

qw:k

qw:k

k , kb

k max

J sb
s

,0

b min

J sb
s

(7.6)
,0

where the water specfic discharges

J w:k J w:k

w w:k w w:k

(7.7)

have been introduced into (7.6). Four approaches for the solution of the mass conservation equations (7.5) and (7.6) have been previously utilized. The solution approaches, hereafter referred to as solution levels, increase in complexity and physical realism and will be briefly summarized. The first level or simplest approach assumes specified time-constant layer thicknesses and void ratios with the left sides of (7.5) and (7.6) being identically zero. Sediment mass flux at all layer interfaces are then identical to the net flux from the bed to the water column.

J s:k J s:k

J sb : k 1, kb 0:k J sb : k kb kb

(7.8)

Bed representations at this level, as exemplified by the RECOVERY model (Boyer, et al., 1994), typically omit the water mass conservation equations. However, it is noted that the water mass conservation is ill posed unless either q1-, the specific discharge at the bottom of the deepest layer or qkb+, the specific discharge at the top of the water column adjacent layer, is specified. If q1- is set to zero, qka+ is then required to exactly cancel the entrainment terms is (7.6). The second level of bed mass conservation representation assumes specified time invariant layer thicknesses. The mass conservation equations (7.5) and (7.6) become

Bk

1
t

1
k s

J s:k

J s:k

k , kb

J sb
s

(7.9)

Bk

k t

qw:k
k

qw:k

k , kb

k max

J sb
s

,0

b min

J sb
s

(7.10)
,0

This system of 2 x kb equations includes kb unknow void ratios, kb unknow internal sediment fluxes, and kb+1 unknow specific discharges and is under determined unless additional information is specified. The constant bed layer thickness option in the WASP5 model (Ambrose, et al., 1993) uses specifed burial velocities to define the internal sediment fluxes
J s:k J s:k wb:k wb:k S k wb:k S k wb:k
s
1

(7.11)
1

Sk

(7.12)
k

where wb is the burial velocity and S is the sediment concentration (mass per unit total volume). Use of the burial velocity eliminates the indetermincy in (7.9) and allowing its solution for the void ratio. In the event that the sediment concentration in the upper most layer becomes negative, the layer is eliminated and the underlying layer become water column adjacent. The left side of the water mass conservation equations (7.10) is now know and the equation is more appropriately written as

qw:k

qw:k

Bk

k t

k , kb
k

k max

J sb
s

,0

b min

J sb
s

(7.13)
,0

The determination of the specific discharges using (7.13) can be viewed is either under determined or physically inconsistent. As shown for the first level approach, the solution of (7.13) is ill posed unless either q1-, the specific discharge at the bottom of the deepest layer or qkb+, the specific discharge at the top of the upper most layer is independently specified. If q1- is specified and the internal specific discharges are determined from (7.13), qka+ is then required to partially cancel the entrainment terms in (7.13). As will be subsequently shown, the specific discharges can be dynamically determined using Darcy's law. However, the specific discharges determined using Darcy's law and the known void ratios are not guaranteed to satisfy (7.13) the level two formulation is dynamically inconsistent with respect to water mass conservation in the sediment bed. The constant bed layer thickness option in the WASP5 ignores this problem entirely by not considering the water mass balance and hence neglecting pore water advection of dissolved contaminants. The third level of bed mass conservation representation assumes specified time invariant layer void ratios. The mass conservation equations (7.5) and (7.6) become

1 1
k

t Bk

1
s

J s:k

J s:k

k , kb

J sb
s

(7.14)

t Bk k

qw:k

qw:k

k , kb

k max

J sb
s

,0

b min

J sb
s

(7.15)
,0

This system of equations exhibits the same under determined nature as (7.9) and (7.10). Specification of internal sediment fluxes or burial velocities allows (7.14) to be solved for the layer thicknesses. Solution of (7.15) for the specific discharges then requires the specification either q1-, the specific discharge at the bottom of the deepest layer or qkb+, the specific discharge at the top of the upper most layer. The variable bed layer thickness option in the WASP5 model (Ambrose, et al., 1993) exemplifies the third level of bed representation. Specifically, the thickness of the water column adjacent layer is allowed to vary in time, while the thicknesses of the underlying layers remain constant. A periodic time variation is specified for the bottom sediment flux in the upper most layer

J s:kb
to N t

0 : to J sb dt : to

to N 1

N 1 t t

t
(7.16)

J s:kb
to

to

N t

where t is the standard water time step and N t is the sediment compaction time. This results in the thickness of the upper most layer periodically returning to its initial value at

time intervals of N t unless the thickness becomes negative due to net resuspension. In that event, the underlying layer becomes the water column adjacent layer. The water mass conservation (7.15) for all but the upper most layer becomes
qk qk q1 : k kb

(7.17)

indicating that all internal specific discharges are equal a specified specific discharge at the bottom of layer 1. Given the solution for the time variation of the water column adjacent thickness and bottom specific discharge, (7.15) can be solved for the specific discharge at the top of the layer. The constant porosity bed option in EFDC is also a level three approach. In EFDC, the internal sediment fluxes are set to zero and the change in thickness of the water column adjacent layer is determined directly using (7.14) while the underlying layers have time invariant thicknesses. As a result, the internal water specific discharges are set to zero and the water entrainment and expulsion in the water column adjacent layer are determined directly from (7.15). The EFDC model is configured to have a user specified maximum number of sediment bed layer. A the start of a simulation, the number of layers containing sediment at a specific horizontal location is specified. Under continued deposition, a new water column layer is created when the thickness of the current layer exceeds a user specified value. If the current water column adjacent layer's index is equal to the maximam number of layers, the bottom two layers are combined and the remaining layers renumbered before addition of the new layer. Under continued resuspension, the layer underlying the current water column adjacent layer becomes the new adjacent layer when all sediment is resuspended form the current layer. The fourth level of bed representation accounts for bed consolidation by allowing the layer void ratios and thicknesses to vary in time. The simplest and most elegant formulations at this level utilize a Lagrangian approach for sediment mass conservation. The Lagrangian approach requires that the sediment mass per unit horizontal area in all layers, except the upper most, be time invariant and without loss of generality, the internal sediment fluxes can be set to zero. Consistent with these requirements (7.5) becomes

Bk
t

k , kb
k

J sb
s

(7.18)

Expanding the left side of the water conservation equation (7.6), and using (7.18) gives

Bk 1
t k k

qw:k

qw:k

k , kb

min

J sb
s

(7.19)

,0

The Lagrangian approach for sediment mass conservation also requires that the number of bed layers vary in time. Under conditions of continued deposition, a new water column adjacent layer would be added when either the thickness, void ratio or mass per unit area of the current water column adjacent layer reaches a predefined value. Under

conditions of continued resuspension, the bed layer immediately under the current water column adjacent layer would become the new water column adjacent layer when the entire sediment mass of the current layer has been resuspended. At the fourth and most realistic level of bed representation, three approaches can be used to represent bed consolidation. Two of the approaches are semi-empirical with the first assuming that the void ratio of a layer decreases with time. A typical relationship which is used for the simple consolidation option in the EFDC model is
m o m

exp

t to

(7.20)

where o is the void ratio at the mean time of deposition, to, m is the ultimate minimum void ratio corresponding to complete consolidation, and is an empirical or experimental constant. Use of (7.20) in the EFDC model involves specifying the depositional void ratio, the ultimate void ratios and the rate constants. The ultimate void ratio can be specified as a function depth below the water column-bed interface. The actual calculation involves using the initial void ratios to determine the deposition time to, after which (7.20) is used to update the void ratios as the simulation progresses. After equation (7.20) is used to calculate the new time level void ratios, equation (7.18) provides the new layer thicknesses. The water conservation equations (7.19) can then be solved using

qw:k

qw:k

Bk 1
t k k

k , kb

min

J sb
s

(7.21)

,0

to determine the water specific discharges, provided that the specific discharge q1-, at the bottom of layer 1 is specified. When this option is specified in the EFDC model, the specific discharge at bottom of the bottom sediment layer is set to zero. Layers are added and deleted in the manner previously described for EFDC's constant porosity option. The SED2D-WES model (Letter et al., 1998) utilizes a similar approach based on a specified time variation of bulk density
s b w bm bo bm

exp

t to

(7.22)

which in turn defines the variation in void ratio. The second semi-empirical approach assumes that the vertical distribution of the bed bulk density or equivalently the, void ratio at any time is given by a self-similar function of vertical position, bed thickness and fixed surface and bottom bulk densities or void ratios. Functionally this equivalent to
V z , BT ,
kb

(7.23)

where V represents the function, z is a vertical coordinate measured upward from the bottom of the lowest layer, and BT is the total thickness of the bed. This approach is used in the original HSTM model (Hayter and Mehta, 1983), the new HSCTM model (Hayter et al., 1998) and is an option in the CE-QUAL-ICM/TOXI model (Dortch, et al., 1998). The determination of the new time level layer thicknesses and void ratios requires an iterative solution of equations (7.18) and (7.23). The solution is completed using (7.21) to determine the water specific discharges. The third and most realistic approach is to dynamically simulate the consolidation of the bed. In the Lagrangian formulation, (7.18) is directly solved for the equivalent sediment thickness
Bk
k

(7.24)
k

and the water conservation equation (7.19) is integrated to determine the void ratio.

t k

qw:k

qw:k

k , kb

min

J sb
s

(7.25)

,0

The specific discharges in (7.25) are determined using the Darcy equation
q K g w u

(7.26)

where K is the hydraulic conductivity and u is the excess pore pressure defined as the difference between the total pore pressure ut, and the hydrostatic pressure uh.
u ut uh

(7.27) and

The total pore pressure is defined as the difference between the total stress effective stress e.
ut
e

(7.28)

The total stress and hydrostatic pressure are given by

zb

pb

g
z

1 1
pb g
w

(7.29)
s

1
zb z

dz

uh

(7.30)

where pb is the water column pressure at the bed zb. Solving for the excess pore pressure using (7.27) through (7.30) gives
zb

s w w

1
z

1 1

(7.31)

dz

The specific discharge (7.26), can alternately be expressed in terms of the effective stress
q K g w
zb s z e w

1 K

1 1

(7.32)
dz

z z

or the void ratio

q K d g w d
zb e z w s

1 K

1 1

(7.33)
dz

z z

where d /d

is a coefficient of compressibility.

For consistency with the Lagrangian representation of sediment mass conservation, a new vertical coordinate , defined by

d dz

1 1

(7.34)

is introduced. The discrete form of (7.34) is

zk 1 zk
k

Bk 1
k k

(7.35)

where D is the equivalent sediment thickness previously defined by (7.24). Introducing (7.34) into (7.26), (7.32), and (7.33) gives
q K g w 1
s e w

(7.36)

K g w 1 q K 1

K 1 K 1

(7.37)

s w

(7.38)

where
1 d g w d
e

(7.39)

is a compressibility length. Three formulations for the solution the consolidation problem can be utilized. The void ratio-excess pore pressure formulation, used in the EFDC model, evaluates the specific discharges at the current time level n, using (7.36) and explicitly integrates (7.25)
n n 1 k n k n k

qw:k

qw:k

k , kb

min

J sb
s

(7.40)

,0

where is the time step, to give the new time level void ratios. The layer thicknesses are then determined by explicit integration of (7.18).
B 1
n 1

B 1
n k

k , kb
k

J sb
s

(7.41)

k n 1 k

k , kb

J sb
s

Constitutive equations required for consolidation prediction generally express the effective stress and hydraulic conductivity as functions of the void ratio. Thus the new time level void ratio is used to determine new time level values of the effective stress and hydraulic conductivity. The new time level excess pore pressures is then given by

s w w

(7.42)

the transformed equivalent of (7.31). The primary advantage of the void ratio-excess pore pressure formulation is the simplicity of its boundary conditions
u ub :
b

(7.43)

u q

uo : or qo :

(7.44)
0

The water column-sediment bed interface boundary condition generally sets ub to zero if the surface water flow is hydrostatic but can incorporate wave induced pore pressures. The bottom boundary conditions allows either the specification of pressure or specific

discharge. The primary disadvantage of this formulation is the stability or positivity criterion imposed on the time step
n n k k n

(7.40)
k

qw:k

qw:k

k , kb
n k

min

J sb
s

,0

k , kb max

J sb
s

(7.41)

,0

In practice, these criteria are readily satisfied if the consolidation time step is identical to the time step of the hydrodynamic model. In the event that these criteria are not met using the hydrodynamic time step, the bed consolidation is sub-cycled using an integer number of time steps, meeting (7.40) and (7.41), per each hydrodynamic time step. Alternately, the consolidation problem can be directly formulated in terms of the effective stress or void ratio. Combining (7.25) and (7.37) using (7.39) gives the effective stress formulation

K
k t k

s w

1
s w

K 1
k b k

(7.42)

K 1

K 1
k

k , kb

min

J sb
s

,0

The continuum equivalent is

1
t e:k

K 1 g
w b k e

g min

K
s w

(7.43)

J sb
s

,0

which is parabolic since is negative. Combining (7.25) and (7.38) using (7.39) gives the void ration formulation

K
k t k

s w

1
s w

K 1
k b k

(7.44)

K 1

K 1
k

k , kb

min

J sb
s

,0

The continuum equivalent is

K
t k s w

K 1
b k b

min

J sb
s

(7.45)
,0

Equation (7.45) is the discrete form of the finite strain consolidation equation first derived by Gibson et al. (1967). Equation (7.45) was used by Cargill (1985) in the formulation of a model for dredge material consolidation and by Le Normant (1998) to represent bed consolidation in a three-dimensional cohesive sediment transport model. The classic linear consolidation equation (Middleton and Wilcock, 1994) omits the second term associated with self weight in (7.45) and introduces a constant consolidation coefficient
Cc 1 K e g w
e

(7.46)

reducing (7.45) to
t

Cc

zz

(7.47)

Equation (7.47) has separable solutions of the form

exp
n n n

Cc t B2 0

(7.48)

z B

which provides some justification for empirical relationship (7.20). The solution of the finite strain consolidation problem in any of its three forms requires constitutive relationships
e e

(7.49) (7.50)

Bear (1979) notes that curve fitting of experimental data typically results in relationships of the form

av

eo

(7.51) (7.52)

Cc ln

e eo

for noncohesive and coheasive soils respectively, where av is the coefficient of compressibility and Cc is the compression index. Graphical presentation of experimental forms of (7.49) and (7.50) are presented in Cargill (1985) and Palermo et al., (1998) which are generally consistent with (7.52) and suggest

ln

K Ko

(7.53)

as a candidate relationship between the void ratio and hydraulic conductivity for cohesive sediment beds. Similarly, a linear relationship
o

Ko

(7.54)

would likely suffice for noncohesive sediment beds.

8.

References

Ackers, P., and W. R. White, 1973: Sediment transport: New approaches and analysis. J. Hyd. Div. ASCE, 99, 2041-2060. Ariathurai, R., and R. B. Krone, 1976: Finite element model for cohesive sediment transport. J. Hyd. Div. ASCE, 102, 323-338. Ambrose, R. B., T. A. Wool, and J. L. Martin, 1993: The water quality analysis and simulation program, WASP5: Part A, model documentation version 5.1. U. S. EPA, Athens Environmental Research Laboratory, 210 pp. Bear, J., 1879: Hydraulics of groundwater, McGraw-Hill, New York. Bagnold, R. A., 1956: The flow of cohesionless grains in fluids. Phil. Trans. Roy. Soc. Lond., Series A, Vol 249, No. 964, 235-297. Blumberg, A. F., B. Galperin, and D. J. O'Connor, 1992: Modeling vertical structure of open-channel flow. J. Hydr. Engr., 118, 1119-1134. Boyer, J. M., S. C. Chapra, C. E. Ruiz, and M. S. Dortch, 1994: RECOVERY, a mathematical model to predict the temporal response of surface water to contaminated

sediment. Tech. Rpt. W-94-4, U. S. Army Engineer Waterways Experiment Station, Vicksburg, MS, 61 pp. Burban, P. Y., W. Lick, and J. Lick, 1989: The flocculation of fine-grained sediments in estuarine waters. J. Geophys. Res., 94, 8323-8330. Burban, P. Y., Y. J. Xu, J. McNeil, and W. Lick, 1990: Settling speeds of flocs in fresh and seawater. J. Geophys. Res., 95, 18,213-18,220. Cargill, K. W., 1985: Mathematical model of the consolidation and desiccation processes in dredge material. U.S. Army Corps of Engineers, Waterways Experiment Station, Technical Report D-85-4. Dortch, M., C. Ruiz, T. Gerald, and R. Hall, 1998: Three-dimensional contaminant transport/fate model. Estuarine and Coastal Modeling, Proceedings of the 5nd International Conference, M. L. Spaulding and A. F. Blumberg, Eds., American Society of Civil Engineers, New York, 75-89. Einstein, H. A., 1950: The bed load function for sediment transport in open channel flows. U.S. Dept. Agric. Tech. Bull., 1026. Galperin, B., L. H. Kantha, S. Hassid, and A. Rosati, 1988: A quasi-equilibrium turbulent energy model for geophysical flows. J. Atmos. Sci., 45, 55-62. Garcia, M., and G. Parker, 1991: Entrainment of bed sediment into suspension. J. Hyd. Engrg., 117, 414-435. Gibbs, R. J., 1985: Estuarine Flocs: their size, settling velocity and density. J. Geophys. Res., 90, 3249-3251. Gibson, R. E., G. L. England, and M. J. L. Hussey, 1967: The theory of one-dimensional consolidation of saturated clays. Geotechnique, 17, 261-273. Hamrick, J. M., 1992: A three-dimensional environmental fluid dynamics computer code: Theoretical and computational aspects. The College of William and Mary, Virginia Institute of Marine Science, Special Report 317, 63 pp. Hamrick, J. M., and T. S. Wu, 1997: Computational design and optimization of the EFDC/HEM3D surface water hydrodynamic and eutrophication models. Next Generation Environmental Models and Computational Methods. G. Delich and M. F. Wheeler, Eds., Society of Industrial and Applied Mathematics, Philadelphia, 143-156. Hayter, E. J., and A. J. Mehta, 1983: Modeling fine sediment transport in estuaries. Report EPA-600/3-83-045, U.S. Environmental Protection Agency. Athens, GA>

Hayter, E.J., M. Bergs, R. Gu, S. McCutcheon, S. J. Smith, and H. J. Whiteley, 1998: HSCTM-2D, a finite element model for depth-averaged hydrodynamics, sediment and contaminant transport. Technical Report, U. S. EPA Environmental Research Laboratory, Athens, GA. Hwang, K.-N, and A. J. Mehta, 1989: Fine sediment erodibility in Lake Okeechobee. Coastal and Oeanographic Enginnering Dept., University of Florida, Report UFL/COEL89/019, Gainsville, FL. Laursen, E., 1958: The total sediment load of streams J. Hyd. Div. ASCE, 84, 1-36. Letter, J. V., L. C. Roig, B. P. Donnell, Wa. A. Thomas, W. H. McAnally, and S. A. Adamec, 1998: A user's manual for SED2D-WES, a generalized computer program for two-dimensional, vertically averaged sediment transport. Version 4.3 Beta Draft Instructional Report, U. S. Army Corps of Engrs., Wtrwy. Exper. Sta., Vicksburg, MS. Le Normant, C., E. Peltier, and C. Teisson, 1998: Three dimensional modelling of cohesive sediment in estuaries. in Physics of Estuaries and Coastal Seas, (J. Dronkers and M. Scheffers, Eds.), Balkema, Rotterdam, pp 65-71. Lick, W., and J. Lick, 1988: Aggregation and disaggregation of fine-grained lake sediments. J Great Lakes Res., 14, 514-523. Mehta, A. J., E. J. Hayter, W. R. Parker, R. B. Krone, A. M. Teeter, 1989: Cohesive sediment transport. I: Process description. J. Hyd. Engrg., 115, 1076-1093. Mehta, A. J., T. M. Parchure, J. G. Dixit, and R. Ariathurai, 1982: Resuspension potential of deposited cohesive sediment beds, in Estuarine Comparisons, V. S. Kennedy, Ed., Academic Press, New York, 348-362. Mehta, A. J., and F. Jiang, 1990: Some field observations on bottom mud motion due to waves. Coastal and Oeanographic Enginnering Dept., University of Florida, Gainsville, FL. Mellor, G. L., and T. Yamada, 1982: Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. Space Phys., 20, 851-875. Meyer-Peter, E. and R. Muller, 1948: Formulas for bed-load transport. Proc. Int. Assoc. Hydr. Struct. Res., Report of Second Meeting, Stockholm, 39-64. Middleton, G. V., and P. R. Wilcock, 1994: Mechanics in the Earth and Environmental Sciences. Cambridge University Press, Cambridge, UK. Nielsen, P., 1992: Coastal bottom boundary layers and sediment transport, World Scientific, Singapore.

Park, K., A. Y. Kuo, J. Shen, and J. M. Hamrick, 1995: A three-dimensional hydrodynamic-eutrophication model (HEM3D): description of water quality and sediment processes submodels. The College of William and Mary, Virginia Institute of Marine Science. Special Report 327, 113 pp. Rahmeyer, W. J., 1999: Lecture notes for CEE5560/6560: Sedimentation Engineering, Dept. of Civil and Environmental Engineering, Utah State University, Logan, Utah. Raukivi, A. J., 1990: Loose boundary hydraulics. 3rd Ed. Pergamon, New York, NY. Ried, I., and L. E. Frostick, 1994: Fluvial sediment transport and deposition. in Sediment Transport and Depositional Processes, K. Pye, ed., Blackwell, Oxford, UK, 89-155. Shrestha, P. A., and G. T. Orlob, 1996: Multiphase distribution of cohesive sediments and heavy metals in estuarine systems. J. Environ. Engrg., 122, 730-740. Smagorinsky, J., 1963: General circulation experiments with the primative equations, Part I: the basic experiment. Mon. Wea. Rev., 91, 99-152. Smith, J. D., and S. R. McLean , 1977: Spatially averaged flow over a wavy bed. J. Geophys. Res., 82, 1735-1746. Smolarkiewicz, P. K., and T. L. Clark, 1986: The multidimensional positive definite advection transport algorithm: further development and applications. J. Comp. Phys., 67, 396-438. Smolarkiewicz, P. K., and W. W. Grabowski, 1990: The multidimensional positive definite advection transport algorithm: nonoscillatory option. J. Comp. Phys., 86, 355375. Spasojevic, M., and F. M. Holly, 1994: Three-dimensional numerical simulation of mobile-bed hydrodynamics. Contract Report HL-94-2, US Army Engineer Waterways Experiment Station, Vicksburg, MS. Stark, T. D., 1996: Program documentation and users guide: PSDDF primary consolidation, secondary compression, and desiccation of dredge fill. Instructional Report EL-96-xx, US Army Engineer Waterways Experiment Station, Vicksburg, MS. Suzuki, K., H. Yamamoto, and A. Kadota, 1998: Mechanism of bed load fluctuations of sand-gravel mixture in a steep slope channel, Proc. of the 11th congress of APD IAHR, Yogyakarta, pp.679-688. Tsai, C. H., S. Iacobellis, and W. Lick, 1987: Floccualtion fo fine-grained lake sediments due to a uniform shear stress. J Great Lakes Res., 13, 135-146.

van Niekerk, A., K. R. Vogel, R. L Slingerland, and J. S. Bridge, 1992: Routing of heterogeneous sedimetns over movable bed: Model development. J. Hyd. Engrg., 118, 246-262. Van Rijin, L. C., 1984a: Sediment transport, Part I: Bed load transport. J. Hyd. Engrg., 110, 1431-1455. Van Rijin, L. C., 1984b: Sediment transport, Part II: Suspended load transport. J. Hyd. Engrg., 110, 1613-1641. Villaret, C., and M. Paulic, 1986: Experiments on the erosion of deposited and placed cohesive sediments in an annular flume and a rocking flume. Coastal and Oeanographic Enginnering Dept., University of Florida, Report UFL/COEL-86/007, Gainsville, FL. Ziegler, C. K., and B. Nesbitt, 1994: Fine-grained sediment transport in Pawtuxet River, Rhode Island. J. Hyd. Engrg., 120, 561-576. Ziegler, C. K., and B. Nesbitt, 1995: Long-term simulation of fine-grained sediment transport in large reservoir. J. Hyd. Engrg., 121, 773-781. Yang, C. T., 1973: Incipient motion and sediment transport. J. Hyd. Div. ASCE, 99, 16791704. Yang, C. T., 1984: Unit stream power equation for gravel. J. Hyd. Engrg., 110, 17831797. Yang, C. T., and A. Molinas, 1982: Sediment transport and unit streams power function. J. Hyd. Div. ASCE, 108, 774-793.

9.

Figures

1e+01 critical shields shear velocity settling velocity 1e+00

1e-01

1e-02

1e-03

1e-04 1e-05 1e-04 1e-03 1e-02 1e-01 1e+00

grain diameter (meters)

Figure 1. Critical Shield's shear velocity and settling velocity as a function of sediment grain size.

8.

Sorptive Contaminant Transport

The transport of a sorptive contaminant in the water column is governed by transport equations for the contaminant dissolved in the water phase, for the contaminant sorbed to material effectively dissolved in the water phase, and for the contaminant sorbed to suspended particles. For the portion of the contaminant dissolved directly in the water phase
t

mx my HCw
z

my HuCw mx my H
i i K aS S i

mx HvCw
i K dS S i i S j

mx my wCw
j K dD D j j D

mx my

Ab H

Cw

(8.1)

Cw
w

i S

i S

mx m y H
j

i j K aD D j

Cw
w

j D

j D

Cw

where Cw is the mass of water dissolved contaminant per unit total volume, S is the mass of contaminant sorbed to sediment class i per mass of sediment, D is the mass of contaminant sorbed to dissolved material j per unit mass of dissolved material, is the porosity, w is the fraction of the water dissolved contaminant available for sorption, Ka is the adsorption rate, Kd is the desorption rate, and is a net linearized decay rate coefficient. The sorption kinetics are based on the Langmuir isotherm (Chapra, 1997) with denoting the saturation sorbed mass per carrier mass. The sediment and dissolved material concentrations, S and D are defined as mass per unit total volume. The transport equation for the portion of material sorbed to a dissolved constituent D is,
t

mx my HD j
z

j D

my HuD j Dj
j D

j D

mx HvD j

j D

mx my wD j
j D j D

j D

mx my

Ab H

j mx my H K sD D j

Cw
w

(8.2)

j mx my H K dD

Dj

j D

The transport equation for the portion of material sorbed to a suspended constituent S is,
t

mx m y HS i

i S

m y HuS i
i S

i S

mx HvS i Ab H
z

i S

mx m y wS i

i S

i mx m y wS S i

mx m y
i S

Si

i S

(8.3)
Si
i S

i mx m y H K aS S i

Cw
w

i S

i mx m y H K dS

Introducing sorbed concentrations defining sorbed mass per unit total volume
j CD

Dj Si

j D

(8.4) (8.5)

i CS

i S

Allows equations (8.1) through (8.3) to be written as

t

mx my HCw
z

my HuCw Cw mx m y H

mx HvCw
i i K dS CS i

mx my wCw
j j K dD CD

mx m y

Ab H

(8.6)

i K aS S i

Cw
w

i S

i S

mx m y H
j

i j K aD D j

Cw
w

j D

j D

Cw

j mx my HCD

j my HuCD

j mx HvCD

j mx my wCD j D j D

mx my

Ab H

j CD

j mx my H K sD D j j mx my H K dD j CD

Cw
w

(8.7)

i mx m y HCS

i m y HuCS

i mx HvCS

z i CS

i mx m y wCS

i i mx m y wS CS

mx m y
i S

Ab H

(8.8)
i CS

i mx m y H K aS S i

Cw
w

i S

i mx m y H K dS

The EFDC sorbed contaminant transport formulation currently employees equilibrium partitioning with the adsorption and desorption terms in (8.7) and (8.8) balancing
j K sD D j

Cw
w

j D

j D

j j K dDCD

(8.9)

i K aS S i

Cw
w

i S

i S

i i K dS CS

(8.10)

Solving (8.9) and (8.10) for the sorbed to water phase concentration ratios gives
j CD Cw j D j Do

f Dj fw 1 P
w

P
j Do

j D

Dj
1

Cw j D

(8.11)

j Do

j j K aD D j K dD

i CS Cw i S i So

f Si fw 1 P
w i So

PSi

Si
1

Cw i S

(8.12)

j PSo

i i K aS S i K dS

where P denotes the partition coefficient, and Po is its linear equilibrium value. For linear equilibrium partitioning, P is set to Po, which in effect approximates ( )-1 terms in (8.11) and (8.12) by unity. Requiring the mass fractions to sum to unity
fw
i

f Si
j

f Dj

(8.13)

gives
fw Cw C
i j CD C i i S

PSi S i
j

PDj D j

f Dj

PDj D j PSi S i
j

PDj D j

(8.14)

f Si

C C

PS PS
i

i S i i S

PDj D j
j

The dissolved concentrations can be alternately expressed by mass per unit volume of the water phase

Cw:w C
j D :w

Cw
j CD

(8.15)

D
with (8.14) becoming
Cw:w C
i

J :w

Dj

1 PS
i S i j

PDj D:J w
J :w

C C C C

j D:w

P D PS
i

j D i i S i S i

(8.16)
PDj D:J w
j

i S

PS PS
i S i

PDj D:J w
j

Which is a generalization of Chapra's (1997) formulation for sorption to dissolved and particulate organic carbon. Adding equations (8.6), (8.7), and (8.8), using the equilibrium partitioning relationships (8.9) and (8.10) gives

mx my HC

1 mx my mx my
i

my HuC
i S

1 mx my
z

mx HvC C

mx my wC
(8.17)

w f C

i S

A mx my b H

mx my H C

the equation for the total concentration, C. The boundary condition at the water columnsediment bed interface, z = 0, is

Ab H
i max J SBS i i S

C
i

i wS f Si C

,0

max

i J SBS i S

Cw ,0
j

j CD

SB i J SBS i S

(8.18)

Cw ,0
j dep

j D

min J
i

i SBS

i S

,0

dep

min

WC i SBB i S

max
i

Cw ,0
j

j CD

SB i J SBB i S

Cw ,0
j dep

j D

dep i

min

WC

Cw max qw , 0
j

j D

Cw min qw , 0
SB j D j D j dep

j CD

WC

Cw qdif
j dep

Cw
j

WC

SB

where JSBS and JSBB are the suspended load and bed load sediment fluxes between the sediment bed and the water column, defined as positive from the bed, s is the sediment density, qw is the water specific discharge due to bed consolidation and groundwater interaction, defined as positive from the bed, and qdif is a diffusion velocity incorporating the effects of molecular diffusion, hydrodynamic dispersion, and biological induced mixing. The subscript SB denotes conditions in the top layer of the sediment bed, while the subscript WC denotes condition in the water column immediately above the bed, with the exception that the specific discharge and diffusion velocity are defined at the water column-bed interface. The subscript, dep, is used to denote the void ratio and porosity of newly depositing sediment. Equation (8.16) indicates that contaminant flux between the bed and water column includes, a flux of suspended sediment sorbed material; fluxes of water dissolved and sorbed to water dissolved material due to the specific discharge of water associated with consolidation and ground water interaction and water entrainment and expulsion associated with both suspended and bed load sediment deposition and resuspension; and a flux of water dissolved and sorbed to water dissolved material due to diffusion like processes. Transport of bed load sediment sorbed material is represented by direct transport between horizontally adjacent top bed layers and is included in the

contaminant mass conservation equations for the sediment bed. The boundary condition at the water free surface is

Ab H

C
i

i wS f Si C

0 : z 1

(8.19)

Using the relationship between the porosity and void ratio (8.20)

1
and (8.5) allows (8.18) to be written as

Ab H
i max J SBS i i CS ,0 Si

C
i

i wS f Si C i J SBS i S i J SBS i S

max

,0

Cw
j

j CD SB j CD j WC

i min J SBS i

i CS ,0 Si

dep

min

,0

Cw
j CD j SB j CD j WC

1
i

max

i J SBB i S

,0

Cw

(8.21)

1
i

dep

min

i SBB i S

,0 1

Cw Cw
j j CD

max qw , 0 min qw , 0

qdif qdif

SB j CD

Cw
j

WC

The sediment concentration can be expressed in terms of the sediment density and void ratio by

Si

Fi 1

i s

(8.22)

where Fi is the fraction of the total sediment volume occupied by each sediment class

i i

Si
i s

Si
i s

(8.23)

Introducing (8.14) and (8.22) into (8.21) gives the final form of the bottom boundary condition

Ab H max
i

C
i

i wS f Si C i F i J SBS ,0 Si

i i SBS S i

, 0 C max

fw
j

f Dj C
SB

min
i

i Fi Ji J SBS f Si , 0 C min depi SBS , 0 Si S dep i J SBB i S i J SBB i S

fw
j j CD j

f Dj C
WC

1
i

max

,0

Cw

(8.24)
SB

1
i

dep

min

,0

Cw
j

j CD WC

max qw , 0 min qw , 0

qdif qdif

fw
j

f Dj C
SB

fw
j

f Dj C
WC

Note that the form of the bed flux associated with bed load transport remains unmodified since a sediment concentration in the water column cannot be readily defined for sediment being transported as bed load. The transport equation (8.17) for the total contaminant concentration in the water column is solved using a fractional step procedure which sequentially treats advection; settling, deposition, and resuspension; pore water advection and diffusion; and reactions. The fractional phase distribution of the contaminant is recalculated between the advection, settling, deposition and resuspension, and pore water advection and diffusion steps using (8.14). The advection step is
HC
n 1/ 4

HC

mx my

my HuC

mx my

mx HvC

wC

(8.25)

with the vertical boundary conditions

wC 0 : z 0, 1

(8.26)

The fractional time level in (8.25) and subsequent equations is used to denote an intermediate result in the fractional step procedure. The spatially discrete from of (8.25) is solved using one of the standard high order, flux limited, advective transport solvers in the EFDC model. The settling, deposition, and resuspension step is

HC

n 1/ 2

HC

n 1/ 4 z i

i wS f SiC

(8.27)

with the boundary conditions

i wS f Si C i i i SBS S i i F i J SBS ,0 Si

max
i

, 0 C max

fw
j

f Dj C
SB

i i Fdep J SBS Ji f i min SBSi S , 0 C min ,0 i S S dep i J SBB i S

(8.28)

fw
j j CD j

j D

C
WC

1
i

max
i J SBB i S

,0

Cw
j CD j

SB

1
i

dep

min

,0

Cw

: z
WC

i wS f SiC

0 : z 1

(8.29)

Integrating (8.27) over a water column layer and using upwind differencing for the settling gives,
n 1/ 2 k k n 1/ 4 k i i wS S i i i wS S i n 1/ 2 k

HC

HC

H
n 1/ 2

f Si Si HC

HC
k 1

n 1/ 2 k 1

(8.30)

f Si Si

n 1/ 2 k

for a layer not adjacent to the bed, and,

HC

n 1/ 2 1 1

HC

n 1/ 4 1 i

i wS S i

f Si Si

n 1/ 2 n C2 2 n 1/ 2 1/ 2

Ji f i max SBSi S , 0 S min

i J SBS f Si ,0 Si

Ji Fi max SBSi , 0 S min

i J SBB i S i J SBB i S i J SBS F i ,0 Si

fw
j

f Dj
sb n 1/ 2

n Csb 1/ 2

fw
j

f Dj
1 n 1/ 2 n Csb 1/ 2 sb n 1/ 2

C1n

1/ 2

(8.31)

1
i

max

,0

fw
j

f Dj f Dj
j 1

1
i

min

,0

fw

C1n

1/ 2

for the first layer adjacent to the bed. Note that (8.31) is also the appropriate form for single layer or depth average application. Since the sediment settling flux is zero at the top of the free surface adjacent layer, (8.27) is integrated downward from the top layer to the bottom layer. The bottom layer equation (8.31) is solved simultaneously with a corresponding equation for the top layer of the sediment bed. The settling fluxes, wSS, and water column-sediment bed fluxes, JSB, in (8.30) and (8.31) are known from the preceding solution for sediment settling, deposition and resuspension. Terms containing the sediment sorbed fraction divided by the sediment concentration in (8.30) and (8.31) are given by

f Si Si
i

PSi PSi S i
j

PDj D j

(8.32)

The diffusion step is given by

HC
n 3/ 4

HC

n 1/ 2 z

Ab H

(8.33)

with boundary conditions

Ab H

max qw , 0 qdif 1
dep

qdif fw

fw
j

f Dj C
SB

(8.34)

min qw , 0

f Dj C
j WC

: z

Ab H

0 : z 1

(8.35)

For the first layer adjacent to the bed

HC
n 3/ 4 1

HC

n 1/ 2 1 1

Ab H f Dj
j j D

n 3/ 4 z

C
1 n 1/ 2 n CSB 3/ 4 SB n 1/ 2

max qw , 0
1

qdif

fw

(8.36)

min qw , 0
1

qdif

fw
j

1
dep 1

C1n

1/ 2

It is noted that the bed concentrations are advanced to the n+3/4 intermediate time level before the advance of the water column concentrations. While for layers not adjacent to the bed,

HC

n 3/ 4 k

HC

n 1/ 2 k 1

Ab H

n 3/ 4 zC k 1

Ab H

n 3/ 4 zC k

(8.37)

The solution is completed by

HC
n 1 k

HC

n 3/ 4 k

HC

n 1 k

(8.38)

an implicit reaction step. Contaminant transport in the sediment bed is represented using the discrete layer formulation developed for bed geomechanical processes. The conservation of mass for the total contaminant concentration in a layer of the sediment bed is given by

BC

BC J
i SBS i

k , kt
i

max

i SBS

f
i

i S

BS

,0

max

F ,0 BS i

fw
j

f Dj
kt

BC

kt

k , kt
i

min

i i i J SBS Fdep J SBS f Si , 0 C min ,0 i Si S dep

fw
j

f Dj C
WC

k , kt
i

i SBB

i SBL

,0 f Dj
j kt

k , kt
i

i J SBB max ,0 i B S

fw ,0 fw

BC f Dj C

kt

(8.39)

k , kt
i

dep

min

i SBB i S

WC

max qw , 0

qdif

min qw , 0 qdif qdif

qdif

1 fw B fw
j

f Dj
j k

BC

k , kt min qw , 0 1 k , kt min qw , 0 qdif

1
kt

f Dj C
WC

1 fw B f Dj
j

f Dj
j k 1

BC

k 1

max qw , 0

1 fw B

BC
k 1

k 1

where

k , kt

0 : k kt 1 : k kt

(8.40)

is used to distinguish processes specific to the top, water column adjacent layer of the bed, kt. Advective fluxes associated with pore water advection in (8.40) are represented in upwind form. In the sediment bed, the actual computational variables for sediment, contaminant, and dissolved material are their concentrations times the thickness of the bed layer. Consistent with this formulation, the fractional phase components in the bed are defined by

fw

BCw BC
j BCD BC

B
i

B P BS i
i S j

PDj BD j
k

f Dj

B
i

PDj BD j PSi BS i
j

(8.41)
PDj BD j
k

f Si

i BCS BC

B
i

PSi BS i PSi BS i
j

PDj BD j
k

The contaminant fluxes associated bed load sediment transport are determined as follows. The net sediment flux from the bed load transport equation
i mx my J SBB x i my QSBLx x i mxQSBLy

(8.42)

is used to evaluate the flux associated with pore water entrainment and expulsion in (8.25) and (8.40). The transport equation for material sorbed to the bed load is
x i my QSBLx i SBL x i mxQSBLy i SBL i mx my J SBB i SBL

(8.43)

Since the contaminant mass per sediment mass in the transport divergence corresponds to conditions in the top layer of the sediment bed, (8.43) can be written as
i my QSBLx

f Si C Si

i mxQSBLy

f Si C Si

i mx my J SBB

i SBL

(8.44)

And solved using an upwind approximation

i mx my J SB i max my QSBLx i SBL i min m y QSBLx C i min m y QSBLx W i min mx QSBLy C i min mx QSBLy S

f Si C Si f C S f Si C Si f C S
i S i i S i

f Si C Si f C S f Si C Si f C S
i S i i S i

i max my QSBLx

(8.45)

i max mx QSBLy

i max mx QSBLy

To evaluate the transport of bed load sorbed material between horizontally adjacent top layers of the sediment bed. Equation (8.39) is solved using a fractional step procedure consistent with that used for the water column transport. Equation (8.41) is used to update the fractional distribution in the bed between the settling, deposition, and resuspension step and the pore water advection and diffusion step. The settling, deposition and resuspension step applies only to the top layer of the bed and is
BC Ji f i max SBS i S , 0 BS
n 1/ 2 kt

BC

n kt n 1/ 2 n 1/ 2 kt

Ji Fi max SBS i , 0 BS

fw
j

j D kt

BC

i i i J SBS Fdep J SBS f Si min , 0 C min ,0 i Si S dep i J SBB i i SBL

n 1/ 2

fw
j

j D

C
WC

(8.46)

,0
n 1/ 2

1
i

Ji max SBB , 0 i B S 1
dep

fw
j

j D kt

BC

n 1/ 2 kt

min

i J SBB i S

n 1/ 2

,0

fw
j

j D

C
WC

This equation is solved simultaneously with equation (8.31) for the bottom layer of the water column. The solution is represented by

BC a11 a12 BC HC
n 1/ 2 kt n 1/ 2 1 1

n kt i i ib

i J SB

n 1/ 2 i SBL

(8.47)
1

a21 a22

HC

n 1/ 4 1 i

i wS S i

f Si Si

n 1/ 2 n C2 2 1/ 2

where the coefficients are given by

a11 1
i

Ji f i max SBS i S , 0 BS 1
i

Ji Fi max SBS i , 0 BS
i J SBB ,0 i B S

n 1/ 2

fw
j n 1/ 2

j D kt

max

fw
j

f Dj
kt n 1/ 2

a12

Ji f i min SBSi S , 0 S 1
i

min
i J SBB i S

i i J SBS Fdep i S dep

,0

fw
j

f
n 1/ 2

j D 1

H a21
i

dep

min

,0

fw
j

f Dj
1 n 1/ 2

Ji f i max SBS i S , 0 BS 1
i

Ji Fi max SBS i , 0 BS
i J SBB ,0 i B S

(8.48)

fw
j n 1/ 2

j D kt

max

fw
j

f Dj
kt n 1/ 2 j D 1

a22

Ji f i min SBSi S , 0 S 1
i

Ji Fi min SBSi , 0 S
i J SBB i S

fw
j n 1/ 2

min

,0

fw
j

f Dj
1

Adding the two equations in (8.47) gives

BC
n 1/ 2 kt 1

HC f Si Si

n 1/ 2 1 n 1/ 2 n C2

BC
1/ 2

n kt 1

HC

n 1/ 4 1 n 1/ 2 i SBL

wS
i

i S

i 1

i J SBS i

(8.49)

This equation verifies the consistency of the water column-sediment bed exchange since the source and sinks on the right side include only settling into the top of the water column layer, and transfer of bed load sediment sorbed contaminant between horizontal sediment bed cells. The pore water advection and diffusion step for the top, water column adjacent, layer is

BC max qw , 0 qdif

n 3/ 4 kt

BC

n 1/ 2 kt n 1/ 2 n 3/ 4 kt

kt

1 fw B 1 fw B 1 H fw

f
j

j D kt n 1/ 2

BC

min qw , 0

qdif

kt

f
j

j D kt n 1/ 2

BC

n 3/ 4 kt

(8.50)

min qw , 0

qdif

kt

f Dj
j 1 n 1/ 2 j D kt 1

HC

n 1/ 2 1

max qw , 0

qdif

kt

1 fw B

f
j

BC

n 1/ 2 kt 1

which is an implicit form. Writing (8.36) in the form

HC

n 3/ 4 1 1

HC

n 1/ 2 1

Ab H f Dj
j

n 3/ 4 z

C
1 n 3/ 4 kt

n 1/ 2

max qw , 0

qdif

kt

1 fw B 1 H fw

BC
SB n 1/ 2

(8.51)

min qw , 0

qdif

kt

f Dj
j 1

HC

n 1/ 2 1

and combining with (4.49) gives

BC

n 3/ 4 kt 1

HC

n 3/ 4 1

BC

n 1/ 2 kt 1

HC f Dj
j

n 1/ 2 1 n 1/ 2

Ab H BC
kt

n 3/ 4 zC 1 n 3/ 4

min qw , 0

qdif

kt

1 fw B 1 fw B

(8.52)

kt n 1/ 2 n 3/ 4 kt 1

max qw , 0

qdif

kt

f Dj
j kt 1

BC

This equation verifies the consistency of the representation of pore water advection and diffusion across water column-sediment bed interface since the source and sink terms on the right side of (8.52) represent fluxes at the top to the water column cell and the bottom of the bed cell. The pore water diffusion and advection step for the remaining bed layers is given by

BC min qw , 0 qdif

n 3/ 4 k

BC

n 1/ 2 k n 1/ 2 n 3/ 4 k 1

1 fw B 1 fw B 1 fw B 1 fw B

f
j

j D k 1 n 1/ 2

BC

max qw , 0

qdif

f
j

j D k n 1/ 2

BC

n 3/ 4 k

(8.53)

min qw , 0

qdif

f Dj
j k n 1/ 2 j D k 1

BC

n 3/ 4 k

max qw , 0

qdif

f
j

BC

n 3/ 4 k 1

For the bottom layer of the bed, k = 1, the bottom, k-, specific discharge and diffusion velocity must be specified as well as the total contaminant concentration, C0. The corresponding thickness of the unresolved layer, k = 0, is set to unity without loss of generality. The system of equations represented by (8.49) and (8.52) is implicit and is solved using a tri-diagonal linear equation solver. It is noted that the n+3/4 time level layer thickness is actually the n+1 time level thickness determined by the solution of (8.23). The specific discharges in (8.49) and (8.52) are given by (8.41) and represent those appearing in (8.23) and guarantee mass conservation for the pore water advection. The bed transport solution is completed by
BC
n 1 k

BC

n 3/ 4 k

BC

n 1 k

(8.54)

an implicit reaction step.

9.

References

Ackers, P., and W. R. White, 1973: Sediment transport: New approaches and analysis. J. Hyd. Div. ASCE, 99, 2041-2060. Ariathurai, R., and R. B. Krone, 1976: Finite element model for cohesive sediment transport. J. Hyd. Div. ASCE, 102, 323-338. Ambrose, R. B., T. A. Wool, and J. L. Martin, 1993: The water quality analysis and simulation program, WASP5: Part A, model documentation version 5.1. U. S. EPA, Athens Environmental Research Laboratory, 210 pp. Bear, J., 1879: Hydraulics of groundwater, McGraw-Hill, New York. Bagnold, R. A., 1956: The flow of cohesionless grains in fluids. Phil. Trans. Roy. Soc. Lond., Series A, Vol 249, No. 964, 235-297. Belleudy, P., 2001: Numerical simulation of sediment mixture deposition, part 1: analysis of flume experiments. J. Hyd. Res., 38, 417-425. Belleudy, P., 2001: Numerical simulation of sediment mixture deposition, part 2: a sensitivity analysis. J. Hyd. Res., 39, 25-31. Blumberg, A. F., B. Galperin, and D. J. O'Connor, 1992: Modeling vertical structure of open-channel flow. J. Hyd. Engr., 118, 1119-1134. Boyer, J. M., S. C. Chapra, C. E. Ruiz, and M. S. Dortch, 1994: RECOVERY, a mathematical model to predict the temporal response of surface water to contaminated sediment. Tech. Rpt. W-94-4, U. S. Army Engineer Waterways Experiment Station, Vicksburg, MS, 61 pp. Burban, P. Y., W. Lick, and J. Lick, 1989: The flocculation of fine-grained sediments in estuarine waters. J. Geophys. Res., 94, 8323-8330. Burban, P. Y., Y. J. Xu, J. McNeil, and W. Lick, 1990: Settling speeds of flocs in fresh and seawater. J. Geophys. Res., 95, 18,213-18,220. Cargill, K. W., 1985: Mathematical model of the consolidation and desiccation processes in dredge material. U.S. Army Corps of Engineers, Waterways Experiment Station, Technical Report D-85-4. Dortch, M., C. Ruiz, T. Gerald, and R. Hall, 1998: Three-dimensional contaminant transport/fate model. Estuarine and Coastal Modeling, Proceedings of the 5nd International Conference, M. L. Spaulding and A. F. Blumberg, Eds., American Society of Civil Engineers, New York, 75-89.

Einstein, H. A., 1950: The bed load function for sediment transport in open channel flows. U.S. Dept. Agric. Tech. Bull., 1026. Fredricks, C., and J. M. Hamrick, 1996: The effect of channel geometry on gravitational circulation in partially mixed estuaries. Buoyancy Effects on Coastal and Estuarine Dynamics, D. Aubrey and C. Fredricks, Eds., AGU, 283-300. Galperin, B., L. H. Kantha, S. Hassid, and A. Rosati, 1988: A quasi-equilibrium turbulent energy model for geophysical flows. J. Atmos. Sci., 45, 55-62. Garcia, M., and G. Parker, 1991: Entrainment of bed sediment into suspension. J. Hyd. Engrg., 117, 414-435. Gibbs, R. J., 1985: Estuarine Flocs: their size, settling velocity and density. J. Geophys. Res., 90, 3249-3251. Gibson, R. E., G. L. England, and M. J. L. Hussey, 1967: The theory of one-dimensional consolidation of saturated clays. Geotechnique, 17, 261-273. Hamrick, J. M., 1992: A three-dimensional environmental fluid dynamics computer code: Theoretical and computational aspects. The College of William and Mary, Virginia Institute of Marine Science, Special Report 317, 63 pp. Hamrick, J. M., 1994: Linking hydrodynamic and biogeochemcial transport models for estuarine and coastal waters. Estuarine and Coastal Modeling, Proceedings of the 3rd International Conference, M. L. Spaulding et al, Eds., American Society of Civil Engineers, New York, 591-608. Hamrick, J. M., and Wm. B. Mills, 2001: Analysis of temperatures in Conowingo Pond as influenced by the Peach Bottom atomic power plant thermal discharge. Environ. Sci. Policy, 3, s197-s209. Hamrick, J. M., and T. S. Wu, 1997: Computational design and optimization of the EFDC/HEM3D surface water hydrodynamic and eutrophication models. Next Generation Environmental Models and Computational Methods. G. Delich and M. F. Wheeler, Eds., Society of Industrial and Applied Mathematics, Philadelphia, 143-156. Hayter, E. J., and A. J. Mehta, 1983: Modeling fine sediment transport in estuaries. Report EPA-600/3-83-045, U.S. Environmental Protection Agency. Athens, GA> Hayter, E.J., M. Bergs, R. Gu, S. McCutcheon, S. J. Smith, and H. J. Whiteley, 1998: HSCTM-2D, a finite element model for depth-averaged hydrodynamics, sediment and contaminant transport. Technical Report, U. S. EPA Environmental Research Laboratory, Athens, GA.

Hwang, K.-N, and A. J. Mehta, 1989: Fine sediment erodibility in Lake Okeechobee. Coastal and Oeanographic Enginnering Dept., University of Florida, Report UFL/COEL89/019, Gainsville, FL. Ji, Z.-G., J. H. Hamrick, and J. Pagenkopf, 2002: Sediment and metals modeling in shallow river, J. Environ. Engrg., 128, 105-119. Jin, K. R., J. M. Hamrick, and T. S. Tisdale, 2000: Application of a three-dimensional hydrodynamic model for Lake Okeechobee, J. Hyd. Engrg., 106, 758-772. Jin, K. R., Z. G. Ji, and J. M. Hamrick, 2002: Modeling winter circulation in Lake Okeechobee, Florida. J. Waterway, Port, Coastal, Ocean Engrg., 128, 114-125. Karim, M. F., and F. M. Holley, Jr., 1986: Armoring and sorting simulation in alluvial rivers. J. Hyd. Engrg., 112, 705-715. Kleinhans, M. G., and L. C. Van Rijin, 2002: Stochastic prediction of sediment transport in sand-gravel bed rivers. J. Hyd. Engrg., 128, 412-425. Kuo, A. Y., J. Shen, and J. M. Hamrick, 1996: The effect of acceleration on bottom shear stress in tidal estuaries. J. Waterway, Port, Coastal, Ocean Engrg., 122, 75-83. Laursen, E., 1958: The total sediment load of streams J. Hyd. Div. ASCE, 84, 1-36. Letter, J. V., L. C. Roig, B. P. Donnell, Wa. A. Thomas, W. H. McAnally, and S. A. Adamec, 1998: A user's manual for SED2D-WES, a generalized computer program for two-dimensional, vertically averaged sediment transport. Version 4.3 Beta Draft Instructional Report, U. S. Army Corps of Engrs., Wtrwy. Exper. Sta., Vicksburg, MS. Le Normant, C., E. Peltier, and C. Teisson, 1998: Three-dimensional modelling of cohesive sediment in estuaries. in Physics of Estuaries and Coastal Seas, (J. Dronkers and M. Scheffers, Eds.), Balkema, Rotterdam, pp 65-71. Lick, W., and J. Lick, 1988: Aggregation and disaggregation of fine-grained lake sediments. J Great Lakes Res., 14, 514-523. Mehta, A. J., E. J. Hayter, W. R. Parker, R. B. Krone, A. M. Teeter, 1989: Cohesive sediment transport. I: Process description. J. Hyd. Engrg., 115, 1076-1093. Mehta, A. J., T. M. Parchure, J. G. Dixit, and R. Ariathurai, 1982: Resuspension potential of deposited cohesive sediment beds, in Estuarine Comparisons, V. S. Kennedy, Ed., Academic Press, New York, 348-362. Mehta, A. J., and F. Jiang, 1990: Some field observations on bottom mud motion due to waves. Coastal and Oeanographic Enginnering Dept., University of Florida, Gainsville, FL.

Mellor, G. L., and T. Yamada, 1982: Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. Space Phys., 20, 851-875. Meyer-Peter, E. and R. Muller, 1948: Formulas for bed-load transport. Proc. Int. Assoc. Hydr. Struct. Res., Report of Second Meeting, Stockholm, 39-64. Middleton, G. V., and P. R. Wilcock, 1994: Mechanics in the Earth and Environmental Sciences. Cambridge University Press, Cambridge, UK. Moustafa, M. Z., and J. M. Hamrick, 2000: Calibration of the wetland hydrodynamic model to the Everglades nutrient removal project. Water Quality and Ecosystem Modeling, 1, 141-167. Nielsen, P., 1992: Coastal bottom boundary layers and sediment transport, World Scientific, Singapore. Park, K., A. Y. Kuo, J. Shen, and J. M. Hamrick, 1995: A three-dimensional hydrodynamic-eutrophication model (HEM3D): description of water quality and sediment processes submodels. The College of William and Mary, Virginia Institute of Marine Science. Special Report 327, 113 pp. Rahmeyer, W. J., 1999: Lecture notes for CEE5560/6560: Sedimentation Engineering, Dept. of Civil and Environmental Engineering, Utah State University, Logan, Utah. Rahuel, J. L., F. M. Holly, Jr., J. P. Chollet, P. J. Belleudy, 1990: Modeling riverbed evolution for bedload sediment mixtures. J. Hyd. Engrg., 115, 1521-1542. Raukivi, A. J., 1990: Loose boundary hydraulics. 3rd Ed. Pergamon, New York, NY. Ried, I., and L. E. Frostick, 1994: Fluvial sediment transport and deposition. in Sediment Transport and Depositional Processes, K. Pye, ed., Blackwell, Oxford, UK, 89-155. Roberts, J., R. Jepson, D. Gotthard, and W. Lick, 1998: Effects of particle size and bulk density on erosion of quartz particles. J. Hyd. Engrg., 124, 1261-1267. Shen, J., J. D. Boon, and A. Y. Kuo, 1999: A modeling study of a tidal intrusion front and its impact on larval dispersion in the James River estuary, Virginia. Estuaries, 22, 681692. Shen, J. and A.Y. Kuo. 1999: Numerical investigation of an estuarine front and its associated topographic eddy. J. Waterway, Port, Coastal Ocean Engrg., 125, 127-135. Shrestha, P. A., and G. T. Orlob, 1996: Multiphase distribution of cohesive sediments and heavy metals in estuarine systems. J. Environ. Engrg., 122, 730-740.

Smagorinsky, J., 1963: General circulation experiments with the primative equations, Part I: the basic experiment. Mon. Wea. Rev., 91, 99-152. Smith, J. D., and S. R. McLean, 1977: Spatially averaged flow over a wavy bed. J. Geophys. Res., 82, 1735-1746. Smolarkiewicz, P. K., and T. L. Clark, 1986: The multidimensional positive definite advection transport algorithm: further development and applications. J. Comp. Phys., 67, 396-438. Smolarkiewicz, P. K., and W. W. Grabowski, 1990: The multidimensional positive definite advection transport algorithm: nonoscillatory option. J. Comp. Phys., 86, 355375. Spasojevic, M., and F. M. Holly, Jr., 1990: 2-D bed evolution in natural watercoursesnew simulation approach. J. Hyd. Engrg., 116, 425-443. Spasojevic, M., and F. M. Holly, Jr., 1994: Three-dimensional numerical simulation of mobile-bed hydrodynamics. Contract Report HL-94-2, US Army Engineer Waterways Experiment Station, Vicksburg, MS. Stark, T. D., 1996: Program documentation and users guide: PSDDF primary consolidation, secondary compression, and desiccation of dredge fill. Instructional Report EL-96-xx, US Army Engineer Waterways Experiment Station, Vicksburg, MS. Styles, R. and S. M. Glenn, 2000: Modeling stratified wave and current boundary layers on the continental shelf. J. Geophys. Res., 105, 24,119-24,139. Suzuki, K., H. Yamamoto, and A. Kadota, 1998: Mechanism of bed load fluctuations of sand-gravel mixture in a steep slope channel, Proc. of the 11th congress of APD IAHR, Yogyakarta, pp.679-688. Tsai, C. H., S. Iacobellis, and W. Lick, 1987: Floccualtion fo fine-grained lake sediments due to a uniform shear stress. J Great Lakes Res., 13, 135-146. van Niekerk, A., K. R. Vogel, R. L Slingerland, and J. S. Bridge, 1992: Routing of heterogeneous sedimetns over movable bed: Model development. J. Hyd. Engrg., 118, 246-262. Van Rijin, L. C., 1984a: Sediment transport, Part I: Bed load transport. J. Hyd. Engrg., 110, 1431-1455. Van Rijin, L. C., 1984b: Sediment transport, Part II: Suspended load transport. J. Hyd. Engrg., 110, 1613-1641.

Villaret, C., and M. Paulic, 1986: Experiments on the erosion of deposited and placed cohesive sediments in an annular flume and a rocking flume. Coastal and Oeanographic Enginnering Dept., University of Florida, Report UFL/COEL-86/007, Gainsville, FL. Vogel, K. R., A. van Niekerk, R. L Slingerland, and J. S. Bridge, 1992: Routing of heterogeneous sedimetns over movable bed: Model verification. J. Hyd. Engrg., 118, 263-279.

Wu, W., S. S. Y. Wang, and Y. Jia, 2000: Nonuniform sediment transport in alluvial rivers. J. Hyd. Res., 38, 427-434. Yang, C. T., 1973: Incipient motion and sediment transport. J. Hyd. Div. ASCE, 99, 16791704. Yang, C. T., 1984: Unit stream power equation for gravel. Journal of Hydraulic Engineering, 110, 1783-1797. Yang, C. T., and A. Molinas, 1982: Sediment transport and unit streams power function. J. Hyd. Div. ASCE, 108, 774-793. Yang, Z., A. Baptista, and J. Darland, 2000: Numerical modeling of flow characteristics in a rotating annular flume. Dyn. Atmos. Oceans, 31, 271-294. Ziegler, C. K., and B. Nesbitt, 1994: Fine-grained sediment transport in Pawtuxet River, Rhode Island. J. Hyd. Engrg., 120, 561-576. Ziegler, C. K., and B. Nesbitt, 1995: Long-term simulation of fine-grained sediment transport in large reservoir. J. Hyd. Engrg., 121, 773-781.