# Numerical Simulation of Flow and Suspended Sediment Transport in the Distributary Channel Networks

Authors: Wei Zhang; Qiong Jia; Xiaowen Chen
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Flow and suspended sediment transport in distributary channel networks play an important role in the evolution of deltas and estuaries, as well as the coastal environment. In this study, a 1D flow and suspended sediment transport model is presented to simulate the hydrodynamics and suspended sediment transport in the distributary channel networks. The governing equations for river flow are the Saint-Venant equations and for suspended sediment transport are the nonequilibrium transport equations. The procedure of solving the governing equations is firstly to get the matrix form of the water level and suspended sediment concentration at all connected junctions by utilizing the transformation of the governing equations of the single channel. Secondly, the water level and suspended sediment concentration at all junctions can be obtained by solving these irregular spare matrix equations. Finally, the water level, discharge, and suspended sediment concentration at each river section can be calculated. The presented 1D flow and suspended sediment transport model has been applied to the Pearl River networks and can reproduce water levels, discharges, and suspended sediment concentration with good accuracy, indicating this that model can be used to simulate the hydrodynamics and suspended sediment concentration in the distributary channel networks.


---

## Body

## 1. Introduction

Distributary channels are very common in the delta area, among which many junctions connect individual channel together to form the river networks system. Channel junctions not only play a crucial role in dividing the river discharge over distributary channels but also govern the distribution of the river sediment transportation in the form of suspended sediment. Since flow and suspended sediment transport in the distributary channels are not only highly related to the evolution of deltas and estuaries but are also related to the environment issue, it is important to understand hydrodynamics and suspended sediment transport in the distributary channels system. The 1D numerical model as an important tool has been widely used to simulate flow and sediment transport over river networks because it can accurately represent complex river network geometries and complex structures like weirs, barrages, and dams [1].In the late 1970’s, Thomas and Prashum [2] formulated the Hydraulic Engineering Center (HEC-6) model in a rectilinear coordinate and described it by using finite-difference schemes. This model was introduced for sediment movement under steady flows in gravel-bed rivers. Then, Karim and Kennedy [3] and Karim [4] developed the Iowa ALLUVIAL (IALLUVIAL) model on the basis of a rectilinear coordinate system and the Saint-Venant flow equations in a quasisteady flow. The obvious difference between IALLUVIAL model and the HEC-6 model is that the HEC-6 solves the differential conservation equation of energy instead of the momentum equation. There are many unsteady flow models (e.g., Chang [5], Molinas and Yang [6], Holly and Rahuel [7], and Runkel and Broshears [8]), which have been developed and applied to river networks and other situations where the unsteady flow prevails over the steady one. Different from rectilinear coordinate, Chang [5] used a curvilinear coordinate system to solve the governing equations of his model FLUVIAL 11. The model accounts for the presence of secondary currents in a curved channel by adjusting the magnitude of the streamwise velocity. Furthermore, Molinas and Yang [6] implemented the theory of minimum stream power to determine the optimum channel width and geometry for a given set of hydraulic and sediment conditions. But the OTIS (one-dimensional transport with inflow and storage) [8] modified the 1D advection diffusion equation with additional terms to account for lateral inflow, first-order decay, sorption of nonconservative solutes, and transient storage of these solutes. Moreover, considering unsteady flow conditions that occur over transcritical flow stream, Papanicolaou et al. [9] developed steep stream sediment transport 1D model (3ST1D) model which is capable of capturing hydraulic jumps and simulating supercritical flows; therefore, it is applicable to flows over step-pool sequences in mountain streams. Han [10] firstly provided a method for nonequilibrium transport of nonuniform suspended load. Van Niekerk et al. [11] simulate deposition and erosion in alluvial channels. Rahuel et al. [12] developed a model to simulate the unsteady water and sediment movement in alluvial rivers. Wu et al. [13] further refined this model to simulate the nonequilibrium transport of nonuniform total load in in the distributary channel networks. Fang et al. [14] presented a one-dimensional numerical model for unsteady flow and nonuniform sediment transport in alluvial rivers. The model calculated the unsteady flow in open channels using the Preissmann implicit method and the Thomas algorithm.Most of the 1D models, however, are widely used to simulate the basic parameters, including the velocity, water level, bed elevation, and suspended sediment transport, in a single channel. Actually, 1D models can also be used to investigate hydrodynamics and sediment transport in distributary channel networks with complex planimetric structure and dynamics [15, 16]. Hence, model representation of network connections and river reach geometry has been heretofore defined by modelers to meet the needs of individual hydraulic and hydrologic models (Liu and Hodges) [17]. So far, river network systems used to be simulated by 1D model, described by Saint-Venant equations, and then solved together with other common 1D rivers. With this kind of models, it takes special experience and much work to divide and simplify the computation domains (Zhang et al.) [18]. What is more, compared with the traditional models, because of their low data, central processor unit (CPU) requirements, and simplicity of use, the presented 1D flow and suspended sediment model is very efficient and can be used in river engineering for long-term predictions of river morphological development. Therefore, the objective of this paper is to (
1
) develop a junction-control method for water level and suspended sediment concentration to simulate unsteady flow and suspended sediment transport over distributary channel networks and to (2) apply this model to the Pearl River networks to check the applicability of this model.


## 2. Mathematical Model and Numerical Approach

### 2.1. Governing Equations

The hydrodynamic module of 1D model is based on the well-known Saint-Venant equations of the mass and momentum conservation. Consider(1)
∂
Q
∂
x
+
B
∂
Z
∂
t
=
q
,
∂
Q
∂
t
+
2
u
∂
Q
∂
x
+
(
g
A
-
B
u
2
)
∂
Z
∂
x
-
u
2
∂
A
∂
x
+
g
n
2
|
u
|
Q
R
4
/
3
=
0
.Suspended sediment model is based on the nonequilibrium transport equation. The continuity equation of suspended load can be expressed as(2)
∂
(
A
S
)
∂
t
+
∂
(
Q
S
)
∂
x
=
-
α
B
ω
(
S
-
S
*
)
.The formula for sediment carrying capacity can be expressed as(3)
S
*
=
K
(
u
3
g
R
ω
)
m
,
where x and t are spatial and temporal axes, Q is flow discharge, A is flow area, q is the lateral flow, g is gravitational acceleration, R is the hydraulic radius, and Z is water surface elevation, B is water surface width, u is water velocity, n is Manning roughness coefficient, S is sediment concentration, α is adjustment coefficient for a cross-section, ω is fall velocity of sediment particles, S
* is averaged sediment carrying capacity, and K and m are the coefficients of sediment carrying capacity.


### 2.2. Solution for Flow Equations

#### 2.2.1. The Finite-Difference Equations

The finite-difference equations are obtained by using Preissmann’s four-point implicit difference scheme (Figure1). Consider
(4)
C
i
Z
i
n
+
1
+
C
i
Z
i
+
1
n
+
1
-
Q
i
n
+
1
+
Q
i
+
1
n
+
1
=
D
i
,
E
i
Q
i
n
+
1
+
G
i
Q
i
+
1
n
+
1
-
F
i
Z
i
n
+
1
+
F
i
Z
i
+
1
n
+
1
=
H
i
,
where C
i, D
i, E
i, G
i, F
i, and H
i are the coefficients of the difference equations as below:
(5)
C
i
=
B
w
i
+
1
/
2
n
Δ
x
i
2
θ
Δ
t
,
B
w
i
+
1
/
2
n
=
0.5
[
B
w
i
n
+
B
w
i
+
1
n
]
,
D
i
=
(
1
-
θ
)
θ
(
Q
i
n
-
Q
i
+
1
n
)
+
C
i
(
Z
i
n
+
Z
i
+
1
n
)
+
q
i
Δ
x
i
θ
,
E
i
=
Δ
x
i
2
θ
Δ
t
-
2
u
i
+
1
/
2
n
+
g
Δ
x
i
2
θ
(
n
2
|
u
|
R
4
/
3
)
i
n
,
F
i
=
(
g
A
-
B
u
2
)
i
+
1
/
2
n
=
g
A
i
+
1
/
2
n
-
B
i
+
1
/
2
n
·
u
i
+
1
/
2
2
n
,
G
i
=
Δ
x
i
2
θ
Δ
t
+
2
u
i
+
1
/
2
n
+
g
Δ
x
i
2
θ
(
n
2
|
u
|
R
4
/
3
)
i
+
1
n
,
H
i
=
Δ
x
i
2
θ
Δ
t
(
Q
i
n
+
Q
i
+
1
n
)
+
2
(
1
-
θ
)
θ
u
i
+
1
/
2
n
(
Q
i
n
-
Q
i
+
1
n
)
+
(
1
-
θ
)
θ
(
g
A
-
B
u
2
)
i
+
1
/
2
n
(
Z
i
n
-
Z
i
+
1
n
)
+
Δ
x
i
θ
(
u
2
∂
A
∂
x
)
i
+
1
/
2
n
,
where θ  (
0.5
≤
θ
≤
1
) is the weighting coefficient and
(6)
(
u
2
∂
A
∂
x
)
i
+
1
/
2
n
=
(
u
i
+
1
/
2
n
)
2
A
i
+
1
n
-
A
i
n
Δ
x
i
.Figure 1
Illustration of the flow on grid and network domains.


