1
CHOESIV
SEDIMENT
TRANSPORT MODULE
ABSTRACT
The simulation of cohesive sediment transport processes is performed solving the 3D

conservative
advection

diffusion equation. For the representation of the processes of flocculation and erosion and
deposition
of bottom sediments refined empirical formulations based on data field was adopted. The
models have been calibrated and tested by simulating tidal flows and suspended sediment transport in
different estuaries. Results show a good agreement between the nume
rical predictions and
corresponding field measurements.
1. INTRODUCTION
The filtering role of estuaries makes them crucial transitional areas that trap
significant quantities of particulate and dissolved matter, originating a wide variety of
physical and
biogeochemical processes for the study of global change phenomena.
Unlike sand, well characterised by its grain size distribution, cohesive sediments are a
complex mixture of different clay minerals (mainly illite, montmorillonite and
kaolinite), organic
matter of diverse origin and nature, a small percentage of sand and
silt and diverse water content.
Hydrodynamics is the most important process involved on sediment transport,
nevertheless the time scales associated to sedimentation are much larger than t
he tidal
time scale, making residual flow a relevant aspect on sediment transport studies. The
complexity of the processes and the time scales involved on sediment transport make
their study a very difficult task for which mathematical modelling can be an
important
contribution.
A large number of mathematical models have been developed following the primary
works of EINSTEIN (1950). One of the earliest studies concerning the effect of the
flow (EINSTEIN and CHIEN, 1955) was undertaken by dividing the water
column
into a bed heavy fluid zone with relatively high sediment concentration, and a light
fluid zone occupying the remainder of the water column and with relatively low
sediment concentration. A velocity distribution was derived for each zone. Using
lab
oratory data, these authors showed that there was a correlation between the von
Karman constant,
K
, and the ratio of the power required to keep the sediment in
suspension to the power required to overcome the friction in the water column. It was
verified t
hat a minimum salinity of 1 gl

1
was required for the onset of flocculation.
The consolidation was quantified by noting the variation of the bed thickness with
time in two different phases and by order of aggregation.
The basic differences between the num
erical modelling approaches arise from the
number of dimensions considered, the numerical technique adopted, the complexity of
the description of bed evolution and the effects considered in the hydrodynamic
forcing e.g. tides, density flows, and wind waves
. Mainly for computer costs, until
2
recently most applications have been restricted to transport processes in one or two

dimension models.
One of the first numerical models for 1

D cohesive sediment transport was developed
by ODD and OWEN (1972). This finit
e difference, cross

sectionally averaged model
neglected the dispersive term in the transport equation and calculated the erosion and
deposition rates based in the formulations proposed by PARTHENAIDES (1965) and
KRONE (1962). In fact, it was a pseudo

2D m
odel since the depth was considered to
be divided into two different layers in which all variables were cross

sectionally
averaged, being the lower layer of constant thickness and thinner than the upper layer
of varying thickness. 1

D models are frequently
used to simulate sediment transport
and the large

scale morphological changes in rivers (DE VRIES
et al
, 1989), in tidal
channels (DYER and EVANS, 1989) and for the simulation of lutocline formation in
estuaries (ROSS and MEHTA, 1989; SMITH and KIRBY, 19
89).
A 2

D vertical integrated model was early presented by O´CONNOR (1971) in order
to obtain a more accurate description of the flow and transport phenomenon. Later
ARIATHURAI and KRONE (1976) developed a 2

D horizontal integrated finite
element model a
dopting triangular elements with a quadratic approximation for the
concentration and a Galerkian weighted residual method. The model utilises the
classical above mentioned relations for erosion and deposition. Continuing
aggregation is accounted for by spe
cifying different settling velocities in each element
for each time step. The authors defend that as the salinity of the water increases
different types of clays minerals become cohesive at different critical salinities.
MULDER and UDINK (1991), presents a
2

DH model applied to the Western Scheldt
estuary (The Netherlands). This model is based on the combined depth

integrated
equations of motion and sediment transport coupled to a wave model for shallow
water waves and producing stationary wave fields. It
solves a spectral action balance
with interpolation between the calculated wave heights and wave periods at different
tidal stages to determine the orbital velocity and bottom shear stress component due to
waves. Classical empirical expressions for sedimen
tation and erosion were used.
Constant values for critical shear stress of erosion and deposition and for settling
velocities were imposed. Consolidation was neglected. For the Gironde estuary
(France) LI
et al
(1994) developed a coupled 2

