# A Coupled Model of the 1D River Network and 3D Estuary Based on Hydrodynamics and Suspended Sediment Simulation

Authors: Wei Zhang; Shujie Lyu; Yuliang Zhu; Xiaowen Chen
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

River networks and estuaries are very common in coastal areas. Runoff from the upper stream interacts with tidal current from open sea in these two systems, leading to a complex hydrodynamics process. Therefore, it is necessary to consider the two systems as a whole to study the flow and suspended sediment transport. Firstly, a 1D model is established in the Pearl River network and a 3D model is applied in its estuary. As sufficient mass exchanges between the river network and its estuary, a strict mathematical relationship of water level at the interfaces can be adopted to couple the 1D model with the 3D model. By doing so, the coupled model does not need to have common nested grids. The river network exchanges the suspended sediment with its estuary by adding the continuity conditions at the interfaces. The coupled model is, respectively, calibrated in the dry season and the wet season. The results demonstrate that the coupled model works excellently in simulating water level and discharge. Although there are more errors in simulating suspended sediment concentration due to some reasons, the coupled model is still good enough to evaluate the suspended sediment transport in river network and estuary systems.


---

## Body

## 1. Introduction

As a link between marine environments and rivers, estuaries are characterized by a variety of complex and complicated processes [1]. Therefore, there has been a growing interest in the field of estuarine and coastal sediment dynamics, including sediment transport, geomorphological process, and the transport of heavy metal by sediment particles [2, 3]. In the researches on hydrodynamics and suspended sediment transport in this region, numerical model is no doubt a highly efficient and low-cost method to describe estuarine flow and sediment transport processes in recent decades. Normally, the 1D model is more suitable than the 3D coastal model to simulate the flow and suspended sediment transport in river network due to the complex river network geometries and complex structures like weirs, barrages, and dams. The 3D model is more useful in evaluating the coastal flow and sediments movement due to the vertical circulations. However, if both the upstream river network geometries and the downstream tidal currents in the estuary are complicated, application of individual 1D or 3D model for the river and the estuarine area will result in either an inadequate representation of these processes (in case of 1D) or a computationally inefficient method (in case of 3D) [4]. Therefore, it is necessary to take upstream river network and downstream estuarine system as a whole in one modelling system, to resolve the relevant physical processes in a consistent way.A lot of efforts have been made recently to develop coupled modelling systems, which were designed for different estuaries [5–7]. Some of the coupled models were used to predict flood effects [8], some to investigate mass fluxes and transformations of nutrients [9], and some to prevent further water pollution [10] and so on. Some researchers developed the “domain decomposition” technique to increase computational efficiency [6]. This technique creates efficient model grids by “stitching” together different partial grids, which enable the high resolution and low resolution areas to be combined and modelled. However, this technique may not be stable due to the existing nested grids. Zhou et al. [10] provide a method that the 1D and 3D schematisations do not need to overlap. The upstream water level of the 3D model is used for the 1D outflow open boundary and the downstream outflow of the 1D model is provided to the upstream inflow of the 1D model [10]. The 1D model and the 3D model run simultaneously and the interface information will be saved for offline computation. Offline coupling requires huge storage space. In order to resolve this, an online coupling method is developed. This method resolves the real-time exchange and interaction between the coastal waters and river network, allowing for accurate and mass conserving modelling of complex coastal water and river network system. It is easy to implement as it only replaced file I/O (read) by communication operations (receive) [11].The Pearl River Estuary (PRE) and its upstream river network form the largest river system in the southern China. Historically a rich and fertile area, it has now become one of the most populated and most developed regions in China. However, as the economy grows in the region, serious problems also occur at the Pearl River network and its estuary, such as salinity intrusion [12], riverbed erosion, and water pollution. Therefore, coupled 1D–3D hydrodynamic modelling, such as horizontal fronts, vertical stratification, and salt intrusion, was developed to simulate the hydrodynamics and salinity in this area [4]. Suspended sediment transport, however, was not involved in this coupled model. Recently, a coupled model has been used to investigate the transport of pollutants and water quality parameters, to provide a scientific basis for the improvement of the water quality of the Pearl River and to prevent further water pollution [10]. In addition, the models were also used for the construction of the nutrient budgets for the entire area and the characterization of the key biogeochemical processes in the transfers and transformations of nutrients [9]. However, previous works on the coupled models fail to provide mathematics derivation on the relations of water level between 1D and 3D models. Moreover, most coupled models focus on the hydrodynamics or water quality in the river network and its estuary, but little attention has been paid to the suspended sediment transport. In fact, the transport of most organic carbon from land to oceans is closely related to river sediment transport. Hence, understanding the estuarine dynamics of suspended sediment is also crucial to resolve estuarine problems [11, 13]. Therefore, this paper provides a method to establish a coupled model of hydrodynamics and suspended sediment transport by strict mathematical derivation, rather than physical interpretation in the river network and its estuary.


## 2. Methods

The model coupled the 1D river network model with the 3D coastal model. It takes into consideration the complex topography of the river network and the vertical gradients of the coastal hydrodynamics. The model is applied in the PRD, which is divided into two parts. The 1D model is set up in the Pearl River network and the 3D model is applied in the PRE to investigate the transport of water flow and suspended sediment in the delta.


### 2.1. One-Dimensional Hydrodynamics and Suspended Sediment Model

Based on the Saint-Venant equations and the nonequilibrium transport equation for the suspended sediment, a numerical model of junction-channel for the river network is established by using the “junction-control” method for water level and suspended sediment concentration. Firstly, the river network is schematized as river channels, junctions, and cross-sections. Secondly, the variables at each junction are calculated by solving the governing equations. Finally, the water levels, discharges, and suspended sediment concentrations at each cross-section in each river channel are calculated. More details about the hydrodynamics difference scheme and transforms of the finite-difference equations can be found in Zhang et al. [12]. The algorithm for suspended sediment concentration in the river network is similar to that for unsteady flow. The details are introduced in the following parts.The governing equation for the nonequilibrium transport of suspended load [14] is expressed as
(1)
∂
(
A
S
)
∂
t
+
∂
(
Q
S
)
∂
x
=
-
α
B
ω
(
S
-
S
*
)
.The riverbed transfiguration equation is(2)
γ
′
B
∂
Z
0
∂
t
=
α
B
ω
(
S
-
S
*
)
-
∂
G
b
∂
x
,
where Q is the discharge, S is the cross-section-averaged suspended sediment concentration, S
* is the cross-section-averaged suspended sediment transport capacity, A is the cross-sectional area, B is the water surface width, H is the mean elevation of the cross-sectional water surface, Z
0 is the mean elevation of the riverbed, U is the cross-section-averaged velocity, ω is the settling velocity of suspended sediment, α is a recovery coefficient relating to the suspended sediment of the saturation, G
b is the bed load sediment transport rate, D is the cross-sectional mean grain diameter of bed load, and γ
′ is the dry density of suspended sediments.For suspended sediment transport equation, the upwind finite difference scheme and Preissmann four-point implicit difference scheme are employed here. According to the direction of flow, the difference equations are written as(3)
α
α
i
S
i
-
1
n
+
1
+
β
β
i
S
i
n
+
1
=
D
D
i
(
if
Q
i
>
0
)
,
(4)
β
β
i
S
i
n
+
1
+
γ
γ
i
S
i
+
1
n
+
1
=
D
D
i
(
if
Q
i
<
0
)
,
where α
α
i, β
β
i, γ
γ
i, and D
D
i are known coefficients of the difference equations.In a tidal river, flow direction changes with ebb and flood current. The flow at a single river channel can be classified by four types (Figure1): (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; (d) depart-flow, Q
s
<
0, Q
e
≥
0; Q
s and Q
e are the discharge at the start and end cross-section of a single channel.Figure 1
The flow direction at a single river channel changes with ebb and flood current in the tidal river. And the direction can be classified by four types (from the left to the right):(
1
) downflow, (
2
) upflow, (
3
) faced-flow, and (
4
) depart-flow.For the downflow, it is single direction flow. If the sediment concentration at the up-boundary node (S
1
n
+
1) is known, the sediment concentration at each node can be obtained by the following formula:
(5)
S
i
n
+
1
=
P
i
+
R
i
S
1
n
+
1
(
i
=
1,2
,
…
,
n
+
1
)
,
where P
i and R
i are the coefficients; they can be gotten by (3) and written as
(6)
P
1
=
0
,
R
1
=
1
,
P
i
=
D
D
i
-
α
α
i
P
i
-
1
β
β
i
,
R
i
=
α
α
i
R
i
-
1
β
β
i
(
i
=
2,3
,
…
,
n
+
1
)
.For the upflow, it is also single direction flow. If the sediment concentration at the down-boundary node (S
n
+
1
n
+
1) is known, the sediment concentration at each node can be obtained by the following formula:
(7)
S
i
n
+
1
=
P
i
+
R
i
S
n
+
1
n
+
1
(
i
=
1,2
,
…
,
n
+
1
)
,
where P
i and R
i are the coefficients; they can be gotten by (4) and written as
(8)
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
=
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 nodes (S
1
n
+
1, S
n
+
1
n
+
1) are known, the sediment concentration at each node can be obtained by formulas (5) and (7).For the depart-flow, it is another double direction flow. Firstly, the position of the stagnant point (Q
k
=
0) must be determined, the 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 in which interboundary is the stagnant point; finally the sediment concentration at each river junction can be obtained by the formulas (5) and (7). The position of the stagnant point is calculated by the following formula:
(9)
Δ
x
m
k
=
|
Q
m
|
|
Q
m
|
+
Q
l
Δ
x
m
l
.Finally, the sediment concentration at the stagnant point is(10)
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
)
.