#### 2.2.2. Supplementary Equation at Junctions

The hydraulic conditions at river networks junctions are reasonably assumed to satisfy two additional conditions: the quality conservation and energy conservation conditions. The quality conservation condition, which is also called discharge connection condition, means that the discharge changes are assumed to be equal to the shortage variation of discharge at a certain junction. Consider(7)
∑
i
=
1
m
Q
i
=
Q
I
n
+
(
1
/
2
)
-
Q
O
n
+
(
1
/
2
)
=
(
Ω
k
n
+
1
-
Ω
k
n
)
Δ
t
,
where k is the junction number, m is the number of incoming (outflow) rivers at the k
th junction, Q
i is the discharge of the i
th river course, Ω
k is the storage of the k
th junction, Q
I is the incoming discharge, and Q
O is the outflow discharge.The equation of energy conservation can be approximated using the following kinematic compatibility condition [19]:
(8)
Z
in
,
k
=
Z
out
,
k
=
H
k
,
where Z
in
,
k and Z
out
,
k represent the water level of any of the incoming and outflowing rivers at the i
th junction, respectively, and H
k is the water level at the i
th junction.Finally, the finite-difference approximation of the continuity equation around the junction can be written as follows:(9)
H
k
n
+
1
-
H
k
n
Δ
t
A
J
=
Q
I
n
+
(
1
/
2
)
-
Q
O
n
+
(
1
/
2
)
,
where A
J is the cross-section area of each junction.The values for water level can be calculated at each junction by successively applying the double sweep Thomas algorithm.For each river section, based on the finite-difference equation (4), the equation of each subsection can be eliminated autocorrelation; then (4), can be rewritten in the following form:
(10)
Q
i
=
α
i
+
β
i
Ζ
i
+
γ
i
Ζ
n
+
1
(
i
=
n
,
…
,
2,1
)
,
(11)
Q
i
+
1
=
ξ
i
+
1
+
ζ
i
+
1
Ζ
i
+
1
+
η
i
+
1
Ζ
1
(
i
=
1,2
,
…
,
n
)
,
where the coefficients in (10) are defined as follows:
(12)
α
n
=
H
n
-
G
n
D
n
E
n
+
G
n
,
β
n
=
G
n
C
n
+
F
n
E
n
+
G
n
,
γ
n
=
G
n
C
n
-
F
n
E
n
+
G
n
,
α
i
=
Y
1
(
H
i
-
G
i
α
i
+
1
)
-
Y
2
(
D
i
-
α
i
+
1
)
Y
1
E
i
+
Y
2
,
β
i
=
Y
2
C
i
+
Y
1
F
i
Y
1
E
i
+
Y
2
,
γ
i
=
γ
i
+
1
(
Y
2
-
Y
1
G
i
)
Y
1
E
i
+
Y
2
,
i
i
i
i
i
(
i
=
n
-
1
,
…
,
1
)
,
where
(13)
Y
1
=
C
i
+
β
i
+
1
,
Y
2
=
G
i
β
i
+
1
+
F
i
.The coefficients in (11) are defined as follows:
(14)
ξ
n
=
E
1
D
1
+
H
1
E
1
+
G
1
,
ζ
n
=
-
C
1
E
1
+
F
1
E
1
+
G
1
,
η
n
=
-
C
1
E
1
-
F
1
E
1
+
G
1
,
ξ
i
+
1
=
Y
2
(
D
i
+
ξ
i
)
+
Y
1
(
H
i
-
E
i
ξ
i
)
Y
2
+
Y
1
G
i
,
ζ
i
+
1
=
-
Y
2
C
i
+
Y
1
F
i
Y
2
+
Y
1
G
i
,
η
i
+
1
=
η
i
(
Y
2
-
Y
1
E
i
)
Y
2
+
Y
1
G
i
,
iiii
(
i
=
2,3
,
…
,
n
)
.
Y
1
=
ζ
i
-
G
i
,
Y
2
=
E
i
ζ
i
-
F
i
.Subsequently, the relationship between the river sections and the junction can be deduced from a single river channel by using the formula (10) and (11) as follows:
(15)
Q
1
=
α
1
+
β
1
Z
1
+
γ
1
Z
n
+
1
,
Q
n
+
1
=
ξ
n
+
1
+
ζ
n
+
1
Z
n
+
1
+
η
n
+
1
Z
1
.Substituting (15) at each junction into (7), water levels at junctions of river networks can be expressed as a matrix form:
(16)
f
i
(
Z
1
,
Z
2
,
…
,
Z
i
,
…
,
Z
M
)
=
0
(
i
=
1
,
2
,
…
,
M
)
,
where M is the total number of junction of the river networks.Finally, the water level and discharge at each river section of the river networks can be calculated.


### 2.3. Solution for Suspended Sediment Transport Equations

#### 2.3.1. The Finite-Difference Equations