DV width

integr
ated hydrodynamic
and sediment transport model. A turbulent

closure model is described to compute
turbulent viscosity and diffusion coefficients. The source and sink terms of the
transport equation are calculated based in the classical formulations.
A th
ree

dimensional approaches are feasible nowadays using small computers and in
situations involving wind action, flow separation, density currents or important
curvature of the flow field. Even in well mixed estuaries there is a vertical
concentration grad
ient of suspended sediment due to exchanges between the bottom
and the water column and due to settling. In absence of density gradients the curvature
of the flow field can create secondary flows of great importance to the suspended load
transport. These s
econdary flows and the vertical concentration gradient justify by
themselves the interest of a 3

D modelling approach. SHENG (1986) discuss several
important aspects concerning the formulation of a comprehensive sediment dispersion
model. Particular attent
ion was given to four process models: bottom boundary layer
model, erosion model, deposition model and flocculation model. O´CONNOR and
3
NICHOLSON (1988) provided a fully

3D mud transport model accounting for
flocculation, deposition, erosion and consolida
tion. A direct comparison of a simple
zero

dimensional model and a complex 3

D model by applying the two schemes to a
field situation is discuss. KATOPODI and RIBBERINK (1992) developed a quasi

3D
model for suspended sediment transport based on an asympt
otic solution of the
advection

diffusion equation for currents and waves. A sensitivity analysis of the
influence of current and wave parameters on the suspended adjustment phenomenon,
shows that the presence of the latter considerably enlarges the suspen
ded load and its
adjustment time and length scales. The model results are compared with 2

DV and
fully 3D numerical solutions. Baroclinic 3

D models for hydrodynamics and coupled
fine sediment transport have been developed for coastal zones (DE KOK
et al
,
1995)
and for estuaries (CANCINO and NEVES, 1994 a, b; CANCINO and NEVES, 1995a,
b). The latter model has been calibrated and successfully tested by simulating tidal
flows and suspended sediment transport in different estuaries (two applications in the
Wes
tern Scheldt and Gironde estuaries are described in Part II of the paper).
3. SUSPENDED SEDIMENT TRANSPORT MODEL
The simulation of cohesive sediment transport processes is performed solving the 3D

advection

diffusion equation, in the same grid used by th
e hydrodynamic model
(CANCINO and NEVES, 1994 a, b; 1995a, b). The governing equation can be written
in the conservative form as:
(
)
(
)
(
)
(
)
C
t
uC
x
vC
y
w
W
C
z
s
x
C
x
y
C
y
z
C
z
x
y
z
(
)
(
)
(
)
(9)
where
C
is the suspended sediment concentration,
t
is time,
x, y
are the horiz
ontal co

ordinates,
z
is
the vertical co

ordinate,
x
,
y
,
z
are the sediment mass diffusion coefficients,
W
S
the sediment fall
velocity and
u, v w
the flow velocity components in
x, y,
z
directions.
The eddy diffus
ivity is assumed to be proportional to the eddy viscosity. The
horizontal one
x
,
y
) depends on the horizontal eddy viscosity which is a function of
the mesh size and of the turbulent energy dissipation rate, and can be calculated
according to Kolmogorov
´s theory (NIHOUL, 1984). The vertical diffusion
coefficient,
z
, is computed using Prandtl mixing length hypothesis (ROBERT and
OUELLET, 1987) corrected in stratified flows using the local Richardson number.
The horizontal transport is solved explicitly,
while the vertical one (including settling)
is solved implicitly, for numerical stability reasons. At each time

step the inversion of
a tridiagonal matrix is performed using the Thomas algorithm (THOMAS, 1949).
Equation (9) stands for conservative proper
ties. This is the case of sediments if we do
not simulate the formation of particles using colloidal matter in the areas of very small
salinity at the entrance of the estuary. In this case the total mass of suspended
sediments can change only due to fluxes
across the estuarine boundaries (open
boundaries, free surface and bottom). The fluxes across the open boundaries and
across the free surface must be imposed using field data. The fluxes across the bottom
interface must be calculated using the information
of concentrations in the water
4
column calculated by the model, the hydrodynamics and the properties of the bottom
sediments. The rate of accumulation of matter is calculated as the divergence of the
fluxes across the surfaces of the computational cell. So
, the boundary conditions have
been imposed in terms of fluxes as well.
In the sigma co