### 2.2. Three-Dimensional Hydrodynamic and Suspended Sediment Model

The 3D hydrodynamic and suspended sediment model is established based on ECOM model, a free-surface, fully nonlinear, primitive-equation estuarine and coastal ocean model. The equations under the orthogonal curve coordinate include water continuity equation and the momentum equations:(11)
∇
V
→
+
∂
W
∂
z
=
0
,
∂
U
∂
t
+
V
→
·
∇
U
+
W
∂
U
∂
z
-
f
V
=
-
1
ρ
0
∂
P
∂
x
+
∂
∂
z
(
K
M
∂
U
∂
z
)
+
∂
∂
x
[
2
A
M
∂
U
∂
x
]
+
∂
∂
y
[
A
M
(
∂
U
∂
y
+
∂
V
∂
x
)
]
,
∂
V
∂
t
+
V
→
·
∇
V
+
W
∂
V
∂
z
+
f
U
=
-
1
ρ
0
∂
P
∂
y
+
∂
∂
z
(
K
M
∂
V
∂
z
)
+
∂
∂
y
[
2
A
M
∂
V
∂
y
]
+
∂
∂
x
[
A
M
(
∂
U
∂
y
+
∂
V
∂
x
)
]
,
ρ
g
=
-
∂
P
∂
z
,
where f is the Coriolis coefficient, ρ
0 is the referenced density, ρ is the local density, g is the gravitational acceleration, K
M and A
M are the turbulent viscosity coefficients of turbulent momentum in vertical and horizontal, respectively, and U, V, and W are the velocity in x, y, and z direction in the Cartesian coordinate system.The diffusion process of salinity is also necessary to consider for the model. The conservation equations are described as(12)
∂
θ
∂
t
+
V
→
·
Δ
θ
+
W
∂
θ
∂
z
=
∂
∂
z
[
K
H
∂
θ
∂
z
]
+
∂
∂
x
[
A
H
∂
θ
∂
x
]
+
∂
∂
y
[
A
H
∂
θ
∂
y
]
,
∂
S
∂
t
+
V
→
·
Δ
S
+
W
∂
S
∂
z
=
∂
∂
z
[
K
H
∂
S
∂
z
]
+
∂
∂
x
[
A
H
∂
S
∂
x
]
+
∂
∂
y
[
A
H
∂
θ
∂
y
]
,
where K
H and A
H are the diffusion coefficients of the vertical and horizontal, respectively, and θ is potential temperature.The transport equation of the suspended sediment is written as(13)
∂
C
k
∂
t
+
∂
U
C
k
∂
x
+
∂
V
C
k
∂
y
+
∂
(
W
-
ω
s
,
k
)
C
k
∂
z
=
∂
∂
x
(
A
H
∂
C
k
∂
x
)
+
∂
∂
y
(
A
H
∂
C
k
∂
y
)
+
∂
∂
z
(
K
H
∂
C
k
∂
z
)
,
where C
k is the sediment concentration and ω
s
,
k is the settling velocity of sediments.The details of ECOM model can refer to Zhu et al. [15].


### 2.3. Coupled 1D and 3D Hydrodynamic and Suspended Sediment Model

