*T T*

done by using the following equation:

For cases where flow is in one direction for example when *v* = 0, transport is given by the

With the aid of Equation (8)–(9), accordingly the determination of bed level changes is now

*ρs*(1 − *po*) · div(*U*3).

 *∂***um** *∂x* + *∂***vm** *∂y*

. (10)

(11)

0 (*U*2 + *V*2)*Vdt*. Next, let us now discuss the

*dt* + √2*DdW*1(*t*),

*dt* + √2*DdW*2(*t*),

*∂h*

∆*h* ≈ −*fdT*

**4. A particle model for sediment transport in shallow waters**

0 (*U*2 + *V*2)*Udt* and **vm** = 1

Lagrangian particle model in the following section.

drift part and a stochastic(diffusive) part [4, 13].

*dX*(*t*) Itoˆ = *U* + *D H ∂H ∂x* + *∂D ∂x*

*dY*(*t*) Itoˆ = *V* + *D H ∂H ∂y* + *∂D*) *∂y*

transport is initiated with zero number of particles.

S*t* =

**4.2. Deposition of sediment particles**

**4.1. Integration of particle movement**

*ρs*(1 − *po*) ·

*T T*

A particle model is a description of a transport process by means of random walk models. Random walk model is defined as the stochastic differential equation that describes the movement of a particle that subsequently undergoes a displacement, which consists of the

In this section, the following 2-dimensional stochastic differential equations is developed:

where the Brownian process *W*1(*t*) and *W*2(*t*) are Gaussian [3], and *D*(*x*, *y*, *t*) is the horizontal dispersion coefficient for sediment transport. Typically, *D* = O(10− 100)*m*2/*s* [14]. Note that *U*(*x*, *y*) and *V*(*x*, *y*) are the flow velocities along the *x* and *y* direction respectively given in *m*/*s*, *H*(*x*, *y*) is the averaged total depth plus relative water levels due to waves, *dW*1(*t*) and *dW*2(*t*) are independent increments of Brownian motions with mean (0, 0)*T* and covariance [*dW*1(*t*)*dW*2(*t*)T] = *Idt* where *I* is an identity matrix ([3, 10]). The simulation of sediment

We associate with each sediment particle a binary state which at any time *t* is given by

0 otherwise (particle is on the sea bed).

1 particle is in suspension

*∂t* = − *fd*

Mass represents concentration of a group of particles at a certain location. A source term is included in our particle model such that the expected number of suspended particles (*enp*) in grid cell *i*, *j* at time *t* is given by

$$emp\_{(i,j,t)} = \frac{\Delta x \cdot \Delta y \cdot \Delta t \cdot \left(\mathcal{U}^2 + V^2\right) \cdot \lambda\_s}{\mathcal{M}\_p} \, , \tag{15}$$

where M*p* is the mass of each particle, ∆*x* and ∆*y* are the width of the grid cells along *x* and *y* directions respectively, ∆*t* is the time step size, and *λs* is the erosion coefficient. For each grid cell *i*, *j* we use the expected number of particles from Equation (15) to determine the actual number of particles to be suspended. This is done by drawing a number from a Poisson distribution function.

## **5. The relationship between the Kolmogorov Forward Partial differential and The Advectiion Diffusion transport model**

In this section the relationship between the Kolmogorov Forward Partial differential corresponding to the 2-dimensional Stochastic differential equations and the well known two dimensional advection diffusion equations for sediment transport is discussed in detail. Since we are interested in the particle being in suspension, we assume that the particle at position (*x*, *y*) at time *t* has expectation of their mass �·� defined by;

$$
\langle \langle m(\mathbf{x}, y, t) \rangle \rangle = f(\mathbf{x}, y, t) \cdot \mathbf{P}(\mathbf{S}\_t = \mathbf{1}).\tag{16}
$$

This is known as mass density of particles per unit area of the grid box. The Kolmogorov Forward Partial differential is known as the Fokker-Planck equation(FPE)[3]. In this section we incorporate the processes such as suspension and sedimentation states of the particles in the model for transport of sediments in shallow water. We let *D* be the diffusion coefficient, P(*St* = 1), is the probability that particle is in suspension, (*St* = 1) denotes a state that a particle is in suspension and (*St* = 0) denotes the state that the particle is deposited on the sea bed or bed of the shallow water. The stochastic process (*Xt*,*Yt*) is a Markov process. The probability density function of the particle position *f*(*x*, *y*, *t*) based on the two dimensional Stochastic Differential equations (SDEs) (11) evolves according to the following Fokker-Planck equation [9].

10.5772/52237

317

http://dx.doi.org/10.5772/52237

The substitution of Equation (20) into equation (19) gives the following Fokker-Planck equation which now have extra terms that model the deposition and erosion of sediments.

· �*m*(*x*, *y*, *t*)�

Modelling of Sediment Transport in Shallow Waters by Stochastic and Partial Differential Equations

1 2 *∂*2

The last two terms in the transport equation (21) represents the process of sediment deposition and erosion of sediments respectively. An average mass �*m*(*x*, *y*, *t*)� per unit area (*kg*/*m*2) of a particle is related to a particle depth averaged concentration *C*(*x*, *y*, *t*) in mass

The concentration of materials is given in *kg*/*m*3, therefore, the expected mass of a particle at position (*x*, *y*) can be related by the concentration in that grid cell(location) as follows:

Substitution of equation (22) into the Fokker-Planck equation (21) leads to the two dimensional Advection diffusion partial differential equation commonly known as Eulerian sediment transport model (4). Consequently, the transport equation (4) which was discussed

Therefore, after having constructed the particle model for sediment transport, it is now necessary to develop the equations that cater for the bed level changes using the particle

In the next section we shall briefly discuss the numerical approximation of our particle model.

Euler scheme is used in the numerical implementation of the particle model. The scheme is convergent in the weak sense with accuracy of order O(∆*t*) and it is half order accurate in the strong sense. Higher order schemes for stochastic differential equations are described in [11]. The discretisations of the hydrodynamic flow models is widely discussed by [15], for example. The particle model (11)-(26) is discretised and uses the following Euler scheme to approximate the numerical solutions. We discretise the two dimensional stochastic differential equations for integrating the movement of the particle in similar way to that as in [9] with the modifications by the inclusion of the sedimentation and deposition parts:

*D*)/*H* +

*∂D ∂x*  ∆*tk* + √

2*D*∆*W*1(*tk*) (23)

in Section 3, is consistent with the particle model for sediment transport (11)-(26).

*U* + ( *∂H ∂x*

**6. Numerical approximation of the particle model**

*X*¯(*tk*+1) = *X*¯(*tk*) +

· �*m*(*x*, *y*, *t*)�

*∂y*2 (2*D* · �*m*(*x*, *y*, *t*)�)

�*m*(*x*, *y*, *t*)� = *H*(*x*, *y*, *t*) · *C*(*x*, *y*, *t*). (22)

− *γ* · �*m*(*x*, *y*, *t*)� + *λs* · E(*U*, *V*). (21)

*∂*�*m*(*x*, *y*, *t*)�

model.

*∂t* = − *∂*

*∂x*

+ 1 2 *∂*2

per unit volume (*kg*/*m*3) ([9, 18]).

*U* +

− *∂ ∂y* *D H ∂H ∂x* + *∂D ∂x*

*V* +

*∂x*2 (2*D* · �*m*(*x*, *y*, *t*)�) +

*D H ∂H ∂y* + *∂D ∂y*

$$\begin{split} \frac{\partial f(\mathbf{x},y,t)}{\partial t} &= -\frac{\partial}{\partial \mathbf{x}} \left( \left[ \mathbf{U} + (\frac{\partial H}{\partial \mathbf{x}} \mathbf{D}) / H + \frac{\partial \mathbf{D}}{\partial \mathbf{x}} \right] \cdot f(\mathbf{x},y,t) \right) \\ &- \frac{\partial}{\partial y} \left( \left[ V + (\frac{\partial H}{\partial y} \mathbf{D}) / H + \frac{\partial \mathbf{D}}{\partial y} \right] \cdot f(\mathbf{x},y,t) \right) \\ &+ \frac{1}{2} \frac{\partial^2}{\partial \mathbf{x}^2} \left( f(\mathbf{x},y,t) \cdot \mathbf{2}D \right) + \frac{1}{2} \frac{\partial^2}{\partial y^2} \left( f(\mathbf{x},y,t) \cdot \mathbf{2}D \right) . \end{split} \tag{17}$$

The resulting sediment transport model (11) is an extension of two dimensional particle model for pollutant dispersion in the shallow waters developed by [9]. The extension in the present model includes the erosion and deposition terms. It is possible to derive the Fokker-Planck equation that describes the probabilistic transport of the sediments from one location to another. The derivation of the Fokker-Planck equation is e done as follows, we first differentiate equation (16) with respect to time *t*, to obtain

$$\frac{\partial}{\partial t} \langle m(\mathbf{x}, y, t) \rangle = \mathbf{P}(\mathcal{S}\_t = 1) \frac{\partial}{\partial t} f(\mathbf{x}, y, t) + f(\mathbf{x}, y, t) \frac{\partial}{\partial t} \mathbf{P}(\mathcal{S}\_t = 1),$$

next with the aid of Equation (12), it follows that,

$$\frac{\partial}{\partial t} \langle m(\mathbf{x}, y, t) \rangle = \mathbf{P}(\mathbf{S}\_t = 1) \frac{\partial}{\partial t} f(\mathbf{x}, y, t) - \gamma f(\mathbf{x}, y, t) \cdot \mathbf{P}(\mathbf{S}\_t = 1). \tag{18}$$

Therefore, we also add the erosion term to Equation (18) and come up with

$$\frac{\partial}{\partial t} \langle m(\mathbf{x}, y, t) \rangle = \mathbf{P}(\mathbf{S}\_{\mathbf{f}} = 1) \frac{\partial}{\partial t} f(\mathbf{x}, y, t) - \gamma f(\mathbf{x}, y, t) \cdot \mathbf{P}(\mathbf{S}\_{\mathbf{f}} = 1) + \lambda\_{\mathbf{s}} \cdot \mathcal{E}(\mathbf{U}, V) \tag{19}$$

next we multiply on both sides of Equation (17) by P(*St* = 1) to obtain

$$\begin{split} \mathbb{P}(\mathcal{S}\_{l} = 1) \frac{\partial}{\partial t} f(\mathbf{x}, y, t) &= -\frac{\partial}{\partial \mathbf{x}} \left( \left[ \mathcal{U} + \frac{D}{H} \frac{\partial H}{\partial \mathbf{x}} + \frac{\partial D}{\partial \mathbf{x}} \right] f(\mathbf{x}, y, t) \cdot \mathbb{P}(\mathcal{S}\_{l} = 1) \right) \\ &- \frac{\partial}{\partial y} \left( \left[ \mathcal{V} + \frac{D}{H} \frac{\partial H}{\partial \mathbf{y}} + \frac{\partial D}{\partial y} \right] f(\mathbf{x}, y, t) \cdot \mathbb{P}(\mathcal{S}\_{l} = 1) \right) \\ &+ \frac{1}{2} \frac{\partial^{2}}{\partial \mathbf{x}^{2}} \left( 2D \cdot f(\mathbf{x}, y, t) \mathbb{P}(\mathcal{S}\_{l} = 1) \right) + \frac{1}{2} \frac{\partial^{2}}{\partial y^{2}} \left( 2D \cdot f(\mathbf{x}, y, t) \cdot \mathbb{P}(\mathcal{S}\_{l} = 1) \right) . \end{split} \tag{20}$$

The substitution of Equation (20) into equation (19) gives the following Fokker-Planck equation which now have extra terms that model the deposition and erosion of sediments.

8 Sediment Transport

Fokker-Planck equation [9].

*∂t* = − *∂*

*∂x*

*∂ ∂t*

*∂ ∂t*

*∂ ∂t*

P(*St* = 1) *∂*

*∂t*