For suspended sediment transport, the upwind finite-difference scheme and Preissmann’s four-point implicit difference scheme are adopted (Figure2). According to the flow direction, the difference equations of the governing equation (2) can be expressed as follows.Figure 2
Illustration of the suspended sediment on grid and network domains.(i) For the downflow,(17)
A
∂
S
∂
t
=
A
¯
i
-
1
/
2
S
i
n
+
1
+
S
i
-
1
n
+
1
-
S
i
n
-
S
i
-
1
n
2
Δ
t
,
Q
∂
S
∂
x
=
Q
¯
i
-
1
/
2
θ
(
S
i
n
+
1
-
S
i
-
1
n
+
1
)
+
(
1
-
θ
)
(
S
i
n
-
S
i
-
1
n
)
Δ
x
i
-
1
,
-
α
B
ω
(
S
-
S
*
)
=
-
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
(
S
i
n
+
1
-
S
*
i
n
+
1
)
-
1
2
α
i
-
1
n
+
1
B
i
-
1
n
+
1
ω
i
-
1
n
+
1
(
S
i
-
1
n
+
1
-
S
*
i
-
1
n
+
1
)
,
where
(18)
A
¯
i
-
1
/
2
=
1
2
(
A
i
n
+
1
+
A
i
-
1
n
+
1
)
,
Q
¯
i
-
1
/
2
=
1
2
(
Q
i
n
+
1
+
Q
i
-
1
n
+
1
)
.The irregular matrix equations are written as(19)
α
α
i
S
i
-
1
n
+
1
+
β
β
i
S
i
n
+
1
=
D
D
i
,
where
(20)
α
α
i
=
A
¯
i
-
1
/
2
2
Δ
t
-
Q
¯
i
-
1
/
2
θ
Δ
x
i
-
1
+
1
2
α
i
-
1
n
+
1
B
i
-
1
n
+
1
ω
i
-
1
n
+
1
,
β
β
i
=
A
¯
i
-
1
/
2
2
Δ
t
+
Q
¯
i
-
1
/
2
θ
Δ
x
i
-
1
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
,
D
D
i
=
A
¯
i
-
(
1
/
2
)
2
Δ
t
(
S
i
n
+
S
i
-
1
n
)
-
Q
¯
i
-
1
/
2
(
1
-
θ
)
Δ
x
i
-
1
(
S
i
n
-
S
i
-
1
n
)
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
S
*
i
n
+
1
+
1
2
α
i
-
1
n
+
1
B
i
-
1
n
+
1
ω
i
-
1
n
+
1
S
*
i
-
1
n
+
1
.(ii) For the upflow,(21)
A
∂
S
∂
t
=
A
¯
i
+
1
/
2
S
i
n
+
1
+
S
i
+
1
n
+
1
-
S
i
n
-
S
i
+
1
n
2
Δ
t
,
Q
∂
S
∂
x
=
Q
¯
i
+
1
/
2
θ
(
S
i
+
1
n
+
1
-
S
i
n
+
1
)
+
(
1
-
θ
)
(
S
i
+
1
n
-
S
i
n
)
Δ
x
i
,
-
α
B
ω
(
S
-
S
*
)
=
-
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
(
S
i
n
+
1
-
S
*
i
n
+
1
)
-
1
2
α
i
+
1
n
+
1
B
i
+
1
n
+
1
ω
i
+
1
n
+
1
(
S
i
+
1
n
+
1
-
S
*
i
+
1
n
+
1
)
,
where
(22)
A
¯
i
+
1
/
2
=
1
2
(
A
i
n
+
1
+
A
i
+
1
n
+
1
)
,
Q
¯
i
+
1
/
2
=
1
2
(
Q
i
n
+
1
+
Q
i
+
1
n
+
1
)
.The irregular matrix equations are written as(23)
β
β
i
S
i
n
+
1
+
r
r
i
S
i
+
1
n
+
1
=
D
D
i
,
where
(24)
β
β
i
=
A
¯
i
+
1
/
2
2
Δ
t
-
Q
¯
i
+
1
/
2
θ
Δ
x
i
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
,
r
r
i
=
A
¯
i
+
1
/
2
2
Δ
t
+
Q
¯
i
+
1
/
2
Δ
x
i
+
1
2
α
i
+
1
n
+
1
B
i
+
1
n
+
1
ω
i
+
1
n
+
1
,
D
D
i
=
A
¯
i
+
1
/
2
2
Δ
t
(
S
i
n
+
S
i
+
1
n
)
-
Q
¯
i
+
1
/
2
(
1
-
θ
)
Δ
x
i
(
S
i
+
1
n
-
S
i
n
)
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
S
*
i
n
+
1
+
1
2
α
i
+
1
n
+
1
B
i
+
1
n
+
1
ω
i
+
1
n
+
1
S
*
i
+
1
n
+
1
,
where α
α
i, β
β
i, γ
γ
i, and D
D
i are as the known coefficients of the equations.


#### 2.3.2. Solution Algorithm for a Single River

In the tidal channel, the flow direction changes with tidal level fluctuation. The flow types of a single channel can be classified as four types (Figure3): (a) downflow, Q
s
≥
0, Q
e
≥
0; (b) upflow, Q
s
<
0, Q
e
<
0; (c) faced-flow, Q
s
≥
0, Q
e
<
0; and (d) depart-flow, Q
s
<
0, Q
e
≥
0.Figure 3
The four types of the flow at a single river channel. The types are as follows: downflow, upflow, faced-flow, and depart-flow, from the left to the right.For the downflow, it is single direction flow. If the sediment concentration at the upboundary junction (S
1
n
+
1) is known, the suspended sediment concentration at each junction can be obtained by the following formula:
(25)
S
i
n
+
1
=
P
i
+
R
i
S
1
n
+
1
(
i
=
1
,
2
,
…
,
n
+
1
)
,
where P
i
,
R
i are the coefficients, and they can be gotten by (19), (25), and written as
(26)
P
1
=
0
,
R
1
=
1
,
P
i
D
D
i
-
α
α
i
P
i
-
1
β
β
i
,
R
i
=
-
α
α
i
R
i
-
1
β
β
i
i
i
i
(
i
=
2,3
,
…
,
n
+
1
)
.For the upflow, it is also single direction flow. If the sediment concentration at the downboundary junction (S
n
+
1
n
+
1) is known, the sediment concentration at each junction can be obtained by the following formula:
(27)
S
i
n
+
1
=
P
i
+
R
i
S
n
+
1
n
+
1
(
i
=
1,2
,
…
,
n
+
1
)
,
where P
i
,
R
i are the coefficients, and they can be gotten by (23), (27), and written as
(28)
P
n
+
1
=
0
,
R
n
+
1
=
1
,
P
i
=
D
D
i
-
γ
γ
i
P
i
+
1
β
β
i
,
R
i
=
-
γ
γ
i
R
i
+
1
β
β
i
i
i
(
i
=
n
,
…
,
2,1
)
.For the faced-flow, it is double direction flow. The single river can be divided into two single direction flow river segments. If the sediment concentrations at the up and down boundary junctions (S
1
n
+
1, S
n
+
1
n
+
1) are known, the suspended sediment concentration at each junction can be obtained by the formulas (25) and (27).For the depart-flow, it is another double direction flow. First, the position of the stagnant point (Q
k
=
0) must be determined, the suspended sediment concentration at the stagnant point (S
1
n
+
1) may be obtained by the transport equation, and then the single river can be divided into two single direction flow river segments which boundary is the stagnant point; finally, the suspended sediment concentration at each junction can be obtained by the formulas (25), (27). The position of the stagnant point is presented as the following formula:
(29)
Δ
x
m
k
=
|
Q
m
|
|
Q
m
|
+
Q
l
Δ
x
m
l
.And the suspended sediment concentration at the stagnant point is(30)
S
k
n
+
1
=
A
k
n
A
k
n
+
1
S
k
n
exp
⁡
(
-
α
k
n
+
1
B
k
n
+
1
ω
k
n
+
1
A
k
n
+
1
Δ
t
)
.This paper also presents a junction-control method for suspended sediment in river networks.Assuming that the erosion and deposition at the junction are small, the suspended sediment mass conservation equation is written as(31)
∑
Q
in
S
in
=
∑
Q
out
S
out
(
S
out
=
S
J
)
,
where S
J is sediment concentration at junctions.Utilizing (25) and (27), the relation formula of sediment concentration between upjunction (
i
=
1
) and downjunction (
i
=
n
+
1
) is obtained, and the control equations for all junctions’ sediment concentration can be written in matrix notation as
(32)
F
j
(
S
1
,
S
2
,
…
,
S
M
)
=
0
j
=
1,2
,
…
,
M
.Solving (31), the suspended sediment concentration at all river junctions can be gotten. Finally, by solving the single river equations (25) and (27), suspended sediment concentration at all channel sections of the distributary channel networks can be calculated.


### 2.4. Modeling Setup

The river networks in the Pearl River Delta (PRD) (Figure4) are selected to check the applicability of the 1D flow and suspended sediment model. The river network in the PRD is considered as one of the world’s most complicated river networks, with the density of channel length per unit at 0.68–1.07 km/km2. The area of the Pearl River networks cover approximately 9750 km2, with a total length of over 1600 km and a coastline extending over 450 km from east to west. The Pearl River networks drain three main tributaries, Dongjiang River, Beijiang River, and Xijiang River, pass through this complex river network, and finally flow into the South China Sea through eight subestuaries that are separately called Humen, Jiaomen, Hongqimen, Hengmen, Modaomen, Jitimen, Hutiaomen, and Yamen.Figure 4
Map of the Pear River networks and sampling stations.The Pearl River network is simplified as 328 single river channels, 13 boundary rivers, 5108 computational river sections with the interval varying from 0.2 to 3.0 km, and 216 river channel junctions. Upstream boundary conditions are given by discharge and suspended sediment concentration at the five upstream boundary junctions (Shizui, Gaoyao, Shijiao, Laoyagang, and Boluo). Tidal level fluctuation and suspended sediment concentration at the downstream outlets are used as downstream boundary conditions. The Manning coefficients in the Pearl rivers are taken to be 0.01–0.05.K and m; the coefficients of sediment carrying capacity are taken to be 0.005 and 0.92. The value of adjustment coefficient for a cross-section α is 0.25 (when deposition) or 1.0 (when erosion). Computations are performed using a time step of 1 minute.