In this paper, the relationship expression of water level between 1D model and 3D model is deduced by strict mathematics method. The 1D and 3D models meet at outlets of river network and estuary, and the hydroinformation is transferred through interfaces, by which the additional hydrodynamic equation of the 1D and 3D coupled numerical model can be derived. The 1D and 3D boundary conditions at interfaces are satisfied through iteration.The water level and discharge at the interfaces of the 1D model and 3D model can be reasonable considered equal:(14)
Z
i
(
1
)
=
Z
i
(
3
)
,
i
=
1,2
,
…
N
,
Q
i
(
1
)
=
Q
i
(
3
)
=
∑
j
=
1
0
∑
l
=
1
m
q
j
,
l
,
i
=
1,2
,
…
N
,
where Z
i
(
1
) and Q
i
(
1
) are water level and discharge at the ith interface of the 1D model, respectively, and Z
i
(
3
) and Q
i
(
3
) stand for water level and discharge at the ith interface of the 3D model, q
j
,
l is the discharge of the grids, m is the number of the grids of the interface, and o is the number of the layers of water depth.Due to the difference of water levels at grids of 3D model, the down-boundary of the 1D model is the averaged water levels of the grids in the 3D model:(15)
Z
i
(
1
)
=
Z
i
(
3
)
=
1
m
z
i
3
(
x
,
y
,
z
)
,
i
=
1,2
,
…
N
.Based on the method of the 1D hydrodynamics model, the governing equation can be deduced to matrix of water level at junctions as below:(16)
A
[
Z
]
=
B
,
where A is the coefficient matrix, Z is the water levels at junctions, and B is constants.However, in a coupled model, the interfaced water levels are not known, which lead toB containing unknown term and constant term. Therefore, the unknown part should be separated from the matrix B:
(17)
A
[
Z
]
=
C
[
H
]
+
D
,
where H is the water levels at the interfaces, C is the corresponding coefficient matrix, and D is the constant term.Gauss elimination method is used to simplify the coefficient matrixA to a diagonal matrix A:
(18)
A
′
[
Z
]
=
C
′
[
H
]
+
D
′
.Then, the water levels at junctions connected to the interfaces can be expressed by the water levels at interfaces of 1D model and 3D model:(19)
Z
i
(
1
)
=
∑
j
=
1
n
c
i
,
j
′
H
j
+
d
i
′′
,
i
=
1
,
…
n
,
where c
i
,
j
′′ and d
i
′′ are coefficient and constant after unitization, respectively.Utilizing the relations between the water level and discharge at junctions [12], the discharge at the beginning and the end of a single river channel can be deduced as below:
(20)
Q
i
,
1
=
α
i
,
1
+
β
i
,
1
H
i
,
1
+
γ
i
,
1
Z
i
,
n
+
1
,
k
k
k
k
k
k
k
k
k
k
i
=
1,2
,
…
,
N
(21)
Q
i
,
n
+
1
=
ξ
i
,
n
+
1
+
ζ
i
,
n
+
1
H
i
,
n
+
1
+
η
i
,
n
+
1
Z
i
,
1
,
k
k
k
k
k
k
k
k
k
k
k
k
k
k
i
=
1,2
,
…
,
N
.Equation (19) is substituted to (20) and (21), respectively; then the discharge at the interfaces of 1D and 3D models can be totally represented by the water level at the interfaces, and ξ
,
ζ
,
η are the coefficients of the equations. Consider the following:
(22)
Q
i
,
1
=
α
i
,
1
+
β
i
,
1
H
i
,
1
+
γ
i
,
1
(
∑
j
=
1
n
c
i
,
j
′′
H
i
,
1
+
d
i
′′
)
,
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
i
=
1,2
,
…
,
N
Q
i
,
n
+
1
=
ξ
i
,
n
+
1
+
ζ
i
,
n
+
1
H
i
,
n
+
1
+
η
i
,
n
+
1
(
∑
j
=
1
n
c
i
,
j
′′
H
i
,
n
+
1
+
d
i
′′
)
,
i
=
1,2
,
…
,
N
.After simplifying (22), the relations between water level and discharge can be established in 1D model. Consider the following:
(23)
Q
i
(
1
)
=
∑
j
=
1
n
p
i
,
j
H
i
+
q
i
,
i
=
1
,
…
,
n
,
and p is the coefficient of equations and q is constant.Finally, taking into account (14), an iterative equation for water level at interface is established as below:
(24)
H
i
(
l
+
1
)
=
1
p
i
,
j
[
Q
i
(
3
)
-
∑
j
=
1
j
≠
i
n
p
i
,
j
H
i
(
l
)
-
q
i
]
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
i
=
1
,
…
,
n
.At the beginning of the calculation, the initial dischargesQ
1
(
0
)
,
Q
2
(
0
)
,
…
,
Q
N
-
1
(
0
)
,
Q
N
(
0
) are used as the upper boundary condition of the 3D hydrodynamics model and can get the water levels H
1
(
0
)
,
H
2
(
0
)
,
…
,
H
N
-
1
(
0
)
,
H
N
(
0
). Then, substituting these water levels into (24) is to get the approximate values H
1
(
1
)
,
H
2
(
1
)
,
…
,
H
m
-
1
(
1
)
,
H
m
(
1
), which can be regarded as the downstream boundaries of 1D hydrodynamics if water levels satisfy the convergence condition MAX
|
H
i
(
1
)
-
H
i
(
0
)
|
≤
ε (ε is a given error). Otherwise, H
1
(
0
)
,
H
2
(
0
)
,
…
,
H
N
-
1
(
0
)
,
H
N
(
0
) have to be substituted into (23) to get discharges Q
1
(
1
)
,
Q
2
(
1
)
,
…
,
Q
N
-
1
(
1
)
,
Q
N
(
1
) for 3D model calculation again. This iterative process stops until the convergence condition is satisfied.As for the suspended sediment model, two additional conditions are needed to be considered. Firstly, the sediment fluxes at the interfaces of the 1D model and 3D model should be equal:(25)
Q
i
(
1
)
S
i
(
1
)
=
Q
i
(
3
)
S
i
(
3
)
,
i
=
1
,
…
,
N
.Secondly, the spatial averaged sediment concentration of 3D model should be closed to that of 1D model:(26)
S
i
(
3
)
(
x
,
y
,
z
)
=
S
i
(
1
)
,
S
i
(
1
)
=
S
i
(
3
)
=
1
m
s
i
(
3
)
(
x
,
y
,
z
)
,
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
i
=
1
,
…
,
N
.


### 2.4. Study Area and Data Processing

The Pearl River Delta (PRD) is located in the northern of the South China Sea. It is a dynamically complex estuarine system, which encompasses Pearl River network and PRE. The Pearl River, a subtropical river, is the second largest Chinese river in terms of annual water discharge (over 3 × 1011 m3). The channel network is formed by a large number of intricately interlaced deltaic branches, including three principal rivers: Xijiang, Beijiang, and Dongjiang, and some small rivers draining the PRD. It is a catchment area of 26,820 km2 and a total length of 1,600 km [16]. Meanwhile, as the distance from one junction to another is 0.68–1.07 km/km2, the Pearl River network is considered as one of the world’s most complicated river networks, which is formed by a large number of intricately interlaced branches [12]. In the river network, the mean suspended sediment concentration is about 0.284 kg/m3 [17], and annual mean sediment load in 1957 to 2006 is 75.3 × 109 kg/yr, with an annual average sediment load of 66.83 Mt/yr in Xijiang alone [18].These fluxes pass through the river network and flow into the PRE through the eight outlets Hutiaomen and Yamen discharging into the Huangmao Bay, Modaomen and Jitimen discharging into open sea directly, and Humen, Jiaomen, Hongqimen, and Hengmen discharging into the Lingding Bay [4]. The PRE has a broad bell shape, in which the Lingding Bay is like an inverted funnel with the narrow neck in the north and wide mouth opening to the south [19]. As the gravitational circulation develops in the PRE, the residual currents are landwards in deep channels and seawards in the western part of the PRE [20]. Normally, the circulation and water properties in the PRE are caused by various forcing mechanisms, including river discharges, tides, monsoon winds, coastal currents, and buoyancy forcing associated with the mixing of freshwater and saline water. The tides in the PRE mainly come from the Pacific oceanic tidal propagation [21], with a mean tidal range between 1.0 and 1.7 m. Thus PRE can be classified as a microtidal estuary. Propagating from offshore towards the estuary, the tidal range increases gradually and reaches a maximum near the river outlets [19]. Salinity in the PRE seldom exceeds 32 psu even during the dry season and is highly stratified in the vertical with a very low salinity surface layer. The distributions of salinity have a significant seasonal variability [7]. Further, the carrying capacity of sediments in the PRE is basically controlled by tidal flow, wave, fluid, and particle features. In the PRE, the very fine-grained silt and clay are contained in the suspended sediments, within which most particles are less than 100 μm in diameter and the median particle diameter d
50 is about 8 μm [22].


### 2.5. The Setup of the Coupled 1D and 3D Model

