Modelling Cohesive Sediment Dynamics in the Marine Environment

1. Introduction

The inflow of fine sediments in the marine environment affects various processes and has significant socioeconomic and environmental impacts, the most straight-forward of which is related to sedimentation in deltaic and coastal areas. However, fine sediments have the ability to affect solar radiation in the water column and to absorb contaminants on their surface [1,2]. Therefore, fine sedimentary plumes also affect primary productivity and the distribution of pollutants in the column and the seabed. The distinction between fine (<63μm) and coarse-grained (>63μm) sediments is rather general and varies with the type of the sedimentary matter, but it is a fact that the dominance of inter-particle cohesion forces over gravitational forces varies inversely with the floc diameter. Thus, the effect of cohesion is much more pronounced on the behaviour of clays (diameter <2μm) than of silts (2-63μm); in fact, the development of cohesion in silty clays is mainly due to the clay fraction in the sediment [1]. Adhesion diversifies the behaviour of cohesive and sandy sediments in the aquatic domain, adding further complexities to modelling efforts for fine sedimentary plumes. Their movement is the result of the combined action of various physical processes and forces (Figure 1).

The processes that control sediment transport in coastal microtidal areas can be divided in the ones that take place in the water column (pelagic processes):

• advection

• dispersion

• coagulation – deflocculation

• settling

• effects of seawater stratification

and the ones that take place at the boundary of the seabed (benthic processes):

• deposition

• erosion

• self-weight consolidation – resuspension

Figure 1.

The main processes that control sediment transport in coastal, deltaic areas.

These processes are analyzed in the following sections in terms of physical interpretation and mathematical parameterization, along with the three basic modelling approaches used for the description of sediments in the aquatic domain. Finally, results from a particle-tracking model applied in a micro-tidal shelf sea are presented and discussed.

2. Sediment transport models

Models describing transport and dispersion of sediments in the aquatic domain can be divided in Eulerian (e.g. [2]), that describe the evolution of suspended matter concentration using finite differences and the mass conservation principle, Lagrangian (e.g. [3]) in which the advection and dispersion of a specific amount of sediment mass is being traced with computational time and in mixed Euler-Lagrange models (e.g. [4]), in which transport is expressed by particles and dispersion through finite differences. The differences between these modelling approaches and their advantages and disadvantages are given further down.

2.1. Eulerian models

The evolution of sediment concentration in an Eulerian Model (hereafter EM), assuming mass conservation and incompressible flow, is written:

$$\frac{\partial c}{\partial t}+\frac{\partial cu}{\partial x}+\frac{\partial cv}{\partial y}+\frac{\partial c\left(w+w_s\right)}{\partial z}=\frac{\partial }{\partial x}\left(K_H\frac{\partial c}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_H\frac{\partial c}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_V\frac{\partial c}{\partial z}\right)+S-L$$ E1

In the equation, c is the mass concentration of suspended sediments, u, v and w are the fluid velocities along the x, y and z directions, w_s the sediment settling velocity, K_H and K_V the turbulent dispersion coefficients in the horizontal and vertical direction and S and L sediment source and loss terms, where applicable. The last three terms of the left-hand side express matter advection, while the first three terms of the right-hand side express matter dispersion. Regarding the source and loss terms, these can be mass exchanges in the bottom boundary, as well as sediment inflow (for example due to riverine discharges). Considering the free surface as the reference for the vertical axis (z=0), the bottom boundary condition reads:

$$K_H{\frac{\partial c}{\partial z}|}_{z=-H}=R_{er}-R_{dep}$$ E2

where R_er και R_dep are the erosion or resuspension and deposition rates, respectively.

2.2. Lagrangian models

Largangian (or particle-tracking) Models (hereafter: LM) are stochastic models that describe suspended matter in the column not as mass concentration in each of the cells of the computational grid, but as passive particles of specific mass that are being advected and dispersed by the flow. Thus the movement of a particle along direction i (i=1, 2 and 3 for directions x, y and z, respectively) is given by formulas of the form:

$$\frac{d{x}_i}{dt}={\overline{u}}_i+u_i{}^{\prime },u_i{}^{\prime }=\xi \sqrt{\frac{6K_i}{dt}},\xi \in \left[\begin{array}{cc}-1,& +1\end{array}\right]$$ E3

In equation 3, the first term expresses the advective part of the particle movement, whose range is defined deterministically by the velocity field, and the second the dispersive part that is defined stochastically using a random number ξ that ranges from -1 to +1; the amplitude of the stochastic part of the motion is determined by the corresponding mass dispersion coefficient Κ_i. It is noted that for the case of the vertical transport (i=3), the deterministic particle velocity is the summation of the vertical flow velocity and the settling velocity of the particle:

$${\overline{u}}_3=w+w_s$$ E4

In a tracer model particles in the domain are traced with simulation time and bear information (like the mass or the characteristic diameter of the particle, its origin etc) that follows the particle throughout its movement, remaining constant or varying with time. Each particle accounts for a specified mass of sediment and, thus, a computational particle represents a cumulus of actual sedimentary particles that are considered to have the same properties and fate.

2.3. Mixed Euler-Lagrange models

Euler-Lagrange Models (hereafter: ELM) are derived from the classic Eulerian advection-dispersion equation, which, applying a simplified temporal discretization is written [4-6]:

$$\frac{\partial c}{\partial t}+u\nabla c=\nabla K\nabla c\Rightarrow \frac{c^n-c^{n-1}}{\Delta t}+{\left[u\nabla c\right]}^{n-1}={\left[\nabla K\nabla c\right]}^n$$ E5

Considering an arbitrary, auxiliary value, c^f, the former equation can be divided into two parts. One that describes advection:

$$\frac{c^f-c^{n-1}}{\Delta t}+{\left[u\nabla c\right]}^{n-1}=0$$ E6

and one expressing dispersion:

$$\frac{c^n-c^f}{\Delta t}={\left[\nabla K\nabla c\right]}^n$$ E7

The advection equation states that concentration remains constant along lines of equal concentration that are defined by:

$${d{x}_i|}_{i=1,2,3}=dt\cdot u_i$$ E8

Thus, these lines are being ‘traced’ backwards in time for each note of the computational grid and the values for the temporal point n-1 are calculated for the corresponding ones of step n. The c^f concentrations for step n are determined using spatial interpolation. Finally, the dispersion equation is solved using finite differences or elements, through which the values of the concentration are updated for the next temporal step.

2.4. Common features and differences of the 3 models

One of the main disadvantages of the EM is that the advection-dispersion equation is solved in all of the grid points of the domain, which, combined with the necessity to adhere to the Courant-Friedrichs-Lewy (CFL) condition, increases the computational load. Contrastingly, a LM performs the necessary calculations only in the locations of particles while an ELM simulates the processes in a decoupled manner. It should be noted that, due to the stochastic nature of dispersion, simulating the same conditions with a random walk model can produce different results regarding the suspended matter in the domain. This effect is lower, the more the particles taken into account by the simulation. However, the number of particles in the domain is a parameter that should be taken into consideration, since increasing the number of particles increases the accuracy of the simulation but also increases computational load. Therefore, it is important to optimize the number of particles selected so that the simulation is accurate and also realizable within reasonable computational time. Even though ELMs present numerical advantages, their most significant disadvantage is that they do not inherently conserve mass, either locally or globally [6] due to tracking errors and non-conservative flow fields; these errors can be eliminated using very fine grids. Using concentrations, in both EM and ELM, provides erroneous results in the location of the sources and does not account for dispersion phenomena in a particular point, unlike LM [5] that describe the processes at the level of the cohesive floc. Furthermore, the resolution of a model that uses concentration to describe suspended mass in the domain is, unavoidably, defined by the discretization of the grid. Contrastingly, a particle-tracking model can describe the particulate matter distribution in scales much smaller than that of the spatial step. The transition from particle distribution to concentrations is easily performed, when and where this is required, from the mass of each particle and the number of particles in each grid cell. Finally, one of the most important and interesting advantages in using LM is the ability to assign information to particles (like origin, kind, state, values of characteristic parameters etc) in the form of indicators and to trace their evolution with time. Namely, knowing the ‘history’ of each sedimentary parcel allows the visualization the trajectory of each particle and of the changes it has undergone, as well as the segregation of particles according to common properties and the investigation of transport patterns related to sedimentary origin.