## 2.1. Governing Equations

The hydrodynamic module of 1D model is based on the well-known Saint-Venant equations of the mass and momentum conservation. Consider(1)
∂
Q
∂
x
+
B
∂
Z
∂
t
=
q
,
∂
Q
∂
t
+
2
u
∂
Q
∂
x
+
(
g
A
-
B
u
2
)
∂
Z
∂
x
-
u
2
∂
A
∂
x
+
g
n
2
|
u
|
Q
R
4
/
3
=
0
.Suspended sediment model is based on the nonequilibrium transport equation. The continuity equation of suspended load can be expressed as(2)
∂
(
A
S
)
∂
t
+
∂
(
Q
S
)
∂
x
=
-
α
B
ω
(
S
-
S
*
)
.The formula for sediment carrying capacity can be expressed as(3)
S
*
=
K
(
u
3
g
R
ω
)
m
,
where x and t are spatial and temporal axes, Q is flow discharge, A is flow area, q is the lateral flow, g is gravitational acceleration, R is the hydraulic radius, and Z is water surface elevation, B is water surface width, u is water velocity, n is Manning roughness coefficient, S is sediment concentration, α is adjustment coefficient for a cross-section, ω is fall velocity of sediment particles, S
* is averaged sediment carrying capacity, and K and m are the coefficients of sediment carrying capacity.


## 2.2. Solution for Flow Equations

### 2.2.1. The Finite-Difference Equations

The finite-difference equations are obtained by using Preissmann’s four-point implicit difference scheme (Figure1). Consider
(4)
C
i
Z
i
n
+
1
+
C
i
Z
i
+
1
n
+
1
-
Q
i
n
+
1
+
Q
i
+
1
n
+
1
=
D
i
,
E
i
Q
i
n
+
1
+
G
i
Q
i
+
1
n
+
1
-
F
i
Z
i
n
+
1
+
F
i
Z
i
+
1
n
+
1
=
H
i
,
where C
i, D
i, E
i, G
i, F
i, and H
i are the coefficients of the difference equations as below:
(5)
C
i
=
B
w
i
+
1
/
2
n
Δ
x
i
2
θ
Δ
t
,
B
w
i
+
1
/
2
n
=
0.5
[
B
w
i
n
+
B
w
i
+
1
n
]
,
D
i
=
(
1
-
θ
)
θ
(
Q
i
n
-
Q
i
+
1
n
)
+
C
i
(
Z
i
n
+
Z
i
+
1
n
)
+
q
i
Δ
x
i
θ
,
E
i
=
Δ
x
i
2
θ
Δ
t
-
2
u
i
+
1
/
2
n
+
g
Δ
x
i
2
θ
(
n
2
|
u
|
R
4
/
3
)
i
n
,
F
i
=
(
g
A
-
B
u
2
)
i
+
1
/
2
n
=
g
A
i
+
1
/
2
n
-
B
i
+
1
/
2
n
·
u
i
+
1
/
2
2
n
,
G
i
=
Δ
x
i
2
θ
Δ
t
+
2
u
i
+
1
/
2
n
+
g
Δ
x
i
2
θ
(
n
2
|
u
|
R
4
/
3
)
i
+
1
n
,
H
i
=
Δ
x
i
2
θ
Δ
t
(
Q
i
n
+
Q
i
+
1
n
)
+
2
(
1
-
θ
)
θ
u
i
+
1
/
2
n
(
Q
i
n
-
Q
i
+
1
n
)
+
(
1
-
θ
)
θ
(
g
A
-
B
u
2
)
i
+
1
/
2
n
(
Z
i
n
-
Z
i
+
1
n
)
+
Δ
x
i
θ
(
u
2
∂
A
∂
x
)
i
+
1
/
2
n
,
where θ  (
0.5
≤
θ
≤
1
) is the weighting coefficient and
(6)
(
u
2
∂
A
∂
x
)
i
+
1
/
2
n
=
(
u
i
+
1
/
2
n
)
2
A
i
+
1
n
-
A
i
n
Δ
x
i
.Figure 1
Illustration of the flow on grid and network domains.


### 2.2.2. Supplementary Equation at Junctions

The hydraulic conditions at river networks junctions are reasonably assumed to satisfy two additional conditions: the quality conservation and energy conservation conditions. The quality conservation condition, which is also called discharge connection condition, means that the discharge changes are assumed to be equal to the shortage variation of discharge at a certain junction. Consider(7)
∑
i
=
1
m
Q
i
=
Q
I
n
+
(
1
/
2
)
-
Q
O
n
+
(
1
/
2
)
=
(
Ω
k
n
+
1
-
Ω
k
n
)
Δ
t
,
where k is the junction number, m is the number of incoming (outflow) rivers at the k
th junction, Q
i is the discharge of the i
th river course, Ω
k is the storage of the k
th junction, Q
I is the incoming discharge, and Q
O is the outflow discharge.The equation of energy conservation can be approximated using the following kinematic compatibility condition [19]:
(8)
Z
in
,
k
=
Z
out
,
k
=
H
k
,
where Z
in
,
k and Z
out
,
k represent the water level of any of the incoming and outflowing rivers at the i
th junction, respectively, and H
k is the water level at the i
th junction.Finally, the finite-difference approximation of the continuity equation around the junction can be written as follows:(9)
H
k
n
+
1
-
H
k
n
Δ
t
A
J
=
Q
I
n
+
(
1
/
2
)
-
Q
O
n
+
(
1
/
2
)
,
where A
J is the cross-section area of each junction.The values for water level can be calculated at each junction by successively applying the double sweep Thomas algorithm.For each river section, based on the finite-difference equation (4), the equation of each subsection can be eliminated autocorrelation; then (4), can be rewritten in the following form:
(10)
Q
i
=
α
i
+
β
i
Ζ
i
+
γ
i
Ζ
n
+
1
(
i
=
n
,
…
,
2,1
)
,
(11)
Q
i
+
1
=
ξ
i
+
1
+
ζ
i
+
1
Ζ
i
+
1
+
η
i
+
1
Ζ
1
(
i
=
1,2
,
…
,
n
)
,
where the coefficients in (10) are defined as follows:
(12)
α
n
=
H
n
-
G
n
D
n
E
n
+
G
n
,
β
n
=
G
n
C
n
+
F
n
E
n
+
G
n
,
γ
n
=
G
n
C
n
-
F
n
E
n
+
G
n
,
α
i
=
Y
1
(
H
i
-
G
i
α
i
+
1
)
-
Y
2
(
D
i
-
α
i
+
1
)
Y
1
E
i
+
Y
2
,
β
i
=
Y
2
C
i
+
Y
1
F
i
Y
1
E
i
+
Y
2
,
γ
i
=
γ
i
+
1
(
Y
2
-
Y
1
G
i
)
Y
1
E
i
+
Y
2
,
i
i
i
i
i
(
i
=
n
-
1
,
…
,
1
)
,
where
(13)
Y
1
=
C
i
+
β
i
+
1
,
Y
2
=
G
i
β
i
+
1
+
F
i
.The coefficients in (11) are defined as follows:
(14)
ξ
n
=
E
1
D
1
+
H
1
E
1
+
G
1
,
ζ
n
=
-
C
1
E
1
+
F
1
E
1
+
G
1
,
η
n
=
-
C
1
E
1
-
F
1
E
1
+
G
1
,
ξ
i
+
1
=
Y
2
(
D
i
+
ξ
i
)
+
Y
1
(
H
i
-
E
i
ξ
i
)
Y
2
+
Y
1
G
i
,
ζ
i
+
1
=
-
Y
2
C
i
+
Y
1
F
i
Y
2
+
Y
1
G
i
,
η
i
+
1
=
η
i
(
Y
2
-
Y
1
E
i
)
Y
2
+
Y
1
G
i
,
iiii
(
i
=
2,3
,
…
,
n
)
.
Y
1
=
ζ
i
-
G
i
,
Y
2
=
E
i
ζ
i
-
F
i
.Subsequently, the relationship between the river sections and the junction can be deduced from a single river channel by using the formula (10) and (11) as follows:
(15)
Q
1
=
α
1
+
β
1
Z
1
+
γ
1
Z
n
+
1
,
Q
n
+
1
=
ξ
n
+
1
+
ζ
n
+
1
Z
n
+
1
+
η
n
+
1
Z
1
.Substituting (15) at each junction into (7), water levels at junctions of river networks can be expressed as a matrix form:
(16)
f
i
(
Z
1
,
Z
2
,
…
,
Z
i
,
…
,
Z
M
)
=
0
(
i
=
1
,
2
,
…
,
M
)
,
where M is the total number of junction of the river networks.Finally, the water level and discharge at each river section of the river networks can be calculated.