The model is applied to the Pearl River network and its estuary. The river network in 1D model is discretized into 347 channels, 13 boundary rivers, and 1,850 cross-sections with the interval varying from 0.2 to 3.0 km and 220 junctions. The length of the river network is about 1,600 km. The upstream boundaries include Shizui at the Tanjiang River, Gaoyao at Xijiang, Shijiao at Beijiang, Laoyagang at the Liuxi River, and Boluo at Dongjiang. The downstream boundaries include Dahu at Humen, Nansha at Jiaomen, Wanqingsha at Hongqili, Hengmen at Hengmen, Denglongshan at Maodaomen, Huangjin at Jitimen, Hutiaomen at Hutiaomen, and Huangchong at Yamen. These eight outlets are also the interfaces of the 1D and 3D model. The 3D estuarine model covers the area from the eight outlets to the −30 m isobaths (Figure2). More than 57,000 orthogonal curvilinear grids are used in the 3D model with spatial resolution from 0.05 km to 1 km and 10 equidistant sigma levels in vertical.Figure 2
Map of the study area including(
1
) the location of the Pearl River network with the names of the three major rivers and (
2
) the major topographic features of the Pearl River Estuary with isobaths in colors. The upstream boundaries of the 1D model are shown on the map, which starts with Shizui, Gaoyao, Shijiao, Laoyagang, and Boluo. The downstream boundaries of the 3D coastal model extend to the −30 m isobaths. Black circles represent major cities. Red triangles and red circles represent the main stations in the Pearl River network and the Pearl River Estuary, respectively.The hourly observed runoff at the upstream stations is the upper boundary conditions of the coupled model. The salinity at those stations is definitely zero. At the open sea part of the coupled model, the tidal levels are evaluated by the harmonic analysis from the most important tidal constituents of the PRE, with the uniform values of temperature and salinity from the observations. Consider the following:(27)
ζ
=
a
0
+
∑
i
=
1
8
a
i
cos
⁡
(
ω
i
t
-
ϕ
i
)
,
where a
0 is the mean water level, a
i is the amplitude of partial tide, and ϕ
i is the phase of partial tide.The simulation is initialized with zero water levels and velocities in the bathymetry of 1998-1999 and runs about 30 days. The results of the model are used as the new initial conditions and the simulation hot-started with the results. The processes are repeated several times until the error between the last two simulations satisfies a given condition. These results are then used as the initial conditions of the coupled model. Two different periods, including July 1998 (wet season) and December 1998 (dry season), are used for model calibration. In both the 1D and 3D model domains, the nonuniform parameters, such as bed roughness parameter and particle diameter, have been adjusted in order to optimize the water level, discharge, and suspended sediment concentration over the upstream river network and its coastal waters. The specific key parameters of the coupled model are shown in Table1.Table 1
Parameters of the modeling.
Settings
Description
1D
Time steps
Based on the scope of model, the time step of the 1D model is 10 s
Roughness coefficient
This coefficient depends on the physical properties of the riverbed. At Beijiang, from upstream boundary to Sanshuin
=  0.026~0.045, from Sanshui to downstream boundary n
=  0.015~0.03; at Xijiang n
=  0.02~0.03; at Dongjiang, from upstream boundary to Shilong n
=  0.025~0.045, from Shilong to downstream boundary n
=  0.015~0.036; at the mouth of the river network n
=  0.01~0.02
Particle diameter
In channel, the median particle diameter of suspended sediment is about 0.01~0.040 mm, and the median particle diameter of bed material is about 0.17~0.44 mm.
Sediment transport capacity
S
*
=
K
(
U
3
/
g
R
ω
)
m,  where m
=
0.92; K
=  0.0012~0.0016.
Recovery coefficient relating to the saturation
When the channels are in the deposition conditions, a value of 0.25 is applied; while, in the erosion conditions, a value of 1.0 is applied.
3D
Time steps
Δ
t
E
≤
1
/
C
t
|
(
1
/
δ
x
2
)
+
(
1
/
δ
y
2
)
|
(
-
1
/
2
), Δ
t
I
≤
1
/
C
T
|
(
1
/
δ
x
2
)
+
(
1
/
δ
y
2
)
|
-
1
/
2;  thus, the internal model is 10 s and the external model is 0.5 s.
Bottom roughness
τ
b
→
=
ρ
0
C
D
|
U
b
→
|
U
b
→, C
D
=
max
⁡
[
κ
2
/
ln
2
⁡
(
z
/
z
0
)
,
0.0025
], where κ
=
0.4, z
0
=
0.15 cm.
Particle diameter
The median particle diameter is 0.003~0.009 mm during the spring tide period.
Sediment settling velocity
ω
s
,
1
=
α
(
C
1
G
)
β, where α
=
2.42, β
=
0.22, and ω
s
,
2
=
υ
/
D
k
[
(
25
+
1.2
D
*
2
)
0.5
-
5
]
1.5


## 2.1. One-Dimensional Hydrodynamics and Suspended Sediment Model

Based on the Saint-Venant equations and the nonequilibrium transport equation for the suspended sediment, a numerical model of junction-channel for the river network is established by using the “junction-control” method for water level and suspended sediment concentration. Firstly, the river network is schematized as river channels, junctions, and cross-sections. Secondly, the variables at each junction are calculated by solving the governing equations. Finally, the water levels, discharges, and suspended sediment concentrations at each cross-section in each river channel are calculated. More details about the hydrodynamics difference scheme and transforms of the finite-difference equations can be found in Zhang et al. [12]. The algorithm for suspended sediment concentration in the river network is similar to that for unsteady flow. The details are introduced in the following parts.The governing equation for the nonequilibrium transport of suspended load [14] is expressed as
(1)
∂
(
A
S
)
∂
t
+
∂
(
Q
S
)
∂
x
=
-
α
B
ω
(
S
-
S
*
)
.The riverbed transfiguration equation is(2)
γ
′
B
∂
Z
0
∂
t
=
α
B
ω
(
S
-
S
*
)
-
∂
G
b
∂
x
,
where Q is the discharge, S is the cross-section-averaged suspended sediment concentration, S
* is the cross-section-averaged suspended sediment transport capacity, A is the cross-sectional area, B is the water surface width, H is the mean elevation of the cross-sectional water surface, Z
0 is the mean elevation of the riverbed, U is the cross-section-averaged velocity, ω is the settling velocity of suspended sediment, α is a recovery coefficient relating to the suspended sediment of the saturation, G
b is the bed load sediment transport rate, D is the cross-sectional mean grain diameter of bed load, and γ
′ is the dry density of suspended sediments.For suspended sediment transport equation, the upwind finite difference scheme and Preissmann four-point implicit difference scheme are employed here. According to the direction of flow, the difference equations are written as(3)
α
α
i
S
i
-
1
n
+
1
+
β
β
i
S
i
n
+
1
=
D
D
i
(
if
Q
i
>
0
)
,
(4)
β
β
i
S
i
n
+
1
+
γ
γ
i
S
i
+
1
n
+
1
=
D
D
i
(
if
Q
i
<
0
)
,
where α
α
i, β
β
i, γ
γ
i, and D
D
i are known coefficients of the difference equations.In a tidal river, flow direction changes with ebb and flood current. The flow at a single river channel can be classified by four types (Figure1): (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; (d) depart-flow, Q
s
<
0, Q
e
≥
0; Q
s and Q
e are the discharge at the start and end cross-section of a single channel.Figure 1
The flow direction at a single river channel changes with ebb and flood current in the tidal river. And the direction can be classified by four types (from the left to the right):(
1
) downflow, (
2
) upflow, (
3
) faced-flow, and (
4
) depart-flow.For the downflow, it is single direction flow. If the sediment concentration at the up-boundary node (S
1
n
+
1) is known, the sediment concentration at each node can be obtained by the following formula:
(5)
S
i
n
+
1
=
P
i
+
R
i
S
1
n
+
1
(
i
=
1,2
,
…
,
n
+
1
)
,
where P
i and R
i are the coefficients; they can be gotten by (3) and written as
(6)
P
1
=
0
,
R
1
=
1
,
P
i
=
D
D
i
-
α
α
i
P
i
-
1
β
β
i
,
R
i
=
α
α
i
R
i
-
1
β
β
i
(
i
=
2,3
,
…
,
n
+
1
)
.For the upflow, it is also single direction flow. If the sediment concentration at the down-boundary node (S
n
+
1
n
+
1) is known, the sediment concentration at each node can be obtained by the following formula:
(7)
S
i
n
+
1
=
P
i
+
R
i
S
n
+
1
n
+
1
(
i
=
1,2
,
…
,
n
+
1
)
,
where P
i and R
i are the coefficients; they can be gotten by (4) and written as
(8)
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
=
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 nodes (S
1
n
+
1, S
n
+
1
n
+
1) are known, the sediment concentration at each node can be obtained by formulas (5) and (7).For the depart-flow, it is another double direction flow. Firstly, the position of the stagnant point (Q
k
=
0) must be determined, the 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 in which interboundary is the stagnant point; finally the sediment concentration at each river junction can be obtained by the formulas (5) and (7). The position of the stagnant point is calculated by the following formula:
(9)
Δ
x
m
k
=
|
Q
m
|
|
Q
m
|
+
Q
l
Δ
x
m
l
.Finally, the sediment concentration at the stagnant point is(10)
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
)
.


## 2.2. Three-Dimensional Hydrodynamic and Suspended Sediment Model