ordinate the vertical velocity is nil at the surface and at the bottom,
and so is the advective transport. The diffusive fluxes across the free surface and
across the
solid boundaries are zero. The exchanges with the bottom are accounted for
by fluxes imposed at bottom interface. Other than those conditions are needed at the
open boundaries. The advective fluxes across the solid boundaries vanishes as a
consequence of t
he normal velocities calculated by the hydrodynamic model and no
diffusive flux is considered at these boundaries.
3.1 The settling velocity
The exchange of sediment particles between the layers of the model can be due to
vertical advection, particle sett
ling or to turbulent diffusion. Vertical flow velocity and
turbulent diffusivity are computed by the hydrodynamic model. The settling velocity
of particles depends on the gravitational forces, and on vertical shear due to the
settling movement.
Gravitation
al forces depend on the density of each individual particle forming flocs
and on how flocs are formed. In fact a floc is formed by terraqueous and biological
particles but also by the water that remains between the individual particles.
The frictional for
ce depends on the form of the floc and on the Reynolds number of
the surrounding flow during settling. For very small particles the flow is laminar and
the ratio between the gravitational and the frictional force varies in the inverse relation
of the diame
ter of the particle. So, settling velocity increases with the size of the
particles. Nevertheless, there is no direct relation between floc size and settling
velocity, since larger flocs can have smaller density.
In a given suspension it should be expecte
d that when the size of the flocs increases
the average settling velocity of the suspended matter should increase, even if the flocs
with the larger settling velocities are not the biggest ones. To increase the size of the
flocs particles must aggregate. T
he probability of aggregation depends on the
probability of the particles to collide which directly increases with the concentration.
Consequently, settling velocity must increase with the concentration of the
suspension, until when the downward movement o
f the particles interferes with the
corresponding upward movement of the water. In this conditions an increase of the
concentration implies a decrease of the settling velocity. This concentration is called
the “hindering settling concentration”.
DYER, (198
6) proposes the following
correlation for the settling velocity:
W
S
=
K
1
C
m
:
C < C
HS
(10a)
5
W
S
= W
0
( 1.0

C
)
m
1
:
C > C
HS
(10b)
where,
W
S
is sediment particles fall velocity (m s

1
);
W
0
is fall velocity of a single particle (m s

1
);
C
HS
is the hin
dering settling concentration (kg m

3
); the constant of proportionality
K
1
(kg m

4
s) depends on
the mineralogy of the mud; exponents
m
and
m
1
depend on particle size and shape.
In the same work is suggested that
K
1
= 6.0 x 10

3
kg m

4
s;
m
= 1.0 and that
m
1
= 4.65
for small particles and one half for large particles.
The flocculation process depends on collision but also on cohesion of particles.
Relation 10 considers the effect of flocculation in differential settling. In the presence
of a growing concentrat
ion the probability of interparticle collision increases, allowing
the production of larger flocs with stronger settling velocity. The probability of
cohesion increases with salinity. Salinities between 1‰ and 2.5‰, are large enough to
allow intensive floc
culation inducing large fall velocities (WOLLAST, 1986).
3.2. Bottom interface module
There is not yet an agreement among the scientific community about the way how
sediment is exchanged between the water column and the bottom. For some
experimentalists e
rosion and deposition are simultaneous processes (e.g. SANFORD
and HALKA, 1993). When resuspension is stronger than deposition there is erosion.
The traditional way of dealing with the problem is following the ideas of EINSTEIN
(1950) and consider that ero
sion and deposition cannot occur simultaneously. In this
case critical shear stresses for erosion and deposition are defined and is assumed that
between those shear stresses there is no net flux at the bottom interface.
In a model both formulations can be
easily included. In this model we preferred the
traditional one because it is much easier to find data in the literature to specify the
parameters needed by the model.
3.3. The erosion model
The erosion algorithm is based on the classical approach of PART
HENIADES,
(1965). Erosion occurs when the ambient shear stress exceeds the threshold of erosion
of the bed particles. The flux of eroded material is measured by the rate of change of
sediments in the bottom (since is assumed that in erosion conditions ther
e is no
deposition):
dM
dt
E
E
E
1
:
>
E
(11a)
dM
dt
E
0
:
<
E
(11b)
6
where the left hand sides of Eqs. 11a and 11b are the erosion fluxes;
is bed shear stress ;
E
is a
critical shear stress for erosion and
E
is the erosion constant (kg m

2
s

1
).
The constant of erosion,
E
(Eq. 11a), depends on the physic