## 2.2.1. The Finite-Difference Equations

The finite-difference equations are obtained by using Preissmann’s four-point implicit difference scheme (Figure1). Consider
(4)
C
i
Z
i
n
+
1
+
C
i
Z
i
+
1
n
+
1
-
Q
i
n
+
1
+
Q
i
+
1
n
+
1
=
D
i
,
E
i
Q
i
n
+
1
+
G
i
Q
i
+
1
n
+
1
-
F
i
Z
i
n
+
1
+
F
i
Z
i
+
1
n
+
1
=
H
i
,
where C
i, D
i, E
i, G
i, F
i, and H
i are the coefficients of the difference equations as below:
(5)
C
i
=
B
w
i
+
1
/
2
n
Δ
x
i
2
θ
Δ
t
,
B
w
i
+
1
/
2
n
=
0.5
[
B
w
i
n
+
B
w
i
+
1
n
]
,
D
i
=
(
1
-
θ
)
θ
(
Q
i
n
-
Q
i
+
1
n
)
+
C
i
(
Z
i
n
+
Z
i
+
1
n
)
+
q
i
Δ
x
i
θ
,
E
i
=
Δ
x
i
2
θ
Δ
t
-
2
u
i
+
1
/
2
n
+
g
Δ
x
i
2
θ
(
n
2
|
u
|
R
4
/
3
)
i
n
,
F
i
=
(
g
A
-
B
u
2
)
i
+
1
/
2
n
=
g
A
i
+
1
/
2
n
-
B
i
+
1
/
2
n
·
u
i
+
1
/
2
2
n
,
G
i
=
Δ
x
i
2
θ
Δ
t
+
2
u
i
+
1
/
2
n
+
g
Δ
x
i
2
θ
(
n
2
|
u
|
R
4
/
3
)
i
+
1
n
,
H
i
=
Δ
x
i
2
θ
Δ
t
(
Q
i
n
+
Q
i
+
1
n
)
+
2
(
1
-
θ
)
θ
u
i
+
1
/
2
n
(
Q
i
n
-
Q
i
+
1
n
)
+
(
1
-
θ
)
θ
(
g
A
-
B
u
2
)
i
+
1
/
2
n
(
Z
i
n
-
Z
i
+
1
n
)
+
Δ
x
i
θ
(
u
2
∂
A
∂
x
)
i
+
1
/
2
n
,
where θ  (
0.5
≤
θ
≤
1
) is the weighting coefficient and
(6)
(
u
2
∂
A
∂
x
)
i
+
1
/
2
n
=
(
u
i
+
1
/
2
n
)
2
A
i
+
1
n
-
A
i
n
Δ
x
i
.Figure 1
Illustration of the flow on grid and network domains.


## 2.2.2. Supplementary Equation at Junctions

The hydraulic conditions at river networks junctions are reasonably assumed to satisfy two additional conditions: the quality conservation and energy conservation conditions. The quality conservation condition, which is also called discharge connection condition, means that the discharge changes are assumed to be equal to the shortage variation of discharge at a certain junction. Consider(7)
∑
i
=
1
m
Q
i
=
Q
I
n
+
(
1
/
2
)
-
Q
O
n
+
(
1
/
2
)
=
(
Ω
k
n
+
1
-
Ω
k
n
)
Δ
t
,
where k is the junction number, m is the number of incoming (outflow) rivers at the k
th junction, Q
i is the discharge of the i
th river course, Ω
k is the storage of the k
th junction, Q
I is the incoming discharge, and Q
O is the outflow discharge.The equation of energy conservation can be approximated using the following kinematic compatibility condition [19]:
(8)
Z
in
,
k
=
Z
out
,
k
=
H
k
,
where Z
in
,
k and Z
out
,
k represent the water level of any of the incoming and outflowing rivers at the i
th junction, respectively, and H
k is the water level at the i
th junction.Finally, the finite-difference approximation of the continuity equation around the junction can be written as follows:(9)
H
k
n
+
1
-
H
k
n
Δ
t
A
J
=
Q
I
n
+
(
1
/
2
)
-
Q
O
n
+
(
1
/
2
)
,
where A
J is the cross-section area of each junction.The values for water level can be calculated at each junction by successively applying the double sweep Thomas algorithm.For each river section, based on the finite-difference equation (4), the equation of each subsection can be eliminated autocorrelation; then (4), can be rewritten in the following form:
(10)
Q
i
=
α
i
+
β
i
Ζ
i
+
γ
i
Ζ
n
+
1
(
i
=
n
,
…
,
2,1
)
,
(11)
Q
i
+
1
=
ξ
i
+
1
+
ζ
i
+
1
Ζ
i
+
1
+
η
i
+
1
Ζ
1
(
i
=
1,2
,
…
,
n
)
,
where the coefficients in (10) are defined as follows:
(12)
α
n
=
H
n
-
G
n
D
n
E
n
+
G
n
,
β
n
=
G
n
C
n
+
F
n
E
n
+
G
n
,
γ
n
=
G
n
C
n
-
F
n
E
n
+
G
n
,
α
i
=
Y
1
(
H
i
-
G
i
α
i
+
1
)
-
Y
2
(
D
i
-
α
i
+
1
)
Y
1
E
i
+
Y
2
,
β
i
=
Y
2
C
i
+
Y
1
F
i
Y
1
E
i
+
Y
2
,
γ
i
=
γ
i
+
1
(
Y
2
-
Y
1
G
i
)
Y
1
E
i
+
Y
2
,
i
i
i
i
i
(
i
=
n
-
1
,
…
,
1
)
,
where
(13)
Y
1
=
C
i
+
β
i
+
1
,
Y
2
=
G
i
β
i
+
1
+
F
i
.The coefficients in (11) are defined as follows:
(14)
ξ
n
=
E
1
D
1
+
H
1
E
1
+
G
1
,
ζ
n
=
-
C
1
E
1
+
F
1
E
1
+
G
1
,
η
n
=
-
C
1
E
1
-
F
1
E
1
+
G
1
,
ξ
i
+
1
=
Y
2
(
D
i
+
ξ
i
)
+
Y
1
(
H
i
-
E
i
ξ
i
)
Y
2
+
Y
1
G
i
,
ζ
i
+
1
=
-
Y
2
C
i
+
Y
1
F
i
Y
2
+
Y
1
G
i
,
η
i
+
1
=
η
i
(
Y
2
-
Y
1
E
i
)
Y
2
+
Y
1
G
i
,
iiii
(
i
=
2,3
,
…
,
n
)
.
Y
1
=
ζ
i
-
G
i
,
Y
2
=
E
i
ζ
i
-
F
i
.Subsequently, the relationship between the river sections and the junction can be deduced from a single river channel by using the formula (10) and (11) as follows:
(15)
Q
1
=
α
1
+
β
1
Z
1
+
γ
1
Z
n
+
1
,
Q
n
+
1
=
ξ
n
+
1
+
ζ
n
+
1
Z
n
+
1
+
η
n
+
1
Z
1
.Substituting (15) at each junction into (7), water levels at junctions of river networks can be expressed as a matrix form:
(16)
f
i
(
Z
1
,
Z
2
,
…
,
Z
i
,
…
,
Z
M
)
=
0
(
i
=
1
,
2
,
…
,
M
)
,
where M is the total number of junction of the river networks.Finally, the water level and discharge at each river section of the river networks can be calculated.


## 2.3. Solution for Suspended Sediment Transport Equations

### 2.3.1. The Finite-Difference Equations