The 3D hydrodynamic and suspended sediment model is established based on ECOM model, a free-surface, fully nonlinear, primitive-equation estuarine and coastal ocean model. The equations under the orthogonal curve coordinate include water continuity equation and the momentum equations:(11)
∇
V
→
+
∂
W
∂
z
=
0
,
∂
U
∂
t
+
V
→
·
∇
U
+
W
∂
U
∂
z
-
f
V
=
-
1
ρ
0
∂
P
∂
x
+
∂
∂
z
(
K
M
∂
U
∂
z
)
+
∂
∂
x
[
2
A
M
∂
U
∂
x
]
+
∂
∂
y
[
A
M
(
∂
U
∂
y
+
∂
V
∂
x
)
]
,
∂
V
∂
t
+
V
→
·
∇
V
+
W
∂
V
∂
z
+
f
U
=
-
1
ρ
0
∂
P
∂
y
+
∂
∂
z
(
K
M
∂
V
∂
z
)
+
∂
∂
y
[
2
A
M
∂
V
∂
y
]
+
∂
∂
x
[
A
M
(
∂
U
∂
y
+
∂
V
∂
x
)
]
,
ρ
g
=
-
∂
P
∂
z
,
where f is the Coriolis coefficient, ρ
0 is the referenced density, ρ is the local density, g is the gravitational acceleration, K
M and A
M are the turbulent viscosity coefficients of turbulent momentum in vertical and horizontal, respectively, and U, V, and W are the velocity in x, y, and z direction in the Cartesian coordinate system.The diffusion process of salinity is also necessary to consider for the model. The conservation equations are described as(12)
∂
θ
∂
t
+
V
→
·
Δ
θ
+
W
∂
θ
∂
z
=
∂
∂
z
[
K
H
∂
θ
∂
z
]
+
∂
∂
x
[
A
H
∂
θ
∂
x
]
+
∂
∂
y
[
A
H
∂
θ
∂
y
]
,
∂
S
∂
t
+
V
→
·
Δ
S
+
W
∂
S
∂
z
=
∂
∂
z
[
K
H
∂
S
∂
z
]
+
∂
∂
x
[
A
H
∂
S
∂
x
]
+
∂
∂
y
[
A
H
∂
θ
∂
y
]
,
where K
H and A
H are the diffusion coefficients of the vertical and horizontal, respectively, and θ is potential temperature.The transport equation of the suspended sediment is written as(13)
∂
C
k
∂
t
+
∂
U
C
k
∂
x
+
∂
V
C
k
∂
y
+
∂
(
W
-
ω
s
,
k
)
C
k
∂
z
=
∂
∂
x
(
A
H
∂
C
k
∂
x
)
+
∂
∂
y
(
A
H
∂
C
k
∂
y
)
+
∂
∂
z
(
K
H
∂
C
k
∂
z
)
,
where C
k is the sediment concentration and ω
s
,
k is the settling velocity of sediments.The details of ECOM model can refer to Zhu et al. [15].


## 2.3. Coupled 1D and 3D Hydrodynamic and Suspended Sediment Model

In this paper, the relationship expression of water level between 1D model and 3D model is deduced by strict mathematics method. The 1D and 3D models meet at outlets of river network and estuary, and the hydroinformation is transferred through interfaces, by which the additional hydrodynamic equation of the 1D and 3D coupled numerical model can be derived. The 1D and 3D boundary conditions at interfaces are satisfied through iteration.The water level and discharge at the interfaces of the 1D model and 3D model can be reasonable considered equal:(14)
Z
i
(
1
)
=
Z
i
(
3
)
,
i
=
1,2
,
…
N
,
Q
i
(
1
)
=
Q
i
(
3
)
=
∑
j
=
1
0
∑
l
=
1
m
q
j
,
l
,
i
=
1,2
,
…
N
,
where Z
i
(
1
) and Q
i
(
1
) are water level and discharge at the ith interface of the 1D model, respectively, and Z
i
(
3
) and Q
i
(
3
) stand for water level and discharge at the ith interface of the 3D model, q
j
,
l is the discharge of the grids, m is the number of the grids of the interface, and o is the number of the layers of water depth.Due to the difference of water levels at grids of 3D model, the down-boundary of the 1D model is the averaged water levels of the grids in the 3D model:(15)
Z
i
(
1
)
=
Z
i
(
3
)
=
1
m
z
i
3
(
x
,
y
,
z
)
,
i
=
1,2
,
…
N
.Based on the method of the 1D hydrodynamics model, the governing equation can be deduced to matrix of water level at junctions as below:(16)
A
[
Z
]
=
B
,
where A is the coefficient matrix, Z is the water levels at junctions, and B is constants.However, in a coupled model, the interfaced water levels are not known, which lead toB containing unknown term and constant term. Therefore, the unknown part should be separated from the matrix B:
(17)
A
[
Z
]
=
C
[
H
]
+
D
,
where H is the water levels at the interfaces, C is the corresponding coefficient matrix, and D is the constant term.Gauss elimination method is used to simplify the coefficient matrixA to a diagonal matrix A:
(18)
A
′
[
Z
]
=
C
′
[
H
]
+
D
′
.Then, the water levels at junctions connected to the interfaces can be expressed by the water levels at interfaces of 1D model and 3D model:(19)
Z
i
(
1
)
=
∑
j
=
1
n
c
i
,
j
′
H
j
+
d
i
′′
,
i
=
1
,
…
n
,
where c
i
,
j
′′ and d
i
′′ are coefficient and constant after unitization, respectively.Utilizing the relations between the water level and discharge at junctions [12], the discharge at the beginning and the end of a single river channel can be deduced as below:
(20)
Q
i
,
1
=
α
i
,
1
+
β
i
,
1
H
i
,
1
+
γ
i
,
1
Z
i
,
n
+
1
,
k
k
k
k
k
k
k
k
k
k
i
=
1,2
,
…
,
N
(21)
Q
i
,
n
+
1
=
ξ
i
,
n
+
1
+
ζ
i
,
n
+
1
H
i
,
n
+
1
+
η
i
,
n
+
1
Z
i
,
1
,
k
k
k
k
k
k
k
k
k
k
k
k
k
k
i
=
1,2
,
…
,
N
.Equation (19) is substituted to (20) and (21), respectively; then the discharge at the interfaces of 1D and 3D models can be totally represented by the water level at the interfaces, and ξ
,
ζ
,
η are the coefficients of the equations. Consider the following:
(22)
Q
i
,
1
=
α
i
,
1
+
β
i
,
1
H
i
,
1
+
γ
i
,
1
(
∑
j
=
1
n
c
i
,
j
′′
H
i
,
1
+
d
i
′′
)
,
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
i
=
1,2
,
…
,
N
Q
i
,
n
+
1
=
ξ
i
,
n
+
1
+
ζ
i
,
n
+
1
H
i
,
n
+
1
+
η
i
,
n
+
1
(
∑
j
=
1
n
c
i
,
j
′′
H
i
,
n
+
1
+
d
i
′′
)
,
i
=
1,2
,
…
,
N
.After simplifying (22), the relations between water level and discharge can be established in 1D model. Consider the following:
(23)
Q
i
(
1
)
=
∑
j
=
1
n
p
i
,
j
H
i
+
q
i
,
i
=
1
,
…
,
n
,
and p is the coefficient of equations and q is constant.Finally, taking into account (14), an iterative equation for water level at interface is established as below:
(24)
H
i
(
l
+
1
)
=
1
p
i
,
j
[
Q
i
(
3
)
-
∑
j
=
1
j
≠
i
n
p
i
,
j
H
i
(
l
)
-
q
i
]
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
i
=
1
,
…
,
n
.At the beginning of the calculation, the initial dischargesQ
1
(
0
)
,
Q
2
(
0
)
,
…
,
Q
N
-
1
(
0
)
,
Q
N
(
0
) are used as the upper boundary condition of the 3D hydrodynamics model and can get the water levels H
1
(
0
)
,
H
2
(
0
)
,
…
,
H
N
-
1
(
0
)
,
H
N
(
0
). Then, substituting these water levels into (24) is to get the approximate values H
1
(
1
)
,
H
2
(
1
)
,
…
,
H
m
-
1
(
1
)
,
H
m
(
1
), which can be regarded as the downstream boundaries of 1D hydrodynamics if water levels satisfy the convergence condition MAX
|
H
i
(
1
)
-
H
i
(
0
)
|
≤
ε (ε is a given error). Otherwise, H
1
(
0
)
,
H
2
(
0
)
,
…
,
H
N
-
1
(
0
)
,
H
N
(
0
) have to be substituted into (23) to get discharges Q
1
(
1
)
,
Q
2
(
1
)
,
…
,
Q
N
-
1
(
1
)
,
Q
N
(
1
) for 3D model calculation again. This iterative process stops until the convergence condition is satisfied.As for the suspended sediment model, two additional conditions are needed to be considered. Firstly, the sediment fluxes at the interfaces of the 1D model and 3D model should be equal:(25)
Q
i
(
1
)
S
i
(
1
)
=
Q
i
(
3
)
S
i
(
3
)
,
i
=
1
,
…
,
N
.Secondly, the spatial averaged sediment concentration of 3D model should be closed to that of 1D model:(26)
S
i
(
3
)
(
x
,
y
,
z
)
=
S
i
(
1
)
,
S
i
(
1
)
=
S
i
(
3
)
=
1
m
s
i
(
3
)
(
x
,
y
,
z
)
,
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
i
=
1
,
…
,
N
.