3. Physical processes for cohesive sediments in the marine environment

3.1. Pelagic processes

The pelagic processes are the ones that take place in the water column and can generally be described as advection and dispersion. However, flocculation and the corresponding changes in the particles’ settling rate and stratification affect advective and dispersive matter transport in the column; these processes are described following.

3.1.1. Turbulent dispersion

It is important to clarify the distinction between diffusion and dispersion of suspensions. Mass diffusion can be molecular in the case of laminar flow, which is the case for example in the stochastic Brownian motion of particles, or turbulent, in the case that the flow in the field presents turbulent eddies. For the movement of particles to be considered diffusive, the existence of a concentration gradient in suspended matter is required. In the opposite case, when the concentration in the field is homogenous, dispersion theoretically has no effect since the lateral expansion of the plume takes place in the same rate in both directions. The mass diffusion coefficient ranges from 10^-9m²/s to 10²m²/s for molecular and turbulent diffusion, respectively. Contrastingly, the case of mass dispersion requires mixing of the flow field due to velocity gradients and the corresponding coefficient ranges from 10^-3m²/s to 10²m²/s. In coastal and marine systems the seawater velocities are variable with depth and along the direction of the flow and, therefore, turbulent dispersion is the main driving mechanism. The scale of turbulent eddies is important to the evolution of a suspended sedimentary plume. In cases that the eddy is much larger than the plume, the movement is essentially advective and the shape of the plume remains unchanged, whereas for intermediate eddies the plume is deformed and its boundaries are stretched; when turbulent scales are much smaller than the plume itself, the main effect is dispersive and leads to evening of the shape of the plume [7].

Richardson [8] studied the evolution of sedimentary plumes from a point, non-constant source considering the segregation between the particles of the plume. Defining standard deviation of particles from their average position as the characteristic length-scale l, he reached a correlation with the dispersion coefficient of the form:

$$K=0.2\cdot l^{\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$$ E9

where l is the length scale of the system. This equation is known as the ‘four-thirds power law’. One of the most widely applied relationships for the calculation of the turbulent dispersion coefficient from the hydrodynamic conditions is the Smagorinski formula [9,10]:

$$K_H=C\cdot dx\cdot dy\sqrt{{\left(\frac{\partial u}{\partial x}\right)}^2+\frac{1}{2}{\left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2}$$ E10

$$K_V=\frac{d{z}^2}{dx\cdot dy}{K}_H$$ E11

C, in equation 10, is a constant that depends on the horizontal discretization step and varies from ~0.01 for steps of a few meters to 0.2 for steps of the order of kilometres, while the subscripts H and V (eq. 10-11) denote horizontal and vertical directions, respectively.

The dispersion coefficients can also be estimated from the corresponding values of turbulent viscosity of seawater (v_ti) using relationships of the form [11]:

$$K_i={\beta }_i{F}_{mi}{v}_{ti}$$ E12

The introduction of coefficients β_i and F_mi is due to the separation of seawater and sediment mixing and to the reduction of vertical mixing, respectively. They are usually taken equal to one, apart from the case of stratified environments, where the coefficient F_mi is estimated as a damping function of the vertical transport; the effect of stratification is analyzed further down.

3.1.2. Flocculation - Settling

A basic distinction of fine-grained from sandy sediments is the development of cohesion that enables the formation of bonds between particles and eventually the formulation of coagulants of larger diameter than their primary constituents. These bonds may be due to attractive van der Waals or electrochemical forces. Van der Waals forces are particularly strong but weaken significantly with distance, which means that they are effective only in the case that the particles are at very low distances. Suspended oblate argillaceous particles are usually negatively charged in their surface and bear positive charge in the edges; the overall electrostatic charge in the surface is negative and therefore repulsive forces exist between particles. Free ions and cations in the marine environment form a cloud in the surface of clay particles, known as (Guoy) double layer [12], that balances the electrostatic charge of the particles and allows them to come closer together, collide and form bonds.

The flocculation process is mainly governed by three mechanisms [13-16]:

1. Random Brownian movement causes collisions between particles and their attachment to a larger coagulant.

2. Turbulence in flows causes particles to collide and aggregate, on one hand, while turbulent shear can induce the breakup of an aggregate, when the shear overcomes the shear strength of the bond between the primary particles.

3. Larger particles settling through a suspension can drift smaller particles, which attach to them, thus forming larger aggregates (differential settling).

It should be noted that the contribution of Brownian motion to flocculation is significant only for small particle diameters (~1-2μm) and results to aggregates of smaller strength, compared to the other two mechanisms [13]. As the diameter of the particle increases, shear becomes the main criterion for the evolution of the process [17].

Other parameters that affect the formation of clayey aggregates, favouring the collision of particles are [18]:

1. Salinity; it has been established that transition to a higher salinity environment leads to increase of the diameter of cohesive suspended matter. The critical salinity value for coagulation, S ₀, related to the clay minerals contained the clay fraction are 0.6, 1 and 2.4 for kaolinite, illite and montmorillonite, respectively [16].

2. The existence of cations increases the attractive forces between particles and, therefore, the aggregate’s strength.

3. The increase of temperature.

4. The existence of organic matter and detritus, bacteria, microphytes and macrophytes. It is essentially a form of biostabilization of the aggregate’s structure by extracellular secretions of the microorganisms.

5. The increase of suspended matter concentration, which increases the frequency of collisions between particles and, therefore, the possibility of adhesion to a larger coagulate.

Dyer [19] established the theoretical background for the evolution of aggregate diameter under the effects of flow shear and suspended concentrations. Concentration increases the diameter of the floc, due to higher availability of matter, while shear stress increases flocculation at low values, favouring inter-particle collisions, and causes floc-breakup and corresponding decrease of the particle diameter for values higher than the strength of the aggregate. Mathematical expressions that can describe this behaviour are exponential formulas with Munk-Anderson type damping functions, of the form [20]:

$$D_{ag}=a_1\frac{1+b_{1}G}{1+b_2{G}^{d_1}}{c}^{d_2}+a_2$$ E13

The constant parameters of equation 13 depend on the concentration of the suspension, c. For high concentrations (c=3.7g/l) these constants have been experimentally estimated at a₁=0.08, a₂=0.86, b₁=29, b₂=125, d₁=4.6 and d₂=-1.1, while for lower concentrations (c=0.8g/l) the same constants read a₁=0.06, a₂=0.88, b₁=33, b₂=138, d₁=3.8 and d₂=-3.5 [21]. The turbulent dissipation parameter G expresses the effects of shear to the diameter evolution and depends on the energy dissipation rate ε and kinematic viscosity v: $$G=\sqrt{\epsilon /\nu }$$

A method widely applied for the parameterization of flocculation is based on the assignment of fractalic dimensions to the cohesive flocs, considering that every aggregate consists of self-similar fractalic entities. Accepting that the linear dimension of the characteristic particle is unitary, it follows that the particle volume is proportional to the number of primary particles. Thus, if n_f is its fractalic dimension, N the number of primary particles that form the aggregate and L its linear dimension, it follows that [15]:

$$n_f=\underset{L\to \infty }{\lim }\frac{\ln \left[N\left(L\right)\right]}{\ln \left[L\right]}$$ E14

Thus, the number of particles in an aggregate is proportional to the diameter of the primary particle, D_p:

$$N~D_p{}^{n_f}$$ E15

The fractalic dimension for cohesive particles ranges from 1.4 for very fragile flocs, like marine snow, to 2.6-2.8 for dense, consolidated beds. It is generally accepted that the fractalic dimension increases from marine snow (1.4-1.8) to estuarine sediments (1.7-2.2) and to oceanic sediments (1.94-2.14) [22,23]. An average fractalic dimension for estuarine and coastal areas is n_f=2 [15].

Applying the fractal theory to the evolution of cohesive aggregate diameter it follows that there is an analogy between the ratio of aggregate and primary particle dimensions and the corresponding volume concentration, $$\varphi$$ [24]:

$$\frac{D_{ag}}{D_p}={\left(\frac{\varphi }{{\varphi }_p}\right)}^{\frac{1}{3-n_f}}$$ E16

The relationship between volume concentration and aggregate mass is written [15]:

$$\varphi =\frac{{\rho }_s-{\rho }_w}{{\rho }_f-{\rho }_w}\frac{c}{{\rho }_s}=f_{s}N{D}_{ag}{}^3$$ E17

ρ_s, ρ_w and ρ_f are solid phase, water and aggregate densities and f_s is a shape constant which, for the case of spherical particles equals to π/6. The increase in aggregate diameter also involves decrease in its density, ρ_f, compared to the primary particles, ρ_p, due to increase in interstitial water [12-14,16,17]. Considering the density difference of solid and liquid phase, the density of the aggregate can be calculated by [24]:

$${\rho }_f={\rho }_w+\left({\rho }_s-{\rho }_w\right){\left(\frac{D_p}{D_{ag}}\right)}^{3-n_f}$$ E18

Based on the fractal theory, the evolution of the aggregate diameter in turbulent flow is [15]:

$$\frac{d{D}_{ag}}{dt}=k_{A}cG{D}_{ag}{}^{4-nf}-k_B{G}^{q+1}{\left(D_{ag}-D_p\right)}^p{D}_{ag}{}^{2q+1}$$ E19

The first term of the right hand-side of equation 19 expresses the effect of coagulation and the second the effect of floc break-up. The corresponding dimensional flocculation and deflocculation parameters (k_A [m²/kg] and k_B [sec^1/2/m²]) are defined as:

$$k_A={k^{\prime }}_A\frac{D_p{}^{n_f-3}}{n_f{f}_s{\rho }_s}{k}_B={k^{\prime }}_B\frac{D_p{}^{-p}}{n_f}{\left(\frac{\mu }{F_y}\right)}^q$$ E20

In equations 18-20, c is the suspended sediment concentration, D_p the diameter of the primary flocs, q and p are constants (q~0.5 and p=3-n_f), μ is the dynamic viscosity of water, F_y is the strength of the sediment (F_y=O{10^-10}) and k_A’ and k_B’ are non-dimensional flocculation and deflocculation coefficients (k_A’=O{10^-1} and k_Β’=O{10^-5}) [15].

The corresponding changes to the settling velocity of the aggregate can be calculated by the Sto- kes velocity, taking into account shape factors α and β (α=β=1 for perfectly spherical particles):

$$w_s=\frac{\alpha }{\beta }\frac{\left({\rho }_f-{\rho }_w\right)g}{18\mu }{D}_{ag}{}^2$$ E21

The shape factors can be simplified considering an irregularity parameter for the flocs, a_irr, (a_irr≤1) equal to the ratio of minimum to maximum aggregate dimensions. Following [25] and after some manipulations:

$$\left\{\begin{array}{l}\alpha =a_{irr}\\ \beta =1-{\lambda }_1\left(1-a_{irr}\right)\end{array}\right\}$$ E22

where λ₁ is a coefficient related to the orientation of the particle during settling. It equals 0.2 or 0.4 when the particle settles with its maximum axis parallel or perpendicular to depth. Considering parallel orientation as the most frequent one, the settling velocity reads:

$$w_s=\frac{a_{irr}}{0.8+0.2a_{irr}}\frac{\left({\rho }_f-{\rho }_w\right)g}{18\mu }{D}_{ag}{}^2$$ E23

The relation between settling velocity, flow shear stress and suspended sediment concentration (c in mg/l) for microflocs (<160μm) and macroflocs (>160μm), defined after experimental measurements from estuarine areas is [26-29]:

$$w_{s_{\mathrm{ma}cro}}=\{\begin{array}{l}0.644+0.000471c+9.36\tau -13.1{\tau }^2,0.04Pa\le \tau \le 0.65Pa\\ 3.960+0.000346c-4.38\tau +1.33{\tau }^2,0.65Pa\tau \le 1.45Pa\\ 1.18+0.000302c-0.49\tau +0.057{\tau }^2,1.45Pa\tau \le 5.00Pa\end{array}$$ E24

$$w_{s_{micro}}=\{\begin{array}{l}0.244+3.25\tau -3.71{\tau }^2,0.04Pa\le \tau \le 0.52Pa\\ 0.65{\tau }^{0.541},0.52Pa\tau \le 10.0Pa\end{array}$$ E25

The predominance of macroflocs versus microflocs can be defined from their concentration quotient, which was found to depend on the overall suspended sediment concentration through a relationship of the form [26]:

$$\frac{c_{macro}}{c_{micro}}=0.815+0.00318c-0.14\cdot {10}^{-6}{c}^2$$ E26

The settling velocity of particles can be distinguished in three ‘phases’, depending on the suspended matter concentration [14]. For low concentrations, C₁=0.1-0.3g/lt, free (Stokes’) settling dominates. For concentrations higher than the limit for free settling and lower than C₂=2-20g/lt settling is highly affected by flocculation and involves increase of the settling velocity due to increase of the particle diameter. The settling velocity can be thus expressed as an exponential function of concentration:

$$w_s=k_1{c}^{n_1},c_1\le c\le c_2$$ E27

where k₁ (~4/3) and n₁ Ο{10^-1} are experimentally determined constants. Taking into account the effect of shear stress, the settling velocity can be determined by the equation [30]:

$$w_s=w_{s,r}\frac{1+aG}{1+b{G}^2}$$ E28

In equation 28, w_s,r is the settling velocity of a single particle in still liquid and a and b (O{1}) are empirical coefficients related to the sediments.

For concentrations higher than C₂ the settling velocities reduce with increasing suspended matter, phenomenon known as hindered settling. Hindered settling essentially expresses the effect of neighbouring particles to the settling velocity of a specific particle, which is important especially for high concentration suspensions. This hindering effect is caused by eddy separation around settling particles, causing reduction to the settling rates of particles within this buoyant flow, but also due to inter-particle collisions and electrochemical interactions. The hindered settling velocity related to the volumetric concentration is:

$$w_s=w_{s,r}{\left(1-\varphi \right)}^n$$ E29

n is a parameter that varies inversely with the particle Reynolds number within the range 2.5<n<5.5. A corresponding expression [15] reads:

$$w_s=w_{s,r}\frac{{\left(1-{\varphi }_{\ast }\right)}^m\left(1-{\varphi }_p\right)}{1+2.5\varphi },{\varphi }_{\ast }=\min \left\{1,\varphi \right\}$$ E30

where ϕ_p is the volumetric concentration of primary particles and m is a constant for the inclusion of possible non-linear effects. The expression that relates the mass concentration in the field to the hindered settling velocity is similar [14]:

$$w_s=w_{so}{\left[1-k_2\left(c-C_2\right)\right]}^{n_2}$$ E31

For highly concentrated suspensions (~80-100g/lt) the settling velocity is practically zero.

A general formulation that can describe the evolution in conditions of flocculation and hindered settling is [14]:

$$w_s=\frac{a\cdot c^n}{{\left(c^2+b^2\right)}^m}$$ E32

The coefficients are determined experimentally; from experimental studies in various estuarine and riverine areas [31] the range of their values is: a=0.01-0.23, b=1.3-25.0, n=0.4-2.8 and m=1.0-2.8. The maximum settling velocity in this case (for concentration equal to C₂) is [14]:

$$w_{so}=a\cdot b^{n-2m}\frac{{\left(\frac{2m}{n}-1\right)}^{m-2n}}{\sqrt{\frac{2m}{n}}},C_2=\frac{b}{\sqrt{\frac{2m}{n}-1}}$$ E33

3.1.3. Effects of stratification

One of the most widely applied methods for the determination of the stability of the water column is the use of the gradient Richardson number [32]:

$$Ri=\frac{N^2}{{\left(\frac{\partial U}{\partial z}\right)}^2}\frac{1}{4},N=\sqrt{-\frac{g}{{\rho }_w}\frac{\partial {\rho }_w}{\partial z}}$$ E34

Ν is the Brunt-Väisälä frequency. In flows with Ri values lower that ¼ Kelvin-Helmholtz instabilities can occur. Based on this criterion, the effect of stratification to mass dispersion, based on the mixing-length theory, is often expressed using Munk-Anderson type velocity and mass dispersion damping functions (F_m) [7]:

$$F_m={\left(1+\beta Ri\right)}^{\alpha }\left\{\alpha \simeq -0.5,\beta \simeq 10\right\}$$ E35

The flux Richardson number, which is the ratio of the buoyant forces to the turbulent kinetic energy production [7], is also frequently used; the mass dispersion coefficient K_m is [33]:

$$K_m=\frac{R_f}{1-R_f}\frac{\overline{\epsilon }}{\overline{N^2}}\stackrel{R_f\ge R_{f,crit}\simeq 0.15}{\to }{K}_m\le 0.2\frac{\overline{\epsilon }}{\overline{N^2}}$$ E36

The overbar in equation 36 indicates time averaging. Regarding the mixing of freshwater plume, expressed as entrainment velocity to the underlying layer (W_e), it can be estimated as [34]:

$$W_e=\frac{e_1{U}_o}{Ri}$$ E37

where U_o is the velocity of the underlying layer and e₁ (O{10^-3}) an empirical coefficient.

Another indicator for the strength of the stratification is the Peclet number (Pe) that expresses the predominance of settling over dispersion (h: width of the mixing zone):

$$Pe=\frac{w_{s}h}{K_m}\ll 1$$ E38

It has recently been established that, even in the case that density increases with depth, the hydrostatic stability of the system is not guaranteed [35], since inhomogeneities in the density field can lead to transporting water masses that can induce a process of matter removal from suspension. The prognosis of such instabilities cannot be described using a Richardson number [33]. Double diffusion is the development of instable density gradients due to differential molecular diffusion of the properties of the fluid, which are capable of changing its density. Given that the thermal diffusion coefficient (K_T≈10^-3m²/s) is two orders of magnitude larger than the salinity coefficient (K_S≈10^-5m²/s), it follows that temperature diffuses faster; this can lead to areas in the flow with higher density than the underlying fluid, thus causing mixing of the water column. Double diffusive instabilities can form as salt fingers, in the case that salinity is the destabilizing parameter (T and S decreasing with depth), or as diffusive layers, when temperature is the destabilizing factor (T and S increasing with depth). Salt fingering instability takes the form of tightly packed blobs of sinking salty and rising fresher water masses near the thermocline, which quickly develops to transport away from the interface. Diffusive layering convection involves a buoyant thermal diffusion flux through the pycnocline that is higher (in terms of density) than the corresponding salinity flux, thus leading to a downward density flux that causes transport from the stratified zone to the homogenous layer. The susceptibility of the column to double diffusive phenomena can be estimated using the stability ratio (R_p), which is the ratio of the stabilizing parameter gradient to the destabilizing one, expressed in terms of density:

$$R_p=\{\begin{array}{l}\alpha \frac{\partial T}{\partial z}/\beta \frac{\partial S}{\partial z}\to \text{salt}\text{fingering}\\ \beta \frac{\partial S}{\partial z}/\alpha \frac{\partial T}{\partial z}\to \text{diffusive}\text{layering}\end{array},\left\{\begin{array}{c}\alpha =-\frac{1}{{\rho }_w}\frac{\partial {\rho }_w}{\partial T}\\ \beta =\frac{1}{{\rho }_w}\frac{\partial {\rho }_w}{\partial S}\end{array}\right\}$$ E39

In equation 39, α is the thermal expansion coefficient and β is the haline contraction coefficient. The gravitational stability of the column increases with the stability ratio; it has been established that double diffusive instabilities develop within the range 1<R_p<10 [36,37], while areas where R_p≥10 are highly stable, with low settling rates of fine sedimentary plumes (~cm/h) [37]. Thus, R_p=10 can be used as a threshold value for the inhibition of double diffusive instabilities and, therefore, for the indication of a highly stable column.

These effects can be taken into account in sediment transport modelling, using damping functions to the vertical advection (F_w) and dispersion (F_Kv) of matter [38]:

$$F_w=\{\begin{array}{l}1,R_{p}10\\ {\left(1+\gamma R_p\right)}^{\delta },R_p\ge 10\end{array}$$ E40

$$F_{Kv}=\{\begin{array}{l}1,R_{p}10\\ \exp \left(-\kappa \cdot R_p\right),R_p\ge 10\end{array}$$ E41

3.2. Benthic processes

Modelling the transport of cohesive sediments in coastal and deltaic areas requires the description of the processes of mass exchange between the column and the seabed through the processes of erosion and deposition. In practice, formulating this module is one of the most challenging tasks in modelling fine sediment transport, mainly due to the lack of general and widely applicable formulae between erosion resistance and bed shear stress. However, even if the main questions regarding mass exchange quantification were answered, the fact that the characteristic parameters vary within the sediment with time and space (horizontally and vertically) still remains [39].

There are three main distinctions regarding fine sediment transport in the boundary of the bed, compared to coarse-grained sediments [40]:

1. Cohesive matter is transported exclusively in suspension, while coarser matter is also transported in semi-contact with the bed, as bed load.

2. Cohesive sediments are not transported as dispersed particles. Flocculation increases settling rates and is responsible their deposition.

3. Cohesive sediment beds undergo self-weight settling and consolidation. In cases of rapid sedimentation the deposited sediments are light and present minimal resistance to shear stress. Cohesive beds can be homogenous or vertically stratified regarding density and strength.

Various modules have been proposed for the exchange processes between seabed and water column, which are mainly based on the transition of the matter to different states. Four characteristic states for a sediment-water mixture are defined in [1]:

1. Horizontally advected mobile suspension.

2. Concentrated benthic suspension (CBS), horizontally static but with vertical mobility.

3. Consolidating soft deposit (fluid mud).

4. Stable, consolidated part of the seabed.

The static suspension is formed by the deposition of the transported suspension, especially under low hydrodynamic conditions, and presents minimal mechanical strength. As deposition progresses, a CBS is formed; CBSs are suspensions with strong interactions with the turbulent flow field through buoyant effects, which however continue to behave as Newtonian fluids [41]. In time, consolidation and the related physicochemical changes lead to strength development and, eventually, to the formation of a fluid mud layer with lower interstitial water, higher shear strength and more stable structure. It must be noted that the transition from one phase to another is gradual and that a change in the flow conditions may lead to phenomena comparable to ones that an increase in concentration would produce. Whether these phenomena take place or not depends on the temporal scales of the processes; for example the time scale of deposition and flocculation define if fluid mud is formed under low flow conditions and the temporal scale of consolidation determines if re-entrainment or erosion phenomena partake during accelerated flow. Regarding the concentration/density (dry and bulk, ρ_d and ρ_b) ranges of the aforementioned suspensions and deposits, the following values are cited [42]:

• Dilute suspension (low concentration C=0.1-10kg/m³)

• CBS (high concentration C=10-100kg/m³)

• Fluid mud (C=ρ_d=100-250kg/m³, ρ_b=1050-1150kg/m³)

• Partially consolidated bed (ρ_d=250-400kg/m³, ρ_b=1150-1250kg/m³)

• Fully consolidated bed (ρ_d=400-550kg/m³, ρ_b=1250-1350kg/m³)

The most straight-forward treatment of the mass-exchange in the boundary of the bed is derived from the balance between settling, deposition and erosion, through conditions of the form [39]:

$${\epsilon }_s\frac{\partial c}{\partial z}-w_{s}c=S_E-S_D$$ E42

In the former relationship, the right hand-side is the balance between turbulent dispersion and settling and the right is the balance of erosion and deposition rates.

3.2.1. Bottom shear stress

The shear stress velocity in the bed can be related to the depth-averaged flow velocity (U) through a Chezy coefficient, Ch [11]:

$$u_{\ast }=U\sqrt{g}/Ch,Ch=18\log \left(12h/k_S\right)$$ E43

where h is the flow depth and k_s an equivalent roughness height. The direct dependence if the shear stress to the flow velocity at a small distance from the bed can also be used [7]:

$${\tau }_b=\rho \cdot C_d\cdot u{(1m)}^2$$ E44

In equation 44, u(1m) is the velocity at a distance of 1m from the bed and C_d the bottom friction coefficient (Ο{10^-3}). This coefficient can be estimated assuming logarithmic velocity profile through equations of the form [2]:

$$C_d=\max \left\{\begin{array}{cc}{\left[\frac{1}{\kappa }\frac{\ln \left(H+z_b\right)}{z_o}\right]}^{-2}& ,0.0025\end{array}\right\}$$ E45

z_o is the roughness length of the bed (Ο{10^-3m}). High density gradients near the bed result to reduction of the apparent roughness with corresponding increase in transport for the same bottom shear stress. Thus, the actual erosion rates for the same flow conditions are lower, compared to the ones calculated not taking into account these buoyancy effects in the bottom boundary conditions. A significant modelling improvement is achieved using damping factors in the bottom boundary conditions; the calculation of the bottom shear can be performed through [43]:

$$u_{*}=F_t\kappa z\frac{\partial U}{\partial z},F_t={\left(1+100Ri\right)}^{\raisebox{1ex}{$-1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$$ E46

where F_t is the velocity damping factor, κ the von Karman constant, ∂U/∂z the velocity gradient at the marginal grid of the bed and Ri the gradient Richardson number. The damping function essentially expresses the ratio of the actual turbulence to the turbulence in the absence of suspended matter. Based on this approach, the gradient at the bed is calculated explicitly from the velocity profile and not using the log-law. Therefore, the ‘corrected’ velocity is:

$$U=\frac{u_{*}}{\kappa z}\ln \left(\frac{z}{\alpha z_o}\right),\frac{1}{F_t}=1-\frac{z}{\alpha }\frac{\partial \alpha }{\partial z}$$ E47

where α is a bottom roughness correction coefficient. The reduction in bottom roughness corresponds to drag reduction, observed both in nature and in the laboratory. Series of numerical experiments [43] suggest that α can be estimated through exponential functions of the form:

$$\alpha =\exp \left[-\left(1+\beta w_s/u_{*}\right)\left(1-\exp \left(-bR{i}^m\right)\right)\right],\left\{\beta =7.7,b=1.25,m=0.6\right\}$$ E48

In coastal areas, where the activity of both waves and currents is significant, their combined effect to the range of the shear stresses can be addressed through equations of the form [2]:

$${\tau }_b=\frac{1}{2}{f}_{cw}{\rho }_w{\left|u_c+u_w\right|}^2$$ E49

In the former relationship u_c and u_w are the current and wave amplitude velocities near the bed and f_cw the wave-current friction factor.

3.2.2. Deposition

The sediment deposition rate (S) is typically defined using the Krone formula [44]:

$$S=\{\begin{array}{l}0,\left|{\tau }_b\right|\ge {\tau }_{cr,dep}\\ c_b\cdot w_s\left(\frac{\left|{\tau }_b\right|}{{\tau }_{cr,dep}}-1\right),\left|{\tau }_b\right|{\tau }_{cr,dep}\end{array}$$ E50

The particulate matter that can be transported by turbulent flow may be expressed by the equilibrium concentration C_e [15,41]:

$$C_e=K_s\frac{{\rho }_b}{g}\frac{{\rho }_s-{\rho }_w}{{\rho }_s}\frac{u_{*}{}^3}{h{W}_s}$$ E51

where Κ_s is a first-order proportionality coefficient (Κ_s~0.7) and W_s is the characteristic (mean) settling velocity. It can be assumed that, in the case that the system is in balance, this velocity equals the entrainment velocity from a CBS, with the sediment suspension to be in equilibrium. At the fluid mud-water interface little or no turbulence is produced. Therefore, as turbulent mixing reduces, C_e is further reduced. This results to a cumulative effect and the complete breakdown of the vertical suspended concentrations’ profile. For cohesive sediments the equilibrium concentration can also be considered as saturation concentration.

The critical shear stress for deposition can be related to the strength of the cohesive floc, F_c (Ο{10^-8Ν}), through relationships of the form [24]:

$${\tau }_{cr,dep}=\frac{1}{5\pi }\frac{F_c}{D^2}\simeq 0.06\frac{F_c}{D^2}$$ E52

Critical shear velocity for deposition can also be expressed with respect to the settling velocity of the particles [44]:

$$u_{\ast cr,dep}=\{\begin{array}{l}0.008,w_s\le 5\cdot {10}^{-5}\\ 0.023,w_s\ge 3\cdot {10}^{-4}\\ 0.094+0.02{\log }_{10}\left(w_s\right),5\cdot {10}^{-5}{w}_{s}3\cdot {10}^{-4}\end{array}$$ E53

3.2.3. Erosion

The sediment that enters suspension during erosion of the seabed is typically quantified using the erosion rate, which is the eroded mass per unit surface and unit time [ΜL^-2Τ^-1]. For homogenous, consolidated beds, the well known Partheniades formula applies [46]:

$$\epsilon ={\epsilon }_o\left(\frac{{\tau }_b}{{\tau }_{cr,er}}-1\right)$$ E54

The critical erosion shear stress (τ_cr,er) and the erosion rate (ε_ο) constant are considered to be constant with time and within the sediment. However, resistance to erosion depends on various parameters, among which sediment composition, porosity and degree of consolidation of the deposit [47]. The seabed may be soft, partially consolidated, with high water content (>100%), or may be a denser, more stable deposit. The way in which erosion takes place also varies with the range of the applied shear stress and the stress history of the deposit. Three erosion processes have been identified [14]:

1. floc-to-floc erosion, or surface erosion of the bed,

2. high concentration fluid mud entrainment and

3. mass erosion of the bed.

A soft bed that usually consists of freshly deposited clayey matter that is under consolidation presents non-uniform characteristics and is therefore vertically stratified. The erosion rate in this case is [14]:

$$\epsilon ={\epsilon }_f\cdot \exp \left[\alpha {\left({\tau }_b-{\tau }_{cr,er}\left(z\right)\right)}^{\beta }\right]$$ E55

In equation 55, ε_f is the floc erosion constant [kg/m²/s] (equals the erosion rate in the case that the shear stress equals the threshold value), τ_cr,er(z) is the erosion resistance profile with depth z [Pa] and β (β~0.5) and α [Pa^-β] are constants. The values of ε_f (10^-4-10^-7 gcm^-2min^-1) and α (5-20mN^-1/2) depend on the sediment-water mixture.

It is clear that one of the most significant parameters in modelling erosion of cohesive bed is the determination of the critical shear threshold. Various experimentally determined formulas have been proposed, among which many relate the erosion threshold to the density of the bed through relationships of the form:

$${\tau }_{cr,er}=\xi \cdot {\rho }_d{}^{\zeta }$$ E56

There is a high variability in the values of the constants involved. We indicatively mention that [48] estimated them at ξ=5.42 10^-6 and ζ=2.28, while [49] defined ξ=1.2 10^-3 and ζ=1.2. It is noted that, alternatively to the use of the sediment dry density, bulk density is also used in self-similar expressions. From insitu measurements in a macrotidal mudflat [50] a positive correlation of the critical erosion shear stress was found with the sediment water content (W), the bulk density (ρ_b) and the mass loss on ignition (LOI):

$${\tau }_{cr,er}=\{\begin{array}{c}-0.0012W+0.51\\ 0.0013{\rho }_b-1.4\\ -0.041LOI+0.65\end{array}$$ E57

It follows that the erosion rate constant and the critical shear threshold are decisive parameters for the quantification of eroding masses. These parameters are usually experimentally determined; however the values proposed in literature present high variability. The main reason for this variability is the dependence of the erosion process on various factors, like:

• Density: the resistance to erosion increases with the density of the bed.

• Floc diameter – fine fraction content: the effect of the percentage of fines in the sediment to the critical erosion stress is strongly non-linear.

• Clay composition: the lithological composition of clay affects the behaviour of the sediment under shear.

• Stress history: the consolidation degree of the sediment is a significant parameter since clayey sediments change their mechanical behaviour in cases of preloading and overconsolidation.

• Organic content: experimental measurements have shown that the erosion threshold is significantly lower for the case of organic-rich beds.

• Biostabilisation: bacteria and benthic diatoms increase the cohesion of the bed through secretion of extracellular polymeric substances, while the existence of macrophytes in the surface of the bed prohibits erosion reducing the flow velocities through their canopy and at the same time stabilizing the sediment through their root system.

• Bioturbation: the existence of micro-fauna in the surface layer of the sediment reduces the resistance to erosion through loosening of the sediment structure during grazing and the production of faecal pellets.

• Measurement method (laboratory or insitu): the measurement method itself introduces uncertainties regarding the accuracy and the general applicability of the values, given that there is high variability for similar areas using different instrumentation. The differences are even higher between insitu and laboratory experiments, with the former to provide significantly higher erodiblity parameters than the latter.

3.2.4. Self-weight consolidation – Resuspension

A fresh mud deposit has the form of a sediment-water mixture with sediment concentrations of the order of a few tens of grams per liter. This mixture evolves in three stages [12,42]:

1. Deposition of aggregates during the first hours forming fluid mud

2. Depletion of the interstitial water in a period of one to two days

3. Depletion of the intra-particle pore water (very slow consolidation)

Typical consolidation time-scales for a deposit with bulk density of 1150-1250kg/m³ are of the order of 1 week, 1 month, 3 months and 0.5-1year for a layer thickness of 0.25m, 0.5m, 1m and 2m, respectively [42].

The resistance of the deposit to resuspension increases with consolidation time, process that can be parameterised considering a simplified vertical structure of the bed. In its simplest form, only two layers are taken into account, a surface, non-consolidated, low-strength layer and a deeper, fully consolidated one; different densities and erodibility parameters are assigned to each layer. This approach can be extended to a polystromatic representation of the vertical composition of the bed and of the consolidation process.

One of the most widely applied parameterizations for the process of consolidation is the Gibson equation [51]:

$$\frac{\partial \varphi }{\partial t}-\frac{{\rho }_s-{\rho }_w}{{\rho }_w}\frac{\partial }{\partial z}\left(k{\varphi }^2\right)-\frac{1}{{\rho }_{w}g}\frac{\partial }{\partial z}\left(k\varphi \frac{\partial {\sigma }^{\prime }}{\partial z}\right)=0$$ E58

Equation 58 is written in Eulerian form and using the solid volume concentration [51], while k and σ’ are the permeability and the effective stress, respectively. Under the reasoning that the critical failure stress for the sediment comes as a result of failure of the bonds between and within the aggregates, a Mohr-Coulomb type criterion applies:

$${\tau }_{cr,er}=c^{\prime }+\tan \phi {\sigma }^{\prime }$$ E59

The first term of equation 59 expresses the effect of cohesion (c’), which is the shear strength in zero effective stress, and φ is the friction angle.

Exponential functions that relate the compactibility of the sediment through parameters like the void ratio, e, to effective stress are also used [52]:

$$e=\left(e_o-e_{\infty }\right)\exp \left(-\lambda {\sigma }^{\prime }\right)+e_{\infty }$$ E60

Another empirical relationship [53] between the critical erosion shear velocity and the yield strength of muddy deposits reads:

$$\left\{\begin{array}{l}{u}_{*cr,er}=0.013{\tau }_y{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.},{\tau }_y=b{\rho }_d{}^{3}3N/m^2\\ u_{*cr,er}=0.009{\tau }_y{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.},{\tau }_y=b{\rho }_d{}^{6}3N/m^2\end{array}\right\}$$ E61

The first of the equations refers to fluid mud and the second to mud that has exceeded the plasticity limit. Accepting that the critical shear stress threshold for resuspension (τ_cr,res) of deposited sediments ranges between the corresponding values for deposition (τ_cr,dep) and erosion (τ_cr,er), at deposition (t_d=0) and after full consolidation (t_d=t_fc), respectively, its value can be estimated through exponential equations of the form [38]:

$${\tau }_{cr,res}={\tau }_{cr,dep}+\left({\tau }_{cr,er}-{\tau }_{cr,dep}\right)\cdot \left[1-e^{-n\cdot t_d}\right],n\simeq 7/t_{fc}[T^{-1}]$$ E62

3.3. The fine sediment transport model

Based on the preceding analysis, the formulation of a sediment transport module involves the mathematical description of all the aforementioned pelagic and benthic processes and the interactions between them. A typical outline of such a module is presented in figure 2, showing the transition of the sedimentary plume to various conditions, the governing processes at each stage and the interactions between the processes. These considerations were taken into account in formulating the Fine Sediment Transport Model (FSTM); the model is outlined in the present work and characteristic results from simulations in Thermaikos Gulf (NW Aegean Sea) with point (rivers, mechanical erosion) and distributed (atmosphere) sources are presented and discussed.

Figure 2.

Schematic representation of a cohesive sediment transport module. Dashed vectors denote transition from one state (ovals) to another and solid vectors indicate effect of a process (rectangles) to another.

FSTM [38] was formulated using the particle-tracking method (LM), thus describing sediments in the marine environment as particles of specific mass and characteristic properties (assigned to each particle in the form of indicators) that are being passively advected and dispersed by the currents. FSTM accounts for all the pelagic and benthic processes that were described in the previous section. More specifically, the position of the particle in the next temporal step is calculated from eq. [3-4], using the dispersion coefficients for the amplitude of the stochastic particle displacement defined in eq. [10-11]. The evolution of the cohesive particle diameter is defined using the fractal method (eq. 19-20), while the changes in the bulk density of the aggregate are calculated from equation 18. These parameters are used to determine the settling velocity of the particle through the modified Stokes law (eq. 23). The effects of stratification to the vertical propagation of the sedimentary plume are taken into account using the double diffusion theory and assigning damping functions (eq. 40-41) to the vertical deterministic and stochastic displacements. The constants γ, δ and κ were estimated at 0.18, -0.9 and 0.03 [38]. Regarding near-bed processes, the benthic shear stress is defined explicitly, from the horizontal velocity profiles (eq. 46), and the critical shear velocity for deposition of a particle is determined from its settling velocity, using equation 53. After the deposition of a particle onto the bed, the critical shear stress for resuspension of the particle is considered to increase with depositional time (eq. 62). After full consolidation (t_d=t_fc) the particle is considered to be part of the seabed, which is assumed homogenous and fully consolidated (eq. 54). Regarding the boundary conditions of FSTM, radiation is applied at the open boundaries of the domain (particles are allowed to cross the boundary and are excluded from following computational cycles) and reflection is applied at the surface (particles are not allowed to escape to the atmosphere), at the bed (in cases that a particle does not deposit it is reflected to its previous vertical position) and at ‘dry’ grid cells (along one or both horizontal directions). The sources of particles may be the topmost grid cell in the vertical sense, in cases of river-borne or aeolian transported particles, or the bottom cell (for physically or mechanically eroded matter).

4. Characteristic applications

4.1. The application domain: Thermaikos Gulf

Thermaikos Gulf (Figure 3), a micro-tidal, elongated shelf in the NW Aegean Sea, is an area of high socio-economic and environmental significance. Numerous anthropogenic activities take place in the vicinity of the gulf, forcing the marine environment with residues from agricultural, urban and industrial residues. Significant marine transport loads also exist, since the port of Thessaloniki is the second largest harbor in Greece and a gateway to the Balkans. Due to its trophic state (mesotrophic, compared to the oligotrophic open Aegean Sea) and its mild-sloped bottom relief, Thermaikos is also one of the key marine sites for aquaculture and fishing (trawling and otherwise) of the Hellenic region. Sources of freshwater and fine sediments by rivers and smaller streams are distributed along the greater part of its northern and western coastlines; four are the main rivers discharging in the area, Axios, Pinios, Aliakmon and Loudias (Figure 3), of which the former two are the largest. It should be noted that the deltaic area of Axios, Aliakmon and Loudias is protected by the Ramsar convention (Wetlands of International Importance) and is included in the EU network of protected areas Natura 2000 (GR1220010).

FSTM was applied in the domain of Thermaikos, covering the aquatic area that extends southwards to longitude 39.6°N and eastwards to latitude 23.5°E (Figure 3). The grid is horizontally curvilinear (dx=dy=1/60°) and vertically Cartesian (dz=2m). The necessary input hydrodynamic and physical parameters (seawater velocities, salinity and temperature fields) were derived by the North Aegean Sea (NAS) model [54]. Time-series of sediment fluxes from the 4 main rivers, used in the simulations, were recorded during the METRO-MED project [55]. Results from three applications are presented further down, for river-borne sediments, aeolian transported dust and mechanically (trawling) eroded matter. It is noted that the hydrodynamics for the former case are after results of NAS with climatological forcing, thus expressing the typical circulation features in the domain.

Figure 3.

Bathymetry and geomorphological parts Thermaikos Gulf and its location in the Mediterranean Sea

4.2. Point source (river-borne sediment) simulation

4.2.1. Transport patterns and particulate matter distribution

One of the major advantages of formulating FSTM using the EL method is the ability to track the history of each particle and to segregate sediment in the domain according to common properties. An example of visualizing such information is given in figure 4; the figure shows the evolution of suspended and surface (top 2m) particles, 30, 120, 210 and 300 days from the beginning of the simulation (corresponding to late January, April, July and October, respectively). The chromatic encoding denotes the riverine origin of the particles. These representations give information regarding the seasonal variability and the transport patterns of sediments by each of the 4 rivers.

It can be noted that, during winter (Figure 4b), Suspended Particulate Matter (SPM) enter the enclosed inner (northmost) part of the gulf, matter that mainly originates from Axios. At the same time, stratification is strong in the northern river-dominated area, leading to high presence of SPM in the surface layer (Figure 4a). The northward expansion of sediment from Axios is probably due to anticyclonic movement of the surficial waters in the gulf during the wintry months [56]. This movement is also responsible for the northward expansion of the Pinios plume, ascertainment that is in accordance with [57]. Matter from Aliakmon is mainly transported towards the south along the western coastline, together with part of the sediment supply of Axios. During the vernal period (Figure 4d) a distinct cyclonic transport pattern is detected in the outer part that forces sediments to move southwards along the west coast, while SPM from Axios are still present in the inner gulf. The same transport pattern is observed for surficial sediments (Figure 4c) with a strong effect of freshwater plumes along the western coastline, effect that keeps sediments in the surface layer at high distances from the outflows. The dispersion of particles is high near the northerly river-outflow system, forcing particles (mainly from Axios) to move eastwards in addition to the prevailing cyclonic movement. The presence of SPM during summer is low (Figure 4f) due to reduced river-borne sediment fluxes; suspended particles appear mainly in the northern part of the outer gulf, with low SPM in the surface layer (Figure 4e). During autumn (Figure 4h), the increase of sediment supplies causes related increase to SPM and a predominant cyclonic transport pattern exists in the northern part of the gulf; a small part of the Axios sediment supply, again, spreads eastwards and towards the inner gulf. Stratification is strong not only for the northern part, but also near the outfalls of Pinios (Figure 4g), with a significant spreading of the surficial plume, which is the most profound of the whole year.

Figure 4.

Variation of surficial and suspended particles at the 30^th (mid-winter) (a, b), 120^th (mid-spring) (c, d), 210^th (mid-summer) (e, f) and 300^th day of the simulation (mid-autumn) (g, h). The chromatic encoding denotes the different riverine origin of the particles.

4.2.2. The contribution of the rivers to the sedimentation of the gulf

The existence of transport and sedimentation patterns related to the riverine origin of SPM and the geomoprhological parts of the gulf becomes evident. The following table (Table 1) presents the percentages of deposition in each region for each of the four rivers, to facilitate the investigation of the spatial distribution of the sediments with respect to its riverine origin.

Part of Thermaikos	River
Axios	Loudias	Aliakmon	Pinios
Inner part	28.2%	19.8%	5.8%	0.0%
Outer part	43.5%	44.7%	61.5%	2.4%
Extended part	28.3%	35.5%	32.7%	97.6%

Table 1.

Sedimentation percentages of the rivers in each part of the gulf

Figure 5.

Simulated sedimentation rate [mm/y] (a) and sediment types after [58] (b).

Sediments from Aliakmon present a constant, throughout the year, transport pattern, with the prevalence of advection against dispersion, and predominant movement along the low-depth areas of the western coast. The majority of the sediment supply of Aliakmon deposits in the outer gulf, while a very small part of the supply reaches the enclosed inner Thermaikos. Contrastingly, Axios appears as the main sediment supplier for the inner gulf, along with a much smaller contribution from the Loudias. Matter from Axios presents a great dispersal in the outer part of the gulf, accumulating sediments along the eastern and western coasts; nearly half of the sediment fluxes of Axios settle in the outer gulf (Table 1) while the remaining half equally supplies the coastal areas of the inner and extended Thermaikos. Similar transport and sedimentation patterns are exhibited by Loudias, mainly due to the proximity of the two river outflows. Particles from Pinios supply the coastal zone along the western coastline, primarily northwards and secondarily southwards from the outflow, remaining almost exclusively in the extended part of the domain. The northward deflection of the Pinios plume is due to the strong bottom slope near the mouth of the river [57] and to the effect of the Sporades eddy, which frequently turns anti-cyclonic [56].

The simulated sedimentation rates (Figure 5a) are in good correlation with corresponding values defined after core-sampling in the area [58]. Furthermore, the areas of high river-borne sediment accumulation coincide with the locations in the gulf where the sediment ranges from mud to sandy mud (Figure 5b), validating the modelling approach. In the coastal zone of the eastern outer Thermaikos, where the benthic material grades from sandy mud to muddy sand, the simulated depositing trends are reducing towards the south.

4.3. Distributed source (aeolian transported dust) simulation

Results of the SKIRON/Eta forecasting system and the Eta/NCEP regional atmospheric model [58] were used for the determination of dust inflow in the domain. The model provides prognoses for total dust deposition over the Mediterranean Sea, parameter that was used after two-linear interpolation to transfer data from the domain of the Mediterranean to the domain of Thermaikos.

The time-series of the total dust mass deposited onto the aquatic surface for the period 13/04-13/06/2005 is depicted in figure 5 b. This specific period was selected on account of an intense dust storm that occurred on April 17, 2005. During this incident thick bands of Saharan dust covered the whole Hellenic area, as depicted in the satellite image of figure 5 a. It can be noted that the dust storm of April 17^th is predicted by the atmospheric model (Figure 5b) and that three more episodes of lower intensity followed in May (7, 12 and 21). The total mass introduced to the domain during the period in question is 0.72 10^-6t, value that exceeds the corresponding annual contribution of the rivers of the area, estimated at 0.65 10^-6t/y [55]. It is noted that the simulation was extended by 1 month (up to 13/07).

4.3.1. Dust distribution

Investigating the temporal evolution of total suspended and deposited dust masses (Figure 6a) it can be noted that the peak-values of the suspended dust masses, as expected, follow the related peaks of dust inflow (Figure 5b), while deposition increases with time. Maximum total suspended dust masses are recorded in mid-May, while almost all of the matter has settled on the bed at the end of the simulation. The average suspended concentrations at the surface (Figure 6b) follow the morphology of the input time-series, which is expected since the surface is the source of the particles, with values ranging from 8 to 16mg/l at peak ‘loading’. The depth-averaged concentrations (Figure 6b) show a small temporal hysteresis in peak values compared to the inflow time-series, of the order of 1-2 days, and range between 0.05 and 1.5mg/l. The maximum suspended concentrations in the field are of the order of 21mg/l and 50mg/l for the water column and the surface layer. The average suspension times of dust in the domain were estimated at 17 days. At the same time, the residence time of dust in deposition increases almost exponentially with time, indicating that deposition prevails over resuspension; near the end of the simulation, the majority of particles settled mainly in the low-depth areas along the coastlines.

Figure 6.

MODIS satellite imagery of the dust storm on April 17, 2005 (a) (the image is a property of NASA) and temporal evolution of the total deposited mass in the domain during the period 13/04-13/06/2005 (b).

Figure 7.

Temporal evolution of total suspended and deposited mass (a) and of suspended depth-averaged and absolute concentrations in the surface layer (b) and horizontal trajectories of randomly selected dust particles (c): the location of entrance in the field is denoted as a dot.

4.3.2. Dust particle trajectories

As in any simulation using a LM, it is possible to track the movement and changes in characteristics of sediments with time using FSTM. Such an example is given in figure 6c as horizontal trajectories of randomly selected dust particles; the location of entrance of the particle in the surface layer of the domain is noted with a dot. Some of these particles escaped the domain through the open boundaries (orange and mauve curves), while others settled on the bed after performing particularly arbitrary movements and presenting very different suspension times. Emphasising on the particle with the maximum residence time in the field (blue curve in figure 6c), it can be noted that it entered the field on 19/04(18:00) and initially moved cyclonically for 4 days, period during which its characteristics remained stable. Following, the particle entered an anti-cyclonic motion for 5 days, at the end of which 28/04(17:30) it started to move cyclonically; this motion lasted for 5 days, during which the particle performed a full circle. Afterward, the particle transported westwards for 3.3 days, period after which its movement became particularly arbitrary. In the vertical sense (data not shown) the particle retained an almost stable settling rate, with an average value of 50mm/hr. The particle finally settled onto the bed 38.5 days after its entrance in the field.

4.4. Mechanically eroded matter

FSTM was also applied for the investigation of mechanically eroded (trawled) masses in Thermaikos, application whose main findings are described briefly in the present chapter; the reader is referred to [60] for further information. The mechanical erosion rates were estimated at 430gr/m², value directly comparable to similar coastal areas, while the sedimentary matter mobilized due to the activity of the benthic trawlers during one trawling season was found significantly higher than the respective contribution of the rivers in the gulf. The sediments, after mobilization by the trawls, continued to move at low distances from the bed, forming Benthic Nefeloid Layers (BNLs), finding that supports the claim that part of the BNLs found in Thermaikos originates from the trawling activity in the area. The redistribution of sediments due to trawling was found generally low, with the matter to redeposit at a small distance from its generation, typically lower than 5km and with average suspension times of 1 to 5 days. The mass export rates from the gulf were also found to be small, especially for the initial part of the trawling period.

5. Conclusions

The present work concentrates on the physical interpretation and modelling approaches used for the main processes that drive sediment transport on coastal and shelf seas. The analysis focused on typical modelling methods and on the physical processes that take place in the water column and near the bed, including the corresponding mass-exchanges. More specifically, advection-dispersion, flocculation and settling mechanisms were analyzed as the processes that take place in the column, along with the effects of stratification to the vertical movement of fine sediments. Regarding near-bed processes, the aspects addressed include shear stress conditions, particle deposition and gradual consolidation and erosion of the bed. In order to accurately describe the fate of cohesive sediments in the marine environment, a numerical tool must generally include all of these processes. However, there are cases in which some of the processes have minimal or no effect; for example, hindered settling does not partake for low concentration environments (<20mg/l) and can therefore be ignored, while stratification becomes significant only in the presence of strong density gradients, typical in estuarine areas. Likewise, wind-generated wave impacts on bottom shear become significant in low-depth (<10m), coastal areas and need to be accounted for only in high resolution regional scale models.

The FSTM model, formulated based on the particle-tracking method and the considerations analyzed in the chapter, was briefly described and indicative results from applications in Thermaikos Gulf, a microtidal shelf sea of the Northern Aegean Sea, regarding point and distributed sedimentary sources were presented and discussed; transport and sedimentation patterns in the gulf were identified for the river-borne sediments of the area. The prevailing cyclonic circulation of the area appears as the main mechanism controlling sedimentation pathways in the domain and forcing higher sediment accumulation along the shallow areas of the western coast. Sediment from Aliakmon is the most characteristic example of this behaviour, with consistent, throughout the year, cyclonic movement and deposition; the greater part of the Aliakmon sediment influx (~⅔) accumulates in the outer gulf, while a significant amount (~⅓) reaches the southern part of the shelf. The geomorphology of Thermaikos also plays a significant part in the phenomena; the location of the Axios delta to the centre of the northern coastline exposes its sediments to stronger currents of alternating directions and, thus, its sediments present much higher dispersion in the field, compared to the other rivers of the area. This behaviour is enhanced by seasonal stratification, which forces sediment to remain suspended longer at periods of high freshwater supplies from the rivers. Nearly half of the sediment from Axios settles in the outer part of the gulf, along the shallow areas of both the western and eastern coasts. The remaining part of its outflow coequally replenishes the sediment of the extended and inner parts. In fact, the mathematic investigation showed that Axios is the main sediment supplier of the enclosed inner Thermaikos. The steep bottom relief near the outflows of Pinios, combined with the hydrodynamics of the area, forces its particles to settle along the western coast close to the delta, with a higher northward spreading. The mass export towards the deep Sporades basin was estimated by the simulation at approximately 10% of the total riverine sediment supply, which corresponds to around 85 10³t/yr; this high sediment accumulation in the shelf is consistent with the characterisation of Thermaikos as a sediment trap.

Results from the three-month simulation of aeolian transported dust during the dust storm reported in mid-April of 2005, showed concentrations of suspended dust between 8 and 16mg/l at maximal loading and corresponding depth-averaged values ranging from 0.05 to 1.5mg/l. The mean suspension times of dust in the water column were around 17 days, while dust deposition prevailed over resuspension. These high concentrations and suspension times reveal the significance of the process of aeolian transport to the suspended masses in the gulf that can exceed the contribution of river-borne matter considerably in cases of dust storms.