For suspended sediment transport, the upwind finite-difference scheme and Preissmann’s four-point implicit difference scheme are adopted (Figure2). According to the flow direction, the difference equations of the governing equation (2) can be expressed as follows.Figure 2
Illustration of the suspended sediment on grid and network domains.(i) For the downflow,(17)
A
∂
S
∂
t
=
A
¯
i
-
1
/
2
S
i
n
+
1
+
S
i
-
1
n
+
1
-
S
i
n
-
S
i
-
1
n
2
Δ
t
,
Q
∂
S
∂
x
=
Q
¯
i
-
1
/
2
θ
(
S
i
n
+
1
-
S
i
-
1
n
+
1
)
+
(
1
-
θ
)
(
S
i
n
-
S
i
-
1
n
)
Δ
x
i
-
1
,
-
α
B
ω
(
S
-
S
*
)
=
-
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
(
S
i
n
+
1
-
S
*
i
n
+
1
)
-
1
2
α
i
-
1
n
+
1
B
i
-
1
n
+
1
ω
i
-
1
n
+
1
(
S
i
-
1
n
+
1
-
S
*
i
-
1
n
+
1
)
,
where
(18)
A
¯
i
-
1
/
2
=
1
2
(
A
i
n
+
1
+
A
i
-
1
n
+
1
)
,
Q
¯
i
-
1
/
2
=
1
2
(
Q
i
n
+
1
+
Q
i
-
1
n
+
1
)
.The irregular matrix equations are written as(19)
α
α
i
S
i
-
1
n
+
1
+
β
β
i
S
i
n
+
1
=
D
D
i
,
where
(20)
α
α
i
=
A
¯
i
-
1
/
2
2
Δ
t
-
Q
¯
i
-
1
/
2
θ
Δ
x
i
-
1
+
1
2
α
i
-
1
n
+
1
B
i
-
1
n
+
1
ω
i
-
1
n
+
1
,
β
β
i
=
A
¯
i
-
1
/
2
2
Δ
t
+
Q
¯
i
-
1
/
2
θ
Δ
x
i
-
1
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
,
D
D
i
=
A
¯
i
-
(
1
/
2
)
2
Δ
t
(
S
i
n
+
S
i
-
1
n
)
-
Q
¯
i
-
1
/
2
(
1
-
θ
)
Δ
x
i
-
1
(
S
i
n
-
S
i
-
1
n
)
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
S
*
i
n
+
1
+
1
2
α
i
-
1
n
+
1
B
i
-
1
n
+
1
ω
i
-
1
n
+
1
S
*
i
-
1
n
+
1
.(ii) For the upflow,(21)
A
∂
S
∂
t
=
A
¯
i
+
1
/
2
S
i
n
+
1
+
S
i
+
1
n
+
1
-
S
i
n
-
S
i
+
1
n
2
Δ
t
,
Q
∂
S
∂
x
=
Q
¯
i
+
1
/
2
θ
(
S
i
+
1
n
+
1
-
S
i
n
+
1
)
+
(
1
-
θ
)
(
S
i
+
1
n
-
S
i
n
)
Δ
x
i
,
-
α
B
ω
(
S
-
S
*
)
=
-
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
(
S
i
n
+
1
-
S
*
i
n
+
1
)
-
1
2
α
i
+
1
n
+
1
B
i
+
1
n
+
1
ω
i
+
1
n
+
1
(
S
i
+
1
n
+
1
-
S
*
i
+
1
n
+
1
)
,
where
(22)
A
¯
i
+
1
/
2
=
1
2
(
A
i
n
+
1
+
A
i
+
1
n
+
1
)
,
Q
¯
i
+
1
/
2
=
1
2
(
Q
i
n
+
1
+
Q
i
+
1
n
+
1
)
.The irregular matrix equations are written as(23)
β
β
i
S
i
n
+
1
+
r
r
i
S
i
+
1
n
+
1
=
D
D
i
,
where
(24)
β
β
i
=
A
¯
i
+
1
/
2
2
Δ
t
-
Q
¯
i
+
1
/
2
θ
Δ
x
i
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
,
r
r
i
=
A
¯
i
+
1
/
2
2
Δ
t
+
Q
¯
i
+
1
/
2
Δ
x
i
+
1
2
α
i
+
1
n
+
1
B
i
+
1
n
+
1
ω
i
+
1
n
+
1
,
D
D
i
=
A
¯
i
+
1
/
2
2
Δ
t
(
S
i
n
+
S
i
+
1
n
)
-
Q
¯
i
+
1
/
2
(
1
-
θ
)
Δ
x
i
(
S
i
+
1
n
-
S
i
n
)
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
S
*
i
n
+
1
+
1
2
α
i
+
1
n
+
1
B
i
+
1
n
+
1
ω
i
+
1
n
+
1
S
*
i
+
1
n
+
1
,
where α
α
i, β
β
i, γ
γ
i, and D
D
i are as the known coefficients of the equations.


### 2.3.2. Solution Algorithm for a Single River

In the tidal channel, the flow direction changes with tidal level fluctuation. The flow types of a single channel can be classified as four types (Figure3): (a) downflow, Q
s
≥
0, Q
e
≥
0; (b) upflow, Q
s
<
0, Q
e
<
0; (c) faced-flow, Q
s
≥
0, Q
e
<
0; and (d) depart-flow, Q
s
<
0, Q
e
≥
0.Figure 3
The four types of the flow at a single river channel. The types are as follows: downflow, upflow, faced-flow, and depart-flow, from the left to the right.For the downflow, it is single direction flow. If the sediment concentration at the upboundary junction (S
1
n
+
1) is known, the suspended sediment concentration at each junction can be obtained by the following formula:
(25)
S
i
n
+
1
=
P
i
+
R
i
S
1
n
+
1
(
i
=
1
,
2
,
…
,
n
+
1
)
,
where P
i
,
R
i are the coefficients, and they can be gotten by (19), (25), and written as
(26)
P
1
=
0
,
R
1
=
1
,
P
i
D
D
i
-
α
α
i
P
i
-
1
β
β
i
,
R
i
=
-
α
α
i
R
i
-
1
β
β
i
i
i
i
(
i
=
2,3
,
…
,
n
+
1
)
.For the upflow, it is also single direction flow. If the sediment concentration at the downboundary junction (S
n
+
1
n
+
1) is known, the sediment concentration at each junction can be obtained by the following formula:
(27)
S
i
n
+
1
=
P
i
+
R
i
S
n
+
1
n
+
1
(
i
=
1,2
,
…
,
n
+
1
)
,
where P
i
,
R
i are the coefficients, and they can be gotten by (23), (27), and written as
(28)
P
n
+
1
=
0
,
R
n
+
1
=
1
,
P
i
=
D
D
i
-
γ
γ
i
P
i
+
1
β
β
i
,
R
i
=
-
γ
γ
i
R
i
+
1
β
β
i
i
i
(
i
=
n
,
…
,
2,1
)
.For the faced-flow, it is double direction flow. The single river can be divided into two single direction flow river segments. If the sediment concentrations at the up and down boundary junctions (S
1
n
+
1, S
n
+
1
n
+
1) are known, the suspended sediment concentration at each junction can be obtained by the formulas (25) and (27).For the depart-flow, it is another double direction flow. First, the position of the stagnant point (Q
k
=
0) must be determined, the suspended sediment concentration at the stagnant point (S
1
n
+
1) may be obtained by the transport equation, and then the single river can be divided into two single direction flow river segments which boundary is the stagnant point; finally, the suspended sediment concentration at each junction can be obtained by the formulas (25), (27). The position of the stagnant point is presented as the following formula:
(29)
Δ
x
m
k
=
|
Q
m
|
|
Q
m
|
+
Q
l
Δ
x
m
l
.And the suspended sediment concentration at the stagnant point is(30)
S
k
n
+
1
=
A
k
n
A
k
n
+
1
S
k
n
exp
⁡
(
-
α
k
n
+
1
B
k
n
+
1
ω
k
n
+
1
A
k
n
+
1
Δ
t
)
.This paper also presents a junction-control method for suspended sediment in river networks.Assuming that the erosion and deposition at the junction are small, the suspended sediment mass conservation equation is written as(31)
∑
Q
in
S
in
=
∑
Q
out
S
out
(
S
out
=
S
J
)
,
where S
J is sediment concentration at junctions.Utilizing (25) and (27), the relation formula of sediment concentration between upjunction (
i
=
1
) and downjunction (
i
=
n
+
1
) is obtained, and the control equations for all junctions’ sediment concentration can be written in matrix notation as
(32)
F
j
(
S
1
,
S
2
,
…
,
S
M
)
=
0
j
=
1,2
,
…
,
M
.Solving (31), the suspended sediment concentration at all river junctions can be gotten. Finally, by solving the single river equations (25) and (27), suspended sediment concentration at all channel sections of the distributary channel networks can be calculated.