## 2.4. Study Area and Data Processing

The Pearl River Delta (PRD) is located in the northern of the South China Sea. It is a dynamically complex estuarine system, which encompasses Pearl River network and PRE. The Pearl River, a subtropical river, is the second largest Chinese river in terms of annual water discharge (over 3 × 1011 m3). The channel network is formed by a large number of intricately interlaced deltaic branches, including three principal rivers: Xijiang, Beijiang, and Dongjiang, and some small rivers draining the PRD. It is a catchment area of 26,820 km2 and a total length of 1,600 km [16]. Meanwhile, as the distance from one junction to another is 0.68–1.07 km/km2, the Pearl River network is considered as one of the world’s most complicated river networks, which is formed by a large number of intricately interlaced branches [12]. In the river network, the mean suspended sediment concentration is about 0.284 kg/m3 [17], and annual mean sediment load in 1957 to 2006 is 75.3 × 109 kg/yr, with an annual average sediment load of 66.83 Mt/yr in Xijiang alone [18].These fluxes pass through the river network and flow into the PRE through the eight outlets Hutiaomen and Yamen discharging into the Huangmao Bay, Modaomen and Jitimen discharging into open sea directly, and Humen, Jiaomen, Hongqimen, and Hengmen discharging into the Lingding Bay [4]. The PRE has a broad bell shape, in which the Lingding Bay is like an inverted funnel with the narrow neck in the north and wide mouth opening to the south [19]. As the gravitational circulation develops in the PRE, the residual currents are landwards in deep channels and seawards in the western part of the PRE [20]. Normally, the circulation and water properties in the PRE are caused by various forcing mechanisms, including river discharges, tides, monsoon winds, coastal currents, and buoyancy forcing associated with the mixing of freshwater and saline water. The tides in the PRE mainly come from the Pacific oceanic tidal propagation [21], with a mean tidal range between 1.0 and 1.7 m. Thus PRE can be classified as a microtidal estuary. Propagating from offshore towards the estuary, the tidal range increases gradually and reaches a maximum near the river outlets [19]. Salinity in the PRE seldom exceeds 32 psu even during the dry season and is highly stratified in the vertical with a very low salinity surface layer. The distributions of salinity have a significant seasonal variability [7]. Further, the carrying capacity of sediments in the PRE is basically controlled by tidal flow, wave, fluid, and particle features. In the PRE, the very fine-grained silt and clay are contained in the suspended sediments, within which most particles are less than 100 μm in diameter and the median particle diameter d
50 is about 8 μm [22].


## 2.5. The Setup of the Coupled 1D and 3D Model

The model is applied to the Pearl River network and its estuary. The river network in 1D model is discretized into 347 channels, 13 boundary rivers, and 1,850 cross-sections with the interval varying from 0.2 to 3.0 km and 220 junctions. The length of the river network is about 1,600 km. The upstream boundaries include Shizui at the Tanjiang River, Gaoyao at Xijiang, Shijiao at Beijiang, Laoyagang at the Liuxi River, and Boluo at Dongjiang. The downstream boundaries include Dahu at Humen, Nansha at Jiaomen, Wanqingsha at Hongqili, Hengmen at Hengmen, Denglongshan at Maodaomen, Huangjin at Jitimen, Hutiaomen at Hutiaomen, and Huangchong at Yamen. These eight outlets are also the interfaces of the 1D and 3D model. The 3D estuarine model covers the area from the eight outlets to the −30 m isobaths (Figure2). More than 57,000 orthogonal curvilinear grids are used in the 3D model with spatial resolution from 0.05 km to 1 km and 10 equidistant sigma levels in vertical.Figure 2
Map of the study area including(
1
) the location of the Pearl River network with the names of the three major rivers and (
2
) the major topographic features of the Pearl River Estuary with isobaths in colors. The upstream boundaries of the 1D model are shown on the map, which starts with Shizui, Gaoyao, Shijiao, Laoyagang, and Boluo. The downstream boundaries of the 3D coastal model extend to the −30 m isobaths. Black circles represent major cities. Red triangles and red circles represent the main stations in the Pearl River network and the Pearl River Estuary, respectively.The hourly observed runoff at the upstream stations is the upper boundary conditions of the coupled model. The salinity at those stations is definitely zero. At the open sea part of the coupled model, the tidal levels are evaluated by the harmonic analysis from the most important tidal constituents of the PRE, with the uniform values of temperature and salinity from the observations. Consider the following:(27)
ζ
=
a
0
+
∑
i
=
1
8
a
i
cos
⁡
(
ω
i
t
-
ϕ
i
)
,
where a
0 is the mean water level, a
i is the amplitude of partial tide, and ϕ
i is the phase of partial tide.The simulation is initialized with zero water levels and velocities in the bathymetry of 1998-1999 and runs about 30 days. The results of the model are used as the new initial conditions and the simulation hot-started with the results. The processes are repeated several times until the error between the last two simulations satisfies a given condition. These results are then used as the initial conditions of the coupled model. Two different periods, including July 1998 (wet season) and December 1998 (dry season), are used for model calibration. In both the 1D and 3D model domains, the nonuniform parameters, such as bed roughness parameter and particle diameter, have been adjusted in order to optimize the water level, discharge, and suspended sediment concentration over the upstream river network and its coastal waters. The specific key parameters of the coupled model are shown in Table1.Table 1
Parameters of the modeling.
Settings
Description
1D
Time steps
Based on the scope of model, the time step of the 1D model is 10 s
Roughness coefficient
This coefficient depends on the physical properties of the riverbed. At Beijiang, from upstream boundary to Sanshuin
=  0.026~0.045, from Sanshui to downstream boundary n
=  0.015~0.03; at Xijiang n
=  0.02~0.03; at Dongjiang, from upstream boundary to Shilong n
=  0.025~0.045, from Shilong to downstream boundary n
=  0.015~0.036; at the mouth of the river network n
=  0.01~0.02
Particle diameter
In channel, the median particle diameter of suspended sediment is about 0.01~0.040 mm, and the median particle diameter of bed material is about 0.17~0.44 mm.
Sediment transport capacity
S
*
=
K
(
U
3
/
g
R
ω
)
m,  where m
=
0.92; K
=  0.0012~0.0016.
Recovery coefficient relating to the saturation
When the channels are in the deposition conditions, a value of 0.25 is applied; while, in the erosion conditions, a value of 1.0 is applied.
3D
Time steps
Δ
t
E
≤
1
/
C
t
|
(
1
/
δ
x
2
)
+
(
1
/
δ
y
2
)
|
(
-
1
/
2
), Δ
t
I
≤
1
/
C
T
|
(
1
/
δ
x
2
)
+
(
1
/
δ
y
2
)
|
-
1
/
2;  thus, the internal model is 10 s and the external model is 0.5 s.
Bottom roughness
τ
b
→
=
ρ
0
C
D
|
U
b
→
|
U
b
→, C
D
=
max
⁡
[
κ
2
/
ln
2
⁡
(
z
/
z
0
)
,
0.0025
], where κ
=
0.4, z
0
=
0.15 cm.
Particle diameter
The median particle diameter is 0.003~0.009 mm during the spring tide period.
Sediment settling velocity
ω
s
,
1
=
α
(
C
1
G
)
β, where α
=
2.42, β
=
0.22, and ω
s
,
2
=
υ
/
D
k
[
(
25
+
1.2
D
*
2
)
0.5
-
5
]
1.5