chemical characteristics
of sediment particles and is a space dependent variable. MU
LDER and UDINK,
(1991) suggest that this parameter is the order of magnitude of 5x10

5
kg m

2
s

1
in the
Western Scheldt.
Critical shear stress for erosion is a function of the cohesiveness of bottom sediments
and of their degree of compaction. The dry den
sity of the bottom sediments (ratio
between the mass of sediment after extraction of the interstitial water at 105
0
C and its
initial volume) can be used as a measure of the degree of compaction of bottom
sediments.
STEPHENS
et al
, (1992) based on the form
ulations proposed by DELO, (1988) use:
E
d
E
A
1
1
(12)
where,
E
is critical erosive stress (kg m

1
s

2
);
d
is dry density of bed sediments (kg m

3
);
A
1
and
E
1
are coefficients depending on mud type.
A
1
=0.0012 m
2
s

2
and
E
1
=1.2 are coeff
icients depending on
mud type.
For this equation to be dimensionally correct
A
1
must be a surface and
E
1
must be 1.
Nevertheless STEPHENS
et al
, (1992) calibrated their model with
A
1
= 0.0012 m
2
and
E
1
= 1.2. In this project we have also used those value
s but we admit that the
value of the coefficient
E
1
is in fact a measure of the error of this correlation plus the
error of the rate of compaction considered in the compaction module.
In fact this relation is a great simplification of reality, since the e
rodibility of a
cohesive bed is also a function of its cohesive nature depending, in a poorly
understood way on clay mineralogy and on the geochemistry and microbiological
processes occurring in the bottom. Some authors argue that it should also depend on
the salinity (HAYTER and METHA,1986). However, no dependency law as yet been
advanced.
3.4. The deposition model
The deposition algorithm, like the erosion algorithm is based on the assumption that
deposition and erosion never occur simultaneously. The a
lgorithm was first proposed
by KRONE (1962) and later on modified by ODD and OWEN, (1972). The
probability that a floc reach the bottom is one if the bottom shear stress is zero; this
probability is zero if the bottom shear stress is larger than the critic
al shear stress for
deposition and has a linear evolution between these two extreme conditions.
dM
dt
C
W
D
B
S
B
D
1
:
<
D
(13a)
7
dM
dt
D
0
:
>
D
(13b)
where
dM
D
/dt
is the mass deposited per unit area per unit time;
D
is the critical stress for
deposition;
(W
S
)
B
is the deposition velocity and
C
B
is evaluated near the sediment

water interface.
Critical shear stress for deposition is assumed a constant value of 0.2 kg m

1
s

2
based
on tidal experiments with natur
al mud from the Western Scheldt (WINTERWERP
et
al
, 1991).
8
5. REFERENCES
ARIATHURAI, R. and R.B. KRONE, 1976. Finite Element Model for Cohesive
Sediment Transport. J. Hydr. Div., ASCE, Vol. 102, No.HY3, p. 323

338.
BECKERS, J.M., 1991. Application of th
e GHER 3D General Circulation Model to
the Western Mediterranean. J. Marine Systems, Elsevier Science Publishers B.V.,
Amsterdam, 1: 315

332.
BLUMBERG, A. F. and G. L. MELLOR. 1983: Diagnostic and prognostic numerical
circulation studies of the South Atla
ntic Bight. J. Geophys. Res., 88, pp. 4579

4592.
CANCINO, L. and R. NEVES, 1994a. Numerical Modelling of Three

Dimensional
Cohesive Sediment Transport in an Estuarine Environment. World Scientific
Publishing,
107