## 2.3.1. The Finite-Difference Equations

For suspended sediment transport, the upwind finite-difference scheme and Preissmann’s four-point implicit difference scheme are adopted (Figure2). According to the flow direction, the difference equations of the governing equation (2) can be expressed as follows.Figure 2
Illustration of the suspended sediment on grid and network domains.(i) For the downflow,(17)
A
∂
S
∂
t
=
A
¯
i
-
1
/
2
S
i
n
+
1
+
S
i
-
1
n
+
1
-
S
i
n
-
S
i
-
1
n
2
Δ
t
,
Q
∂
S
∂
x
=
Q
¯
i
-
1
/
2
θ
(
S
i
n
+
1
-
S
i
-
1
n
+
1
)
+
(
1
-
θ
)
(
S
i
n
-
S
i
-
1
n
)
Δ
x
i
-
1
,
-
α
B
ω
(
S
-
S
*
)
=
-
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
(
S
i
n
+
1
-
S
*
i
n
+
1
)
-
1
2
α
i
-
1
n
+
1
B
i
-
1
n
+
1
ω
i
-
1
n
+
1
(
S
i
-
1
n
+
1
-
S
*
i
-
1
n
+
1
)
,
where
(18)
A
¯
i
-
1
/
2
=
1
2
(
A
i
n
+
1
+
A
i
-
1
n
+
1
)
,
Q
¯
i
-
1
/
2
=
1
2
(
Q
i
n
+
1
+
Q
i
-
1
n
+
1
)
.The irregular matrix equations are written as(19)
α
α
i
S
i
-
1
n
+
1
+
β
β
i
S
i
n
+
1
=
D
D
i
,
where
(20)
α
α
i
=
A
¯
i
-
1
/
2
2
Δ
t
-
Q
¯
i
-
1
/
2
θ
Δ
x
i
-
1
+
1
2
α
i
-
1
n
+
1
B
i
-
1
n
+
1
ω
i
-
1
n
+
1
,
β
β
i
=
A
¯
i
-
1
/
2
2
Δ
t
+
Q
¯
i
-
1
/
2
θ
Δ
x
i
-
1
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
,
D
D
i
=
A
¯
i
-
(
1
/
2
)
2
Δ
t
(
S
i
n
+
S
i
-
1
n
)
-
Q
¯
i
-
1
/
2
(
1
-
θ
)
Δ
x
i
-
1
(
S
i
n
-
S
i
-
1
n
)
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
S
*
i
n
+
1
+
1
2
α
i
-
1
n
+
1
B
i
-
1
n
+
1
ω
i
-
1
n
+
1
S
*
i
-
1
n
+
1
.(ii) For the upflow,(21)
A
∂
S
∂
t
=
A
¯
i
+
1
/
2
S
i
n
+
1
+
S
i
+
1
n
+
1
-
S
i
n
-
S
i
+
1
n
2
Δ
t
,
Q
∂
S
∂
x
=
Q
¯
i
+
1
/
2
θ
(
S
i
+
1
n
+
1
-
S
i
n
+
1
)
+
(
1
-
θ
)
(
S
i
+
1
n
-
S
i
n
)
Δ
x
i
,
-
α
B
ω
(
S
-
S
*
)
=
-
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
(
S
i
n
+
1
-
S
*
i
n
+
1
)
-
1
2
α
i
+
1
n
+
1
B
i
+
1
n
+
1
ω
i
+
1
n
+
1
(
S
i
+
1
n
+
1
-
S
*
i
+
1
n
+
1
)
,
where
(22)
A
¯
i
+
1
/
2
=
1
2
(
A
i
n
+
1
+
A
i
+
1
n
+
1
)
,
Q
¯
i
+
1
/
2
=
1
2
(
Q
i
n
+
1
+
Q
i
+
1
n
+
1
)
.The irregular matrix equations are written as(23)
β
β
i
S
i
n
+
1
+
r
r
i
S
i
+
1
n
+
1
=
D
D
i
,
where
(24)
β
β
i
=
A
¯
i
+
1
/
2
2
Δ
t
-
Q
¯
i
+
1
/
2
θ
Δ
x
i
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
,
r
r
i
=
A
¯
i
+
1
/
2
2
Δ
t
+
Q
¯
i
+
1
/
2
Δ
x
i
+
1
2
α
i
+
1
n
+
1
B
i
+
1
n
+
1
ω
i
+
1
n
+
1
,
D
D
i
=
A
¯
i
+
1
/
2
2
Δ
t
(
S
i
n
+
S
i
+
1
n
)
-
Q
¯
i
+
1
/
2
(
1
-
θ
)
Δ
x
i
(
S
i
+
1
n
-
S
i
n
)
+
1
2
α
i
n
+
1
B
i
n
+
1
ω
i
n
+
1
S
*
i
n
+
1
+
1
2
α
i
+
1
n
+
1
B
i
+
1
n
+
1
ω
i
+
1
n
+
1
S
*
i
+
1
n
+
1
,
where α
α
i, β
β
i, γ
γ
i, and D
D
i are as the known coefficients of the equations.


## 2.3.2. Solution Algorithm for a Single River

In the tidal channel, the flow direction changes with tidal level fluctuation. The flow types of a single channel can be classified as four types (Figure3): (a) downflow, Q
s
≥
0, Q
e
≥
0; (b) upflow, Q
s
<
0, Q
e
<
0; (c) faced-flow, Q
s
≥
0, Q
e
<
0; and (d) depart-flow, Q
s
<
0, Q
e
≥
0.Figure 3
The four types of the flow at a single river channel. The types are as follows: downflow, upflow, faced-flow, and depart-flow, from the left to the right.For the downflow, it is single direction flow. If the sediment concentration at the upboundary junction (S
1
n
+
1) is known, the suspended sediment concentration at each junction can be obtained by the following formula:
(25)
S
i
n
+
1
=
P
i
+
R
i
S
1
n
+
1
(
i
=
1
,
2
,
…
,
n
+
1
)
,
where P
i
,
R
i are the coefficients, and they can be gotten by (19), (25), and written as
(26)
P
1
=
0
,
R
1
=
1
,
P
i
D
D
i
-
α
α
i
P
i
-
1
β
β
i
,
R
i
=
-
α
α
i
R
i
-
1
β
β
i
i
i
i
(
i
=
2,3
,
…
,
n
+
1
)
.For the upflow, it is also single direction flow. If the sediment concentration at the downboundary junction (S
n
+
1
n
+
1) is known, the sediment concentration at each junction can be obtained by the following formula:
(27)
S
i
n
+
1
=
P
i
+
R
i
S
n
+
1
n
+
1
(
i
=
1,2
,
…
,
n
+
1
)
,
where P
i
,
R
i are the coefficients, and they can be gotten by (23), (27), and written as
(28)
P
n
+
1
=
0
,
R
n
+
1
=
1
,
P
i
=
D
D
i
-
γ
γ
i
P
i
+
1
β
β
i
,
R
i
=
-
γ
γ
i
R
i
+
1
β
β
i
i
i
(
i
=
n
,
…
,
2,1
)
.For the faced-flow, it is double direction flow. The single river can be divided into two single direction flow river segments. If the sediment concentrations at the up and down boundary junctions (S
1
n
+
1, S
n
+
1
n
+
1) are known, the suspended sediment concentration at each junction can be obtained by the formulas (25) and (27).For the depart-flow, it is another double direction flow. First, the position of the stagnant point (Q
k
=
0) must be determined, the suspended sediment concentration at the stagnant point (S
1
n
+
1) may be obtained by the transport equation, and then the single river can be divided into two single direction flow river segments which boundary is the stagnant point; finally, the suspended sediment concentration at each junction can be obtained by the formulas (25), (27). The position of the stagnant point is presented as the following formula:
(29)
Δ
x
m
k
=
|
Q
m
|
|
Q
m
|
+
Q
l
Δ
x
m
l
.And the suspended sediment concentration at the stagnant point is(30)
S
k
n
+
1
=
A
k
n
A
k
n
+
1
S
k
n
exp
⁡
(
-
α
k
n
+
1
B
k
n
+
1
ω
k
n
+
1
A
k
n
+
1
Δ
t
)
.This paper also presents a junction-control method for suspended sediment in river networks.Assuming that the erosion and deposition at the junction are small, the suspended sediment mass conservation equation is written as(31)
∑
Q
in
S
in
=
∑
Q
out
S
out
(
S
out
=
S
J
)
,
where S
J is sediment concentration at junctions.Utilizing (25) and (27), the relation formula of sediment concentration between upjunction (
i
=
1
) and downjunction (
i
=
n
+
1
) is obtained, and the control equations for all junctions’ sediment concentration can be written in matrix notation as
(32)
F
j
(
S
1
,
S
2
,
…
,
S
M
)
=
0
j
=
1,2
,
…
,
M
.Solving (31), the suspended sediment concentration at all river junctions can be gotten. Finally, by solving the single river equations (25) and (27), suspended sediment concentration at all channel sections of the distributary channel networks can be calculated.