*U* + ( *∂H ∂x*

first differentiate equation (16) with respect to time *t*, to obtain

�*m*(*x*, *y*, *t*)� = P(*St* = 1) *∂*

�*m*(*x*, *y*, *t*)� = P(*St* = 1) *∂*

*∂x U* + *D H ∂H ∂x* + *∂D ∂x*

next with the aid of Equation (12), it follows that,

�*m*(*x*, *y*, *t*)� = P(*St* = 1) *∂*

*f*(*x*, *y*, *t*) = − *∂*

+ 1 2 *∂*2

− *∂ ∂y V* + *D H ∂H ∂y* + *∂D ∂y*

− *∂ ∂y*

*D*)/*H* +

*V* + ( *∂H ∂y*

> + 1 2 *∂*2

*∂D ∂x*

*D*)/*H* +

The resulting sediment transport model (11) is an extension of two dimensional particle model for pollutant dispersion in the shallow waters developed by [9]. The extension in the present model includes the erosion and deposition terms. It is possible to derive the Fokker-Planck equation that describes the probabilistic transport of the sediments from one location to another. The derivation of the Fokker-Planck equation is e done as follows, we

*∂t*

*∂t*

Therefore, we also add the erosion term to Equation (18) and come up with

*∂t*

next we multiply on both sides of Equation (17) by P(*St* = 1) to obtain

*∂x*2 (2*D* · *f*(*x*, *y*, *t*)P(*St* = 1)) +

· *f*(*x*, *y*, *t*)

*∂D ∂y*

*∂x*2 (*f*(*x*, *y*, *t*) · 2*D*) +

*f*(*x*, *y*, *t*) + *f*(*x*, *y*, *t*) *∂*

· *f*(*x*, *y*, *t*)

*∂y*2 (*f*(*x*, *y*, *t*) · 2*D*). (17)

1 2 *∂*2

*∂t*

*f*(*x*, *y*, *t*) − *γ f*(*x*, *y*, *t*) · P(*St* = 1) + *λs* · E(*U*, *V*) (19)

*f*(*x*, *y*, *t*) · P(*St* = 1)

*f*(*x*, *y*, *t*) · P(*St* = 1)

1 2 *∂*2 P(*St* = 1),

*f*(*x*, *y*, *t*) − *γ f*(*x*, *y*, *t*) · P(*St* = 1). (18)

*∂y*2 (2*D* · *f*(*x*, *y*, *t*) · P(*St* = 1)). (20)

*∂ f*(*x*, *y*, *t*)

This is known as mass density of particles per unit area of the grid box. The Kolmogorov Forward Partial differential is known as the Fokker-Planck equation(FPE)[3]. In this section we incorporate the processes such as suspension and sedimentation states of the particles in the model for transport of sediments in shallow water. We let *D* be the diffusion coefficient, P(*St* = 1), is the probability that particle is in suspension, (*St* = 1) denotes a state that a particle is in suspension and (*St* = 0) denotes the state that the particle is deposited on the sea bed or bed of the shallow water. The stochastic process (*Xt*,*Yt*) is a Markov process. The probability density function of the particle position *f*(*x*, *y*, *t*) based on the two dimensional Stochastic Differential equations (SDEs) (11) evolves according to the following

$$\begin{split} \frac{\partial \langle m(\mathbf{x},y,t) \rangle}{\partial t} &= -\frac{\partial}{\partial \mathbf{x}} \left( \left[ \mathcal{U} + \frac{D}{H} \frac{\partial H}{\partial \mathbf{x}} + \frac{\partial D}{\partial \mathbf{x}} \right] \cdot \langle m(\mathbf{x},y,t) \rangle \right) \\ &- \frac{\partial}{\partial y} \left( \left[ V + \frac{D}{H} \frac{\partial H}{\partial y} + \frac{\partial D}{\partial y} \right] \cdot \langle m(\mathbf{x},y,t) \rangle \right) \\ &+ \frac{1}{2} \frac{\partial^2}{\partial \mathbf{x}^2} \left( 2D \cdot \langle m(\mathbf{x},y,t) \rangle \right) + \frac{1}{2} \frac{\partial^2}{\partial y^2} \left( 2D \cdot \langle m(\mathbf{x},y,t) \rangle \right) \\ &- \gamma \cdot \langle m(\mathbf{x},y,t) \rangle + \lambda\_s \cdot \mathcal{E}(\mathcal{U}\_t V). \end{split} \tag{21}$$

The last two terms in the transport equation (21) represents the process of sediment deposition and erosion of sediments respectively. An average mass �*m*(*x*, *y*, *t*)� per unit area (*kg*/*m*2) of a particle is related to a particle depth averaged concentration *C*(*x*, *y*, *t*) in mass per unit volume (*kg*/*m*3) ([9, 18]).

The concentration of materials is given in *kg*/*m*3, therefore, the expected mass of a particle at position (*x*, *y*) can be related by the concentration in that grid cell(location) as follows:

$$
\langle \mathcal{m}(\mathbf{x}, \mathbf{y}, t) \rangle = H(\mathbf{x}, \mathbf{y}, t) \cdot \mathbb{C}(\mathbf{x}, \mathbf{y}, t). \tag{22}
$$

Substitution of equation (22) into the Fokker-Planck equation (21) leads to the two dimensional Advection diffusion partial differential equation commonly known as Eulerian sediment transport model (4). Consequently, the transport equation (4) which was discussed in Section 3, is consistent with the particle model for sediment transport (11)-(26).

Therefore, after having constructed the particle model for sediment transport, it is now necessary to develop the equations that cater for the bed level changes using the particle model.

In the next section we shall briefly discuss the numerical approximation of our particle model.

## **6. Numerical approximation of the particle model**

Euler scheme is used in the numerical implementation of the particle model. The scheme is convergent in the weak sense with accuracy of order O(∆*t*) and it is half order accurate in the strong sense. Higher order schemes for stochastic differential equations are described in [11]. The discretisations of the hydrodynamic flow models is widely discussed by [15], for example. The particle model (11)-(26) is discretised and uses the following Euler scheme to approximate the numerical solutions. We discretise the two dimensional stochastic differential equations for integrating the movement of the particle in similar way to that as in [9] with the modifications by the inclusion of the sedimentation and deposition parts:

$$\dot{X}(t\_{k+1}) = X(t\_k) + \left[ \mathcal{U} + (\frac{\partial H}{\partial \mathbf{x}} D)/H + \frac{\partial D}{\partial \mathbf{x}} \right] \Delta t\_k + \sqrt{2D} \Delta \mathcal{W}\_1(t\_k) \tag{23}$$

$$\bar{Y}(t\_{k+1}) = \bar{Y}(t\_k) + \left[V + (\frac{\partial H}{\partial y}D)/H + \frac{\partial D}{\partial y}\right] \Delta t\_k + \sqrt{2D} \Delta W\_2(t\_k) \tag{24}$$

$$\mathbf{P}\_{k+1}(\mathbf{S}\_t = 1) = (1 - \gamma(\mathbf{x}, y, t)\Delta\_k)\mathbf{P}\_k(\mathbf{S}\_t = 1). \tag{25}$$

10.5772/52237

319

http://dx.doi.org/10.5772/52237

(26)

• If a closed boundary is crossed during the diffusive part of integration, the step size halving procedure described above is maintained with the modification that in addition to the position, the white noise process is also restored to its state prior to the invalidated integration step. The process of halving the time step and continuing integration with the reduced step size is repeated until the full ∆*t* time step has been integrated without

Modelling of Sediment Transport in Shallow Waters by Stochastic and Partial Differential Equations

In addition, this chapter have also considered the pick up of sediment particles at the inflows

We now consider particle flux at open boundaries. This flux is the difference between particles flowing into and out of the domain. The number of particles flowing out should not be controlled as it is a natural consequence of the movement of a particle. As soon as a particle crosses an open boundary it is considered gone and further integration is no longer possible as no data outside the domain is given. For this very reason, however, we do need to explicitly model the particles flowing in. We determine the expected number of particles

where, *γ* is the deposition coefficient. The actual number of particles added at the domain boundary at each iteration is obtained by drawing a value from the Poisson distribution

To determine the actual number of particles to be suspended in a grid cell *i*, *j* we draw a number from a Poisson distribution function with mean *enp*(*i*,*j*,*t*) determined by Equation (15). We assume that particles are infinitely many on the sea bed. However, the

Before implementing the particle model (11)-(26), we first required to show the consistence between the Fokker-Planck equation and its the Eulerian transport model. This will be

Comparing with a simplified form of Equation (4), where in this case we assume for the local

balance (see Section 3.1). The approximation of the change in mass with respect to time, in

1 ∆*x*∆*y*

each grid cell *i*, *j* in the this particle model is determined by the following equations:

*∂t* ≈ *de* − *se*. We also assume that the deposition and erosion processes

M*p*. (27)

*γ*M*p* inflow parallel to y-axis

*γ*M*p* inflow parallel to x-axis,

crossing a boundary.

this will be discussed in the next section.

**6.2. Particle flux at open boundaries**

*enp*(*i*,*j*,*t*) =

parameterised by the above expectation value.

  ∆*y*·∆*t*·*V*·(*U*2+*V*2)·*λs*

∆*x*·∆*t*·*U*·(*U*2+*V*2)·*λs*

particle that is suspended is not the same as the one that is deposited.

**7. Determination of bedlevel changes using particle models**

*∂m ∂t* = ∆N*p* ∆*t*

entering the domain as follows:

described in the next Section.

change in mass *∂m*

Where *X*¯(*tk*+1) and *Y*¯(*tk*+1) are the numerical approximations of *X*(*t*) and *Y*(*t*) respectively, while *X*¯(*t*0) = *X*(*t*0) = *x*0 and *Y*¯(*t*0) = *Y*(*t*0) = *y*0 are initial locations of a particle. In addition to Eqns. (23)-(25), we also use Eqns. (15)-(26) to make the simulation of sediment transport complete. There are several schemes that can be used for simulation process for instance, Euler, Heun, Milstein scheme and Runge kutta methods. Much detailed work on numerical methods for stochastic models can be found in (e,g., [11]).

Numerical schemes such as the Euler scheme often show very poor convergence behaviour. This implies that, in order to get accurate results, small time steps are needed thus requiring much computation. Another problem with the Euler (or any other numerical scheme) is its undesirable behaviour in the vicinity of boundaries; a time step that is too large may result in particles unintentionally crossing boundaries. Therefore, the treatment of boundary condition for particle Models is often done in the following section as follows.

#### **6.1. Boundaries**

One problem with numerical integration of particle positions arises in the vicinity of boundaries. Given the current location, (*X*(*t*),*Y*(*t*)), we may find that the new location, (*X*(*t* + ∆*t*),*Y*(*t* + ∆*t*)), is on the other side of a boundary, i.e. the particle has crossed a boundary. Depending on the type of boundary this may be physically impossible. We consider two types of boundaries. The first type, closed boundaries, represents boundaries intrinsic to the domain such as banks, sea bed, and coastal lines. The second type of boundaries are open boundaries, which arise from the modeller's decision to artificially limit the domain because particles are not expected to reach any further or simply because no domain information is available at those locations. It is clear that it is undesirable to have particles cross the first type of boundary, whereas for the second type it is quite natural. Based on this classification, we apply the following rules to particles crossing borders during integration:

• If a closed boundary is crossed during the diffusive part of integration, the step size halving procedure described above is maintained with the modification that in addition to the position, the white noise process is also restored to its state prior to the invalidated integration step. The process of halving the time step and continuing integration with the reduced step size is repeated until the full ∆*t* time step has been integrated without crossing a boundary.

In addition, this chapter have also considered the pick up of sediment particles at the inflows this will be discussed in the next section.

## **6.2. Particle flux at open boundaries**

10 Sediment Transport

**6.1. Boundaries**

integration:

*Y*¯(*tk*+1) = *Y*¯(*tk*) +

numerical methods for stochastic models can be found in (e,g., [11]).

condition for particle Models is often done in the following section as follows.

*V* + ( *∂H ∂y*

*D*)/*H* +

Where *X*¯(*tk*+1) and *Y*¯(*tk*+1) are the numerical approximations of *X*(*t*) and *Y*(*t*) respectively, while *X*¯(*t*0) = *X*(*t*0) = *x*0 and *Y*¯(*t*0) = *Y*(*t*0) = *y*0 are initial locations of a particle. In addition to Eqns. (23)-(25), we also use Eqns. (15)-(26) to make the simulation of sediment transport complete. There are several schemes that can be used for simulation process for instance, Euler, Heun, Milstein scheme and Runge kutta methods. Much detailed work on

Numerical schemes such as the Euler scheme often show very poor convergence behaviour. This implies that, in order to get accurate results, small time steps are needed thus requiring much computation. Another problem with the Euler (or any other numerical scheme) is its undesirable behaviour in the vicinity of boundaries; a time step that is too large may result in particles unintentionally crossing boundaries. Therefore, the treatment of boundary

One problem with numerical integration of particle positions arises in the vicinity of boundaries. Given the current location, (*X*(*t*),*Y*(*t*)), we may find that the new location, (*X*(*t* + ∆*t*),*Y*(*t* + ∆*t*)), is on the other side of a boundary, i.e. the particle has crossed a boundary. Depending on the type of boundary this may be physically impossible. We consider two types of boundaries. The first type, closed boundaries, represents boundaries intrinsic to the domain such as banks, sea bed, and coastal lines. The second type of boundaries are open boundaries, which arise from the modeller's decision to artificially limit the domain because particles are not expected to reach any further or simply because no domain information is available at those locations. It is clear that it is undesirable to have particles cross the first type of boundary, whereas for the second type it is quite natural. Based on this classification, we apply the following rules to particles crossing borders during

• In case an open boundary is crossed by a particle, the particle remains in the sea but is

• In case a closed boundary is crossed by a particle during the drift step of integration, the step taken is cancelled and the time step halved until the boundary is no longer crossed. However, because of the halving, say *n* times, the integration time is reduced to 2−*n*∆*t*, leaving a remaining (1 − 2−*n*)∆*t* integration time, which, at a constant step size, requires at least another 2*n* − 1 steps in order to complete the full time-step ∆*t*. Note that at each of these steps it may be needed to further reduce the step size. This further reduction applies only to the current time step, leaving the step size of following sub-steps unaffected. This

now outside the scope of the model and is therefore removed.

method effectively models shear along the coastline.

*∂D ∂y*

P*k*+1(*St* = 1)=(1 − *γ*(*x*, *y*, *t*)∆*k*)P*k*(*St* = 1). (25)

∆*tk* + √

2*D*∆*W*2(*tk*) (24)

We now consider particle flux at open boundaries. This flux is the difference between particles flowing into and out of the domain. The number of particles flowing out should not be controlled as it is a natural consequence of the movement of a particle. As soon as a particle crosses an open boundary it is considered gone and further integration is no longer possible as no data outside the domain is given. For this very reason, however, we do need to explicitly model the particles flowing in. We determine the expected number of particles entering the domain as follows:

$$emp\_{(i,j,t)} = \begin{cases} \frac{\Delta y \cdot \Delta t \cdot V \cdot (\mathcal{U}^2 + V^2) \cdot \lambda\_i}{\gamma \mathcal{M}\_p} & \text{ifflow parallel to y-axis} \\\\ \frac{\Delta x \cdot \Delta t \cdot \mathcal{U} \cdot (\mathcal{U}^2 + V^2) \cdot \lambda\_i}{\gamma \mathcal{M}\_p} & \text{ifflow parallel to x-axis} \end{cases} \tag{26}$$

where, *γ* is the deposition coefficient. The actual number of particles added at the domain boundary at each iteration is obtained by drawing a value from the Poisson distribution parameterised by the above expectation value.

To determine the actual number of particles to be suspended in a grid cell *i*, *j* we draw a number from a Poisson distribution function with mean *enp*(*i*,*j*,*t*) determined by Equation (15). We assume that particles are infinitely many on the sea bed. However, the particle that is suspended is not the same as the one that is deposited.

Before implementing the particle model (11)-(26), we first required to show the consistence between the Fokker-Planck equation and its the Eulerian transport model. This will be described in the next Section.

## **7. Determination of bedlevel changes using particle models**

Comparing with a simplified form of Equation (4), where in this case we assume for the local change in mass *∂m ∂t* ≈ *de* − *se*. We also assume that the deposition and erosion processes balance (see Section 3.1). The approximation of the change in mass with respect to time, in each grid cell *i*, *j* in the this particle model is determined by the following equations:

$$\frac{\partial m}{\partial t} = \frac{\Delta \mathcal{N}\_p}{\Delta t} \frac{1}{\Delta \mathbf{x} \Delta y} \mathcal{M}\_p. \tag{27}$$

Using equation (5) and the fact that *∂m ∂t* ≈ *de* − *se*, the equation for the bed level change can be derived as follows,

$$\frac{\partial h}{\partial t} = \frac{\Delta \mathcal{N}\_p}{\Delta t} \frac{1}{\rho\_s (1 - po) \Delta x \Delta y} \mathcal{M}\_p.$$

Where, M*p* is the mass of each particle, *ρs* and *po* denote respectively the density of an individual grain particle and the bed porosity. From [14], we find *ρs* = 2650*kg*/*m*3 and *po* ≈ 0.5. ∆N*p* is the difference between the number of deposited and suspended particles at each iteration in each grid cell *i*, *j*. Hence the cumulated (integrated) change of the level of the sea ∆*h* in *m* for all time steps is determined by the following equation:

$$
\Delta h = \int\_0^1 \frac{\partial h}{\partial t} dt
$$

$$
\Delta h \approx \sum \frac{\Delta \mathcal{N}\_p}{\rho\_s (1 - p o) \Delta x \Delta y} \mathcal{M}\_p. \tag{28}
$$

10.5772/52237

321

http://dx.doi.org/10.5772/52237

x 104

along *y* is evaluated using central differences wherever possible and forward or backward differences on the boundaries. Note that the average change of the bedforms in each grid cell are estimated by using Equation (10) in Section 3.1. The results obtained are eventually compared with those obtained when the estimations of the changes of bedforms are carried out using Equation (28) with the aid of the two dimensional SDEs or sometimes in this chapter we call it as the Lagrangian particle model. This stochastic model is discussed in Section 7. Note that both results due to Equation (10) and Equation (28) are similar interms

Modelling of Sediment Transport in Shallow Waters by Stochastic and Partial Differential Equations

x 104

−1.5 −1 −0.5 0 0.5 1 1.5 x 10−4

(a) (b)

**Figure 1.** Change of bed level in *m*/*year* for a two dimensional channel (a) is due to the particle model, while (b) the result is

−1 −0.5 0 0.5 1

x−axis

Rate of bedlevel

of shapes, (see Fig. 2).

−2 0 2 4 x 10−4

computed by using the Eulerian approach.

∆ h [m/year]

∆ h [m/year]

−1 −0.5 0 0.5 1

−1

y x

0

the occurrence of erosion of sediments.

2000 4000 6000−4 −2 0 2 4 x 10−4

computed by using the Eulerian approach.

0

1 x 104 −2 −1 0 1 2 3 x 10−4

(a) (b)

**Figure 2.** Change of bed level in *m*/*year* for a two dimensional channel (a) is due to the particle model, while (b) the result is

The Figures 1 and 2 represent the results obtained by solving the same problem by using two approaches. Fig. 1(a) and Fig. 2(a) are due to a simplified particle model with very small effect of the diffusion. On the other hand, in Figures 1(b) and Fig. 2(b) are due to a simplified Eulerian model, no diffusion is considered. We should expected to get deposition at the retardation of the flow and erosion at the acceleration of the flow. Resulting in the net migration of the channels in the direction of (*U*3) as in [19]. The positive sign on the colorbar in the figures imply that deposition is taking place while the negative sign implying

Some results, in Fig. 3(a,b), represent the local depth change in two selected cells. Part (a) shows a steady deposition in the grid cell at the location (*x*, *y*)=(5*km*, 2*km*) whereas another grid cell in the location (*x*, *y*)=(−2*km*, 0.8*km*), part (b), there is also a steady deposition. The diffusion coefficient in the test case for the particle model is 0.00001*m*2/*s*, *γ* = 0.00013*s*−1.

x

More information about the effect of parameters on the sea bed level changes can be found in [16], for example.

#### **7.1. Primary input of the model**

The initial field is defined as

$$depth = h\_0 + h\_1 \cdot \exp(-((\text{x} - 0.0)^2)/(2 \cdot (wd)^2)),$$

where the initial amplitude of the disturbance is *h*1 = 10.0, width of disturbance (*wd*) = 2000*m*, the tidal period *T* = 720 minutes. Constant for sediment transport rate *K* = 0.16*kgs*3/*m*6, porosity *p*0 = 0.5, *ρ* = 2600[*kg*/*m*3], density of sediment, horizontal domain length of *x* goes from −10000*m* to 10000*m*. The horizontal domain length of *y* goes from 2*m* to 4500*m*.