121.
CANCINO, L. and R. NEVES, 1994b. 3D

Nu
merical Modelling of Cohesive
Suspended Sediment in the Western Scheldt Estuary (The Netherlands).
Netherlands
Journal of Aquatic Ecology, 28 (3

4), 337

345, 1994b.
CANCINO, L. and R. NEVES, 1995a.
Biogeochemistry of the Maximum Turbidity
Zone in Estuaries
. Final Report on Hydrodynamic and Transport Modelling
Contribution for CEC

Environment Programme, Inst. Superior Técnico, Lisbon, 80 p.
CANCINO, L. and R. NEVES, 1995b. Three

Dimensional Model System for
Baroclinic Estuarine Dynamics and Suspended Sedimen
t Transport in a Mesotidal
Estuary. In Computer Modelling of Seas and Coastal Regions. Computational
Mechanics Publications, 353

360.
DELLEERSNIJDER, E. and J.

M. BECKERS, 1992. On the Use of the sigma

coordinate System in Regions of Large Bathymetric Var
iations. Journal of Marine
Systems, 3, 381

390.
DELO, E.A., 1988. Estuarine Muds Manual. Report No. SR 164 , Hydraulics
Research, Wallingford, UK, 64 pp.
DE KOK, J.M., R. SALDEN, I.D.M. ROZENDAAL, P. BLOKLAND, J. LANDER,
1995. Transport Path of Suspended
Matter Along the Dutch Coast. In Modelling of
Seas and Coastal Regions. Computational Mechanics Publications, 000

000.
DE VRIES, M., G.J. KLAASSEN and N. STRIKSMA, 1989. On the Use of Movable
Bed Models for Rivers Problems: State of the Art, Symp. River
Sedimentation,
Beijing, China.
DYER, K.R. and E.M. EVANS, 1989. Dynamics of Turbidity Maximum in a
Homogeneous tidal Channel. J. Coastal Research, Special Issue No.5, Fort
Lauderdale, Florida, p. 23

30.
DYER, K.R., 1986. Coastal and Estuarine Sediment Dyn
amics. Wiley

Interscience,
New York, 342 pp.
9
EINSTEIN, H. A., 1950. The Bedload Function for Sediment Transportation in Open
Channel Flows. Soil Cons. Serv. U. S. Dept. Agric. Tech. Bull., No. 1026, 78 pp.
EINSTEIN, H. A. and N. CHIEN, 1955. Effects of He
avy Sediment Concentration
Near Bed on Velocity and Sediment Distribution. MRD Series Rep.8, Inst. Engrg.
Research, Univ. of Calif., Berkeley.
KATOPODI, I. and J. S. RIBBERINK, 1992. Quasi

3D Modelling of Suspended
Sediment Transport by Currents and Waves
. Coastal Engineering, 18, 83

110.
KRONE, R.B., 1962. Flume Studies of the Transport in Estuarine Shoaling Processes.
Hydr. Eng. Lab., Univ. of Berkeley, California, USA, 110 pp.
LARDNER, R.W. and Y. SONG, 1991. An Algoritm for Three

Dimensional
Convecti
on and Diffusion with Very Different Horizontal and Vertical Length Scales.
Int. Journal for Num. Meth. in Eng., 32: 1303

1319.
LI, Z.H, K.D. NGUYEN, J.C. BRUN

COTTAN and J.M. MARTIN, 1994. Numerical
Simulation of the Turbidity Maximum Transport in the G
ironde Estuary (France),
Oceanologica Acta.
MARTIN, J.

M. and H.L. WINDOM, 1991. Present and Future Roles of Ocean
Margins in Regulating Marine Biogeochemical Cycles of Trace Elements. In:
Montoura, J:

M:. Martin, R. Wollast (Eds). Ocean Margin Processes
in Global
Change, Wiley