## 2.4. Modeling Setup

The river networks in the Pearl River Delta (PRD) (Figure4) are selected to check the applicability of the 1D flow and suspended sediment model. The river network in the PRD is considered as one of the world’s most complicated river networks, with the density of channel length per unit at 0.68–1.07 km/km2. The area of the Pearl River networks cover approximately 9750 km2, with a total length of over 1600 km and a coastline extending over 450 km from east to west. The Pearl River networks drain three main tributaries, Dongjiang River, Beijiang River, and Xijiang River, pass through this complex river network, and finally flow into the South China Sea through eight subestuaries that are separately called Humen, Jiaomen, Hongqimen, Hengmen, Modaomen, Jitimen, Hutiaomen, and Yamen.Figure 4
Map of the Pear River networks and sampling stations.The Pearl River network is simplified as 328 single river channels, 13 boundary rivers, 5108 computational river sections with the interval varying from 0.2 to 3.0 km, and 216 river channel junctions. Upstream boundary conditions are given by discharge and suspended sediment concentration at the five upstream boundary junctions (Shizui, Gaoyao, Shijiao, Laoyagang, and Boluo). Tidal level fluctuation and suspended sediment concentration at the downstream outlets are used as downstream boundary conditions. The Manning coefficients in the Pearl rivers are taken to be 0.01–0.05.K and m; the coefficients of sediment carrying capacity are taken to be 0.005 and 0.92. The value of adjustment coefficient for a cross-section α is 0.25 (when deposition) or 1.0 (when erosion). Computations are performed using a time step of 1 minute.


## 3. Results and Discussion

The distributary channel networks model has been extensively calibrated with respect to the field measurements to optimize the representation of the water level, discharge, and suspended sediment concentration over the study area. The field measurement data in 1998 and 2001 are used for calibration and validation, respectively. The Nash-Sutcliffe model efficiency (ME) which is the ratio of model error to variability in observational data is employed in this paper to evaluate the performance of the distributary channel networks model. Consider(33)
ME
=
1
-
∑
(
D
-
M
)
2
∑
(
D
-
D
¯
)
2
,whereD are the observational data, D
¯ is the mean of the observational data, and M is the corresponding model estimate. In the squaring of the error, a good fit receives a high score and a poor fit receives a low score. Performance levels are categorized as follows: ME
>
0.65 is excellent, 0.5–0.65 is very good, 0.2–0.5 is good, and ME
<
0.2 is poor [20]. The calibration and verification of water level, discharge, and suspended sediment concentration are shown in Figures 5 to 7.(a) The calibration of water levels from July 18 to July 19, 1998. (b) The validation of water levels from February 7 to February 16, 2001.
(a)
(b)The calibration result shows that the majority of the ME values for water levels over the Pearl River networks is over 0.65, indicating that the numerical model of the distributary channel networks can accurately simulate water level in the river networks (Figure5(a)). The water level is further verified by comparing the model results and measured value from February 7 to February 16, 2001 (Figure 5(b)). The comparison revealed that the simulated water level and distribution over the main branches matched excellently with the observed value and the ME value of almost all stations, which are greater than 0.65.In terms of water discharge (Figure6(a)), the modeled flow discharges follow the observations reasonably well, and the ME values are over 0.65 at 8 stations, showing an excellent skill. It can be noticed that there are two stations with efficiencies below 0.65, which are probably due to the let-out discharges of the reservoir nearby during the simulation period. The verification results also show a very good agreement in the discharge distribution over the Pearl River networks (Figure 6(b)). Based on the agreement in water level and discharge between the simulated and the measured throughout the distributary river networks, it is reasonable to conclude that the distributary channel networks model can simulate hydrodynamics process well in the river networks system.(a) The calibration of discharges from July 18 to July 19, 1998. (b) The validation of discharges from February 7 to February 16, 2001.
(a)
(b)(a) The calibration of suspended sediment concentration from July 18 to July 19, 1998. (b) The validation of suspended sediment concentration from February 7 to 16, 2001.
(a)
(b)Figures7(a) and 7(b) show the results of calibration and validation of the suspended sediment concentration in the Pearl River networks. It can be noticed that the simulated suspended sediment concentration in the Pearl River network still seem to be in a good agreement with the observed data, although the simulated results do not perform quite as well as water level and the discharge. The calibration results show that among the 10 stations observed, Sanshui station has ME value between 0.5 and 0.65, implying a very good skill, and nine stations have scores between 0.2 and 0.5, showing a good skill. The verification results show two stations with scores less than 0.2 accounting for 20% of the total number of stations with poor skills. Generally speaking, most of the stations have ME value for simulation sediment transport between 0.2 and 0.5 for both the wet and dry seasons.There are some potential reasons for errors between the simulated and observed suspended sediment concentration. Firstly, the bathymetry of the Pearl River network changes dramatically due to the intensive sand excavation during the last three decades. Therefore, the asynchronous of hydrographic and bathymetry measurement may directly lead to the inaccuracy of the simulation results. Actually, the field hydrographic measurements took place in 1998 and 2001, but the bathymetry survey happened in 1999 in our model. The asynchronous of hydrographic and bathymetry measurement should be one of the main reasons for the simulation errors. Secondly, large domain simulation needs more data and parameters, especially in suspended sediment transport simulation. However, some parameters of the suspended sediment, such as the diameter of the sediment, are regarded as relatively uniform distribution over the Pearl River network due to the limit of the data, which definitely lead to errors when simulating sediment transport. In general, the modeled results can be considered acceptable for simulating the flow and suspended sediment transport in the distributary channel networks.


## 4. Conclusion

A 1D hydrodynamics and suspended sediment model is developed for simulating the flow and sediment transport in the distributary channel networks. Firstly, by utilizing the finite-difference equations of the single river, as well as the mass and energy conservation at junctions, the equations of water level at junctions can be transformed to the irregular spare matrix equations. Water level and discharge at each channel section of a certain channel then can be obtained by the governing equations of the single channel. Assuming that the the erosion and deposition at the junctions are small, the nonequilibrium transport equation for suspended sediment can also be written in a matrix form. Solving the irregular spare matrix equations, suspended sediment concentration can also be obtained at all junctions. Finally, the suspended sediment concentration at each river section in river networks can be calculated.The distributary channel networks model was applied to the Pearl River networks, one of the most complicated river networks in the world, to check the performance of the model. The model was extensively calibrated and validated against field measurements to provide an accurate representation of water level and discharge, as well as suspended sediment transport in the Pearl River networks. The calibrated and validated results show that the distributary channel networks model display an excellent performance in simulating water level and discharge. Although the performance in simulating suspended sediment concentration is not as excellent as that of water level and discharge, the model is still good enough to simulate the suspended sediment transport in the distributary channel networks model.



---

*Source: 948731-2014-04-13.xml*