Diffusion constant *D* = 10*m*2/*s*, starting time *Tstart* = 0*s*, number of seconds in a year *Tyear* = 365 · 24 · 60 · 60*s*/*year*. Final time *s Tstop* = 100 · *Tyear*, *M* = 50 is the number of grid points across the channel, *dt* = 1.0 · *year* is the time-step *s*. Tidal mean discharge per unit width *m*2/*s* is given by *q*0 = *h*0 · *u*0, tidal amplitude discharge per unit width *m*2/*s* is given by *q*1 = *h*0 · *u*1. We then determine the changes in the bedforms as described in Section 3.1. The rate of the changes in the level of the floor of the sea is described by the divergence of transport in each grid cell. For comparison purposes, the divergence of transport is we computed by using the original data from the hydrodynamic model as described in this section.

Note that the two dimensional channel in a Matlab code is used whereby in the routine, we first compute the flow fields *U*, *V* as well as its depths in each grid cell. Then compute average load transport and finally compute the divergence using a finite difference methods in a Matlab code. Where the divergence of the *U*-vector field along *x* and the *V*-vector field along *y* is evaluated using central differences wherever possible and forward or backward differences on the boundaries. Note that the average change of the bedforms in each grid cell are estimated by using Equation (10) in Section 3.1. The results obtained are eventually compared with those obtained when the estimations of the changes of bedforms are carried out using Equation (28) with the aid of the two dimensional SDEs or sometimes in this chapter we call it as the Lagrangian particle model. This stochastic model is discussed in Section 7. Note that both results due to Equation (10) and Equation (28) are similar interms of shapes, (see Fig. 2).

12 Sediment Transport

be derived as follows,

in [16], for example.

to 4500*m*.

section.

**7.1. Primary input of the model**

The initial field is defined as

Using equation (5) and the fact that *∂m*

*∂h ∂t* = ∆N*p* ∆*t*

the sea ∆*h* in *m* for all time steps is determined by the following equation:

 *T* 0

*depth* = *h*0 + *h*1 · exp(−((*x* − 0.0).

*∂h ∂t dt*

*ρs*(1 − *po*)∆*x*∆*y*

More information about the effect of parameters on the sea bed level changes can be found

where the initial amplitude of the disturbance is *h*1 = 10.0, width of disturbance (*wd*) = 2000*m*, the tidal period *T* = 720 minutes. Constant for sediment transport rate *K* = 0.16*kgs*3/*m*6, porosity *p*0 = 0.5, *ρ* = 2600[*kg*/*m*3], density of sediment, horizontal domain length of *x* goes from −10000*m* to 10000*m*. The horizontal domain length of *y* goes from 2*m*

Diffusion constant *D* = 10*m*2/*s*, starting time *Tstart* = 0*s*, number of seconds in a year *Tyear* = 365 · 24 · 60 · 60*s*/*year*. Final time *s Tstop* = 100 · *Tyear*, *M* = 50 is the number of grid points across the channel, *dt* = 1.0 · *year* is the time-step *s*. Tidal mean discharge per unit width *m*2/*s* is given by *q*0 = *h*0 · *u*0, tidal amplitude discharge per unit width *m*2/*s* is given by *q*1 = *h*0 · *u*1. We then determine the changes in the bedforms as described in Section 3.1. The rate of the changes in the level of the floor of the sea is described by the divergence of transport in each grid cell. For comparison purposes, the divergence of transport is we computed by using the original data from the hydrodynamic model as described in this

Note that the two dimensional channel in a Matlab code is used whereby in the routine, we first compute the flow fields *U*, *V* as well as its depths in each grid cell. Then compute average load transport and finally compute the divergence using a finite difference methods in a Matlab code. Where the divergence of the *U*-vector field along *x* and the *V*-vector field

∆*h* ≈ ∑ ∆N*p*

∆*h* =

*∂t* ≈ *de* − *se*, the equation for the bed level change can

M*p*. (28)

M*p*.

2)/(2 · (*wd*)2)),

1 *ρs*(1 − *po*)∆*x*∆*y*

Where, M*p* is the mass of each particle, *ρs* and *po* denote respectively the density of an individual grain particle and the bed porosity. From [14], we find *ρs* = 2650*kg*/*m*3 and *po* ≈ 0.5. ∆N*p* is the difference between the number of deposited and suspended particles at each iteration in each grid cell *i*, *j*. Hence the cumulated (integrated) change of the level of

**Figure 1.** Change of bed level in *m*/*year* for a two dimensional channel (a) is due to the particle model, while (b) the result is computed by using the Eulerian approach.

**Figure 2.** Change of bed level in *m*/*year* for a two dimensional channel (a) is due to the particle model, while (b) the result is computed by using the Eulerian approach.

The Figures 1 and 2 represent the results obtained by solving the same problem by using two approaches. Fig. 1(a) and Fig. 2(a) are due to a simplified particle model with very small effect of the diffusion. On the other hand, in Figures 1(b) and Fig. 2(b) are due to a simplified Eulerian model, no diffusion is considered. We should expected to get deposition at the retardation of the flow and erosion at the acceleration of the flow. Resulting in the net migration of the channels in the direction of (*U*3) as in [19]. The positive sign on the colorbar in the figures imply that deposition is taking place while the negative sign implying the occurrence of erosion of sediments.