## 3. Results

The coupled model uses the field measurements in July 1998 (wet season) and December 1998 (dry season) at eighteen stations (Figure2) as calibration. The model skill score (SS) is employed to evaluate quantitatively the model performance of water level, discharge, and suspended sediment concentration against the measured time series:
(28)
SS
=
1
-
∑
i
=
1
N
(
X
mod
⁡
-
X
obs
)
2
∑
i
=
1
N
(
X
mod
⁡
-
X
-
obs
)
2
.This method is widely utilized by numerical modeling system [12, 23] and the performance levels are categorized as follows: SS
>
0.65 is excellent, 0.5–0.65 is very good, and 0.2–0.5 is good; if SS
<
0.2, it means a poor fit [24].The model skill scores for water levels are fairly high (Figures3(a) and 4(a)), indicating that the coupled model can accurately simulate water level in river network. Among the 18 stations, most of them performed excellently with only 4 stations scoring between 0.5 and 0.65 in the wet season (Figure 3(a)). In the dry season, 14 stations score more than 0.65 and the rest score between 0.5 and 0.65 (Figure 4(a)). In terms of the discharge, although the skill scores are not as excellent as that of water level, the model skill scores still show a very good skill of simulating discharge. In the wet season, the skill scores for most stations are greater than 0.5, implying a very good skill, although 4 stations score between 0.2 and 0.5 (Figure 3(b)). The coupled model also performed well in simulating discharge in dry season. Among the 18 stations observed, the skill scores of 13 stations are greater than 0.5 (Figure 4(b)), also showing a very good skill. Based on the agreement between the simulated data and the measured data throughout the river network, it can be concluded that the coupled model can simulate hydrodynamics process well in the river network and its estuary. Although the coupled model results do not perform quite as well as water level and the discharge, the simulated suspended sediment concentration in the Pearl River network also seems to be in a good agreement with the observed data. Figures 3(c) and 4(c) show that most of stations get their skill scores between 0.2 and 0.5 for both the wet and dry season. In the dry season two stations score below 0.2 and in the wet season three. There are further two stations in the wet season and one in the dry season that score greater than 0.5. Generally speaking, the simulated suspended sediment concentrations do not match the observed excellently but still satisfy the demand.The coupled model uses the field measurements from July 18, 1998, to July 19, 1998 (wet season), at eighteen stations as calibration. The model skill score is employed to evaluate quantitatively the model performance of water levels (a), discharges (b), and suspended sediment concentration (c) against the measured time series. The circles with the different colors represent the different performances of the model. The high scores of the water levels and the discharges show the good performance of the model, which include the contribution of the calibration of the suspended sediment concentration.
(a)
(b)
(c)The calibration of the coupled model results with the measured time series at eighteen stations is shown, from December 19, 1998, to December 20, 1998 (dry season). The model skill score is used to evaluate the model performance of water levels (a), discharges (b), and suspended sediment concentration (c) quantitatively. The circles with the different colors represent the different performances of the model. With the high scores of the water levels and the discharges, the performances indicate the success of the model. The skill scores of the suspended sediment concentration are also good enough.
(a)
(b)
(c)Figure5 shows the results of the skill scores of the velocity at different layers in the PRE. In general, the skill scores of the velocity range from 0.4 to 0.6, showing a good or a very good performance. What is more, the skill scores of the velocity at the surface layer are a little bit better than those at the bottom. As for the suspended sediments, although the skill scores of the suspended sediment concentration (Figure 6) are relatively smaller compared with the skill scores velocities (Figure 5), the skill scores of all the stations are above 0.2, indicating that the model does a good job to simulate the suspended sediment concentration in the PRE. The skill scores of suspended sediment for some stations are even greater than 0.5, showing a very good performance. The PRE is a sedimentary system characterized by intricate depositional structures in space and time [11]. Due to the intensive anthropogenic activity, the bathymetry of the entire Pearl River network has changed dramatically. In particular, because of the large-scale sand excavation, the impacts of the anthropogenic activity cannot be neglected in the investigation of the erosion and deposition in the river network and its estuary. Figure 7 quantifies the impacts of the sand excavation in the bathymetry evolution. It compares the sediment erosion between the sand excavation and the natural conditions from Sixianjiao to Chaolianzhoutou in 1998. Figure 7 provides an insight into the amount of sediment in reaches under the sand excavation and the natural erosion and deposition. It is clear that the amount of the sand excavation is close to that of sediment erosion in Reach 1 and Reach 2 (Figure 8(a)) and larger than the amount of sediment erosion in Reach 3 and Reach 4. In Reach 4, the amount of sand excavation is almost 2.5 times that of the natural erosion. Therefore, the amount of sand excavation has to be taken into account in the comparison between the observed and the modelled sediment erosion and deposition in the river network. Figure 8 shows the thickness of simulated sediment erosion and deposition in the aforementioned reaches in 1998. The results show that erosion generally occurs in these channels in a year, although the river channel sedimentation plays a dominant role in the wet season. There is some discrepancy between the simulated and the measured river bed (Figure 8(b)), and the results still show that the model can be used to evaluate sediment transport in this region.Figure 5
The calibration of velocity between the model results and observations at 9 stations in the Pearl River Estuary is shown. The calibration of the velocity at the surface (blue), the middle (yellow), and the bottom (green) is deployed. The result shows a good performance with the skill scores of the velocity ranging from 0.4 to 0.6.Figure 6
The comparison of suspended sediment concentration for calibration at stations in the Pearl River Estuary is displayed. The skill scores of the surface (blue), the middle (yellow), and the bottom (green) show the good performance of the model with the skill scores of all the stations being above 0.2.Figure 7
From Sixianjiao to Chaolianzhoutou, the annual amount of erosion and deposition in the sand excavation (up) and under the natural conditions (down) is provided in 1998. The comparison quantifies the impacts of the sand excavation in the bathymetry evolution and indicates that the amount of sand excavation has to be taken into account in the comparison between the observed and the modelled sediment erosion and deposition.The map of the part of the Pearl River network (a) includes the locations of the main river reaches and stations. The comparison of the observed and simulated thickness of erosion or deposition (b) in 1998 is exhibited, which shows that the performance of the model can be used to evaluate sediment transport in this region.
(a)
(b)


## 4. Discussion