Interscience, (45

689
MULDER, H.P.J. and C. UDINK, 1991. Modelling of Cohesive Sediment Transport.
A case study: The Western Scheldt Estuary. In: B.L. Edge, Ed., Proceedings of the
22nd International Conference on Coastal Engineerin
g. American Society of Civil
Engineers, New York, p. 3012

3023.
NIHOUL, J.C.J., 1984. A Three

dimensional General Marine Circulation Model in a
Remote Sensing Perspective. Annales Geophysicae, 2

4: 433

442.
O´CONNOR, B.A., 1971. Mathematical Model for Sed
iment Distribution. Proc. 14th
IAHR Conf., Paper D23, Paris, France.
O´CONNOR, B.A. and J.A. NICHOLSON, 1988. Mud Transport Modelling. In
Physical Processes in Estuaries, J. Dronkers and W. van Leussen (Eds.), Springer

Verlag, 532

544.
ODD, N.V.M. and M
.W. OWEN, 1972. A Two

Layer Model of Mud Transport in the
Thames Estuary. Proceedings, Institution of Civil Engineers, London, p. 195

202.
PARTHENIADES, E., 1965. Erosion and Deposition of Cohesive Soils. Journal of
the Hydr. Div., ASCE, vol. 91, No. HY1:
105

139.
PHILLIPS, N.A., 1957. A Coordinate System Having Some Spacial Advantages for
Numerical Forecasting.
J. Meteorol., 14: 184

186 .
10
ROBERT, J.L. and Y. OUELLET, 1987. A Three

Dimensional Finite Element Model
for the Study of Steady and Non

Steady
Natural Flows. In: J.C.J. Nihoul and B. M
Jamart, Ed., Elsevier Oceanography Series, 45: 359

372.
ROSS, M.A. and A.J. MEHTA, 1989. On the Mechanics of Lutoclines and Fluid
Mud. J. Coastal Research, 5 (51

62).
SANTOS, A.J and R. NEVES, 1991. Radiative Arti
ficial Boundaries in Ocean
Barotropic Models. In: A. S. Arcilla, Ed., Proceedings of the 2nd Int. Conf. on
Computer Modelling in Ocean Engineering, Barcelona, p. 373

383.
SANTOS, A.J., 1995. Modelo 3

D para Escoamentos de Superfície Livre. Thesis for
the d
egree of Doctor of Philosophy, I.S.T., Portugal.
SHENG, Y.P., 1986. Modeling Bottom Boundary Layer and Cohesive Sediment
Dynamics in Estuarine and Coastal Waters. In: A. J. Mehta, Ed., Estuarine Cohesive
Sediment Dynamics, Springer, Berlin, 14: 360

400.
SMITH, T.J. and R. KIRBY, 1989. Generation, Stabilization and Dissipation of
Layered Fine Sediment Suspensions. J. Coastal Research, 5 (63

73).
STANFORD L.P., and J.P. HALKA, 1993. Assessing the Paradigm of Mutually
Exclusive Erosion and Deposition of Mud
, with Exemples from Upper Chesapeake
Bay, Marine Geology, vol. 114, (37

57).
STEPHENS, J.A., R.J. UNCLES, M.L. BARTON and F. FITZPATRICK., 1992. Bulk
Properties of Intertidal Sediments in a Muddy Macrotidal Estuary. Mar. Geol. 103: 15
pp.
THOMAS, L.H., 1
994. Elliptic Problems in Linear Difference Equations Over a
Network. Watson Sci. Comput. Lab. Rept., Columbia University, New York. p. 373

383.
WOLLAST, R., 1986. The Scheldt Estuary. In: W. Salomon, B.L. Bayne, E.K.
Duursma and U. Forstner, Ed., Pollutio
n of the North Sea. An Assessment. Springer

Verlag, p. 183

193.
11
a)
b)
 h
z =
z =
h
ref.l evel
z
x
z =
z =
 h
c)
d)
x = L
3
3
x = 0
H
L
x = L
3
3
x = 0
Figure 1.
12
Fig. 1. Representation of an idealised vert
ical profile a) in a real domain, b) in a Cartesian model, c) in
a model using the sigma co

ordinate and d) in the sigma plan (adapted from Cancino and Neves,
1994a).
13
a)
b)
c)
d)
Figure 2.
14
Figure 2. Iso


lines and ideal
ised thermocline in (a) real space and (b) in the sigma co

ordinate. Iso


lines and idealised thermocline in (c) real space and (b) in the two

fold sigma co

ordinate, (adapted
from Beckers, 1991).
15
Σχόλια 0
Συνδεθείτε για να κοινοποιήσετε σχόλιο