Some results, in Fig. 3(a,b), represent the local depth change in two selected cells. Part (a) shows a steady deposition in the grid cell at the location (*x*, *y*)=(5*km*, 2*km*) whereas another grid cell in the location (*x*, *y*)=(−2*km*, 0.8*km*), part (b), there is also a steady deposition. The diffusion coefficient in the test case for the particle model is 0.00001*m*2/*s*, *γ* = 0.00013*s*−1.

10.5772/52237

323

http://dx.doi.org/10.5772/52237

= *Se* − *Sd*, (31)

However, the difficulty with the heavy (coarser) material is that the parameter *γ* becomes larger (settling velocity) and this makes the equation (29) 'stiff'. The maximum allowed time step thus becomes very small and since we want to make long simulation this is not very convenient. Therefore, the solution is to make a first order approximation and assume that the two source terms *Se* and *Sd* are much bigger than the right side terms of Equation (29).

+

*CH* = *Q* = *Q*0 + *δQ*1 + O(*δ*2),

*∂*(*HVC*) *∂y*

Modelling of Sediment Transport in Shallow Waters by Stochastic and Partial Differential Equations

*Se* − *Sd* = 0, (32)

*∂y* ) = −(*γQ*1). (33)

*∂y* = 0, (34)

*∂t* . (35)

*∂*(*HUC*) *∂x*

into Equation (31) to get the following equations such that 0-order is given by :

*∂*(*UQ*0) *∂x*

+

*∂*(*UQ*) *∂x*

with (32) as an initial condition for every time step separately and then set

In other words, in Equation (33), we are looking for a source term (*γQ*1) that balances the advective transport of the 0 -order solution. Although this is easy in finite difference approach, however, in the particle model, Equation (33) can be approached as follows. The

+

*γQ*1 ≈ *∂Q*

In other words, Equation (33) should be more or less balanced. If we now omit the source term and measures the rate of changes, these should be approximately equal. Therefore for the heavy particles, between the beginning and end of the integration time loop should do

*∂*(*VQ*)

*∂*(*VQ*0)

now we substitute a Taylor series expansion of the depth averaged concentration *CH*

Therefore, we can write equation (29):

while the 1st order is:

the following

best we can do so far is to solve:

(i) First we remove all deposited particles.

*δ*

 *∂*(*CH*) *∂t* +

> *∂*(*Q*0) *∂t* +

> > *∂Q ∂t* +

(ii) Followed by generating particles according to (see Equation 32)

**Figure 3.** Cumulative local changes in depth by using data of 90 days for two selected grid cells in the ideal two dimensional channel domain light particles

The parameters for the test case experiments using ideal domain with the PM are found in (Table 1.)

**Table 1.** Parameters used by particle model for the sediment transport in the test case

#### **8. Transport of heavy particles in shallow water**

As mentioned earlier, sediment transport is a complex process determined by various properties of the sediment materials. Let us in this section consider that we deal with heavy particles. Heavy particles unlike the lighter ones, tend to attain the equilibrium much faster than the tidal cycle. The description of the derivation of the equations for the heavy particles have followed similar lines using the same equations as those for the lighter particles.

$$\frac{\partial(\text{CH})}{\partial t} + \frac{\partial(H\text{UC})}{\partial \text{x}} + \frac{\partial(HV\text{C})}{\partial y} = \mathcal{S}\_{\varepsilon} - \mathcal{S}\_{d} \tag{29}$$

where in this case *Se* = *λsUk* is erosion term, −*Sd* = *γCH*, is the deposition term where *λs*, *γ*, *k* are constants. In addition, an equation for the change of the bedlevel is presented:

$$\frac{\partial h}{\partial t} = \frac{1}{\rho\_s (1 - po)} \cdot (\mathcal{S}\_\varepsilon - \mathcal{S}\_d). \tag{30}$$

However, the difficulty with the heavy (coarser) material is that the parameter *γ* becomes larger (settling velocity) and this makes the equation (29) 'stiff'. The maximum allowed time step thus becomes very small and since we want to make long simulation this is not very convenient. Therefore, the solution is to make a first order approximation and assume that the two source terms *Se* and *Sd* are much bigger than the right side terms of Equation (29). Therefore, we can write equation (29):

$$\delta \left( \frac{\partial (\mathcal{C}H)}{\partial t} + \frac{\partial (H\mathcal{UC})}{\partial x} + \frac{\partial (HV\mathcal{C})}{\partial y} \right) = \mathcal{S}\_{\varepsilon} - \mathcal{S}\_{d\prime} \tag{31}$$

now we substitute a Taylor series expansion of the depth averaged concentration *CH*

$$\mathcal{C}H = Q = Q\_0 + \delta Q\_1 + \mathcal{O}(\delta^2)\_\prime$$

into Equation (31) to get the following equations such that 0-order is given by :

$$S\_{\varepsilon} - S\_d = 0,\tag{32}$$

while the 1st order is:

14 Sediment Transport

∆ h [m/year]

channel domain light particles

(Table 1.)

presented:

10 x 10−4

Jan '99 Feb '99 Mar '99 Apr '99 −2

time (a) Location at (*x*, *y*) = 5000, 2000)

initial location *m* (1800, 400) *f* = *λs*

**Table 1.** Parameters used by particle model for the sediment transport in the test case

**8. Transport of heavy particles in shallow water**

*∂*(*CH*) *∂t* +

*∂h*

*∂t* = 1

Jan '99 Feb '99 Mar '99 Apr '99 −0.01

(b) Location at (*x*, *y*) = −5000, 800)

time

*γ kg* · *s*2 · *m*−4 O(10−4)

*∂y* = *Se* − *Sd* (29)

*ρs*(1 − *po*) · (*Se* − *Sd*). (30)

(x,y)=(−5000m, 800m)

0

0.01

0.02

∆ h [m/year]

**Figure 3.** Cumulative local changes in depth by using data of 90 days for two selected grid cells in the ideal two dimensional

The parameters for the test case experiments using ideal domain with the PM are found in

**Constant Unit Value Constant Unit Value** porosity - 0.50 *γ s*−1 0.00013 grid offset *m* (−10000, 2) *λs kg*·*s*·*m*−4 0.001 grid size - 22 × 101 *D m*2 · *s*−1 0.000010 cell size *m* 200 × 200 *δ* - 0 sand density *kgm*−3 2650 M*p kg* 3000

As mentioned earlier, sediment transport is a complex process determined by various properties of the sediment materials. Let us in this section consider that we deal with heavy particles. Heavy particles unlike the lighter ones, tend to attain the equilibrium much faster than the tidal cycle. The description of the derivation of the equations for the heavy particles have followed similar lines using the same equations as those for the lighter particles.

+

where in this case *Se* = *λsUk* is erosion term, −*Sd* = *γCH*, is the deposition term where *λs*, *γ*, *k* are constants. In addition, an equation for the change of the bedlevel is

*∂*(*HVC*)

*∂*(*HUC*) *∂x*

0.03

0.04

$$\frac{\partial(Q\_0)}{\partial t} + \frac{\partial(UQ\_0)}{\partial x} + \frac{\partial(VQ\_0)}{\partial y}) = -(\gamma Q\_1). \tag{33}$$

In other words, in Equation (33), we are looking for a source term (*γQ*1) that balances the advective transport of the 0 -order solution. Although this is easy in finite difference approach, however, in the particle model, Equation (33) can be approached as follows. The best we can do so far is to solve:

$$\frac{\partial \mathcal{Q}}{\partial t} + \frac{\partial (\mathcal{U}\mathcal{Q})}{\partial \mathbf{x}} + \frac{\partial (V\mathcal{Q})}{\partial y} = \mathbf{0},\tag{34}$$

with (32) as an initial condition for every time step separately and then set

$$
\gamma Q\_1 \approx \frac{\partial Q}{\partial t}.\tag{35}
$$

In other words, Equation (33) should be more or less balanced. If we now omit the source term and measures the rate of changes, these should be approximately equal. Therefore for the heavy particles, between the beginning and end of the integration time loop should do the following

For now we have assumed that

$$\vert \frac{\partial Q\_0}{\partial t} \vert \ll \vert \gamma Q\_1 \vert.$$

Note that

$$\text{transport vector } q\_{\text{x}} = \text{LUC} = \frac{\lambda\_{\text{s}}}{\gamma} \mathcal{U} (\mathcal{U}^2 + V^2) \tag{36}$$

10.5772/52237

325

http://dx.doi.org/10.5772/52237

**9. Discussion and conclusions**

a type of PDEs.

called point spread function(PSF).

In this chapter we have developed a two dimensional stochastic differential equations. This type of Mathematical model is known as Particle model in some cases. In particular we have developed a two dimensional particle model for sediment transport in which we have added two more equations to model erosion and deposition processes. The description of Transport of materials in shallow waters or in atmosphere can be described by Eulerian approach. This is a deterministic approach in which Partial Differential Equations are used. Numerical methods are usually implemented to approximate the solutions of the PDEs. In this case, one can be faced with computation problems if the dimensional of the model is high. An in most cases if you modify the model to include in more processes it may not be possible to get a closed solution. An alternative approach called the Lagrangian approach can be also used to describe the transport of sediments[35]. This is a probabilistic approach that uses the transition probabilities to derive the state of the transported material. The transition probabilities or density functions is the solution of the Fokker -Planck equation which is also

Modelling of Sediment Transport in Shallow Waters by Stochastic and Partial Differential Equations

The crucial part is to show that FPE is consistent with ADEs, once that is done then one can derive the underlying SDEs. The derived SDEs can be used as a particle model for the

In this chapter we have derived the Fokker-Planck equation and included the deposition and erosion terms based on the developed particle model. Moreover, we have also shown that by interpreting the Eulerian sediment transport as the Fokker Planck equation with the additional terms, it becomes possible to derived the underlying particle model that is consistent with the Eulerian sediment transport model. Furthermore, the results of sediment transport due to particle model has been compared to that obtained by computing the transport using Eulerian model in their simplest form. We have got some results for the changes in the bed level for the data of 90 days for an idealized two dimension domain. We have also used our model to test the prediction of bed-level changes by using the approach of parallel computations of bedforms using the real data of the Dutch North sea[5]. Therefore, at least for now we can say that we have solved the set of mathematical equations called particle model for sediment transport. These equations have given us reasonable results for the sea bedlevel changes. Nevertheless, the determination of the morphological changes is a complicated process that depends on several factors such as waves, the size of sand, mass and density[36]. The particle model in this work has been simplified, what we can say is that the results are reasonable. But more factors will have to be taken into account. For instance, in the particle model we have considered that all particles have equal mass while in reality each particle has different mass. For better predictions of the complex behaviour sediment trasnsport that vary with time and space a feedback among **water motion, sediment transport and bottom changes**, we recommend the coupling of both the hydrodynamic and the transport models. Moreover, make sure the grid mesh are of the same form for the hydrodynamic and particle model. Sometimes you may also be required to change the number of particles into concentration. That can be done by using a function

simulation of the transport of materials in the shallow or atmosphere.

$$\text{transport vector } q\_{\mathcal{Y}} = VHC = \frac{\lambda\_s}{\gamma} V(\mathcal{U}^2 + V^2), \tag{37}$$

since *λsU*2 − *γCH* = 0, the equation for the rate change of the bedlevel due to the transport coarse material is given by

$$\frac{\partial h}{\partial t} = -\left(\frac{\partial q\_x}{\partial x} + \frac{\partial q\_y}{\partial y}\right) \cdot \frac{1}{\rho\_s (1 - po)}\,,\tag{38}$$

where *fd* = *λs γ* is the pick up function which depends into the characteristics of the sediment/sand materials. For example in case of a mixture of larger sand of volume=*l* 3 will have a mass=*ρsl* 3. While a sand of double size whose volume is 8*l* 3 will have mass=*ρs*8*l* 3. Therefore, the two particle require different value of the drag force particle *fd*. But in this chapter we assume that all particles have the same mass, *l*=length.

**Figure 4.** Cumulative local changes in depth by using data of 90 days for two selected grid cells in the ideal two dimensional channel domain for heavy particles

## **9. Discussion and conclusions**

16 Sediment Transport

Note that

For now we have assumed that

coarse material is given by

−0.14 −0.12 −0.1 −0.08 −0.06 −0.04 −0.02 0

channel domain for heavy particles

change in h meters

where *fd* = *λs*

have a mass=*ρsl*

(iii) Store net change in number of particles (concentration).

(v) Compute differences over the previous step using (iii).

*∂h ∂t* = −

chapter we assume that all particles have the same mass, *l*=length.

For a fixed grid cell

Jan '99 Feb '99 Mar '99 Apr '99 −0.16

(a) Local erosion for heavy particles

time (days)

(iv) Next we do one time step integration of the particle using Equation (34).

transport vector *qx* = *UHC* = *λs*

transport vector *qy* = *VHC* = *λs*

 *∂qx ∂x* + *∂qy ∂y*

sediment/sand materials. For example in case of a mixture of larger sand of volume=*l*

Therefore, the two particle require different value of the drag force particle *fd*. But in this

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

change in h meters

**Figure 4.** Cumulative local changes in depth by using data of 90 days for two selected grid cells in the ideal two dimensional

3. While a sand of double size whose volume is 8*l*

*∂t* |≪|*γQ*1|.

since *λsU*2 − *γCH* = 0, the equation for the rate change of the bedlevel due to the transport

*γ*

*γ*

· 1 *ρs*(1 − *po*)

*γ* is the pick up function which depends into the characteristics of the

Jan '99 Feb '99 Mar '99 Apr '99 0

(b) Local deposition for heavy particles

time (days)

For a fixed grid cell

*U*(*U*2 + *V*2) (36)

*V*(*U*2 + *V*2), (37)

, (38)

3 will have mass=*ρs*8*l*

3 will

3.

In this chapter we have developed a two dimensional stochastic differential equations. This type of Mathematical model is known as Particle model in some cases. In particular we have developed a two dimensional particle model for sediment transport in which we have added two more equations to model erosion and deposition processes. The description of Transport of materials in shallow waters or in atmosphere can be described by Eulerian approach. This is a deterministic approach in which Partial Differential Equations are used. Numerical methods are usually implemented to approximate the solutions of the PDEs. In this case, one can be faced with computation problems if the dimensional of the model is high. An in most cases if you modify the model to include in more processes it may not be possible to get a closed solution. An alternative approach called the Lagrangian approach can be also used to describe the transport of sediments[35]. This is a probabilistic approach that uses the transition probabilities to derive the state of the transported material. The transition probabilities or density functions is the solution of the Fokker -Planck equation which is also a type of PDEs.

The crucial part is to show that FPE is consistent with ADEs, once that is done then one can derive the underlying SDEs. The derived SDEs can be used as a particle model for the simulation of the transport of materials in the shallow or atmosphere.

In this chapter we have derived the Fokker-Planck equation and included the deposition and erosion terms based on the developed particle model. Moreover, we have also shown that by interpreting the Eulerian sediment transport as the Fokker Planck equation with the additional terms, it becomes possible to derived the underlying particle model that is consistent with the Eulerian sediment transport model. Furthermore, the results of sediment transport due to particle model has been compared to that obtained by computing the transport using Eulerian model in their simplest form. We have got some results for the changes in the bed level for the data of 90 days for an idealized two dimension domain. We have also used our model to test the prediction of bed-level changes by using the approach of parallel computations of bedforms using the real data of the Dutch North sea[5]. Therefore, at least for now we can say that we have solved the set of mathematical equations called particle model for sediment transport. These equations have given us reasonable results for the sea bedlevel changes. Nevertheless, the determination of the morphological changes is a complicated process that depends on several factors such as waves, the size of sand, mass and density[36]. The particle model in this work has been simplified, what we can say is that the results are reasonable. But more factors will have to be taken into account. For instance, in the particle model we have considered that all particles have equal mass while in reality each particle has different mass. For better predictions of the complex behaviour sediment trasnsport that vary with time and space a feedback among **water motion, sediment transport and bottom changes**, we recommend the coupling of both the hydrodynamic and the transport models. Moreover, make sure the grid mesh are of the same form for the hydrodynamic and particle model. Sometimes you may also be required to change the number of particles into concentration. That can be done by using a function called point spread function(PSF).

## **Author details**

Charles Wilson Mahera and Anton Mtega Narsis

1 University of Dar-es-salaam, College of Natural and Applied Sciences, Department of Mathematics, Dar-es-salaam, Tanzania

10.5772/52237

327

http://dx.doi.org/10.5772/52237

[13] A. Matheja (2000). Sediment transport using the random walk method on top of cobra

Modelling of Sediment Transport in Shallow Waters by Stochastic and Partial Differential Equations

[14] H. M. Schuttelaars and H. E. de Swart (1996). An idealized long-term morhpodynamic

[15] G. S. Stelling (1983). Communications on construction of computational methods for shallow water flow problems *PhD Dissertation*, Delft University of Technology, The

[16] K. de Jong (1998). Tidally Averaged Transport Models *PhD Dissertation*, Delft University

[17] A. F. B. Tompson and L. W. Grelhar(1990). Numerical simulation of the solute transport in three-dimensional,randomly heterogeneous porous media *Water Resources Research*,

[18] H. F. P. van den Boogaard, M. J. J . Hoogkamer, and A. W. Heemink (1993). Parameter identification in particle models *Stochastic Hydrology and hydraulics*, Vol. 7(2), :109-130.

[19] J. van den Kreeke, S. E. Hoogewoning, and M.Verlaan,H. F. P. van den Boogaard, M. J. J . Hoogkamer, and A. W. Heemink (1993). Analytical model for the morphodynamics of a trench in the presence of tidal currents *Continental Shelf Research*, Vol. 22(11-12),

[21] Amoudry, L. O., and A. J. Souza, A. J. (2011) .Deterministic coastal morphological and sediment transport modelling: A review and discussion", Rev. Geophys., 49, RG2002,

[22] Bhuiyan, F. 2011 Stochastic and Deterministic Methods of Computing Graded Bedload Transport, Sediment Transport - Flow and Morphological Processes", Prof. Faruk

[23] Brown, J. M., & Davies, A. G. (2009) Methods for medium term prediction of the net sediment transport by waves and currents in complex coastal regions", Cont. Shelf Res.,

[24] de Vriend, H. J. (1987)2DH mathematical modelling of morphological evolutions in

[25] Dewals, B, Rulot, F., Erpicum, S., Archambeau, P. & Pirotton, M. (2011) *Advanced Topics in Sediment Transport Modelling: Non-alluvial Beds and Hyperconcentrated Flows, Sediment*

[26] Frey, P. & Church, M. (2009) How River Beds Move, *Science*, Vol 325(5947), pp. 1509 -

*Transport*, Dr. Silvia Susana Ginsberg (Ed.), ISBN: 978-953-307-189-3.

[20] L. C. van Rijn (1993). Principles of sediment transport in rivers, estuaries and

*Proc. of the 4th of Inter. Confer. on Hydroinformatics, Iowa City, USA*, 1–7

model of a tidal embayment *Eur. J.Mech., B/Fluids* Vol. 15(1), :55–80

Netherlands

:1811-1820.

of Technology, The Netherlands.

Vol. 26(10) :2541–2562.

doi:10.1029/2010RG000341.

29, pp.1502-1514.

1510.

Bhuiyan (Ed.), ISBN: 978-953-307-374-3.

shallow water, *Coastal Eng.*, Vol 11, pp. 1-27.

Department of General Studies, Dar es Salaam Institute of Technology, Dar es Salaam, Tanzania

## **References**

[13] A. Matheja (2000). Sediment transport using the random walk method on top of cobra *Proc. of the 4th of Inter. Confer. on Hydroinformatics, Iowa City, USA*, 1–7

18 Sediment Transport

Tanzania

**References**

**Author details**

Charles Wilson Mahera and Anton Mtega Narsis

*methods in fluids*, Vol, 25, No. 0, pp. 1373–1384.

WIT Trans. on Ecol. and the Env., Vol 80, Algarve

*and Coastal Waters*, 6th Ed,Academic Press, New York.

˝

*Natural Sciences*, Springer-Verlag, Berlin.

*Hydrology and hydraulics*, Vol. 4, :161-174

Mathematics, Dar-es-salaam, Tanzania

Ed, Springer, New York,.

*Computing*,Vol 68 717U728.

Scientific, Singapore.

Berlin Heidelberg.

Scientific, Singapore.

York.

1 University of Dar-es-salaam, College of Natural and Applied Sciences, Department of

Department of General Studies, Dar es Salaam Institute of Technology, Dar es Salaam,

[1] P. S. Addison, Bo Qu, A. Nisbet, and G. Pender.(1997). A non-Fickian particle-tracking diffusion model based on fractional Brownian motion. *International Journal for numerical*

[2] B.Øksendal (2005). *Stochastic Differential Equations. An introduction with applications*, 6th

[3] Arnold. L. (1974). *Stochastic differential equations: Theory and applications*, Wiley, London.

[4] R. W. Barber & N. P. Volakos (2005). Modelling depth-integrated pollution dispersion in the gulf of thermaikos using a lagrangian particle technique, *Proceedings of Inter. Confer. on Water Resources Management*, pp. 173–184, ISBN, Portugal, April and 2005,

[5] W. M. Charles, E.van den Berg, H. X. Lin, A. W. Heemink, and M. Verlaan (2008). Parallel and distributed simulation of sediment dynamics in shallow water using particle decomposition approach. *International Journal of Parallel and Distributed*

[6] H. B. Fischer, E.J. List, R.C.Y. Koh, J. Imberger, & N.H. Brooks (1979). *Mixing in Inland*

[7] C. W. Gardiner, (2004). *Handbook of Stochastic Methods for physics, Chemistry and the*

[8] J. Fredsøe and R. Deigaard, (1992). *Mechanics of Coastal Sediment transport*, World

[9] A. W. Heemink (1990). Stochastic modeling of dispersion in shallow water *Stochastic*

[10] A. H. Jazwinski, (1970). *Stochastic Processes and Filtering Theory*, Academic Press, New

[11] P. E. Kloeden and E. Platen, (1999). *Numerical solutions of Stochastic Differential equations. Application of Mathematics,Stochastic Modelling and applied probability*,Springer-Verlag,

[12] J. Fredsøe and R. Deigaard, (1992). *Mechanics of Coastal Sediment transport*, World

[27] Harris, C. K., Sherwood, C. R., Signell, R. P., Bever, A. J., & Warner, J. C. (2008) Sediment dispersal in the northwestern Adriatic Sea, *J. Geophys. Res.*, 113, C11S03, doi:10.1029/2006JC003868.

**Chapter 13**

**Numerical Modeling Tidal Circulation and**

Additional information is available at the end of the chapter

**Embayment**

Yu-Hai Wang

**1. Introduction**

http://dx.doi.org/10.5772/51760

grate, merge or perish by siltation.

**1.1. Research question**

**Morphodynamics in a Dumbbell-Shaped Coastal**

Channel-shoal (ridge) system is a common morphological feature in wide, shallow coastal bays and estuarine mouth where tidal flow is relatively stronger. Sediment transport and morphological evolution is complex within such a system as constrained by interacting tidal

As the deeper tidal channels are usually utilized as navigation courses or water supply source for coal & nuclear power plants or other engineering purposes, it is vitally important to maintain the stability of these tidal channels, that is, they should not be allowed to mi‐

This chapter chooses the dumbbell-shaped Qinzhou Bay as the study site to investigate the sediment transport process and resultant morphological evolution of the channel-shoal sys‐ tem within the Bay using numerical simulations under the status quo situation. This is bene‐ ficial to the planned large-scale coastal engineering projects that might exert a profound

Qinzhou Harbor is one of the most important sea harbors connecting southwestern China's inland and the southeast Asian countries, it delivered up to 47.162 million ton cargos in 2011. The harbor has been using waterways as its navigational channel. Meanwhile, coal & nuclear power plants, industrial development zones, recreational parks, land reclamations and many other coastal projects have been or planned to be constructed around the Qinzhou

> © 2013 Wang; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

© 2013 Wang; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

distribution, and reproduction in any medium, provided the original work is properly cited.

Bay coast, all of which compete for the limited shoreline and water area resources.

force, river current, sediment source and characteristics, shoreline configuration, etc.

long-term influence upon the stability of the channel-shoal system.

## **Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment**

Yu-Hai Wang

20 Sediment Transport

p. 182.

doi:10.1029/2006JC003868.

328 Sediment Transport Processes and Their Modelling Applications

*Meteorological Society*, 84pp, Boston.

*Hydraulic Engineering, ASCE, 14pgs*

through Øresund, *J. Mar. Syst.*, Vol 56, pp. 179-193.

[27] Harris, C. K., Sherwood, C. R., Signell, R. P., Bever, A. J., & Warner, J. C. (2008) Sediment dispersal in the northwestern Adriatic Sea, *J. Geophys. Res.*, 113, C11S03,

[28] Hu, K., Ding, P., Wang, Z., and Yang, S. (2009) A 2D/3D hydrodynamic and sediment transport model for the Yangtze Estuary, China, *J. Mar. Syst.*, Vol 77, pp. 114-136.

[29] Lumborg, U. (2005) Modelling the deposition, erosion, and flux of cohesive sediment

[30] Man, C. (2007) *Stochastic Modelling of Suspended Sediment Transport in Regular and Extreme Flow Environments*, PhD Dissertation, State University of New York at Buffalo,

[31] Németh, A. A. & Hulscher, S. J. M. H., (2003) Finite amplitude sand waves in shallow seas, *IAHR-symposium on River, Coastal and Estuarine Morphodynamics*, Barcelona,

[32] Rodean, H.C. (1996) Stochastic Lagrangian Models of Turbulence Diffusion, *American,*

[33] Struiksma, N., Olesen, K. W., Flokstra, C., & de Vriend, H. J. (1985) Bed deformation

[34] Papanicolauo, A. N., Elhakeem, M., Krallis, G., Prakash, S., & Edinger, J., (2008) Sediment Transport Modeling Review-Current and Future Developments, *Journalof*

[35] Hunter, J. R,Crais, P.D and Phillips H. E.(1993) On the use of random walk models

[36] H. M. Schuttelaars and H. E. de Swart (2000) Multiple morphodynamic equilibria in

with spatially variable diffusivity.*Computation Physics*, 106:366–376, 1993.

Spanje, editors Sánchez-Arcilla, S. en Bateman, A., pp. 435-444.

in curved alluvial channels, *J. Hydraul. Res.*, Vol 23(1), pp. 57-79.

tidal embayments.*J. Geophysical Research*, 105(c10):14,105–14,118.

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/51760

## **1. Introduction**

Channel-shoal (ridge) system is a common morphological feature in wide, shallow coastal bays and estuarine mouth where tidal flow is relatively stronger. Sediment transport and morphological evolution is complex within such a system as constrained by interacting tidal force, river current, sediment source and characteristics, shoreline configuration, etc.

As the deeper tidal channels are usually utilized as navigation courses or water supply source for coal & nuclear power plants or other engineering purposes, it is vitally important to maintain the stability of these tidal channels, that is, they should not be allowed to mi‐ grate, merge or perish by siltation.

This chapter chooses the dumbbell-shaped Qinzhou Bay as the study site to investigate the sediment transport process and resultant morphological evolution of the channel-shoal sys‐ tem within the Bay using numerical simulations under the status quo situation. This is bene‐ ficial to the planned large-scale coastal engineering projects that might exert a profound long-term influence upon the stability of the channel-shoal system.

#### **1.1. Research question**

Qinzhou Harbor is one of the most important sea harbors connecting southwestern China's inland and the southeast Asian countries, it delivered up to 47.162 million ton cargos in 2011. The harbor has been using waterways as its navigational channel. Meanwhile, coal & nuclear power plants, industrial development zones, recreational parks, land reclamations and many other coastal projects have been or planned to be constructed around the Qinzhou Bay coast, all of which compete for the limited shoreline and water area resources.

For the optimum planning of the coastal engineering projects regarding size, location and the sustainable regional economic, social and environmental development, it is critically ur‐ gent to know how the stability of the channel-shoal system under the present coastal config‐ uration and bathymetry. This might be answered by investigating the sediment transport processes and morphodynamics of the channel-shoal system.

sand ridges/shoals between them (Figure 1). The Yingling tidal inlet is 10.1 km long and 1.1~3.5 km wide with a water depth of 5~20 m; it is a rocky inlet with a total of 71 vari‐

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment

http://dx.doi.org/10.5772/51760

331

The spring tide and middle tide in Qinzhou Bay are diurnal throughout most of the year but become irregularly diurnal in March and September each year, while the neap tide is usual‐ ly semi-diurnal throughout the whole year. The mean tidal range is 2.51m, and the maxi‐ mum tidal range is 5.27m. The flood tides last longer than the ebb tides in spring, middle and neap tides. It is 13 hours plus 14 minutes, 11 hours plus 18 minutes for spring tide; 14 hours plus 36 minutes and 10 hours plus 7 minutes for middle tide, and 6 hours plus 33 mi‐

Tidal flow in the outer bay demonstrates as reciprocating flow in parallel to the major water‐ ways; the mean flood velocity, ebb velocity of spring tide is 0.37m/s and 0.51 m/s, respec‐ tively, the mean flood velocity, ebb velocity of middle tide is 0.33m/s and 0.38 m/s, respectively, the mean flood velocity, ebb velocity of neap tide is 0.22m/s and 0.18 m/s; while flow velocity in the Yingling inlet becomes significantly larger, the mean flood veloci‐ ty, ebb velocity of spring tide is 0.67m/s and 0.90 m/s, the mean flood velocity, ebb velocity of middle tide is 0.57m/s and 0.68 m/s, the mean flood velocity, ebb velocity of neap tide is 0.42m/s and 0.33m/s, respectively, and the maximum flood flow velocity reaches up to 1.40

The Qinzhou Bay is influenced by subtropical monsoon and the waves within the Bay are

The waves in winter season (October~Apirl) prevail in N-NE direction while they prevail in S-SW direction in summer season (May-September) and the stronger waves propagate in SSW, SSE directions; the mean wave height is 0.52m and mean wave period is 3.1s [2].

Suspended sediment concentration within Qinzhou Bay water is generally low. In summer

, occurred in spring ebb tide, the maximum middle tidal concentration

in middle tide and 0.013 kg/m3

, among which, it is 0.035 kg/m3 in

in neap tide; the maximum concen‐

ous-sized islands and 72 narrow, small waterways.

nutes and 5 hours plus 40 minutes for neap tide, respectively [1].

m/s and the maximum ebb flow velocity is up to 1.32 m/s [1].

2009, the mean full tidal concentration is 0.022 kg/m3

mainly wind-driven with some surge waves traveled from the open sea.

**2. Dynamics**

**2.1. Tidal regime**

**2.2. Wave climate**

**3. Sediment characteristics**

**3.1. Suspended sediment**

spring tide, 0.020 kg/m3

tration is 0.081 kg/m3

## **1.2. Site description**

Qinzhou Bay is located on Guangxi Province's coast facing South China Sea (Figure 1).

**Figure 1.** Sketch map of Qinzhou Bay coastal configuration (the axis is local coordinate, unit: m).

Appearing as a dumbbell shape it consists of three parts: the inner enclosed bay, also known as Maoweihai, the outer bay and the Yingling tidal inlet connecting two bays. Two rivers, i.e. Maolingjiang and Qinjiang flow into the inner bay, delivering annually 27.73×109 m3 water and 86.4×103 t sediment; the outer bay is a show, trumpet-shaped bay, its area is nearly 2.55×108 m2 with a mean depth of 4.67m (calculated by mean sea level). A complex channel-shoal system is present within the outer bay, consisting of three domi‐ nant tidal waterways, i.e. the east waterway, middle waterway and west waterway and sand ridges/shoals between them (Figure 1). The Yingling tidal inlet is 10.1 km long and 1.1~3.5 km wide with a water depth of 5~20 m; it is a rocky inlet with a total of 71 vari‐ ous-sized islands and 72 narrow, small waterways.

## **2. Dynamics**

For the optimum planning of the coastal engineering projects regarding size, location and the sustainable regional economic, social and environmental development, it is critically ur‐ gent to know how the stability of the channel-shoal system under the present coastal config‐ uration and bathymetry. This might be answered by investigating the sediment transport

Qinzhou Bay is located on Guangxi Province's coast facing South China Sea (Figure 1).

**Figure 1.** Sketch map of Qinzhou Bay coastal configuration (the axis is local coordinate, unit: m).

Appearing as a dumbbell shape it consists of three parts: the inner enclosed bay, also known as Maoweihai, the outer bay and the Yingling tidal inlet connecting two bays. Two rivers, i.e. Maolingjiang and Qinjiang flow into the inner bay, delivering annually 27.73×109 m3 water and 86.4×103 t sediment; the outer bay is a show, trumpet-shaped bay, its area is nearly 2.55×108 m2 with a mean depth of 4.67m (calculated by mean sea level). A complex channel-shoal system is present within the outer bay, consisting of three domi‐ nant tidal waterways, i.e. the east waterway, middle waterway and west waterway and

processes and morphodynamics of the channel-shoal system.

330 Sediment Transport Processes and Their Modelling Applications

**1.2. Site description**

## **2.1. Tidal regime**

The spring tide and middle tide in Qinzhou Bay are diurnal throughout most of the year but become irregularly diurnal in March and September each year, while the neap tide is usual‐ ly semi-diurnal throughout the whole year. The mean tidal range is 2.51m, and the maxi‐ mum tidal range is 5.27m. The flood tides last longer than the ebb tides in spring, middle and neap tides. It is 13 hours plus 14 minutes, 11 hours plus 18 minutes for spring tide; 14 hours plus 36 minutes and 10 hours plus 7 minutes for middle tide, and 6 hours plus 33 mi‐ nutes and 5 hours plus 40 minutes for neap tide, respectively [1].

Tidal flow in the outer bay demonstrates as reciprocating flow in parallel to the major water‐ ways; the mean flood velocity, ebb velocity of spring tide is 0.37m/s and 0.51 m/s, respec‐ tively, the mean flood velocity, ebb velocity of middle tide is 0.33m/s and 0.38 m/s, respectively, the mean flood velocity, ebb velocity of neap tide is 0.22m/s and 0.18 m/s; while flow velocity in the Yingling inlet becomes significantly larger, the mean flood veloci‐ ty, ebb velocity of spring tide is 0.67m/s and 0.90 m/s, the mean flood velocity, ebb velocity of middle tide is 0.57m/s and 0.68 m/s, the mean flood velocity, ebb velocity of neap tide is 0.42m/s and 0.33m/s, respectively, and the maximum flood flow velocity reaches up to 1.40 m/s and the maximum ebb flow velocity is up to 1.32 m/s [1].

#### **2.2. Wave climate**

The Qinzhou Bay is influenced by subtropical monsoon and the waves within the Bay are mainly wind-driven with some surge waves traveled from the open sea.

The waves in winter season (October~Apirl) prevail in N-NE direction while they prevail in S-SW direction in summer season (May-September) and the stronger waves propagate in SSW, SSE directions; the mean wave height is 0.52m and mean wave period is 3.1s [2].

## **3. Sediment characteristics**

#### **3.1. Suspended sediment**

Suspended sediment concentration within Qinzhou Bay water is generally low. In summer 2009, the mean full tidal concentration is 0.022 kg/m3 , among which, it is 0.035 kg/m3 in spring tide, 0.020 kg/m3 in middle tide and 0.013 kg/m3 in neap tide; the maximum concen‐ tration is 0.081 kg/m3 , occurred in spring ebb tide, the maximum middle tidal concentration is 0.034 kg/m3 , occurred also in ebb tidal period, the maximum neap tidal concentration is 0.025 kg/m3 , occurred also in ebb tidal period [1].

The gain sizes vary remarkably in 0.0027~1.099mm (Figure 2). The spatial mean median di‐ ameter (D50) in the inner bay is 0.334mm; it is 0.356mm in the deep channel but becomes 0.0041mm in the shallow parts in the Yingling inlet; in the outer bay, the spatial mean D50 is about 0.298mm with an overall deposition pattern: coarser in waterways but finer in shoals,

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment

http://dx.doi.org/10.5772/51760

333

The Qinzhou Bay has been a drowned rocky valley by the last sea level transgression since 7,000-8,000 year before present [3]. Therefore, the huge amount of sand deposits in the outer bay has come from the deposits by paleo-Maolingjiang river and paleo-Qinjiang river, they have been reformed into the contemporary channel-shoal geomorphology by tidal dynamics. At the present day, sediments delivered by these two rivers are deposited within the inner bay with limited amount of fine particles transported into the outer bay and open sea; meanwhile, limited amount of sediment eroded from the adjacent slopes by storm rains also enter the outer bay. Generally speaking, sediment from the open sea into the Qinzhou Bay is very limited.

A 3D unstructured grid, finite-volume coastal ocean model (called FVCOM) has been devel‐ oped in the Marine Ecosystem Dynamics Modeling Laboratory led by Dr. C. Chen at the University of Massachusetts–Dartmouth (UMASS-D) in collaboration with Dr. R. Beardsley at the Woods Hole Oceanographic Institute. FVCOM is a three-dimensional (3D) primitive equation ocean model, consisting of momentum, continuity, sediment, temperature, salinity, and density equations and is closed physically and mathematically using the Mellor and Ya‐ mada level-2.5 turbulent closure submodel; the irregular bottom slope is represented using a σ-coordinate transformation, and the horizontal grids comprise unstructured triangular cells; the finite-volume method used in the model combines the advantages of a finite-ele‐ ment method for geometric flexibility and a finite-difference method for simple discrete computation; current, sediment, temperature, and salinity in the model are computed in the integral form of the equations, which provides a better representation of the conservative

laws for mass, momentum, and heat in the coastal region with complex geometry [4].

The governing equations consist of the following momentum, continuity, temperature, sal‐

0 *uvw xyz* ¶¶¶ ++ =

¶¶¶ (1)

and coarser in the western part than in the eastern part of the bay [1].

**3.3. Sediment source**

**4. Numerical model**

**4.1. The primitive equations**

inity, and density equations:

(1) continuity equation

The medium diameters of suspended sediment vary within 0.0067~0.0152mm with a mean value of 0.0101mm. The suspended sediment is mainly clayey silt with 30.8% clay particles, 53.4% silt particles and 15.8% sand particles (Table 1). The sorting index is 1.90 [1].

**Table 1.** Grading of suspended sediment.

#### **3.2. Bottom sediment**

The bottom sediments in Qinzhou Bay mainly consist of gravel, coarse sand, medium sand, fine sand, silty clay and clayey silt, etc.

**Figure 2.** Bottom sediment grain size distribution in Qinzhou Bay.

The gain sizes vary remarkably in 0.0027~1.099mm (Figure 2). The spatial mean median di‐ ameter (D50) in the inner bay is 0.334mm; it is 0.356mm in the deep channel but becomes 0.0041mm in the shallow parts in the Yingling inlet; in the outer bay, the spatial mean D50 is about 0.298mm with an overall deposition pattern: coarser in waterways but finer in shoals, and coarser in the western part than in the eastern part of the bay [1].

#### **3.3. Sediment source**

is 0.034 kg/m3

**Table 1.** Grading of suspended sediment.

fine sand, silty clay and clayey silt, etc.

**Figure 2.** Bottom sediment grain size distribution in Qinzhou Bay.

**3.2. Bottom sediment**

, occurred also in ebb tidal period [1].

332 Sediment Transport Processes and Their Modelling Applications

0.025 kg/m3

, occurred also in ebb tidal period, the maximum neap tidal concentration is

The medium diameters of suspended sediment vary within 0.0067~0.0152mm with a mean value of 0.0101mm. The suspended sediment is mainly clayey silt with 30.8% clay particles,

**Grading sand silt clay Sorting index**

The bottom sediments in Qinzhou Bay mainly consist of gravel, coarse sand, medium sand,

percentage% 15.8 53.4 30.8 1.90

53.4% silt particles and 15.8% sand particles (Table 1). The sorting index is 1.90 [1].

The Qinzhou Bay has been a drowned rocky valley by the last sea level transgression since 7,000-8,000 year before present [3]. Therefore, the huge amount of sand deposits in the outer bay has come from the deposits by paleo-Maolingjiang river and paleo-Qinjiang river, they have been reformed into the contemporary channel-shoal geomorphology by tidal dynamics. At the present day, sediments delivered by these two rivers are deposited within the inner bay with limited amount of fine particles transported into the outer bay and open sea; meanwhile, limited amount of sediment eroded from the adjacent slopes by storm rains also enter the outer bay. Generally speaking, sediment from the open sea into the Qinzhou Bay is very limited.

## **4. Numerical model**

A 3D unstructured grid, finite-volume coastal ocean model (called FVCOM) has been devel‐ oped in the Marine Ecosystem Dynamics Modeling Laboratory led by Dr. C. Chen at the University of Massachusetts–Dartmouth (UMASS-D) in collaboration with Dr. R. Beardsley at the Woods Hole Oceanographic Institute. FVCOM is a three-dimensional (3D) primitive equation ocean model, consisting of momentum, continuity, sediment, temperature, salinity, and density equations and is closed physically and mathematically using the Mellor and Ya‐ mada level-2.5 turbulent closure submodel; the irregular bottom slope is represented using a σ-coordinate transformation, and the horizontal grids comprise unstructured triangular cells; the finite-volume method used in the model combines the advantages of a finite-ele‐ ment method for geometric flexibility and a finite-difference method for simple discrete computation; current, sediment, temperature, and salinity in the model are computed in the integral form of the equations, which provides a better representation of the conservative laws for mass, momentum, and heat in the coastal region with complex geometry [4].

#### **4.1. The primitive equations**

The governing equations consist of the following momentum, continuity, temperature, sal‐ inity, and density equations:

(1) continuity equation

$$\frac{\partial u}{\partial \mathbf{x}} + \frac{\partial v}{\partial \mathbf{y}} + \frac{\partial w}{\partial \mathbf{z}} = \mathbf{0} \tag{1}$$

#### (2) *x*-direction momentum equation

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial \mathbf{x}} + \nu \frac{\partial u}{\partial \mathbf{y}} + \nu \frac{\partial u}{\partial z} - f\nu = -\frac{1}{\rho\_0} \frac{\partial P}{\partial \mathbf{x}} + \frac{\partial}{\partial z} (K\_m \frac{\partial u}{\partial z}) + F\_u \tag{2}$$

tions, the latter is solved numerically using a simple combined explicit and implicit scheme, in which the local change of the current is integrated using the first-order accuracy upwind scheme; the advection terms are computed explicitly by a second-order accuracy Runge–

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment

http://dx.doi.org/10.5772/51760

335

FVCOM adopts the Community Numerical Modeling System to simulate erosion, trans‐ port, deposition and the fate of sediments in the coastal ocean developed by experts from USGS [5]. The sediment-transport algorithms are implemented for an unlimited number of user-defined noncohesive/cohesive sediment classes. Each class has attributes of grain diameter, density, settling velocity, critical stress threshold for erosion, and erodibility constant. These properties are used to determine bulk properties of each bed layer. Sus‐ pended-sediment transport in the water column is computed with the same advection-dif‐ fusion algorithm used for all passive tracers and an additional algorithm for vertical

( ) ( )( )( ) *iii i i i i i*

here *Ci* is the concentration of the *i th* sediment class, *As*,*h* is the horizontal sediment diffusivi‐

*C uC vC C C C C AAK t x y z x xy yz z*

> , 0, *i s h C K z t*

, , *i s h i i C K E Dz D*

*i i b bi*

*E tE P F*

here *D* is water depth. While the erosion flux of the *i th* sediment class is computed as:

0 (1 ) ( 1) *b*

*ci*

t

t

*t*

here *ζ*is surface elevation. At the bottom boundary, the vertical sediment flux is specified by:

+++ = + +

,, ,

¶¶¶ ¶¶¶

¶¶ ¶ ¶ ¶ ¶¶ ¶¶ ¶ (8)

z

*s h s h s h*

is the settling velocity of the *i th* sediment class

¶ = = ¶ (9)

¶ = - =- ¶ (10)

=D - - (11)

Kutta time-stepping scheme [4].

settling that is not limited by the CFL criterion [5].

w w

At the top boundary, the vertical diffusive flux is set to be zero:

Suspended sediment transport equation is:

ty, *Ks*,*h* is the vertical sediment diffusivity, *ωi*

¶ ¶ ¶ ¶-

**4.3. Sediment computation**

*4.3.1. Suspended sediment*

given by the user.

#### (3) *y*-direction momentum equation

$$\frac{\partial \mathbf{v}}{\partial t} + \boldsymbol{\mu} \frac{\partial \mathbf{v}}{\partial \mathbf{x}} + \boldsymbol{\nu} \frac{\partial \mathbf{v}}{\partial \mathbf{y}} + \boldsymbol{\nu} \frac{\partial \mathbf{v}}{\partial \mathbf{z}} + f \mathbf{u} = -\frac{1}{\rho\_0} \frac{\partial P}{\partial \mathbf{y}} + \frac{\partial}{\partial \mathbf{z}} (K\_\mathbf{u} \frac{\partial \mathbf{v}}{\partial \mathbf{z}}) + F\_\mathbf{v} \tag{3}$$

#### (4) temperature equation

$$\frac{\partial T}{\partial t} + \mu \frac{\partial T}{\partial \mathbf{x}} + \nu \frac{\partial T}{\partial \mathbf{y}} + \nu \frac{\partial T}{\partial \mathbf{z}} = \frac{\partial}{\partial \mathbf{z}} (K\_h \frac{\partial T}{\partial \mathbf{z}}) + F\_T \tag{4}$$

(5) salinity equation

$$\frac{\partial \mathbf{\hat{s}}}{\partial t} + \mu \frac{\partial \mathbf{\hat{s}}}{\partial \mathbf{x}} + \nu \frac{\partial \mathbf{\hat{s}}}{\partial \mathbf{y}} + \nu \frac{\partial \mathbf{\hat{s}}}{\partial z} = \frac{\partial}{\partial z} (K\_h \frac{\partial \mathbf{\hat{s}}}{\partial z}) + F\_s \tag{5}$$

(6) pressure equation

$$\frac{\partial P}{\partial t} = -\rho \mathbf{g} \tag{6}$$

$$
\rho = \rho(T, s) \tag{7}
$$

where *x, y,* and *z* are the east, north, and vertical axis of the Cartesian coordinate; *u,* v, and *w* are the *x, y, z* velocity components; *T* is the potential temperature; *s* is the salinity; *P* is the pressure; *f* is the Coriolis parameter; *g* is the gravitational acceleration; *Km* is the vertical ed‐ dy viscosity coefficient; and *Kh* is the thermal vertical eddy diffusion coefficient. Here *Fu*, *F*y, *F*T, and *Fs* represent the horizontal momentum, thermal, and salt diffusion terms.

#### **4.2. Numerical solutions**

The momentum and continuity equations are solved using a 'model splitting' method [4], that is, the current is divided into external and internal modes that can be computed using two distinct time steps. The external mode is used to solve the 2D vertically integrated mo‐ mentum and continuity equations while the internal mode is computed for the 3D equa‐ tions, the latter is solved numerically using a simple combined explicit and implicit scheme, in which the local change of the current is integrated using the first-order accuracy upwind scheme; the advection terms are computed explicitly by a second-order accuracy Runge– Kutta time-stepping scheme [4].

#### **4.3. Sediment computation**

(2) *x*-direction momentum equation

334 Sediment Transport Processes and Their Modelling Applications

(3) *y*-direction momentum equation

(4) temperature equation

(5) salinity equation

(6) pressure equation

**4.2.  Here *Fu*, *F*y,

The momentum and continuity equations are solved using a 'model splitting' method [4], that is, the current is divided into external and internal modes that can be computed using two distinct time steps. The external mode is used to solve the 2D vertically integrated mo‐ mentum and continuity equations while the internal mode is computed for the 3D equa‐

r

¶ ¶ ¶ ¶ ¶¶ ¶ + + + + =- + +

¶ ¶ ¶ ¶ ¶¶ ¶ + + + - =- + +

1 ( ) *m u*

1 ( ) *m v*

( ) *h T*

( ) *h s*

¶ ¶ ¶ ¶¶ ¶ (5)

¶ = - ¶ (6)

= (,) *T s* (7)

¶ ¶ ¶ ¶¶ ¶ (4)

¶ ¶ ¶ ¶ ¶¶ ¶ (2)

¶ ¶ ¶ ¶ ¶¶ ¶ (3)

FVCOM adopts the Community Numerical Modeling System to simulate erosion, trans‐ port, deposition and the fate of sediments in the coastal ocean developed by experts from USGS [5]. The sediment-transport algorithms are implemented for an unlimited number of user-defined noncohesive/cohesive sediment classes. Each class has attributes of grain diameter, density, settling velocity, critical stress threshold for erosion, and erodibility constant. These properties are used to determine bulk properties of each bed layer. Sus‐ pended-sediment transport in the water column is computed with the same advection-dif‐ fusion algorithm used for all passive tracers and an additional algorithm for vertical settling that is not limited by the CFL criterion [5].

#### *4.3.1. Suspended sediment*

Suspended sediment transport equation is:

$$\frac{\partial C\_{i}}{\partial t} + \frac{\partial u C\_{i}}{\partial \mathbf{x}} + \frac{\partial \mathbf{v} C\_{i}}{\partial \mathbf{y}} + \frac{\partial (\alpha - \alpha\_{i}) C\_{i}}{\partial \mathbf{z}} = \frac{\partial}{\partial \mathbf{x}} (A\_{s,h} \frac{\partial C\_{i}}{\partial \mathbf{x}}) + \frac{\partial}{\partial \mathbf{y}} (A\_{s,h} \frac{\partial C\_{i}}{\partial \mathbf{y}}) + \frac{\partial}{\partial \mathbf{z}} (K\_{s,h} \frac{\partial C\_{i}}{\partial \mathbf{z}}) \tag{8}$$

here *Ci* is the concentration of the *i th* sediment class, *As*,*h* is the horizontal sediment diffusivi‐ ty, *Ks*,*h* is the vertical sediment diffusivity, *ωi* is the settling velocity of the *i th* sediment class given by the user.

At the top boundary, the vertical diffusive flux is set to be zero:

$$K\_{s,h} \frac{\partial C\_i}{\partial t} = 0, z = \zeta \tag{9}$$

here *ζ*is surface elevation. At the bottom boundary, the vertical sediment flux is specified by:

$$K\_{s,h} \frac{\partial C\_i}{\partial t} = E\_i - D\_i,\\ \mathbf{z} = -D \tag{10}$$

here *D* is water depth. While the erosion flux of the *i th* sediment class is computed as:

$$E\_i = \Delta t E\_{0i} (1 - P\_b) F\_{bi} (\frac{\tau\_b}{\tau\_{cl}} - 1) \tag{11}$$

where *E*0*i* is an empirical constant set by the user as the erosion rate of the *i th* sediment class, *Pb*is the porosity of the bed,*Fbi* is the fraction of the *i th* sediment class, *τb*is the bottom flow shear stress, *τci* is the critical stress for the incipient motion of the *i th* sediment class.

The bed layers are modified at each time step to account for erosion and deposition (Figure 3) and track stratigraphy; at the beginning of each time step, an active-layer thickness *z a* is

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment

where *τsf* is bottom skin-friction stress; *τce* ¯is the critical stress for erosion; and the overbar

indicates this is averaged over all sediment classes; D50 is the median grain diameter of sur‐

Each sediment class can be transported by suspended-load and/or bedload; suspended-load mass is exchanged vertically between the water column and the top bed layer; mass of each sediment class available for transport is limited to the mass available in the active layer; bed‐ load mass is exchanged horizontally between the top layers of the bed; mass of each sedi‐

If continuous deposition results in a top layer thicker than a user-defined threshold, a new lay‐ er is provided to begin accumulation of depositing mass; the bottom two layers are then com‐ bined to conserve the number of layers; after erosion and deposition have been calculated, the active-layer thickness is recalculated and bed layers are readjusted to accommodate it [5].

*k D* ë û (12)

http://dx.doi.org/10.5772/51760

337

*a* max ( ) , 0 1 *sf ce* 0 2 50 *z k* = -+ é ù t tr

face sediment; and *k 1* and *k 2* are empirical constants, 0.007 and 6.0, respectively.

ment class available for transport is limited to the mass available in the top layer [5].

Figure 4 shows the computation domain consisting of unstructured triangular grids.

It has 16466 nodes and 29722 elements in each horizontal layer and 7 sigma-levels in the ver‐ tical. The horizontal grid resolution varies from 2,000 m at the open boundary to 350m in the

calculated based on the relation of Harris and Wiberg [7]:

**5. Model setup**

**5.1. Model domain**

**Figure 4.** Unstructured grids of Qinzhou Bay.

When the bottom shear stress passes the critical stress erosion occurs on the bed. The sedi‐ ment concentration profile in the water body is determined by horizontal convection, diffu‐ sion, vertical diffusion, settling and bottom erosion flux [5].

#### *4.3.2. Bedload transport computation*

Bedload transport rate is computed using established empirical formula, i.e. the Meyer-Peter and Muller formula or using formulas that the modeler considers appropriate, for example, the theoretical-based formula [6].

#### *4.3.3. Sediment bed*

The sediment bed consists of a constant number layers, and each layer is initialized with a thickness, sediment-class distribution, porosity, and age, the mass of each sediment class can be determined from these values and the grain density; the bed evolving properties include bulk properties of the surface layer (active-layer thickness, mean grain diameter, mean den‐ sity, mean settling velocity, mean critical stress for erosion) [5].

**Figure 3.** Distribution of vertical layers in bed model (from [5]).

The bed layers are modified at each time step to account for erosion and deposition (Figure 3) and track stratigraphy; at the beginning of each time step, an active-layer thickness *z a* is calculated based on the relation of Harris and Wiberg [7]:

$$z\_a = \max\left[k\_1(\tau\_{sf} - \overline{\tau\_{ca}})\rho\_0, 0\right] + k\_2 D\_{90} \tag{12}$$

where *τsf* is bottom skin-friction stress; *τce* ¯is the critical stress for erosion; and the overbar indicates this is averaged over all sediment classes; D50 is the median grain diameter of sur‐ face sediment; and *k 1* and *k 2* are empirical constants, 0.007 and 6.0, respectively.

Each sediment class can be transported by suspended-load and/or bedload; suspended-load mass is exchanged vertically between the water column and the top bed layer; mass of each sediment class available for transport is limited to the mass available in the active layer; bed‐ load mass is exchanged horizontally between the top layers of the bed; mass of each sedi‐ ment class available for transport is limited to the mass available in the top layer [5].

If continuous deposition results in a top layer thicker than a user-defined threshold, a new lay‐ er is provided to begin accumulation of depositing mass; the bottom two layers are then com‐ bined to conserve the number of layers; after erosion and deposition have been calculated, the active-layer thickness is recalculated and bed layers are readjusted to accommodate it [5].

## **5. Model setup**

where *E*0*i* is an empirical constant set by the user as the erosion rate of the *i th* sediment class,

When the bottom shear stress passes the critical stress erosion occurs on the bed. The sedi‐ ment concentration profile in the water body is determined by horizontal convection, diffu‐

Bedload transport rate is computed using established empirical formula, i.e. the Meyer-Peter and Muller formula or using formulas that the modeler considers appropriate, for example,

The sediment bed consists of a constant number layers, and each layer is initialized with a thickness, sediment-class distribution, porosity, and age, the mass of each sediment class can be determined from these values and the grain density; the bed evolving properties include bulk properties of the surface layer (active-layer thickness, mean grain diameter, mean den‐

is the critical stress for the incipient motion of the *i th* sediment class.

is the fraction of the *i th* sediment class, *τb*is the bottom flow

*Pb*is the porosity of the bed,*Fbi*

336 Sediment Transport Processes and Their Modelling Applications

*4.3.2. Bedload transport computation*

the theoretical-based formula [6].

sion, vertical diffusion, settling and bottom erosion flux [5].

sity, mean settling velocity, mean critical stress for erosion) [5].

**Figure 3.** Distribution of vertical layers in bed model (from [5]).

shear stress, *τci*

*4.3.3. Sediment bed*

## **5.1. Model domain**

Figure 4 shows the computation domain consisting of unstructured triangular grids.

**Figure 4.** Unstructured grids of Qinzhou Bay.

It has 16466 nodes and 29722 elements in each horizontal layer and 7 sigma-levels in the ver‐ tical. The horizontal grid resolution varies from 2,000 m at the open boundary to 350m in the channel-shoal region to 150 m in the Yingling inlet, especially, down-to-30-m elements are interpolated around islands and in the estuarine channels.

The present study firstly performs the calibrations of water level, flow velocity& direction

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment

http://dx.doi.org/10.5772/51760

339

The Calibrations for water levels at three tidal gauges are shown in Figure 6.

and suspended-sediment concentrations.

**Figure 6.** Tidal water level calibrations (the ordinate unit of is meter).

The flow calibrations are shown in Figure 7. For limited space only calibrations for spring

**6.2. Flow calibrations**

tide are show here.

**6.1. Tidal water level**

#### **5.2. Boundary conditions**

The open boundary uses observed water level as input condition, for this purpose halfmonth water levels at each open boundary grid are interpolated from two tidal station, i.e. the Beihai Station and the Bailongwei Station (Figure 4). River boundaries use annuallymean discharges as input conditions.

Suspended-sediment concentrations at the open boundary grids are interpolated from ad‐ jacent observation sites, while those at river input grids are annually-mean suspendedsediment discharges.

The bathymetry in the computation domain consists of a local bathymetric survey in the Qinzhou Bay in 2008 and sea maps surveys in 2004 and 1997 to supplement other parts.

The time step for both of external mode and internal mode is 1 s.

## **6. Calibrations**

TCZC [1] measured half month water level at three temporary tidal gauges, namely, Guo‐ zishan, Shabadun and Wulei, also measured 26-hour spring, middle and neap tidal flow ve‐ locity & direction and suspended-sediment concentrations at 10 sites (Figure 5).

**Figure 5.** Tidal gauges (red square) and flow observations sites (green triangular) in Qinzhou Bay.

The present study firstly performs the calibrations of water level, flow velocity& direction and suspended-sediment concentrations.

#### **6.1. Tidal water level**

channel-shoal region to 150 m in the Yingling inlet, especially, down-to-30-m elements are

The open boundary uses observed water level as input condition, for this purpose halfmonth water levels at each open boundary grid are interpolated from two tidal station, i.e. the Beihai Station and the Bailongwei Station (Figure 4). River boundaries use annually-

Suspended-sediment concentrations at the open boundary grids are interpolated from ad‐ jacent observation sites, while those at river input grids are annually-mean suspended-

The bathymetry in the computation domain consists of a local bathymetric survey in the Qinzhou Bay in 2008 and sea maps surveys in 2004 and 1997 to supplement other parts.

TCZC [1] measured half month water level at three temporary tidal gauges, namely, Guo‐ zishan, Shabadun and Wulei, also measured 26-hour spring, middle and neap tidal flow ve‐

locity & direction and suspended-sediment concentrations at 10 sites (Figure 5).

**Figure 5.** Tidal gauges (red square) and flow observations sites (green triangular) in Qinzhou Bay.

interpolated around islands and in the estuarine channels.

338 Sediment Transport Processes and Their Modelling Applications

The time step for both of external mode and internal mode is 1 s.

**5.2. Boundary conditions**

sediment discharges.

**6. Calibrations**

mean discharges as input conditions.

The Calibrations for water levels at three tidal gauges are shown in Figure 6.

**Figure 6.** Tidal water level calibrations (the ordinate unit of is meter).

#### **6.2. Flow calibrations**

The flow calibrations are shown in Figure 7. For limited space only calibrations for spring tide are show here.

**Figure 7.** Flow calibrations.

#### **6.3. Sediment calibration**

Sediment calibrations for selected sites are shown on Figure 8. It needs to explain that the overall sediment calibrations are satisfactory, but results for some sites are not satisfactory enough due to observation and computation errors, ship activities and dredging during the observation period, etc.

**Figure 8.** Suspended-sediment calibrations at selected sites for spring tide (the ordinate unit is kg/m3).

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment

http://dx.doi.org/10.5772/51760

341

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment http://dx.doi.org/10.5772/51760 341

**Figure 8.** Suspended-sediment calibrations at selected sites for spring tide (the ordinate unit is kg/m3).

**Figure 7.** Flow calibrations.

**6.3. Sediment calibration**

340 Sediment Transport Processes and Their Modelling Applications

the observation period, etc.

Sediment calibrations for selected sites are shown on Figure 8. It needs to explain that the overall sediment calibrations are satisfactory, but results for some sites are not satisfactory enough due to observation and computation errors, ship activities and dredging during

#### **6.4. Morphological calibration**

Based upon a local bathymetric survey conducted in 2008 summer and a sea map sur‐ veyed in 2004, statistics shows that the total siltation amount at this area is about 1.52 million m3 with a spatial average value of 0.592m (the deposition volume is divided by the deposition area, the same hereinafter), the total erosion amount is nearly 3.35 million m3 with a spatial average value of 1.104m (the erosion volume is divided by the erosion area, the same hereinafter), the net eroded sediment amounts to 1.83 million m3 (Figure 9).

**Figure 10.** The computed 2004~2008 erosion & deposition distribution map.

ern direction to the South China Sea (Figure 11).

As for neap tide and middle tide, tidal water floods into the Qinzhou Bay in the northeast‐ ern direction while it floods into the Bay nearly in the northern direction during spring tide. As constrained by the shoreline and channel-shoal geomorphology tidal water propagates in the northwestern direction into the Yingling inlet and further into the inner bay, where it flows anticlockwise till stack water; then tidal water rushes into Yingling inlet and diverges among east, middle and west waterways; finally it leaves the Qinzhou Bay in the southwest‐

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment

http://dx.doi.org/10.5772/51760

343

Generally speaking, the tidal flow field in the Qinzhou Bay is characteristics of reciprocating flow in parallel to the major waterways, large scale eddies occur during flow reversal periods.

The computation results show that the ebb-mean velocities of spring tide and middle tide are all larger than the flood-mean velocities; while flood-mean velocity becomes larger than

ebb-mean velocity during neap tide in the overall flow field of Qinzhou Bay.

**7. Computation results**

**7.1. Tidal flow field**

**Figure 9.** Erosion & deposition distribution map in a part of Qinzhou Bay.

The present study computes the morphological evolution using the 2004 bathymetry as the initial bathymetry, the computed four-year accumulative erosion & siltation distribution is shown in Figure 10. The computed total siltation amount is 1.244 million m3 with a spatial average value of 0.461m, amounting 81.89% and 77.87% of surveyed quantities, respectively; the computed total erosion amount is 2.97 million m3 with a spatial average value of 1.073m, amounting 88.53% and 97.19% of surveyed quantities, respectively; the computed net erod‐ ed sediment amounts to 1.72 million m3 , amounting 94.03% of surveyed quantities.

In view of the discrepancies of computed results vs. surveyed quantities the present mor‐ phological calibration is quite satisfactory. This lays down very good basis for further mor‐ phodynamic study.

**Figure 10.** The computed 2004~2008 erosion & deposition distribution map.

## **7. Computation results**

#### **7.1. Tidal flow field**

**6.4. Morphological calibration**

342 Sediment Transport Processes and Their Modelling Applications

**Figure 9.** Erosion & deposition distribution map in a part of Qinzhou Bay.

ed sediment amounts to 1.72 million m3

phodynamic study.

million m3

Based upon a local bathymetric survey conducted in 2008 summer and a sea map sur‐ veyed in 2004, statistics shows that the total siltation amount at this area is about 1.52

deposition area, the same hereinafter), the total erosion amount is nearly 3.35 million m3 with a spatial average value of 1.104m (the erosion volume is divided by the erosion area, the same hereinafter), the net eroded sediment amounts to 1.83 million m3 (Figure 9).

The present study computes the morphological evolution using the 2004 bathymetry as the initial bathymetry, the computed four-year accumulative erosion & siltation distribution is

average value of 0.461m, amounting 81.89% and 77.87% of surveyed quantities, respectively; the computed total erosion amount is 2.97 million m3 with a spatial average value of 1.073m, amounting 88.53% and 97.19% of surveyed quantities, respectively; the computed net erod‐

In view of the discrepancies of computed results vs. surveyed quantities the present mor‐ phological calibration is quite satisfactory. This lays down very good basis for further mor‐

, amounting 94.03% of surveyed quantities.

with a spatial

shown in Figure 10. The computed total siltation amount is 1.244 million m3

with a spatial average value of 0.592m (the deposition volume is divided by the

As for neap tide and middle tide, tidal water floods into the Qinzhou Bay in the northeast‐ ern direction while it floods into the Bay nearly in the northern direction during spring tide. As constrained by the shoreline and channel-shoal geomorphology tidal water propagates in the northwestern direction into the Yingling inlet and further into the inner bay, where it flows anticlockwise till stack water; then tidal water rushes into Yingling inlet and diverges among east, middle and west waterways; finally it leaves the Qinzhou Bay in the southwest‐ ern direction to the South China Sea (Figure 11).

Generally speaking, the tidal flow field in the Qinzhou Bay is characteristics of reciprocating flow in parallel to the major waterways, large scale eddies occur during flow reversal periods.

The computation results show that the ebb-mean velocities of spring tide and middle tide are all larger than the flood-mean velocities; while flood-mean velocity becomes larger than ebb-mean velocity during neap tide in the overall flow field of Qinzhou Bay.

**Figure 11.** a) Flood peak flow of spring tide in Qinzhou Bay, b) Ebb peak flow of spring tide in Qinzhou Bay.

**Figure 12.** a) Upper-layer flood peak flow of spring tide in Qinzhou Bay, b) Middle-layer flood peak flow of spring tide

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment

http://dx.doi.org/10.5772/51760

345

in Qinzhou Bay, c) Lower-layer flood peak flow of spring tide in Qinzhou Bay.

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment http://dx.doi.org/10.5772/51760 345

**Figure 12.** a) Upper-layer flood peak flow of spring tide in Qinzhou Bay, b) Middle-layer flood peak flow of spring tide in Qinzhou Bay, c) Lower-layer flood peak flow of spring tide in Qinzhou Bay.

**Figure 11.** a) Flood peak flow of spring tide in Qinzhou Bay, b) Ebb peak flow of spring tide in Qinzhou Bay.

344 Sediment Transport Processes and Their Modelling Applications

Among three major tidal channels in the outer bay, the flood-mean and ebb-mean velocities of spring, middle, neap tides in the middle channel are all larger than those in the east and west channels; though flood-mean velocity in spring tide in the west channel is somewhat smaller than that in the east channel the flood-mean velocities in middle & neap tides in the west channel are all larger than those in the east channel; as for ebb-mean velocity, they are all larger in the west channel than those in the east channel. These data demonstrates that the west channel is the dominant channel for tidal water flowing into and out the Qinzhou Bay, the middle channel comes second and the east channel is third; ebb tide dominates in the west channel and middle channel but flood tide dominates in the east channel.

suspended in situ in the Bay. As a result, the spatial & temporal variations of suspended-

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment

Generally speaking, suspended-sediment concentrations are larger in major channels than those on shoals and intertidal zones (Figure 14), decrease from bottom to top. The majority of sediment delivered by Maolingjiang and Qinjiang rivers is deposited within the inner bay

**Figure 14.** a) Suspended-sediment concentration field at flood peak of spring tide, b) Suspended-sediment concentra‐

The computation results show that suspended-sediment concentration at spring tidal flood

at middle tidal ebb stack, respectively, the mean middle-tidal concentration is 0.026 kg/m3

The west channel has been the major channel to transport sediment from outer bay thought Yingling inlet to inner bay and vice versa; the middle channel comes second and the east

at neap tidal ebb peak and 0.012 kg/m3

; sediment concentration is 0.031 kg/m3

at spring tidal flood stack, 0.034 kg/m3 at spring tidal ebb

at neap tidal flood peak, 0.013 kg/m3

at middle tidal flood peak,

http://dx.doi.org/10.5772/51760

347

at neap tidal ebb stack, re‐

. Generally speaking, suspend‐

;

at neap tidal

at spring tidal ebb stack, respectively, the mean spring-tidal concen‐

at middle tidal flood stack, 0.031 kg/m3 at middle tidal ebb peak and 0.018 kg/m3

tion field at ebb peak of middle tide.

, 0.021 kg/m3

spectively, The mean neap-tidal concentration is 0.013 kg/m3

ed-sediment concentration is relatively low in the Qinzhou Bay.

the sediment concentration is 0.013 kg/m3

peak is 0.037 kg/m3

peak and 0.023 kg/m3

tration is 0.029 kg/m3

flood stack, 0.013kg/m3

0.023 kg/m3

with limited amount of finer particles transported into the outer bay and deep water.

sediment concentrations are in accordance with the processes of tidal flows.

Generally speaking, flow velocity at the upper water layer is the largest and decreases from top to bottom (Figure 12). The depth-averaged residual flow is shown in Figure 13. Various-sized residual eddies occur in the Qinzhou Bay. Mean residual flow velocity in the Yingling inlet is around 0.15m/s with largest velocity of 0.489m/s, it is generally below 0.05m/s in other parts.

**Figure 13.** Depth-averaged residual flow in Qinzhou Bay.

#### **7.2. Sediment transport**

The suspended-sediment sources include those delivered by Maolingjiang river and Qinjing river and limited amount transported from the open sea, but the majority is eroded and re‐ suspended in situ in the Bay. As a result, the spatial & temporal variations of suspendedsediment concentrations are in accordance with the processes of tidal flows.

Among three major tidal channels in the outer bay, the flood-mean and ebb-mean velocities of spring, middle, neap tides in the middle channel are all larger than those in the east and west channels; though flood-mean velocity in spring tide in the west channel is somewhat smaller than that in the east channel the flood-mean velocities in middle & neap tides in the west channel are all larger than those in the east channel; as for ebb-mean velocity, they are all larger in the west channel than those in the east channel. These data demonstrates that the west channel is the dominant channel for tidal water flowing into and out the Qinzhou Bay, the middle channel comes second and the east channel is third; ebb tide dominates in

Generally speaking, flow velocity at the upper water layer is the largest and decreases from top to bottom (Figure 12). The depth-averaged residual flow is shown in Figure 13. Various-sized residual eddies occur in the Qinzhou Bay. Mean residual flow velocity in the Yingling inlet is around 0.15m/s with largest velocity of 0.489m/s, it is generally below 0.05m/s in other parts.

The suspended-sediment sources include those delivered by Maolingjiang river and Qinjing river and limited amount transported from the open sea, but the majority is eroded and re‐

the west channel and middle channel but flood tide dominates in the east channel.

346 Sediment Transport Processes and Their Modelling Applications

**Figure 13.** Depth-averaged residual flow in Qinzhou Bay.

**7.2. Sediment transport**

Generally speaking, suspended-sediment concentrations are larger in major channels than those on shoals and intertidal zones (Figure 14), decrease from bottom to top. The majority of sediment delivered by Maolingjiang and Qinjiang rivers is deposited within the inner bay with limited amount of finer particles transported into the outer bay and deep water.

**Figure 14.** a) Suspended-sediment concentration field at flood peak of spring tide, b) Suspended-sediment concentra‐ tion field at ebb peak of middle tide.

The computation results show that suspended-sediment concentration at spring tidal flood peak is 0.037 kg/m3 , 0.021 kg/m3 at spring tidal flood stack, 0.034 kg/m3 at spring tidal ebb peak and 0.023 kg/m3 at spring tidal ebb stack, respectively, the mean spring-tidal concen‐ tration is 0.029 kg/m3 ; sediment concentration is 0.031 kg/m3 at middle tidal flood peak, 0.023 kg/m3 at middle tidal flood stack, 0.031 kg/m3 at middle tidal ebb peak and 0.018 kg/m3 at middle tidal ebb stack, respectively, the mean middle-tidal concentration is 0.026 kg/m3 ; the sediment concentration is 0.013 kg/m3 at neap tidal flood peak, 0.013 kg/m3 at neap tidal flood stack, 0.013kg/m3 at neap tidal ebb peak and 0.012 kg/m3 at neap tidal ebb stack, re‐ spectively, The mean neap-tidal concentration is 0.013 kg/m3 . Generally speaking, suspend‐ ed-sediment concentration is relatively low in the Qinzhou Bay.

The west channel has been the major channel to transport sediment from outer bay thought Yingling inlet to inner bay and vice versa; the middle channel comes second and the east channel contributes the least. The sediment discharges at ebb tide are all larger than those at flood tide in the three channels, demonstrating net sediment transport into the open deep water. The sediment discharges at spring tide in three channels are all larger than those at middle & neap tide, and the dividing ratio of west channel, middle channel and east channel is nearly 5:2:1 for spring flood tide and 7:5:1 for spring ebb tide, 4:2:1 for middle flood tide and 8:4:1 for middle ebb tide, respectively.

The sediment transport pattern is normally accordance with tidal flow asymmetry in three channels, that is, ebb flow strength and discharge are superior to flood flow strength and discharge.

## **7.3. Morphological evolution**

### *7.3.1. 2009-year erosion and deposition*

Due to lack of data on the deposit thickness distribution in Qinzhou Bay and considering that rock is exposed locally within the deep channel in Yingling inlet by strong tidal flows [8], the present study assumes the initial deposit thickness is 0.3m within the deep channel in Yingling inlet and 20m in other parts of the Qinzhou Bay. The morphological computa‐ tion starts from year 2008.

The computed 2009 annual erosion & deposition distribution map is shown in Figure 15. It can be observed that erosions mainly occur within channels including three major chan‐ nels in the outer bay, deep-water channel in the Yingling inlet and those in the inner bay while depositions occur at shoal & ridge area and at the end of channels. This asserts that tide flow is really the dominant force for maintaining and reforming the channel-shoal morphology in the Qinzhou Bay.

Generally speaking, eroded sediments exceed deposited sediment for the whole Qinzhou Bay with net erosion nearly up to 10.288 million m3 . Except for the inner bay where net dep‐ osition occurs with a quantity of 3.190 million m3 net erosions all occur in the Yingling inlet and the outer bay, they are 2.999 million m3 and 3.503 million m3 , respectively. Due to lack sediment supply, the offshore slope outside the Qinzhou Bay is also subjected to net erosion of 6.976 million m3 , where erosion mainly occurs at the middle and southeastern part while deposition occurs at the southwestern part.

**Figure 15.** 2009-year annual erosion & deposition distribution map.

and 4.161 million m3

hou Bay is still subjected to net erosion of 8.331 million m3

erosion is roughly 1,488,779.676 m3

*7.3.2. 2012-year erosion and deposition*

erosion nearly up to 10.469 million m3

erosion is roughly 2,031,819.599 m3

tion of 2.832 million m3

are 0.809 million m3

0.623m and the maximum erosion is -1.285m.

The total deposition in the east channel in the outer bay is roughly 499,259.028 m3

The total deposition in the west channel in the outer bay is roughly 240,510.501 m3

mean deposition is 0.110m, the spatial mean erosion is -0.213m, the maximum deposition is

The computed 2012-year annual erosion & deposition distribution map is shown in Figure 16. Eroded sediments still exceed deposited sediment for the whole Qinzhou Bay with net

, and the net erosion is about 989,520.6 m3

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment

, net erosions still occur in the Yingling inlet and the outer bay, they

.

, and the net erosion is about 1,791,309.0 m3

. The inner bay continues to accommodate net deposi‐

, respectively; the offshore slope outside the Qinz‐

, the total

, the total

; the spatial

; the spatial

http://dx.doi.org/10.5772/51760

349

The total deposition in the west channel (bounded by -5m bathymetric contour, the same for middle channel and east channel) in the outer bay is roughly 490,887.486 m3 , the total ero‐ sion is roughly 2,257,125.612 m3 , and the net erosion is about 1,766,238.1 m3 ; the spatial mean deposition is 0.097m, the spatial mean erosion is -0.232m, the maximum deposition is 0.355m and the maximum erosion is -1.244m.

The total deposition in the middle channel in the outer bay is roughly 9,830.569 m3 , the total erosion is roughly 551,595.451 m3 , and the net erosion is about 453,285.88 m3 ; the spatial mean deposition is 0.073m, the spatial mean erosion is -0.394m, the maximum deposition is 0.157m and the maximum erosion is -1.785m. The erosion mainly occurs at the channel mouth connecting with the Yingling inlet.

**Figure 15.** 2009-year annual erosion & deposition distribution map.

The total deposition in the east channel in the outer bay is roughly 499,259.028 m3 , the total erosion is roughly 1,488,779.676 m3 , and the net erosion is about 989,520.6 m3 ; the spatial mean deposition is 0.110m, the spatial mean erosion is -0.213m, the maximum deposition is 0.623m and the maximum erosion is -1.285m.

## *7.3.2. 2012-year erosion and deposition*

channel contributes the least. The sediment discharges at ebb tide are all larger than those at flood tide in the three channels, demonstrating net sediment transport into the open deep water. The sediment discharges at spring tide in three channels are all larger than those at middle & neap tide, and the dividing ratio of west channel, middle channel and east channel is nearly 5:2:1 for spring flood tide and 7:5:1 for spring ebb tide, 4:2:1 for middle flood tide

The sediment transport pattern is normally accordance with tidal flow asymmetry in three channels, that is, ebb flow strength and discharge are superior to flood flow

Due to lack of data on the deposit thickness distribution in Qinzhou Bay and considering that rock is exposed locally within the deep channel in Yingling inlet by strong tidal flows [8], the present study assumes the initial deposit thickness is 0.3m within the deep channel in Yingling inlet and 20m in other parts of the Qinzhou Bay. The morphological computa‐

The computed 2009 annual erosion & deposition distribution map is shown in Figure 15. It can be observed that erosions mainly occur within channels including three major chan‐ nels in the outer bay, deep-water channel in the Yingling inlet and those in the inner bay while depositions occur at shoal & ridge area and at the end of channels. This asserts that tide flow is really the dominant force for maintaining and reforming the channel-shoal

Generally speaking, eroded sediments exceed deposited sediment for the whole Qinzhou

sediment supply, the offshore slope outside the Qinzhou Bay is also subjected to net erosion

The total deposition in the west channel (bounded by -5m bathymetric contour, the same for

deposition is 0.097m, the spatial mean erosion is -0.232m, the maximum deposition is

mean deposition is 0.073m, the spatial mean erosion is -0.394m, the maximum deposition is 0.157m and the maximum erosion is -1.785m. The erosion mainly occurs at the channel

middle channel and east channel) in the outer bay is roughly 490,887.486 m3

The total deposition in the middle channel in the outer bay is roughly 9,830.569 m3

and 3.503 million m3

, where erosion mainly occurs at the middle and southeastern part while

, and the net erosion is about 1,766,238.1 m3

, and the net erosion is about 453,285.88 m3

. Except for the inner bay where net dep‐

, respectively. Due to lack

, the total ero‐

, the total

; the spatial

; the spatial mean

net erosions all occur in the Yingling inlet

and 8:4:1 for middle ebb tide, respectively.

348 Sediment Transport Processes and Their Modelling Applications

strength and discharge.

**7.3. Morphological evolution**

tion starts from year 2008.

morphology in the Qinzhou Bay.

of 6.976 million m3

Bay with net erosion nearly up to 10.288 million m3

osition occurs with a quantity of 3.190 million m3

and the outer bay, they are 2.999 million m3

deposition occurs at the southwestern part.

0.355m and the maximum erosion is -1.244m.

mouth connecting with the Yingling inlet.

sion is roughly 2,257,125.612 m3

erosion is roughly 551,595.451 m3

*7.3.1. 2009-year erosion and deposition*

The computed 2012-year annual erosion & deposition distribution map is shown in Figure 16. Eroded sediments still exceed deposited sediment for the whole Qinzhou Bay with net erosion nearly up to 10.469 million m3 . The inner bay continues to accommodate net deposi‐ tion of 2.832 million m3 , net erosions still occur in the Yingling inlet and the outer bay, they are 0.809 million m3 and 4.161 million m3 , respectively; the offshore slope outside the Qinz‐ hou Bay is still subjected to net erosion of 8.331 million m3 .

The total deposition in the west channel in the outer bay is roughly 240,510.501 m3 , the total erosion is roughly 2,031,819.599 m3 , and the net erosion is about 1,791,309.0 m3 ; the spatial mean deposition is 0.063m, the spatial mean erosion is -0.185m, the maximum deposition is 0.232m and the maximum erosion is -0.889m.

channel has increased substantially, leading almost to whole-channel erosion, and tidal

Numerical Modeling Tidal Circulation and Morphodynamics in a Dumbbell-Shaped Coastal Embayment