The results show that the model skill score for sediment concentration is lower than that for water level and discharge. There are some potential reasons for the relative larger discrepancy between the simulated and observed value. First of all, more parameters in sediment simulation need to be selected than that of water level and discharge simulation. Normally, the roughness coefficient is the most important parameter in hydrodynamics simulation. However, large domain simulation in suspended sediment transport needs more data and parameters, such as the particle diameter and sediment settling rate. These parameters have large influences on the suspended sediment concentration simulation. However, due to the lack of data, some of these parameters are treated as relatively uniform distribution over the Pearl River network. This leads to error when simulating sediment transport. Secondly, the bathymetry of the Pearl River network changes dramatically due to the intensive sand excavation during the last three decades. Therefore, the asynchrony of hydrographic and bathymetry measurement may directly lead to the inaccuracy of the simulation results, especially in suspended sediment transport simulation. In fact, in our model, the field measurements were taken in 1998, but the bathymetry survey was done in 1999. The asynchrony of hydrographic and bathymetry measurement should be another reason for the simulation errors. Finally, some complex impacts, such as wind and discharge let-out by nearby reservoir, are not considered in this model, which probably affect the results of the coupled model.This model was then used to simulate the flow division in the PRD, because flow division at some important junctions is important for the flow and sediment transport in distributary channel networks. Makou and Sanshui Stations are located at the apex of the delta. Their water discharge enters the Xijiang and the Beijiang, respectively. Table2 shows the diversion ratio of flow at the two main junctions in the Pearl River network: (a) Makou and Sanshui and (b) Tianhe and Nanhua. The diversion ratios of flow between Makou and Sanshui Stations are 75.7% and 24.3%, respectively, showing the disequilibrium of the flow distribution in the waterway network. At the same time, the diversion ratios of flows at Tianhe and Nanhua are relative equal with the 54.2% in Tianhe and 45.8% in Nanhua. Although the flow flux displays an obvious seasonal variation, the diversion ration of flow is relative stabile. The flow flux at Makou Station is 1185 × 108 m3 in the wet season and 429 × 108 m3 in the dry season. The diversion ratio of flow is 74.5% in the wet season and 79.3% in the dry season. The diversion ratios of flows at the Tianhe and Nanhua Stations are also relatively stable in the wet season and the dry season. It can be observed that diversion ratio of the suspended sediment transport is similar to the flow ratio (Table 3). More sediment enters the Xijiang River as the flow distributes at the Makou and Sanshui station. Although the amount of sediment is much lower in dry season than that in the wet season, the majority of sediment still enters Xijiang River. Table 4 shows the annual flow flux of the eight outlets of the PRD. It can be noticed that Modaomen, with the annual water flux of 623 × 108 m3, is the largest contributor of the flow flux to the sea, accounting for 24.8% of the total flow flux from the delta to the sea. Humen comes second, accounting for the 21.4% of the total flow flux to the sea. Flow flux from Jiaomen outlet is 515 × 108 m3, accounting for 20.5%. Flow fluxes from the rest of the outlets are relatively small, which are 205 × 108 m3, 301 × 108 m3, 64 × 108 m3, 79 × 108 m3, and 187 × 108 m3, respectively. The annual sediment fluxes from outlets to the sea are in general consistent with the individual flow ratio. Modaomen is still the most important outlet to transport sediment from delta to the sea, accounting for nearly 30% of the total sediment flux to the sea. Following Modaomen, there are Humen and Jiaomen, both of which transport more than 40% sediment from PRD to the sea. The changes in sediment transport in the Pearl River delta and its estuary have great impacts on the water quality in this region. Therefore further study should be conducted to improve the water quality of the river. Nowadays, the application of nanotechnology for the remediation of contaminants may give promising results in future [25, 26]. An insight of the flow and sediment transportation mechanism in the river and estuary system has been obtained through the paper, and a reliable basis to control water pollution of the Pearl River Delta can be proposed combined with the present studies.Table 2
The water flux and the diversion ratio at the two main joints in the Pearl River network: (a) Makou and Sanshui and (b) Tianhe and Nanhua (unit: ×108 m3).
Bathymetry
1999
Station
Makou
Sanshui
Tianhe
Nanhua
Annual water flux
2132
1521
Wet
1185 [73.4%]
406 [78.3%]
598 [72.5%]
515 [74.0%]
Diversion ratio
74.5%
25.5%
53.7%
46.3%
Dry
429 [26.6%]
112 [21.7%]
227 [27.5%]
181 [26.0%]
Diversion ratio
79.3%
20.7%
55.6%
44.4%
Annual
1614 [100%]
518 [100%]
825 [100%]
696 [100%]
Diversion ratio
75.7%
24.3%
54.2%
45.8%Table 3
The suspended sediment flux and the diversion ratio at the two main joints in the Pearl River network: (a) Makou and Sanshui and (b) Tianhe and Nanhua (unit: ×104 t).
Bathymetry
1999
Station
Makou
Sanshui
Tianhe
Nanhua
Annual sediment flux
1494
1052
Wet
1198 [96.8%]
253 [98.8%]
553 [96.3%]
461 [96.4%]
Diversion ratio
82.6%
17.4%
54.5%
45.5%
Dry
40 [3.2%]
3 [1.2%]
21 [3.7%]
17 [3.6%]
Diversion ratio
93.0%
7.0%
55.3%
44.7%
Annual
1238 [100%]
256 [100%]
574 [100%]
478 [100%]
Diversion ratio
82.9%
17.1%
54.6%
45.4%Table 4
The simulated flow flux and suspended sediment flux at the eight outlets in 1999.
Outlet
Humen
Jiaomen
Hongqimen
Hengmen
Maodaomen
Jitimen
Hutiaomen
Yamen
Total
Flow flux(unit: ×108 m3)
Annual
537
514
205
301
623
64
79
187
2510
[21.4%]
[20.5%]
[8.2%]
[12.0%]
[24.8%]
[2.5%]
[3.1%]
[7.5%]
[100%]
Wet
401
381
156
223
447
46
59
123
1836
[21.8%]
[20.8%]
[8.5%]
[12.1%]
[24.4%]
[2.5%]
[3.2%]
[6.7%]
[100%]
Dry
136
133
49
78
176
18
20
64
674
[20.2%]
[19.7%]
[7.3%]
[11.5%]
[26.1%]
[2.7%]
[3.0%]
[9.5%]
[100%]
Suspended sedimentflux (unit: ×104 t)
Annual
208
199
90
115
289
17
24
41
983
[21.2%]
[20.2%]
[9.2%]
[11.7%]
[29.4%]
[1.7%]
[2.4%]
[4.2%]
[100%]
Wet
157
185
85
111
258
15
22
26
859
[18.3%]
[21.5%]
[9.9%]
[12.9%]
[30.0%]
[1.8%]
[2.6%]
[3.0%]
[100%]
Dry
51
14
5
4
31
2
2
15
124
[41.1%]
[11.3%]
[4.1%]
[3.2%]
[25%]
[1.6%]
[1.6%]
[12.1%]
[100%]


## 5. Conclusions

In deltaic estuaries, the exchanges of the mass flux between river network and estuary require the researchers to consider these two systems as a whole. Therefore, this paper provides a mathematical method to connect the 1D river networks with the 3D estuarine model to study the flow and suspended sediment transport in this region. Based on the classic Saint-Venant equations and the nonequilibrium transport equation, 1D flow and suspended sediment model is established. Based on the orthogonal curve coordinate and the ECOM model, 3D estuarine model is established to describe the complex horizontal and vertical variability of the flow and suspended sediment transport. The relationship expression of water level at interfaces of 1D model and 3D model is strictly deduced by utilizing the relations between the water level and discharge at junctions. The 1D and 3D hydrodynamics models are then coupled by iterative computations. Finally, the 1D and 3D suspended sediment models are also coupled by adding two additional conditions. The advantage of this coupled method is that the schematisations do not overlap at the interface of the two systems, which can improve the stability of modelling.The coupled model has been applied successful in the Pearl River network and its estuary, which is a complicated system. The simulated results are compared with the field measurements to evaluate the accuracy of the coupled model. The results show that the coupled model is capable of simulating water levels, discharges, and the suspended sediment concentration in the river network, as well as the velocity and suspended sediment at the shallow estuary. Generally speaking, the coupled model is elaborate enough to simulate water levels, discharges, and velocity in the study area. Due to the limit of the data and the simplification of the model, although the simulated suspended sediment is not as good as water level and discharge, it is still acceptable to investigate the suspended sediment transport in river network and estuary system.



---

*Source: 798579-2014-08-12.xml*