# Quantification of Multiple Climate Change and Human Activity Impact Factors on Flood Regimes in the Pearl River Delta of China

Authors: Yihan Tang; Shufeng Xi; Xiaohong Chen; Yanqing Lian
Journal: Advances in Meteorology
 (2016)

---

## Abstract

Coastal flood regimes have been irreversibly altered by both climate change and human activities. This paper aims to quantify the impacts of multiple factors on delta flood. The Pearl River Delta (PRD), with dense river network and population, is one of the most developed coastal areas in China. The recorded extreme water level (m.s.l.) in flood season has been heavily interfered with by varied income flood flow, sea-level rise, and dredged riverbeds. A methodology, composed of a numerical model and the indexR, has been developed to quantify the impacts of these driving factors in the the PRD. Results show that the flood level varied 4.29%–53.49% from the change of fluvial discharge, 3.35%–38.73% from riverbed dredging, and 0.12%–16.81% from sea-level rise. The variation of flood flow apparently takes the most effect and sea-level rise the least. In particular, dense river network intensifies the impact of income flood change and sea-level rise. Findings from this study help understand the causes of the the PRD flood regimes and provide theoretical support for flood protection in the delta region.


---

## Body

## 1. Introduction

Under the impacts of climate change and intensive human activities, flood regimes in deltas have been irreversibly changed [1–4]. Changes in income flood flow of deltas have been detected across the world in the last decade [5, 6]. Sea-level rise has also affected flood regimes in deltas all over the world [7]. At the same time, fluvial geomorphologies of many deltas have been severely altered by human interventions [8, 9]. Human activities like the construction of municipal water infrastructure and floodplain engineering, along with climate changes such as inflow variation and sea-level rise, are threatening the safety of coastal habitats all over the world [10, 11]. Under the integrated impacts of climate change and human activities, magnitude and frequency of coastal floods have changed, and deltas are more exposed to flood risk than ever [12, 13].To better understand the cause and effect of delta flood variation, a number of previous works have intended to quantify the fluvial impacts of different causes. Impact of sea-level rise has been investigated the most. Numerical model was built to predict the return period of flood level in tidal reaches under the effect of sea-level rise [14]. Percentage of flood submerged area was calculated by generating different sea-level rise scenarios [15]. Flood vulnerable index under various sea-level rise scenarios was estimated, and the differences among estimations were taken as the quantitative impact of climate change [16]. Quantitative effect of storm surge under the impact of sea-level rise was also studied [17]. In particular, impacts of human activities like urbanization were studied [18]. The impact of land-cover changes was considered in a hydrological model [19]. Furthermore, effects of sediment variation as well as riverbed dredging were studied [20, 21].Most previous works attended to the influence of a single factor. Somehow, floods in coastal area are simultaneously affected by multiple factors, which are also mutually affected [22]. Quantitative cause analysis for multiple factors has yet to be done. In this study, therefore, we intend to quantify the impacts of multiple factors on coastal flood. Effects of different factors are also to be compared. A multilinear regression model is proposed to estimate the on-site flood level. Index R is also defined to calculate the impacts of different factors. This research offers a simple and efficient approach to quantify impacts of climate change and human activities on coastal flood.


## 2. Study Area and Data

### 2.1. Study Area

Pearl River Delta (21°30′~23°40′N, 109°40′~117°20′E), covering an area of 39,380 km2, is the third largest delta in China (Figure 1). The altitude of terrain is higher in its northwest than the southeast. River channels in this area are all tidal rivers; that is, the in-channel flood level is affected by both the flood flow from upstream and the tidal waves at outlets. West River (WR) and North River (NR) are the upstream rivers of the the PRD. In the past 20 years, the 6 greatest floods in WR and NR all happened in June and July (Table 1). High tides in the South China Sea mainly happen in August and November within a year. Affected by both the upstream flood and downstream high tides, flood season in the the PRD lasts from June to August.Table 1
Big floods in the past 20 years in Pearl River Delta.
Year
Flood river
Peak time
Gaugingstation
Peak volume (m3/s)
1994
West River
Jun. 19
Wuzhou
41,100
1998
West River
Jun. 28
Wuzhou
51,600
2005
West River
Jun. 22
Wuzhou
53,900
North River
Jun. 24
Shijiao
13,500
2006
North River
Jul. 17
Shijiao
17,500
2008
West River
Jun. 15
Wuzhou
46,000Figure 1
Map of Pearl River Delta and locations of gauging stations.River network in the the PRD area is extremely complex. Inflow from WR and NR runs into the the PRD through Makou and Sanshui. Makou and Sanshui are connected by Sixianjiao Channel. Through dense river network around Jiangmen and Zhuyin, flood flow runs into the South China Sea through eight outlets.Cities as Guangzhou and Foshan locate alongside the dense river networks, and metropolises like Macau and Zhuhai gather around the outlets. Rapid economic development and dense population have made the PRD extremely vulnerable to flood risk. Ever since the “reform and open up” policy in 1979, GDP growth in the PRD has stayed around 8% per year and the local population has annually increased by more than 10%. By 2011, the PRD already had a population of 120 million and a GDP of 502.5 billion dollars. For the use of municipal construction, large scale of sand in riverbed was dredged from the mid 1980s to the early 21st century. Large-scale sand-dredging has caused irreversible change in the regional river channels.


### 2.2. Dataset

Hydrological data of fifteen gauge stations in the PRD were used (Table2). Recorded flood flows from upstream river basins and flood level in tidal reaches (m.s.l.) were both applied. Topographic statistics of the PRD river network from 1980 to 2006 were also applied. All the flow data was provided by the Hydrology Bureau of Guangdong Province, China (HBGPC), and all the stage data came from the Chinese Hydrological Almanacs. All the topographic statistics were obtained from the Pearl River Commission. The homogeneity and reliability of all applied data were strictly examined and controlled before release.Table 2
Dataset of hydrological series.
Section
Number
Abbr.
Station
Elevation from the Pearl Delta datum (m)
Period
Flow
Stage
Upstream
1
wz
Wuzhou
1956–2011
Daily; Monthly Maximum
None
2
sj
Shijiao
Income
3
mk
Makou
-
0.003
1959–2011
Daily; Monthly Maximum
Monthly Maximum
4
ss
Sanshui
0.011
Network
5
th
Tianhe
0.059
1958–2011
None
6
jm
Jiangmen
0.067
7
zy
Zhuyin
0.004
8
rq
Rongqi
-
0.033
9
nh
Nanhua
0.018
1956–2011
10
sd
Sanduo
-
0.079
1959–2011
Outlets
11
dls
Denglongshan
0.022
1959–2011
12
ssk
Sanshakou
0.034
13
nh
Nansha
0.007
14
wqsx
Wanqinshaxi
-
0.088
15
hm
Hengmen
-
0.069
Note: number is the station number shown in Figure1.


## 2.1. Study Area

Pearl River Delta (21°30′~23°40′N, 109°40′~117°20′E), covering an area of 39,380 km2, is the third largest delta in China (Figure 1). The altitude of terrain is higher in its northwest than the southeast. River channels in this area are all tidal rivers; that is, the in-channel flood level is affected by both the flood flow from upstream and the tidal waves at outlets. West River (WR) and North River (NR) are the upstream rivers of the the PRD. In the past 20 years, the 6 greatest floods in WR and NR all happened in June and July (Table 1). High tides in the South China Sea mainly happen in August and November within a year. Affected by both the upstream flood and downstream high tides, flood season in the the PRD lasts from June to August.Table 1
Big floods in the past 20 years in Pearl River Delta.
Year
Flood river
Peak time
Gaugingstation
Peak volume (m3/s)
1994
West River
Jun. 19
Wuzhou
41,100
1998
West River
Jun. 28
Wuzhou
51,600
2005
West River
Jun. 22
Wuzhou
53,900
North River
Jun. 24
Shijiao
13,500
2006
North River
Jul. 17
Shijiao
17,500
2008
West River
Jun. 15
Wuzhou
46,000Figure 1
Map of Pearl River Delta and locations of gauging stations.River network in the the PRD area is extremely complex. Inflow from WR and NR runs into the the PRD through Makou and Sanshui. Makou and Sanshui are connected by Sixianjiao Channel. Through dense river network around Jiangmen and Zhuyin, flood flow runs into the South China Sea through eight outlets.Cities as Guangzhou and Foshan locate alongside the dense river networks, and metropolises like Macau and Zhuhai gather around the outlets. Rapid economic development and dense population have made the PRD extremely vulnerable to flood risk. Ever since the “reform and open up” policy in 1979, GDP growth in the PRD has stayed around 8% per year and the local population has annually increased by more than 10%. By 2011, the PRD already had a population of 120 million and a GDP of 502.5 billion dollars. For the use of municipal construction, large scale of sand in riverbed was dredged from the mid 1980s to the early 21st century. Large-scale sand-dredging has caused irreversible change in the regional river channels.


## 2.2. Dataset

Hydrological data of fifteen gauge stations in the PRD were used (Table2). Recorded flood flows from upstream river basins and flood level in tidal reaches (m.s.l.) were both applied. Topographic statistics of the PRD river network from 1980 to 2006 were also applied. All the flow data was provided by the Hydrology Bureau of Guangdong Province, China (HBGPC), and all the stage data came from the Chinese Hydrological Almanacs. All the topographic statistics were obtained from the Pearl River Commission. The homogeneity and reliability of all applied data were strictly examined and controlled before release.Table 2
Dataset of hydrological series.
Section
Number
Abbr.
Station
Elevation from the Pearl Delta datum (m)
Period
Flow
Stage
Upstream
1
wz
Wuzhou
1956–2011
Daily; Monthly Maximum
None
2
sj
Shijiao
Income
3
mk
Makou
-
0.003
1959–2011
Daily; Monthly Maximum
Monthly Maximum
4
ss
Sanshui
0.011
Network
5
th
Tianhe
0.059
1958–2011
None
6
jm
Jiangmen
0.067
7
zy
Zhuyin
0.004
8
rq
Rongqi
-
0.033
9
nh
Nanhua
0.018
1956–2011
10
sd
Sanduo
-
0.079
1959–2011
Outlets
11
dls
Denglongshan
0.022
1959–2011
12
ssk
Sanshakou
0.034
13
nh
Nansha
0.007
14
wqsx
Wanqinshaxi
-
0.088
15
hm
Hengmen
-
0.069
Note: number is the station number shown in Figure1.


## 3. Methodology for Quantitative Cause Analysis

“In the PRD, varied income flood flow, sea-level rise and sand mining are the most important driving factors for flood regime variation” [14, 23–26]. In order to quantify the impacts of these three factors, this study proposed a methodology composed of a multilinear regression model and an index R. The regression model was based on the one proposed by Huang et al. [14], and it was used to estimate the on-site flood levels. The index R, improved from Jiang’s theory in 2011 [31], was defined to calculate the impacts from different driving factors.


### 3.1. Multilinear Regression Model

There are mainly four approaches to study the influence of multiple factors, that is, mathematical model, field observation, qualitative theoretical reasoning, and hydrodynamic model [27, 28]. However, qualitative theoretical reasoning cannot provide quantitative results. Field observation costs great, and it is beyond the financial reach of this research. At the same time, large amounts of terrain data from field observation are required to build a precise hydrodynamic model for a delta with complex river network. Therefore, three approaches were denied, and mathematical model was chosen for this research.


#### 3.1.1. Modification of Huang’s Model

In 2004, Huang et al. proposed a multilinear regression model to estimate the backwater effect in the PRD under different scenarios of income discharge [14]. The annual extreme water level in any gauge station was supposed to be a dependent variable of the income discharge and the annual maximum tide level at outlets. And the on-site annual maximum water level can be represented as(1)
Z
i
=
f
Q
,
Z
=
C
1
Q
+
C
2
Z
+
α
Q
Z
+
C
0
,where Z
i is the annual extreme water level of a gauge station; Z is the annual maximum tide level at the outlet; Q is the river discharge from WR detected at Makou Station; α and β are the parameters; C
1 and C
2 are the constants. Parameters and constants together represent the regional topographic feature. Given the rise in sea level as Δ
Z, the variability of annual maximum water level in the gauge station can be therefore estimated by the following equation:(2)
Δ
Z
i
=
f
Q
,
Z
′
-
f
Q
,
Z
=
C
1
Q
+
C
2
Z
2
+
Δ
Z
+
α
Q
Z
+
Δ
Z
+
C
0
-
C
1
Q
+
C
2
Z
+
α
Q
Z
+
C
0
=
C
2
+
α
Q
Δ
Z
.Despite its reliable application in the PRD, there are some underlying limitations in Huang’s model. In (1), Z
i is only affected by tide at one outlet. However, flood flow runs towards sea from multiple outlets, and sea-level rise acts differently on tidal levels at difference outlets, even in the same delta. In this case, the impact from sea-level rise should be represented by tidal levels at multiple outlets. Likewise, Huang’s equation takes only the impact from one income flow source into consideration. Somehow, discharge at a gauge station usually comes from more than one upstream river. Impacts from different inflows should therefore be taken into consideration.In this way,C
1
Q
+
C
2
Z has been developed into ∑
i
=
1
n
C
i
Q
j
+
∑
j
=
n
+
1
n
+
m
C
j
Z
j, where n and m are the numbers of inflow resources and outlets and C
i, C
j, the parameters related to the local topographic features. α
Q
Z in (1) has been modified into α
∏
i
=
1
n
Q
i
∏
j
=
n
+
1
n
+
m
Z
j, and α also relates with the local topography. At the same time, the time scale in Huang’s model is year. Somehow, in our research, it has become month, that is, June, July, and August in the flood season. Therefore, the time scale of inputs and outcomes has been narrowed down. As a result, Huang’s model has been modified into the following equation:(3)
Z
=
f
Q
i
,
Z
j
=
∑
i
=
1
n
C
i
Q
i
+
∑
j
=
n
+
1
n
+
m
C
j
Z
j
+
α
∏
i
=
1
n
Q
i
∏
j
=
1
m
Z
j
+
C
0
,where Z is the monthly extreme flood level of the studied gauge station; Q
i is the income flood volume; Z
j is the monthly extreme tide level at estuary of the PRD; α, C
0, C
i, and C
j are the constants representing topographic features.


#### 3.1.2. Model Adjustment

Based on (3), regression models for all the stations in the PRD have been built (Table 3). All the the PRD stations have first been divided into three groups according to their locations. Flood from upstream runs through Makou and Sanshui Station into delta, and therefore they are income stations. Flood volumes in income stations have been taken as Q
i. Stations located within 5 km to the river mouth in delta are outlet stations. Tide levels of these stations have been applied in models as Z
j. Other stations in the PRD are network stations. Flood levels in these stations are Z, which is the dependent variable of Q
i, Z
j, and the parameters α, C
0, C
i, and C
j. For all the equations in Table 3, subscriptions of Q and Z are abbreviations of corresponding gauge stations, and the abbreviations for all gauge stations are listed in Table 2.Table 3
Regression model of the PRD stations.
Group
Station
Equation
Income stations
Makou
Z
mk
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
α
Q
mk
Q
ss
Z
dls
+
C
0
Sanshui
Z
ss
=
C
1
Q
s
s
+
C
2
Q
m
k
+
C
1
Z
s
s
k
+
C
2
Z
ns
+
C
3
Z
wqsx
+
α
Q
s
s
Q
m
k
Z
s
s
k
Z
n
s
Z
w
q
s
x
+
C
0
Network stations
Tianhe/Jiangmen/Zhuyin
Z
th/jm/zy
=
C
1
Q
m
k
+
C
2
Q
ss
+
C
3
Z
d
l
s
+
α
Q
m
k
Q
ss
Z
d
l
s
+
C
0
Nanhua/Rongqi
Z
nh/rq
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
wqsx
+
C
4
Z
hm
+
α
Q
mk
Q
ss
Z
wqsx
Z
hm
+
C
0
Sanduo
Z
sd
⁡
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
ssk
+
C
4
Z
ns
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
ssk
Z
ns
Z
wqsx
Z
hm
+
C
0
Outlet stations
Denglongshan
Z
dls
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
ssk
+
C
4
Z
ns
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
ssk
Z
ns
Z
wqsx
Z
hm
+
C
0
Nansha
Z
ns
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ssk
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
dls
Z
ssk
Z
wqsx
Z
hm
+
C
0
Sanshakou
Z
ssk
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ns
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
dls
Z
ns
Z
wqsx
Z
hm
+
C
0
Wanqinshaxi
Z
wqsx
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ns
+
C
5
Z
ssk
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
dls
Z
ns
Z
ssk
Z
hm
+
C
0
Hengmen
Z
hm
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ns
+
C
5
Z
ssk
+
C
6
Z
wqsx
+
α
Q
mk
Q
ss
Z
dls
Z
ns
Z
ssk
Z
wqsx
+
C
0As for the income stations, the model has been adapted. Since income stations are connected by Sixianjiao Channel, both flood flows and flood levels in these two income stations are mutually affected. In addition, flood flow from Makou goes into the South China Sea through Denglongshan and Sanshui through the other four outlet stations. Therefore, flood level in Makou is also affected by the tide level in Denglongshan. Likewise, flood level in Sanshui is related to flood flows from both income stations as well as and the tide levels in Sanshakou, Nansha, Wanqinshaxi, and Hengmen.The model of Zhuyin Station is taken as an example for network stations. Due to the connection between Makou and Sanshui, flood from both income stations, taken asQ
m
k and Q
s
s, would affect the flood level in Zhuyin. And since flood in Zhuyin runs into South China Sea through Denglongshan, tide level in Denglongshan Z
d
l
s takes effect on Z
z
y, that is, the flood level in Zhuyin.The same as the network stations, flood levels in all the outlet stations are affected by the flood flow from both Makou and Sanshui. Since flood levels in all the outlet stations are simultaneously affected by the tidal wave in the South China Sea, flood levels in outlet station in the the PRD are mutually affected. Therefore, to estimateZ for each outlet station, the flows in both Makou and Sanshui are taken as Q
j and the water level in other outlet stations is taken as Z
j.


#### 3.1.3. Model Validation

For each station,Q
i in income stations, Z
j in outlet stations, and the on-site extreme level Z of the same month within the flood season in a year are identified as a sample. For each year, there are three samples, that is, the one in June, July, and August. Take Zhuyin, for example; Q
1 and Q
2 are the peak flood flow in Makou and Sanshui in June 1959, Z
1 is the extreme tide level in Denglongshan in June 1959, and Z is the flood level in Zhuyin. Therefore, Q
1, Q
2, Z
1, and Z together are taken as a sample. In the sample, Q
i and Z
j are instances and Z is the target.Cross-validation (CV) makes full use of data information [29], and the validation and accuracy of regression model have been tested by the 10-fold cross-validation method [30]. For each gauge station, we randomly selected 90% of 159 samples to estimate the coefficients and tested the accuracy of the at-site model on the remaining 10% of data. Mean Squared Error (MSE) was applied to evaluate the accuracy of each test. After the procedure was repeated for 10 times on different data partition, a series of 10 MSEs was achieved. The average of this MSE series, that is, CV mean, was used to measure the accuracy of the regression model. The standard deviation, that is, CV variance, represented the robustness of the model. The lower the CV mean is, the more accurate the model is. Likewise, the smaller the CV variance is, the more robust the model is. Since Huang’s model in 2004 was successfully applied and shared, cross-validation results of Huang’s model were taken as the benchmark. Regression model with a lower CV mean and variance than the benchmark was proved eligible. The whole CV test was performed in MATLAB.


### 3.2. Changing RatioR

#### 3.2.1. Jiang’s Idea

IndexR defined by Jiang et al. [31] was used to quantify the impacts of different driving factors on monthly flow of a small river basin in China. The local hydrological cycle was severely disturbed by human activities and the time series was no longer consistent. They split the time series into two according to the timing of severe human interference, that is, one before the severe human interventions and one afterwards. The previous period was anatural period, and the latter was ahuman-affected period. First they calculated the model constants with data fromnatural period, and then they estimated monthly flow in thehuman-affected period with these constants. The estimated monthly flow represented the monthly flow in thehuman-affected period if none of the severe human activities happened. The gap between the estimated and recorded flow in thehuman-affected period was the impact of human activities on monthly flow.


#### 3.2.2. Division of Time Series

Based on Jiang’s theory [31], index R was defined to show the impacts of three driving factors. The difference between the factor-affected flood level and the untouched one represents the impact of a certain factor. To estimate the untouched series, hydrological series need to be divided first. Data length should be unified so that results in different stations can be compared. Since the available data varied among gauge stations (Table 2), we have only applied data from 1959 to 2011. In order to lower the uncertainties in analysis, the time series have only been divided into two periods. To find the dividing time for all three driving factors, we have studied the temporal pattern of these factors.Changes in upstream flow have been detected. Wuzhou, located downstream of West River, is more than 300 km away from Makou. Shijiao, located downstream of North River, is 50 km away from Sanshui. Mann-Kendall test was applied in detecting changes in recorded flow from 1959 to 2011 of these four stations (Table4). Flood flow increased in Wuzhou and decreased in Shijiao, yet neither of them reached the significance of 90%. Flood flow in both Makou and Sanshui has increased at a 95% significant level. It can be concluded that flood flows in upstream river basins do not have significant trend. As for the abruption in series, flood flow in Wuzhou jumps in 1991, and that in Sanshui shows no significant abruption. Flood flow in Makou stays consistent, while that in Sanshui took a jump in 1992. Therefore, there was abruption in West River flow in the early 1990s, probably due to dam construction upstream [32].Table 4
Mann-Kendall test results of income flood flow.
Station
West River
North River
Wuzhou
Makou
Shijiao
Sanshui
Mann-Kendall statistics
1.55
1.97
−1.12
2.03
Year of abruption
1991
—
—
1992
Note: Mann-Kendall statistics is known asZ. Positive Z indicates an increasing trend just as a negative one refers to a decreasing trend. Moreover, Z
>
1.65, Z
>
1.96, and Z
>
2.57 show tendencies at the 0.01, 0.05, and 0.1 significance levels.The inconsistent change of flow between Wuzhou and Makou is partly due to dyke construction and variations in precipitation. In the early 1990s, dykes were built in large scale around Wuzhou [32], and the flood flow from Wuzhou to Makou is locked within the channel. Furthermore, precipitation in flood season increased faster on the catchment area of Makou than Wuzhou within the last 50 years [33]. Therefore, flood flow in Makou has increased more severely than that in Wuzhou.The discrepancies between Wuzhou and Makou, Shijiao and Sanshui are mainly due to severe sand-dredging in the PRD. Table5 shows the variation of channel geometries in both Sanshui and Makou. Apparently, the average channel depth in Sanshui increased 0.5 m more than that in Makou from 1985 to 2006. Channel volume in Sanshui therefore increased by 69% during 1985–1999 and 36% during 1999–2006. Meanwhile, the percentage in Makou is only 18% and 29%. From 1980 to 2003, discharge area over both 2 m and 6 m in Sanshui increased three times that in Makou. Due to different extent of riverbed alteration, flood flows from Makou to Sanshui through Sixianjiao Channel. As a result, flow partition ratio between 1959 and 2011 in flood seasons increased in Sanshui (Figure 2). The increase became sharp from the early 1990s. Since mean flow in WR triples that of NR (Table 6), increase of flow partition ratio for more than 5% in Sanshui has significantly increased the on-site flow.Table 5
Variation of channel geometries in Sanshui and Makou.
Geometry index
Unit
Period
Makou
Sanshui
Maximum cutting depth
m
1980–2003
+7.11
+8.07
Average depth
1985–19991999–2006
0.802.00
2.801.50
Channel volume
1985–1999
+18.00
+69.00
%
1999–2006
+29.00
+36.00
Change of discharge area
>2 m
1980–2003
+19.30
+99.50
>6 m
1980–2003
+15.20
+40.00
Note: “+” means increase.Table 6
Average monthly flow in West River and North River in flood season (units: m3).
River
Station
Jun.
Jul.
Aug.
West River
Wuzhou
25,222
24,971
20,043
Makou
27,658
26,780
20,564
North River
Shijiao
7,406
4,844
3,501
Sanshui
7,229
6,917
4,893Figure 2
Tendency of flow partition ration of Sanshui.MK test was also carried out on the average tide level (m.s.l.) of the outlet stations. Results of MK in Sanshakou are 2.15, 1.75 for Nansha, 1.83 for Wanqinshaxi, 3.10 for Hengmen, and 0.08 for Denglongshan. All results are positive, which indicate increasing tendency in all outlets. According to the previous study, sea level at South China Sea has kept rising since the early 1960s [7, 14, 25]. Among all the outlet stations, increase in Denglongshan has not reached a significance of 1%. Due to the wandering rim of continental shelf, sea-level rise in South China Sea has caused different extent of increase at different outlets. m.s.l. in Denglongshan increased at an annual speed of 1.5–2.0 mm around 1990 and 3.0 mm in the early 21st century [14, 25]. However, m.s.l.s at the other four outlets have increased 3–5 mm annually since 1992 [7]. Also, abruption of the tide levels in all outlet stations was also tested (Table 7). Abruption took place in three out of five outlet stations in the late 1980s. Besides, abruption took place in Sanshakou in 1992 and in Wanqinshaxi in 1964. To conclude, abruption of flood tidal level in outlet stations mainly took place from 1985 to 1992.Table 7
Abruption of flood water levels around Pearl River Delta.
Year of abruption
Stations
1964
—
Wanqinshaxi
1986
Hengmen
Nansha
1989
Denglongshan
—
1991
Zhuyin
—
1992
Rongqi
Sanshui SanduoSanshakou
Note: all the abruption events have reached a 95% significant level.In the PRD, large-scale sand-dredging started around 1985 and was not officially banned until 2008. From 1985 to 1990, the speed of sand excavation was 7~10 × 108 m3 per year, while the natural speed of sand siltation in the PRD is only 800~1000 × 104 m3 per year. According to Kong [25], “with the speed of natural deposits, the dredged sands in riverbed cannot be returned within 100 years.” The volume of sand excavating from 1995 to 2008 was one over ten thousand comparing to the 1990s [34], and sand mining in the PRD has been thoroughly banned since 2008.Time series have been divided according to the change in upstream flow, sea-level rise, and local history of sand-dredging. According to the previous discussion, all three driving factors took irreversible alteration during 1985–1995. Thus, the years before 1985 are theprechange period and the years after 1996 theprochange period. As a result, the time series, one piece of data for each month in flood season, that is, three samples per year, has been divided and applied in the quantitative cause analysis.


#### 3.2.3. Calculation ofR

In order to calculate indexR, time series in each gauge station have been first divided into three periods. The first period lasts from 1959 to 1984, and data recorded in this period is the sample forprechange period. The second period lasts from 1985 to 1995 and all the driving factors took abrupt change. The third period lasts from 1996 to 2011, and the samples recorded in this period represent theprochange period.Parametersα, C
0, C
i, and C
j in each model should then be estimated, for both thepre- andprochange period. The parameter set (α
,
C
0
,
C
i
,
C
j) for each model in theprechange period is θ, and (α
′
,
C
0
′
,
C
i
′
,
C
j
′) inprochange period is θ
′. Nonlinear least-squares fit has been applied to estimate the parameters and it has also been carried out in MATLAB. Q
i
- and Z
j
-, mean of Q
i and Z
j in theprechange period, as well as Q
i
′
- and Z
j
′
-, mean in theprochange period, have then been calculated. Z
p
r
e
-, the average on-site flood level in theprechange period, has been estimated with Q
i
- and Z
j
- as the input and θ as the parameter.IndexR has been defined to quantify the impacts of driving factors. According to Huang et al. [14], the quantified impact of sea-level rise can be calculated by (2). To calculate the index R for sea-level rise, we suppose the sea-level rise is the only changing variable in inputs. In other words, sand-dredging does not occur; the topographic feature and the upstream flow of the PRD stay the same as they were during theprechange period. To eliminate the dimensions and differences among stations, we defined index R of sea-level rise by dividing the absolute value of variability in flood level, that is, Δ
Z, by Z
p
r
e
-. R
sea-level  rise can be calculated as(4)
R
sea-level rise
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro  sea-level  rise
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
,where Z
pro sea-level rise
- is the “sea-level rise affected” average on-site flood level during theprochange period, and it has been estimated, with Q
i
- and Z
j
′
- as the input and θ as the parameter.R
sand-dredging and R
upstream have been defined the same way as R
sea-level rise: (5)
R
sand-dredging
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro  sand-dredging
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
,where Z
pro sand-dredging
- is the “sand-dredging affected” average on-site flood level during theprochange period, and it has been estimated, with Q
i
- and Z
j
- as the input and θ
′ as the parameter. Consider(6)
R
upstream
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro  upstream
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
,where Z
pro upstream
- is the “upstream flood variation affected” average on-site flood level during theprochange period, and it has been estimated with Q
i
′
- and Z
j
- as the input and θ as the parameter.Z
pro upstream
-, Z
pro sand-dredging
-, and Z
pro sea-level rise
- together are called Z
pro
-, the average on-site flood level during theprochange period. Therefore, index R for any of the three driving factors can be calculated as(7)
R
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
.There are certain limitations of calculating the quantitative impact of each factor with indexR. Since the numerator in (7) is the absolute value of variability, the value of index R can only work as a comparative percentage, without showing how the impact factors have taken effect on the target.


## 3.1. Multilinear Regression Model

There are mainly four approaches to study the influence of multiple factors, that is, mathematical model, field observation, qualitative theoretical reasoning, and hydrodynamic model [27, 28]. However, qualitative theoretical reasoning cannot provide quantitative results. Field observation costs great, and it is beyond the financial reach of this research. At the same time, large amounts of terrain data from field observation are required to build a precise hydrodynamic model for a delta with complex river network. Therefore, three approaches were denied, and mathematical model was chosen for this research.


### 3.1.1. Modification of Huang’s Model

In 2004, Huang et al. proposed a multilinear regression model to estimate the backwater effect in the PRD under different scenarios of income discharge [14]. The annual extreme water level in any gauge station was supposed to be a dependent variable of the income discharge and the annual maximum tide level at outlets. And the on-site annual maximum water level can be represented as(1)
Z
i
=
f
Q
,
Z
=
C
1
Q
+
C
2
Z
+
α
Q
Z
+
C
0
,where Z
i is the annual extreme water level of a gauge station; Z is the annual maximum tide level at the outlet; Q is the river discharge from WR detected at Makou Station; α and β are the parameters; C
1 and C
2 are the constants. Parameters and constants together represent the regional topographic feature. Given the rise in sea level as Δ
Z, the variability of annual maximum water level in the gauge station can be therefore estimated by the following equation:(2)
Δ
Z
i
=
f
Q
,
Z
′
-
f
Q
,
Z
=
C
1
Q
+
C
2
Z
2
+
Δ
Z
+
α
Q
Z
+
Δ
Z
+
C
0
-
C
1
Q
+
C
2
Z
+
α
Q
Z
+
C
0
=
C
2
+
α
Q
Δ
Z
.Despite its reliable application in the PRD, there are some underlying limitations in Huang’s model. In (1), Z
i is only affected by tide at one outlet. However, flood flow runs towards sea from multiple outlets, and sea-level rise acts differently on tidal levels at difference outlets, even in the same delta. In this case, the impact from sea-level rise should be represented by tidal levels at multiple outlets. Likewise, Huang’s equation takes only the impact from one income flow source into consideration. Somehow, discharge at a gauge station usually comes from more than one upstream river. Impacts from different inflows should therefore be taken into consideration.In this way,C
1
Q
+
C
2
Z has been developed into ∑
i
=
1
n
C
i
Q
j
+
∑
j
=
n
+
1
n
+
m
C
j
Z
j, where n and m are the numbers of inflow resources and outlets and C
i, C
j, the parameters related to the local topographic features. α
Q
Z in (1) has been modified into α
∏
i
=
1
n
Q
i
∏
j
=
n
+
1
n
+
m
Z
j, and α also relates with the local topography. At the same time, the time scale in Huang’s model is year. Somehow, in our research, it has become month, that is, June, July, and August in the flood season. Therefore, the time scale of inputs and outcomes has been narrowed down. As a result, Huang’s model has been modified into the following equation:(3)
Z
=
f
Q
i
,
Z
j
=
∑
i
=
1
n
C
i
Q
i
+
∑
j
=
n
+
1
n
+
m
C
j
Z
j
+
α
∏
i
=
1
n
Q
i
∏
j
=
1
m
Z
j
+
C
0
,where Z is the monthly extreme flood level of the studied gauge station; Q
i is the income flood volume; Z
j is the monthly extreme tide level at estuary of the PRD; α, C
0, C
i, and C
j are the constants representing topographic features.


### 3.1.2. Model Adjustment

Based on (3), regression models for all the stations in the PRD have been built (Table 3). All the the PRD stations have first been divided into three groups according to their locations. Flood from upstream runs through Makou and Sanshui Station into delta, and therefore they are income stations. Flood volumes in income stations have been taken as Q
i. Stations located within 5 km to the river mouth in delta are outlet stations. Tide levels of these stations have been applied in models as Z
j. Other stations in the PRD are network stations. Flood levels in these stations are Z, which is the dependent variable of Q
i, Z
j, and the parameters α, C
0, C
i, and C
j. For all the equations in Table 3, subscriptions of Q and Z are abbreviations of corresponding gauge stations, and the abbreviations for all gauge stations are listed in Table 2.Table 3
Regression model of the PRD stations.
Group
Station
Equation
Income stations
Makou
Z
mk
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
α
Q
mk
Q
ss
Z
dls
+
C
0
Sanshui
Z
ss
=
C
1
Q
s
s
+
C
2
Q
m
k
+
C
1
Z
s
s
k
+
C
2
Z
ns
+
C
3
Z
wqsx
+
α
Q
s
s
Q
m
k
Z
s
s
k
Z
n
s
Z
w
q
s
x
+
C
0
Network stations
Tianhe/Jiangmen/Zhuyin
Z
th/jm/zy
=
C
1
Q
m
k
+
C
2
Q
ss
+
C
3
Z
d
l
s
+
α
Q
m
k
Q
ss
Z
d
l
s
+
C
0
Nanhua/Rongqi
Z
nh/rq
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
wqsx
+
C
4
Z
hm
+
α
Q
mk
Q
ss
Z
wqsx
Z
hm
+
C
0
Sanduo
Z
sd
⁡
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
ssk
+
C
4
Z
ns
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
ssk
Z
ns
Z
wqsx
Z
hm
+
C
0
Outlet stations
Denglongshan
Z
dls
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
ssk
+
C
4
Z
ns
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
ssk
Z
ns
Z
wqsx
Z
hm
+
C
0
Nansha
Z
ns
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ssk
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
dls
Z
ssk
Z
wqsx
Z
hm
+
C
0
Sanshakou
Z
ssk
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ns
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
dls
Z
ns
Z
wqsx
Z
hm
+
C
0
Wanqinshaxi
Z
wqsx
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ns
+
C
5
Z
ssk
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
dls
Z
ns
Z
ssk
Z
hm
+
C
0
Hengmen
Z
hm
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ns
+
C
5
Z
ssk
+
C
6
Z
wqsx
+
α
Q
mk
Q
ss
Z
dls
Z
ns
Z
ssk
Z
wqsx
+
C
0As for the income stations, the model has been adapted. Since income stations are connected by Sixianjiao Channel, both flood flows and flood levels in these two income stations are mutually affected. In addition, flood flow from Makou goes into the South China Sea through Denglongshan and Sanshui through the other four outlet stations. Therefore, flood level in Makou is also affected by the tide level in Denglongshan. Likewise, flood level in Sanshui is related to flood flows from both income stations as well as and the tide levels in Sanshakou, Nansha, Wanqinshaxi, and Hengmen.The model of Zhuyin Station is taken as an example for network stations. Due to the connection between Makou and Sanshui, flood from both income stations, taken asQ
m
k and Q
s
s, would affect the flood level in Zhuyin. And since flood in Zhuyin runs into South China Sea through Denglongshan, tide level in Denglongshan Z
d
l
s takes effect on Z
z
y, that is, the flood level in Zhuyin.The same as the network stations, flood levels in all the outlet stations are affected by the flood flow from both Makou and Sanshui. Since flood levels in all the outlet stations are simultaneously affected by the tidal wave in the South China Sea, flood levels in outlet station in the the PRD are mutually affected. Therefore, to estimateZ for each outlet station, the flows in both Makou and Sanshui are taken as Q
j and the water level in other outlet stations is taken as Z
j.


### 3.1.3. Model Validation

For each station,Q
i in income stations, Z
j in outlet stations, and the on-site extreme level Z of the same month within the flood season in a year are identified as a sample. For each year, there are three samples, that is, the one in June, July, and August. Take Zhuyin, for example; Q
1 and Q
2 are the peak flood flow in Makou and Sanshui in June 1959, Z
1 is the extreme tide level in Denglongshan in June 1959, and Z is the flood level in Zhuyin. Therefore, Q
1, Q
2, Z
1, and Z together are taken as a sample. In the sample, Q
i and Z
j are instances and Z is the target.Cross-validation (CV) makes full use of data information [29], and the validation and accuracy of regression model have been tested by the 10-fold cross-validation method [30]. For each gauge station, we randomly selected 90% of 159 samples to estimate the coefficients and tested the accuracy of the at-site model on the remaining 10% of data. Mean Squared Error (MSE) was applied to evaluate the accuracy of each test. After the procedure was repeated for 10 times on different data partition, a series of 10 MSEs was achieved. The average of this MSE series, that is, CV mean, was used to measure the accuracy of the regression model. The standard deviation, that is, CV variance, represented the robustness of the model. The lower the CV mean is, the more accurate the model is. Likewise, the smaller the CV variance is, the more robust the model is. Since Huang’s model in 2004 was successfully applied and shared, cross-validation results of Huang’s model were taken as the benchmark. Regression model with a lower CV mean and variance than the benchmark was proved eligible. The whole CV test was performed in MATLAB.


## 3.1.1. Modification of Huang’s Model

In 2004, Huang et al. proposed a multilinear regression model to estimate the backwater effect in the PRD under different scenarios of income discharge [14]. The annual extreme water level in any gauge station was supposed to be a dependent variable of the income discharge and the annual maximum tide level at outlets. And the on-site annual maximum water level can be represented as(1)
Z
i
=
f
Q
,
Z
=
C
1
Q
+
C
2
Z
+
α
Q
Z
+
C
0
,where Z
i is the annual extreme water level of a gauge station; Z is the annual maximum tide level at the outlet; Q is the river discharge from WR detected at Makou Station; α and β are the parameters; C
1 and C
2 are the constants. Parameters and constants together represent the regional topographic feature. Given the rise in sea level as Δ
Z, the variability of annual maximum water level in the gauge station can be therefore estimated by the following equation:(2)
Δ
Z
i
=
f
Q
,
Z
′
-
f
Q
,
Z
=
C
1
Q
+
C
2
Z
2
+
Δ
Z
+
α
Q
Z
+
Δ
Z
+
C
0
-
C
1
Q
+
C
2
Z
+
α
Q
Z
+
C
0
=
C
2
+
α
Q
Δ
Z
.Despite its reliable application in the PRD, there are some underlying limitations in Huang’s model. In (1), Z
i is only affected by tide at one outlet. However, flood flow runs towards sea from multiple outlets, and sea-level rise acts differently on tidal levels at difference outlets, even in the same delta. In this case, the impact from sea-level rise should be represented by tidal levels at multiple outlets. Likewise, Huang’s equation takes only the impact from one income flow source into consideration. Somehow, discharge at a gauge station usually comes from more than one upstream river. Impacts from different inflows should therefore be taken into consideration.In this way,C
1
Q
+
C
2
Z has been developed into ∑
i
=
1
n
C
i
Q
j
+
∑
j
=
n
+
1
n
+
m
C
j
Z
j, where n and m are the numbers of inflow resources and outlets and C
i, C
j, the parameters related to the local topographic features. α
Q
Z in (1) has been modified into α
∏
i
=
1
n
Q
i
∏
j
=
n
+
1
n
+
m
Z
j, and α also relates with the local topography. At the same time, the time scale in Huang’s model is year. Somehow, in our research, it has become month, that is, June, July, and August in the flood season. Therefore, the time scale of inputs and outcomes has been narrowed down. As a result, Huang’s model has been modified into the following equation:(3)
Z
=
f
Q
i
,
Z
j
=
∑
i
=
1
n
C
i
Q
i
+
∑
j
=
n
+
1
n
+
m
C
j
Z
j
+
α
∏
i
=
1
n
Q
i
∏
j
=
1
m
Z
j
+
C
0
,where Z is the monthly extreme flood level of the studied gauge station; Q
i is the income flood volume; Z
j is the monthly extreme tide level at estuary of the PRD; α, C
0, C
i, and C
j are the constants representing topographic features.


## 3.1.2. Model Adjustment

Based on (3), regression models for all the stations in the PRD have been built (Table 3). All the the PRD stations have first been divided into three groups according to their locations. Flood from upstream runs through Makou and Sanshui Station into delta, and therefore they are income stations. Flood volumes in income stations have been taken as Q
i. Stations located within 5 km to the river mouth in delta are outlet stations. Tide levels of these stations have been applied in models as Z
j. Other stations in the PRD are network stations. Flood levels in these stations are Z, which is the dependent variable of Q
i, Z
j, and the parameters α, C
0, C
i, and C
j. For all the equations in Table 3, subscriptions of Q and Z are abbreviations of corresponding gauge stations, and the abbreviations for all gauge stations are listed in Table 2.Table 3
Regression model of the PRD stations.
Group
Station
Equation
Income stations
Makou
Z
mk
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
α
Q
mk
Q
ss
Z
dls
+
C
0
Sanshui
Z
ss
=
C
1
Q
s
s
+
C
2
Q
m
k
+
C
1
Z
s
s
k
+
C
2
Z
ns
+
C
3
Z
wqsx
+
α
Q
s
s
Q
m
k
Z
s
s
k
Z
n
s
Z
w
q
s
x
+
C
0
Network stations
Tianhe/Jiangmen/Zhuyin
Z
th/jm/zy
=
C
1
Q
m
k
+
C
2
Q
ss
+
C
3
Z
d
l
s
+
α
Q
m
k
Q
ss
Z
d
l
s
+
C
0
Nanhua/Rongqi
Z
nh/rq
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
wqsx
+
C
4
Z
hm
+
α
Q
mk
Q
ss
Z
wqsx
Z
hm
+
C
0
Sanduo
Z
sd
⁡
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
ssk
+
C
4
Z
ns
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
ssk
Z
ns
Z
wqsx
Z
hm
+
C
0
Outlet stations
Denglongshan
Z
dls
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
ssk
+
C
4
Z
ns
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
ssk
Z
ns
Z
wqsx
Z
hm
+
C
0
Nansha
Z
ns
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ssk
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
dls
Z
ssk
Z
wqsx
Z
hm
+
C
0
Sanshakou
Z
ssk
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ns
+
C
5
Z
wqsx
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
dls
Z
ns
Z
wqsx
Z
hm
+
C
0
Wanqinshaxi
Z
wqsx
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ns
+
C
5
Z
ssk
+
C
6
Z
hm
+
α
Q
mk
Q
ss
Z
dls
Z
ns
Z
ssk
Z
hm
+
C
0
Hengmen
Z
hm
=
C
1
Q
mk
+
C
2
Q
ss
+
C
3
Z
dls
+
C
4
Z
ns
+
C
5
Z
ssk
+
C
6
Z
wqsx
+
α
Q
mk
Q
ss
Z
dls
Z
ns
Z
ssk
Z
wqsx
+
C
0As for the income stations, the model has been adapted. Since income stations are connected by Sixianjiao Channel, both flood flows and flood levels in these two income stations are mutually affected. In addition, flood flow from Makou goes into the South China Sea through Denglongshan and Sanshui through the other four outlet stations. Therefore, flood level in Makou is also affected by the tide level in Denglongshan. Likewise, flood level in Sanshui is related to flood flows from both income stations as well as and the tide levels in Sanshakou, Nansha, Wanqinshaxi, and Hengmen.The model of Zhuyin Station is taken as an example for network stations. Due to the connection between Makou and Sanshui, flood from both income stations, taken asQ
m
k and Q
s
s, would affect the flood level in Zhuyin. And since flood in Zhuyin runs into South China Sea through Denglongshan, tide level in Denglongshan Z
d
l
s takes effect on Z
z
y, that is, the flood level in Zhuyin.The same as the network stations, flood levels in all the outlet stations are affected by the flood flow from both Makou and Sanshui. Since flood levels in all the outlet stations are simultaneously affected by the tidal wave in the South China Sea, flood levels in outlet station in the the PRD are mutually affected. Therefore, to estimateZ for each outlet station, the flows in both Makou and Sanshui are taken as Q
j and the water level in other outlet stations is taken as Z
j.


## 3.1.3. Model Validation

For each station,Q
i in income stations, Z
j in outlet stations, and the on-site extreme level Z of the same month within the flood season in a year are identified as a sample. For each year, there are three samples, that is, the one in June, July, and August. Take Zhuyin, for example; Q
1 and Q
2 are the peak flood flow in Makou and Sanshui in June 1959, Z
1 is the extreme tide level in Denglongshan in June 1959, and Z is the flood level in Zhuyin. Therefore, Q
1, Q
2, Z
1, and Z together are taken as a sample. In the sample, Q
i and Z
j are instances and Z is the target.Cross-validation (CV) makes full use of data information [29], and the validation and accuracy of regression model have been tested by the 10-fold cross-validation method [30]. For each gauge station, we randomly selected 90% of 159 samples to estimate the coefficients and tested the accuracy of the at-site model on the remaining 10% of data. Mean Squared Error (MSE) was applied to evaluate the accuracy of each test. After the procedure was repeated for 10 times on different data partition, a series of 10 MSEs was achieved. The average of this MSE series, that is, CV mean, was used to measure the accuracy of the regression model. The standard deviation, that is, CV variance, represented the robustness of the model. The lower the CV mean is, the more accurate the model is. Likewise, the smaller the CV variance is, the more robust the model is. Since Huang’s model in 2004 was successfully applied and shared, cross-validation results of Huang’s model were taken as the benchmark. Regression model with a lower CV mean and variance than the benchmark was proved eligible. The whole CV test was performed in MATLAB.


## 3.2. Changing RatioR

### 3.2.1. Jiang’s Idea

IndexR defined by Jiang et al. [31] was used to quantify the impacts of different driving factors on monthly flow of a small river basin in China. The local hydrological cycle was severely disturbed by human activities and the time series was no longer consistent. They split the time series into two according to the timing of severe human interference, that is, one before the severe human interventions and one afterwards. The previous period was anatural period, and the latter was ahuman-affected period. First they calculated the model constants with data fromnatural period, and then they estimated monthly flow in thehuman-affected period with these constants. The estimated monthly flow represented the monthly flow in thehuman-affected period if none of the severe human activities happened. The gap between the estimated and recorded flow in thehuman-affected period was the impact of human activities on monthly flow.


### 3.2.2. Division of Time Series

Based on Jiang’s theory [31], index R was defined to show the impacts of three driving factors. The difference between the factor-affected flood level and the untouched one represents the impact of a certain factor. To estimate the untouched series, hydrological series need to be divided first. Data length should be unified so that results in different stations can be compared. Since the available data varied among gauge stations (Table 2), we have only applied data from 1959 to 2011. In order to lower the uncertainties in analysis, the time series have only been divided into two periods. To find the dividing time for all three driving factors, we have studied the temporal pattern of these factors.Changes in upstream flow have been detected. Wuzhou, located downstream of West River, is more than 300 km away from Makou. Shijiao, located downstream of North River, is 50 km away from Sanshui. Mann-Kendall test was applied in detecting changes in recorded flow from 1959 to 2011 of these four stations (Table4). Flood flow increased in Wuzhou and decreased in Shijiao, yet neither of them reached the significance of 90%. Flood flow in both Makou and Sanshui has increased at a 95% significant level. It can be concluded that flood flows in upstream river basins do not have significant trend. As for the abruption in series, flood flow in Wuzhou jumps in 1991, and that in Sanshui shows no significant abruption. Flood flow in Makou stays consistent, while that in Sanshui took a jump in 1992. Therefore, there was abruption in West River flow in the early 1990s, probably due to dam construction upstream [32].Table 4
Mann-Kendall test results of income flood flow.
Station
West River
North River
Wuzhou
Makou
Shijiao
Sanshui
Mann-Kendall statistics
1.55
1.97
−1.12
2.03
Year of abruption
1991
—
—
1992
Note: Mann-Kendall statistics is known asZ. Positive Z indicates an increasing trend just as a negative one refers to a decreasing trend. Moreover, Z
>
1.65, Z
>
1.96, and Z
>
2.57 show tendencies at the 0.01, 0.05, and 0.1 significance levels.The inconsistent change of flow between Wuzhou and Makou is partly due to dyke construction and variations in precipitation. In the early 1990s, dykes were built in large scale around Wuzhou [32], and the flood flow from Wuzhou to Makou is locked within the channel. Furthermore, precipitation in flood season increased faster on the catchment area of Makou than Wuzhou within the last 50 years [33]. Therefore, flood flow in Makou has increased more severely than that in Wuzhou.The discrepancies between Wuzhou and Makou, Shijiao and Sanshui are mainly due to severe sand-dredging in the PRD. Table5 shows the variation of channel geometries in both Sanshui and Makou. Apparently, the average channel depth in Sanshui increased 0.5 m more than that in Makou from 1985 to 2006. Channel volume in Sanshui therefore increased by 69% during 1985–1999 and 36% during 1999–2006. Meanwhile, the percentage in Makou is only 18% and 29%. From 1980 to 2003, discharge area over both 2 m and 6 m in Sanshui increased three times that in Makou. Due to different extent of riverbed alteration, flood flows from Makou to Sanshui through Sixianjiao Channel. As a result, flow partition ratio between 1959 and 2011 in flood seasons increased in Sanshui (Figure 2). The increase became sharp from the early 1990s. Since mean flow in WR triples that of NR (Table 6), increase of flow partition ratio for more than 5% in Sanshui has significantly increased the on-site flow.Table 5
Variation of channel geometries in Sanshui and Makou.
Geometry index
Unit
Period
Makou
Sanshui
Maximum cutting depth
m
1980–2003
+7.11
+8.07
Average depth
1985–19991999–2006
0.802.00
2.801.50
Channel volume
1985–1999
+18.00
+69.00
%
1999–2006
+29.00
+36.00
Change of discharge area
>2 m
1980–2003
+19.30
+99.50
>6 m
1980–2003
+15.20
+40.00
Note: “+” means increase.Table 6
Average monthly flow in West River and North River in flood season (units: m3).
River
Station
Jun.
Jul.
Aug.
West River
Wuzhou
25,222
24,971
20,043
Makou
27,658
26,780
20,564
North River
Shijiao
7,406
4,844
3,501
Sanshui
7,229
6,917
4,893Figure 2
Tendency of flow partition ration of Sanshui.MK test was also carried out on the average tide level (m.s.l.) of the outlet stations. Results of MK in Sanshakou are 2.15, 1.75 for Nansha, 1.83 for Wanqinshaxi, 3.10 for Hengmen, and 0.08 for Denglongshan. All results are positive, which indicate increasing tendency in all outlets. According to the previous study, sea level at South China Sea has kept rising since the early 1960s [7, 14, 25]. Among all the outlet stations, increase in Denglongshan has not reached a significance of 1%. Due to the wandering rim of continental shelf, sea-level rise in South China Sea has caused different extent of increase at different outlets. m.s.l. in Denglongshan increased at an annual speed of 1.5–2.0 mm around 1990 and 3.0 mm in the early 21st century [14, 25]. However, m.s.l.s at the other four outlets have increased 3–5 mm annually since 1992 [7]. Also, abruption of the tide levels in all outlet stations was also tested (Table 7). Abruption took place in three out of five outlet stations in the late 1980s. Besides, abruption took place in Sanshakou in 1992 and in Wanqinshaxi in 1964. To conclude, abruption of flood tidal level in outlet stations mainly took place from 1985 to 1992.Table 7
Abruption of flood water levels around Pearl River Delta.
Year of abruption
Stations
1964
—
Wanqinshaxi
1986
Hengmen
Nansha
1989
Denglongshan
—
1991
Zhuyin
—
1992
Rongqi
Sanshui SanduoSanshakou
Note: all the abruption events have reached a 95% significant level.In the PRD, large-scale sand-dredging started around 1985 and was not officially banned until 2008. From 1985 to 1990, the speed of sand excavation was 7~10 × 108 m3 per year, while the natural speed of sand siltation in the PRD is only 800~1000 × 104 m3 per year. According to Kong [25], “with the speed of natural deposits, the dredged sands in riverbed cannot be returned within 100 years.” The volume of sand excavating from 1995 to 2008 was one over ten thousand comparing to the 1990s [34], and sand mining in the PRD has been thoroughly banned since 2008.Time series have been divided according to the change in upstream flow, sea-level rise, and local history of sand-dredging. According to the previous discussion, all three driving factors took irreversible alteration during 1985–1995. Thus, the years before 1985 are theprechange period and the years after 1996 theprochange period. As a result, the time series, one piece of data for each month in flood season, that is, three samples per year, has been divided and applied in the quantitative cause analysis.


### 3.2.3. Calculation ofR

In order to calculate indexR, time series in each gauge station have been first divided into three periods. The first period lasts from 1959 to 1984, and data recorded in this period is the sample forprechange period. The second period lasts from 1985 to 1995 and all the driving factors took abrupt change. The third period lasts from 1996 to 2011, and the samples recorded in this period represent theprochange period.Parametersα, C
0, C
i, and C
j in each model should then be estimated, for both thepre- andprochange period. The parameter set (α
,
C
0
,
C
i
,
C
j) for each model in theprechange period is θ, and (α
′
,
C
0
′
,
C
i
′
,
C
j
′) inprochange period is θ
′. Nonlinear least-squares fit has been applied to estimate the parameters and it has also been carried out in MATLAB. Q
i
- and Z
j
-, mean of Q
i and Z
j in theprechange period, as well as Q
i
′
- and Z
j
′
-, mean in theprochange period, have then been calculated. Z
p
r
e
-, the average on-site flood level in theprechange period, has been estimated with Q
i
- and Z
j
- as the input and θ as the parameter.IndexR has been defined to quantify the impacts of driving factors. According to Huang et al. [14], the quantified impact of sea-level rise can be calculated by (2). To calculate the index R for sea-level rise, we suppose the sea-level rise is the only changing variable in inputs. In other words, sand-dredging does not occur; the topographic feature and the upstream flow of the PRD stay the same as they were during theprechange period. To eliminate the dimensions and differences among stations, we defined index R of sea-level rise by dividing the absolute value of variability in flood level, that is, Δ
Z, by Z
p
r
e
-. R
sea-level  rise can be calculated as(4)
R
sea-level rise
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro  sea-level  rise
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
,where Z
pro sea-level rise
- is the “sea-level rise affected” average on-site flood level during theprochange period, and it has been estimated, with Q
i
- and Z
j
′
- as the input and θ as the parameter.R
sand-dredging and R
upstream have been defined the same way as R
sea-level rise: (5)
R
sand-dredging
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro  sand-dredging
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
,where Z
pro sand-dredging
- is the “sand-dredging affected” average on-site flood level during theprochange period, and it has been estimated, with Q
i
- and Z
j
- as the input and θ
′ as the parameter. Consider(6)
R
upstream
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro  upstream
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
,where Z
pro upstream
- is the “upstream flood variation affected” average on-site flood level during theprochange period, and it has been estimated with Q
i
′
- and Z
j
- as the input and θ as the parameter.Z
pro upstream
-, Z
pro sand-dredging
-, and Z
pro sea-level rise
- together are called Z
pro
-, the average on-site flood level during theprochange period. Therefore, index R for any of the three driving factors can be calculated as(7)
R
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
.There are certain limitations of calculating the quantitative impact of each factor with indexR. Since the numerator in (7) is the absolute value of variability, the value of index R can only work as a comparative percentage, without showing how the impact factors have taken effect on the target.


## 3.2.1. Jiang’s Idea

IndexR defined by Jiang et al. [31] was used to quantify the impacts of different driving factors on monthly flow of a small river basin in China. The local hydrological cycle was severely disturbed by human activities and the time series was no longer consistent. They split the time series into two according to the timing of severe human interference, that is, one before the severe human interventions and one afterwards. The previous period was anatural period, and the latter was ahuman-affected period. First they calculated the model constants with data fromnatural period, and then they estimated monthly flow in thehuman-affected period with these constants. The estimated monthly flow represented the monthly flow in thehuman-affected period if none of the severe human activities happened. The gap between the estimated and recorded flow in thehuman-affected period was the impact of human activities on monthly flow.


## 3.2.2. Division of Time Series

Based on Jiang’s theory [31], index R was defined to show the impacts of three driving factors. The difference between the factor-affected flood level and the untouched one represents the impact of a certain factor. To estimate the untouched series, hydrological series need to be divided first. Data length should be unified so that results in different stations can be compared. Since the available data varied among gauge stations (Table 2), we have only applied data from 1959 to 2011. In order to lower the uncertainties in analysis, the time series have only been divided into two periods. To find the dividing time for all three driving factors, we have studied the temporal pattern of these factors.Changes in upstream flow have been detected. Wuzhou, located downstream of West River, is more than 300 km away from Makou. Shijiao, located downstream of North River, is 50 km away from Sanshui. Mann-Kendall test was applied in detecting changes in recorded flow from 1959 to 2011 of these four stations (Table4). Flood flow increased in Wuzhou and decreased in Shijiao, yet neither of them reached the significance of 90%. Flood flow in both Makou and Sanshui has increased at a 95% significant level. It can be concluded that flood flows in upstream river basins do not have significant trend. As for the abruption in series, flood flow in Wuzhou jumps in 1991, and that in Sanshui shows no significant abruption. Flood flow in Makou stays consistent, while that in Sanshui took a jump in 1992. Therefore, there was abruption in West River flow in the early 1990s, probably due to dam construction upstream [32].Table 4
Mann-Kendall test results of income flood flow.
Station
West River
North River
Wuzhou
Makou
Shijiao
Sanshui
Mann-Kendall statistics
1.55
1.97
−1.12
2.03
Year of abruption
1991
—
—
1992
Note: Mann-Kendall statistics is known asZ. Positive Z indicates an increasing trend just as a negative one refers to a decreasing trend. Moreover, Z
>
1.65, Z
>
1.96, and Z
>
2.57 show tendencies at the 0.01, 0.05, and 0.1 significance levels.The inconsistent change of flow between Wuzhou and Makou is partly due to dyke construction and variations in precipitation. In the early 1990s, dykes were built in large scale around Wuzhou [32], and the flood flow from Wuzhou to Makou is locked within the channel. Furthermore, precipitation in flood season increased faster on the catchment area of Makou than Wuzhou within the last 50 years [33]. Therefore, flood flow in Makou has increased more severely than that in Wuzhou.The discrepancies between Wuzhou and Makou, Shijiao and Sanshui are mainly due to severe sand-dredging in the PRD. Table5 shows the variation of channel geometries in both Sanshui and Makou. Apparently, the average channel depth in Sanshui increased 0.5 m more than that in Makou from 1985 to 2006. Channel volume in Sanshui therefore increased by 69% during 1985–1999 and 36% during 1999–2006. Meanwhile, the percentage in Makou is only 18% and 29%. From 1980 to 2003, discharge area over both 2 m and 6 m in Sanshui increased three times that in Makou. Due to different extent of riverbed alteration, flood flows from Makou to Sanshui through Sixianjiao Channel. As a result, flow partition ratio between 1959 and 2011 in flood seasons increased in Sanshui (Figure 2). The increase became sharp from the early 1990s. Since mean flow in WR triples that of NR (Table 6), increase of flow partition ratio for more than 5% in Sanshui has significantly increased the on-site flow.Table 5
Variation of channel geometries in Sanshui and Makou.
Geometry index
Unit
Period
Makou
Sanshui
Maximum cutting depth
m
1980–2003
+7.11
+8.07
Average depth
1985–19991999–2006
0.802.00
2.801.50
Channel volume
1985–1999
+18.00
+69.00
%
1999–2006
+29.00
+36.00
Change of discharge area
>2 m
1980–2003
+19.30
+99.50
>6 m
1980–2003
+15.20
+40.00
Note: “+” means increase.Table 6
Average monthly flow in West River and North River in flood season (units: m3).
River
Station
Jun.
Jul.
Aug.
West River
Wuzhou
25,222
24,971
20,043
Makou
27,658
26,780
20,564
North River
Shijiao
7,406
4,844
3,501
Sanshui
7,229
6,917
4,893Figure 2
Tendency of flow partition ration of Sanshui.MK test was also carried out on the average tide level (m.s.l.) of the outlet stations. Results of MK in Sanshakou are 2.15, 1.75 for Nansha, 1.83 for Wanqinshaxi, 3.10 for Hengmen, and 0.08 for Denglongshan. All results are positive, which indicate increasing tendency in all outlets. According to the previous study, sea level at South China Sea has kept rising since the early 1960s [7, 14, 25]. Among all the outlet stations, increase in Denglongshan has not reached a significance of 1%. Due to the wandering rim of continental shelf, sea-level rise in South China Sea has caused different extent of increase at different outlets. m.s.l. in Denglongshan increased at an annual speed of 1.5–2.0 mm around 1990 and 3.0 mm in the early 21st century [14, 25]. However, m.s.l.s at the other four outlets have increased 3–5 mm annually since 1992 [7]. Also, abruption of the tide levels in all outlet stations was also tested (Table 7). Abruption took place in three out of five outlet stations in the late 1980s. Besides, abruption took place in Sanshakou in 1992 and in Wanqinshaxi in 1964. To conclude, abruption of flood tidal level in outlet stations mainly took place from 1985 to 1992.Table 7
Abruption of flood water levels around Pearl River Delta.
Year of abruption
Stations
1964
—
Wanqinshaxi
1986
Hengmen
Nansha
1989
Denglongshan
—
1991
Zhuyin
—
1992
Rongqi
Sanshui SanduoSanshakou
Note: all the abruption events have reached a 95% significant level.In the PRD, large-scale sand-dredging started around 1985 and was not officially banned until 2008. From 1985 to 1990, the speed of sand excavation was 7~10 × 108 m3 per year, while the natural speed of sand siltation in the PRD is only 800~1000 × 104 m3 per year. According to Kong [25], “with the speed of natural deposits, the dredged sands in riverbed cannot be returned within 100 years.” The volume of sand excavating from 1995 to 2008 was one over ten thousand comparing to the 1990s [34], and sand mining in the PRD has been thoroughly banned since 2008.Time series have been divided according to the change in upstream flow, sea-level rise, and local history of sand-dredging. According to the previous discussion, all three driving factors took irreversible alteration during 1985–1995. Thus, the years before 1985 are theprechange period and the years after 1996 theprochange period. As a result, the time series, one piece of data for each month in flood season, that is, three samples per year, has been divided and applied in the quantitative cause analysis.


## 3.2.3. Calculation ofR

In order to calculate indexR, time series in each gauge station have been first divided into three periods. The first period lasts from 1959 to 1984, and data recorded in this period is the sample forprechange period. The second period lasts from 1985 to 1995 and all the driving factors took abrupt change. The third period lasts from 1996 to 2011, and the samples recorded in this period represent theprochange period.Parametersα, C
0, C
i, and C
j in each model should then be estimated, for both thepre- andprochange period. The parameter set (α
,
C
0
,
C
i
,
C
j) for each model in theprechange period is θ, and (α
′
,
C
0
′
,
C
i
′
,
C
j
′) inprochange period is θ
′. Nonlinear least-squares fit has been applied to estimate the parameters and it has also been carried out in MATLAB. Q
i
- and Z
j
-, mean of Q
i and Z
j in theprechange period, as well as Q
i
′
- and Z
j
′
-, mean in theprochange period, have then been calculated. Z
p
r
e
-, the average on-site flood level in theprechange period, has been estimated with Q
i
- and Z
j
- as the input and θ as the parameter.IndexR has been defined to quantify the impacts of driving factors. According to Huang et al. [14], the quantified impact of sea-level rise can be calculated by (2). To calculate the index R for sea-level rise, we suppose the sea-level rise is the only changing variable in inputs. In other words, sand-dredging does not occur; the topographic feature and the upstream flow of the PRD stay the same as they were during theprechange period. To eliminate the dimensions and differences among stations, we defined index R of sea-level rise by dividing the absolute value of variability in flood level, that is, Δ
Z, by Z
p
r
e
-. R
sea-level  rise can be calculated as(4)
R
sea-level rise
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro  sea-level  rise
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
,where Z
pro sea-level rise
- is the “sea-level rise affected” average on-site flood level during theprochange period, and it has been estimated, with Q
i
- and Z
j
′
- as the input and θ as the parameter.R
sand-dredging and R
upstream have been defined the same way as R
sea-level rise: (5)
R
sand-dredging
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro  sand-dredging
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
,where Z
pro sand-dredging
- is the “sand-dredging affected” average on-site flood level during theprochange period, and it has been estimated, with Q
i
- and Z
j
- as the input and θ
′ as the parameter. Consider(6)
R
upstream
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro  upstream
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
,where Z
pro upstream
- is the “upstream flood variation affected” average on-site flood level during theprochange period, and it has been estimated with Q
i
′
- and Z
j
- as the input and θ as the parameter.Z
pro upstream
-, Z
pro sand-dredging
-, and Z
pro sea-level rise
- together are called Z
pro
-, the average on-site flood level during theprochange period. Therefore, index R for any of the three driving factors can be calculated as(7)
R
=
Δ
Z
Z
p
r
e
-
×
100
%
=
Z
pro
-
-
Z
p
r
e
-
Z
p
r
e
-
×
100
%
.There are certain limitations of calculating the quantitative impact of each factor with indexR. Since the numerator in (7) is the absolute value of variability, the value of index R can only work as a comparative percentage, without showing how the impact factors have taken effect on the target.


## 4. Results and Discussions

### 4.1. Quantitative Impacts of Driving Factors on Flood Level in the PRD

The proposed hydrological model for each station in the PRD has first been proved valid through the cross-validation test. In the CV test, the modified models for all the stations had lower CV mean and variance than Huang’s model. Based on the model and index defined in Section3, we have achieved the quantified impacts of three driving factors on flood level in the PRD. Nanhua is taken as the example to explain the results of R. In Nanhua, R
upstream is 17.26%, R
sand-dredging is 15.35%, and R
sea-level rise is 2.61%. Variation of flood volume in both income stations has increased fromprechange toprochange period, and it has caused the flood level varied by 17.26% in Nanhua (Table 8). Due to the large-scale sand-dredging, the maximum cut depth of the upstream and downstream channel of Nanhua has been cut down by at least 0.5 m fromprechange toprochange period (Figure 3). It has altered flood level in Nanhua by 15.35%. And the sea-level rise, taking effect through tidal waves in Denglongshan and Hengmen, has changed flood level in Nanhua by 2.61%. Owing to the definition of R as well as the complex working mechanism in the PRD, whether these changes caused by driving factor has actually pushed up or lowered down the flood level is unknown. It can be told from the results that the change in income flood takes the most effect and then the sand-dredging and sea-level rise take the last.Table 8
Changing ratioR caused by three driving factors (unit: %).
Station
R
upstream
R
sand-dredging
R
sea-level rise
Makou
40.68
30.48
0.19
Sanshui
53.49
38.73
0.12
Tianhe
18.50
12.74
2.47
Nanhua
17.26
15.35
2.61
Jiangmen
15.11
10.13
5.79
Zhuyin
12.57
9.91
10.44
Sanduo
22.46
23.71
1.47
Rongqi
14.52
13.01
7.71
Denglongshan
5.15
8.73
12.46
Hengmen
5.22
3.35
16.81
Sanshakou
4.29
5.12
13.03
Nansha
5.18
3.37
8.91
Wanqinshaxi
5.14
4.11
14.43Figure 3
Change of maximum cut depth in Pearl River Delta from 1985 to 2006.We have calculatedR
upstream, R
sand-dredging, and R
sea-level rise for all the 13 stations in the PRD (Table 8). Effects due to upstream flood variation are negatively correlated with the distance between the location of the station and the outlets (Table 9). Changes of upstream flood flows have taken effects on the flood levels in income stations by Rs over 40%. It has changed the flood levels in network stations by 10% and 25% and the outlet stations between 4% and 5.5%. R
upstream peaks around income station Makou by 53.49% and comes to the minimum around outlet station Sanshakou by 4.29%. Located downstream of Sanshui, flood level in Sanduo has been affected the most by the income flood variation among all the network stations. Meanwhile, Tianhe and Nanhua have been similarly affected. R
upstream in Tianhe is higher than Jiangmen, and that in Jiangmen is higher than Zhuyin. It is not hard to discover that the further a network station is away from the income stations, the less the flood flow change will affect the on-site flood level. In general, variation in flood flow affects the most flood levels of income stations, then the network stations, and then the outlet stations. In other words, R
upstream decreases from income stations to outlet stations along the flow direction.Table 9
Distance among stations in West River Delta (units: km).
A
B
Distance between A and B
Makou
Denglongshan
123.8
Tianhe
67.8
Jiangmen
54.8
Zhuyin
19.8Riverbed alteration has taken effects on the flood levels by variations over 30% in both income stations, between 9 and 25% in network stations, and lower than 9% in outlet stations. The impact level of sand-dredging,R
sand-dredging, also decreases from income stations to outlet stations. R
sand-dredging in Sanshui is higher than that in Makou. It is consistent with the fact that riverbed alteration in Sanshui, from 1984 to 1995, was more severe than that in Makou (Figure 3 and Table 5). For R
sand-dredging of the network stations, Sanduo is the highest. R
sand-dredging in Tianhe, Nanhua, Jiangmen, and Rongqi are above 10%. R
sand-dredging in Zhuyin are the lowest, slightly below 10%. According to the riverbed alteration between 1985 and 2006 (Figure 3), riverbed in Sanduo has been cut down the most, then Nanhua, and then the other network stations. Despite the slight differences between R
sand-dredging in Tianhe, Jiangmen, and Zhuyin, R
sand-dredging of network stations are consistent with the extent of riverbed alteration. Among all the outlet stations, since the maximum cut depth varies the most in Denglongshan, flood level changes the most in Denglongshan under the impact of riverbed alteration.Sea-level rise has influenced the flood level in the PRD stations by a changing ration between 0.1% and 20%. It has hardly affected the flood level in income stations sinceR
sea-level rise are only 0.12% and 0.19%. Among network stations, R
sea-level rise tops in Zhuyin by 10.44%. Meanwhile, R
sea-level rise in other network stations are all below 6%. R
sea-level rise for Tianhe and Nanhua are close. In addition, R
sea-level rise in Tianhe is lower than that in Jiangmen, and that in Jiangmen is lower than Zhuyin at the same time. R
sea-level rise in outlet stations are all above 8%, much higher than that in income stations or most network stations. Similar to R
upstream, the value of R
sea-level rise seems to be correlated with the distance between stations and delta income, only in an opposite way. In another way, R
sea-level rise decreases from outlet stations to income stations in the the PRD.To compareR
upstream, R
sand-dredging, and R
sea-level rise of the income stations, it is easy to find that R
upstream is the highest and R
sea-level rise the lowest. It can be deducted that flood levels in income stations are influenced the most by the flood flow variation and the least by the sea-level rise. It is almost true for all the network stations as well, except for stations Sanduo and Zhuyin. In Sanduo, sand-dredging has caused the most change in on-site flood level. This is probably because that riverbed has been the most severely cut around Sanduo (Figure 3). In Zhuyin, somehow, the sea-level rise has taken the most effect. The cause will be talked about in Section 4.3. As for all the outlet stations, R
upstream and R
sand-dredging of outlet stations are close, and the difference between these two is all below 3.6%. At the same time, impact from sea-level rise ranks the first among the three driving factors.


### 4.2. Reliability of Quantitative Analysis Results

The efficiency of the quantitative cause analysis results can be generally verified by its consistency with the changing environment. To better tell the spatial patterns of impacts from different driving factors, the results have been interpolated on the whole delta with Kriging method and the results are shown in Figure4. Contour lines in Figure 4(a) show the spatial distribution of flood level change caused by variation of income flood flow. Figure 4(b) shows the spatial distribution of sand-dredging impact, and Figure 4(c) shows the spatial pattern of R
sea-level rise. As discussed in Section 3.2.1, variation of flood flow in Sanshui is greater than that in Makou. R
upstream in Sanshui is greater in these two stations as well. Meanwhile, R
upstream decreases when the station stands further from the upstream. R
sea-level rise also decreases when the station is further from the outlet. The impacts of the varied flood flow and sea-level rise are also changing mainly along the “northwest-southeast” direction, because river channels in the PRD are most developed along the “northwest-southeast” direction, and river discharge runs easily towards estuary by gravity force [23]. At the same time, R
sand-dredging of all network stations are in accordance with the extent of sand-dredging between the 1980s and 1990s (Figure 3). It is higher around Sanduo among network stations and lower around outlet stations. As is shown in Figure 3, channel varies in the greatest extent around Sanduo. Dredging of the riverbeds around Rongqi, Hengmen, Wanqinshaxi, Nansha, and Sanshakou is lower. Besides, by comparing the impacts of three driving factors, we found that, despite Sanduo, flood levels of all the income and network stations were most affected by the upstream flood flow. Previous works have examined the impact upstream flood flow had on network stations like Tianhe, and they argued that upstream flood flow is the most important factor for the PRD flood, despite the severe dredging of the riverbed and obvious sea-level rise [14, 23, 25, 34]. Our result here has verified this argument in network stations as well as the income stations.Figure 4
Spatial pattern ofR.
(a)
R
upstream
(b)
R
sand-dredging
(c)
R
sea-level riseCorrelation coefficients (CCs) among stations have also been calculated to double-check the reliability of our quantitative analysis result. CCs among stations that are not physically connected were not considered. For instance, there are no river channels connecting Tianhe and Rongqi. Thus, flood flow or flood level in Tianhe has no influence on the flood flow or level in Rongqi. These two stations are defined to be not physically connected. Except for the ones from 1985 to 1995, flood level series from 1959 to 2011 of all gauge stations have been extracted to calculate CC. Applied flood level series of each station consists of extreme water level in flood season for 42 years. CCs among stations are thus calculated by the following equation:(8)
C
o
v
X
,
Y
=
E
X
Y
-
E
X
E
Y
,where X is the time series of one station and Y the other; E, the expectation of a time series; Cov, the correlation coefficient of two time series.Higher value of CC indicates better correlation of hydrological time series in two stations. As is shown in Table10, flood level in Tianhe is related to Nanhua the most, and it is also true for Nanhua the other way around. Flood levels in these two stations take almost the same impact from flood variation, sand-dredging, and sea-level rise (Table 8). CCs among stations are positively connected with the distance among stations (Tables 9 and 10). It is the same as the spatial pattern of R
upstream and R
sea-level rise. Furthermore, flood levels in all the network stations are both connected with Sanshui and Makou, but CCs of stations with Sanshui are higher. It is true in Figure 4(a) that R
upstream is higher in Sanshui than Makou. What is more, as for the connection between network stations and outlet stations, CC between Zhuyin and Denglongshan is 2 times higher than any CC between stations downstream of Makou and other outlet stations. Last but not the least, Rongqi is more related with Wanqinshaxi than Nansha (Table 10) and it is the same in Figure 4(c) that R
sea-level rise of Rongqi and Wanqinshaxi is closer.Table 10
Correlation coefficient (CC) among physically relevant stations.
Tianhe
Jiangmen
Zhuyin
Nanhua
Rongqi
Sanduo
Makou
0.834
0.801
0.548
0.834
0.739
0.804
Sanshui
0.966
0.941
0.667
0.964
0.878
0.938
Nanhua
0.9997
0.983
0.681
—
0.939
—
Rongqi
—
—
—
0.939
—
—
Tianhe
—
0.983
0.684
0.9997
—
—
Jiangmen
0.983
—
0.696
0.983
—
—
Zhuyin
0.684
0.696
—
0.681
—
—
Sanshakou
—
—
—
—
—
0.383
Nansha
—
—
—
—
—
0.423
Wanqinshaxi
—
—
—
0.427
0.461
0.397
Hengmen
—
—
—
0.223
0.360
—
Denglongshan
0.241
0.265
0.723
—
—
—


### 4.3. Impact of Network Density

Obviously,R
upstream and R
sea-level rise change from income stations to outlet stations along the flow direction. However, the changing rate slows down around Jiangmen and Zhuyin. There is 35 km between Jiangmen and Zhuyin, almost 3 times the distance between Jiangmen and Tianhe (Table 9). As discussed above, the impacts on flood level by inflow flood change and rising sea level are positively correlated with spatial distances among stations. In this way, the difference of R
upstream or R
sea-level rise between Jiangmen and Zhuyin is supposed to be almost 3 times that between Jiangmen and Tianhe. The difference of R
upstream between Jiangmen and Tianhe is 3.39%. Somehow, the one between Jiangmen and Zhuyin is 2.54% rather than 10.17%. Meanwhile, the difference of R
sea-level rise between Jiangmen and Tianhe is 2.61%, and the one between Jiangmen and Zhuyin is less than that by being only 0.22%.Density of river network is inferred to take up the major responsibility. As is shown in Figure1 and mentioned in Section 2.1, river networks around Jiangmen are significantly denser than other places in the the PRD. There are more river nodes within the same area. At a node, flood flows towards diverse directions, or flood flow from multiple branches joins in the same channel. Due to the terrain features in the PRD, river channels have developed the most along the “northwest-southeast” direction [23], and the velocity of flood flow slows down in river channels that have developed otherwise. Flood from upstream is thus held up among dense river networks. As a result, variation of income flood takes more impact on the dense river network.Meanwhile, sea-level rise takes effect on the flood level through the fluctuation of water body. Its effect tracks back up the river way and jerks up the in-channel flood within dense river network when the flow velocity drops. In this way, the impact of sea-level rise on flood level has also been intensified in dense river network. To conclude, dense river network in delta area intensifies the impact on flood level from both the income flood change and sea-level rise, which explains the irregular pattern ofR
upstream and R
sea-level rise around Jiangmen and Zhuyin.


## 4.1. Quantitative Impacts of Driving Factors on Flood Level in the PRD

The proposed hydrological model for each station in the PRD has first been proved valid through the cross-validation test. In the CV test, the modified models for all the stations had lower CV mean and variance than Huang’s model. Based on the model and index defined in Section3, we have achieved the quantified impacts of three driving factors on flood level in the PRD. Nanhua is taken as the example to explain the results of R. In Nanhua, R
upstream is 17.26%, R
sand-dredging is 15.35%, and R
sea-level rise is 2.61%. Variation of flood volume in both income stations has increased fromprechange toprochange period, and it has caused the flood level varied by 17.26% in Nanhua (Table 8). Due to the large-scale sand-dredging, the maximum cut depth of the upstream and downstream channel of Nanhua has been cut down by at least 0.5 m fromprechange toprochange period (Figure 3). It has altered flood level in Nanhua by 15.35%. And the sea-level rise, taking effect through tidal waves in Denglongshan and Hengmen, has changed flood level in Nanhua by 2.61%. Owing to the definition of R as well as the complex working mechanism in the PRD, whether these changes caused by driving factor has actually pushed up or lowered down the flood level is unknown. It can be told from the results that the change in income flood takes the most effect and then the sand-dredging and sea-level rise take the last.Table 8
Changing ratioR caused by three driving factors (unit: %).
Station
R
upstream
R
sand-dredging
R
sea-level rise
Makou
40.68
30.48
0.19
Sanshui
53.49
38.73
0.12
Tianhe
18.50
12.74
2.47
Nanhua
17.26
15.35
2.61
Jiangmen
15.11
10.13
5.79
Zhuyin
12.57
9.91
10.44
Sanduo
22.46
23.71
1.47
Rongqi
14.52
13.01
7.71
Denglongshan
5.15
8.73
12.46
Hengmen
5.22
3.35
16.81
Sanshakou
4.29
5.12
13.03
Nansha
5.18
3.37
8.91
Wanqinshaxi
5.14
4.11
14.43Figure 3
Change of maximum cut depth in Pearl River Delta from 1985 to 2006.We have calculatedR
upstream, R
sand-dredging, and R
sea-level rise for all the 13 stations in the PRD (Table 8). Effects due to upstream flood variation are negatively correlated with the distance between the location of the station and the outlets (Table 9). Changes of upstream flood flows have taken effects on the flood levels in income stations by Rs over 40%. It has changed the flood levels in network stations by 10% and 25% and the outlet stations between 4% and 5.5%. R
upstream peaks around income station Makou by 53.49% and comes to the minimum around outlet station Sanshakou by 4.29%. Located downstream of Sanshui, flood level in Sanduo has been affected the most by the income flood variation among all the network stations. Meanwhile, Tianhe and Nanhua have been similarly affected. R
upstream in Tianhe is higher than Jiangmen, and that in Jiangmen is higher than Zhuyin. It is not hard to discover that the further a network station is away from the income stations, the less the flood flow change will affect the on-site flood level. In general, variation in flood flow affects the most flood levels of income stations, then the network stations, and then the outlet stations. In other words, R
upstream decreases from income stations to outlet stations along the flow direction.Table 9
Distance among stations in West River Delta (units: km).
A
B
Distance between A and B
Makou
Denglongshan
123.8
Tianhe
67.8
Jiangmen
54.8
Zhuyin
19.8Riverbed alteration has taken effects on the flood levels by variations over 30% in both income stations, between 9 and 25% in network stations, and lower than 9% in outlet stations. The impact level of sand-dredging,R
sand-dredging, also decreases from income stations to outlet stations. R
sand-dredging in Sanshui is higher than that in Makou. It is consistent with the fact that riverbed alteration in Sanshui, from 1984 to 1995, was more severe than that in Makou (Figure 3 and Table 5). For R
sand-dredging of the network stations, Sanduo is the highest. R
sand-dredging in Tianhe, Nanhua, Jiangmen, and Rongqi are above 10%. R
sand-dredging in Zhuyin are the lowest, slightly below 10%. According to the riverbed alteration between 1985 and 2006 (Figure 3), riverbed in Sanduo has been cut down the most, then Nanhua, and then the other network stations. Despite the slight differences between R
sand-dredging in Tianhe, Jiangmen, and Zhuyin, R
sand-dredging of network stations are consistent with the extent of riverbed alteration. Among all the outlet stations, since the maximum cut depth varies the most in Denglongshan, flood level changes the most in Denglongshan under the impact of riverbed alteration.Sea-level rise has influenced the flood level in the PRD stations by a changing ration between 0.1% and 20%. It has hardly affected the flood level in income stations sinceR
sea-level rise are only 0.12% and 0.19%. Among network stations, R
sea-level rise tops in Zhuyin by 10.44%. Meanwhile, R
sea-level rise in other network stations are all below 6%. R
sea-level rise for Tianhe and Nanhua are close. In addition, R
sea-level rise in Tianhe is lower than that in Jiangmen, and that in Jiangmen is lower than Zhuyin at the same time. R
sea-level rise in outlet stations are all above 8%, much higher than that in income stations or most network stations. Similar to R
upstream, the value of R
sea-level rise seems to be correlated with the distance between stations and delta income, only in an opposite way. In another way, R
sea-level rise decreases from outlet stations to income stations in the the PRD.To compareR
upstream, R
sand-dredging, and R
sea-level rise of the income stations, it is easy to find that R
upstream is the highest and R
sea-level rise the lowest. It can be deducted that flood levels in income stations are influenced the most by the flood flow variation and the least by the sea-level rise. It is almost true for all the network stations as well, except for stations Sanduo and Zhuyin. In Sanduo, sand-dredging has caused the most change in on-site flood level. This is probably because that riverbed has been the most severely cut around Sanduo (Figure 3). In Zhuyin, somehow, the sea-level rise has taken the most effect. The cause will be talked about in Section 4.3. As for all the outlet stations, R
upstream and R
sand-dredging of outlet stations are close, and the difference between these two is all below 3.6%. At the same time, impact from sea-level rise ranks the first among the three driving factors.


## 4.2. Reliability of Quantitative Analysis Results

The efficiency of the quantitative cause analysis results can be generally verified by its consistency with the changing environment. To better tell the spatial patterns of impacts from different driving factors, the results have been interpolated on the whole delta with Kriging method and the results are shown in Figure4. Contour lines in Figure 4(a) show the spatial distribution of flood level change caused by variation of income flood flow. Figure 4(b) shows the spatial distribution of sand-dredging impact, and Figure 4(c) shows the spatial pattern of R
sea-level rise. As discussed in Section 3.2.1, variation of flood flow in Sanshui is greater than that in Makou. R
upstream in Sanshui is greater in these two stations as well. Meanwhile, R
upstream decreases when the station stands further from the upstream. R
sea-level rise also decreases when the station is further from the outlet. The impacts of the varied flood flow and sea-level rise are also changing mainly along the “northwest-southeast” direction, because river channels in the PRD are most developed along the “northwest-southeast” direction, and river discharge runs easily towards estuary by gravity force [23]. At the same time, R
sand-dredging of all network stations are in accordance with the extent of sand-dredging between the 1980s and 1990s (Figure 3). It is higher around Sanduo among network stations and lower around outlet stations. As is shown in Figure 3, channel varies in the greatest extent around Sanduo. Dredging of the riverbeds around Rongqi, Hengmen, Wanqinshaxi, Nansha, and Sanshakou is lower. Besides, by comparing the impacts of three driving factors, we found that, despite Sanduo, flood levels of all the income and network stations were most affected by the upstream flood flow. Previous works have examined the impact upstream flood flow had on network stations like Tianhe, and they argued that upstream flood flow is the most important factor for the PRD flood, despite the severe dredging of the riverbed and obvious sea-level rise [14, 23, 25, 34]. Our result here has verified this argument in network stations as well as the income stations.Figure 4
Spatial pattern ofR.
(a)
R
upstream
(b)
R
sand-dredging
(c)
R
sea-level riseCorrelation coefficients (CCs) among stations have also been calculated to double-check the reliability of our quantitative analysis result. CCs among stations that are not physically connected were not considered. For instance, there are no river channels connecting Tianhe and Rongqi. Thus, flood flow or flood level in Tianhe has no influence on the flood flow or level in Rongqi. These two stations are defined to be not physically connected. Except for the ones from 1985 to 1995, flood level series from 1959 to 2011 of all gauge stations have been extracted to calculate CC. Applied flood level series of each station consists of extreme water level in flood season for 42 years. CCs among stations are thus calculated by the following equation:(8)
C
o
v
X
,
Y
=
E
X
Y
-
E
X
E
Y
,where X is the time series of one station and Y the other; E, the expectation of a time series; Cov, the correlation coefficient of two time series.Higher value of CC indicates better correlation of hydrological time series in two stations. As is shown in Table10, flood level in Tianhe is related to Nanhua the most, and it is also true for Nanhua the other way around. Flood levels in these two stations take almost the same impact from flood variation, sand-dredging, and sea-level rise (Table 8). CCs among stations are positively connected with the distance among stations (Tables 9 and 10). It is the same as the spatial pattern of R
upstream and R
sea-level rise. Furthermore, flood levels in all the network stations are both connected with Sanshui and Makou, but CCs of stations with Sanshui are higher. It is true in Figure 4(a) that R
upstream is higher in Sanshui than Makou. What is more, as for the connection between network stations and outlet stations, CC between Zhuyin and Denglongshan is 2 times higher than any CC between stations downstream of Makou and other outlet stations. Last but not the least, Rongqi is more related with Wanqinshaxi than Nansha (Table 10) and it is the same in Figure 4(c) that R
sea-level rise of Rongqi and Wanqinshaxi is closer.Table 10
Correlation coefficient (CC) among physically relevant stations.
Tianhe
Jiangmen
Zhuyin
Nanhua
Rongqi
Sanduo
Makou
0.834
0.801
0.548
0.834
0.739
0.804
Sanshui
0.966
0.941
0.667
0.964
0.878
0.938
Nanhua
0.9997
0.983
0.681
—
0.939
—
Rongqi
—
—
—
0.939
—
—
Tianhe
—
0.983
0.684
0.9997
—
—
Jiangmen
0.983
—
0.696
0.983
—
—
Zhuyin
0.684
0.696
—
0.681
—
—
Sanshakou
—
—
—
—
—
0.383
Nansha
—
—
—
—
—
0.423
Wanqinshaxi
—
—
—
0.427
0.461
0.397
Hengmen
—
—
—
0.223
0.360
—
Denglongshan
0.241
0.265
0.723
—
—
—


## 4.3. Impact of Network Density

Obviously,R
upstream and R
sea-level rise change from income stations to outlet stations along the flow direction. However, the changing rate slows down around Jiangmen and Zhuyin. There is 35 km between Jiangmen and Zhuyin, almost 3 times the distance between Jiangmen and Tianhe (Table 9). As discussed above, the impacts on flood level by inflow flood change and rising sea level are positively correlated with spatial distances among stations. In this way, the difference of R
upstream or R
sea-level rise between Jiangmen and Zhuyin is supposed to be almost 3 times that between Jiangmen and Tianhe. The difference of R
upstream between Jiangmen and Tianhe is 3.39%. Somehow, the one between Jiangmen and Zhuyin is 2.54% rather than 10.17%. Meanwhile, the difference of R
sea-level rise between Jiangmen and Tianhe is 2.61%, and the one between Jiangmen and Zhuyin is less than that by being only 0.22%.Density of river network is inferred to take up the major responsibility. As is shown in Figure1 and mentioned in Section 2.1, river networks around Jiangmen are significantly denser than other places in the the PRD. There are more river nodes within the same area. At a node, flood flows towards diverse directions, or flood flow from multiple branches joins in the same channel. Due to the terrain features in the PRD, river channels have developed the most along the “northwest-southeast” direction [23], and the velocity of flood flow slows down in river channels that have developed otherwise. Flood from upstream is thus held up among dense river networks. As a result, variation of income flood takes more impact on the dense river network.Meanwhile, sea-level rise takes effect on the flood level through the fluctuation of water body. Its effect tracks back up the river way and jerks up the in-channel flood within dense river network when the flow velocity drops. In this way, the impact of sea-level rise on flood level has also been intensified in dense river network. To conclude, dense river network in delta area intensifies the impact on flood level from both the income flood change and sea-level rise, which explains the irregular pattern ofR
upstream and R
sea-level rise around Jiangmen and Zhuyin.


## 5. Conclusions

Flood levels in the Pearl River Delta of China are typically affected by three major factors, that is, changing flood flow, riverbed alteration caused by sand-dredging, and sea-level rise. This study has adopted the multiple variant regression method proposed by Huang et al. [14] and used and index R inspired by Jiang et al. [31] to quantify impacts of these three factors on flood levels in the PRD. It was found that the impact of flood flow variation on flood level decreased from income stations to outlet stations. The impact of sea-level rise on the flood level decreased for station moving away from the tidal station on coast. Flood level was positively related to the extent of sand-dredging between the 1980s and 1990s as a result of the deepening of channel bottom. Spatial pattern of quantitative impacts from three driving factors is consistent with both the changing pattern of driving factors and the correlations among stations in the PRD (Figure 4 and Table 10). Among income and network stations, flood flow variation causes the most influence on on-site flood level, then the dredging of the riverbed, and at last the sea-level rise (Table 8). High channel density can strengthen the impact from factor variation resulting from both change in upstream flood flow and sea-level rise around delta outlets.Result from quantitative cause analysis in this study would help for further understanding of the connection between cause and effect and for better misunderstanding of the causes of floods in the PRD which will further help in developing flood protection strategies. Recent works have focused on the delta flood threat from sea-level rise. However, results in the PRD show that flood variation is the greatest threat to delta flood. In addition, the results of this research have also revealed that even severe dredging of the riverbed resulting from human interference takes more effect than the sea-level rise. It may indicate that more attention should be put on the effect from changing income flood flow and urbanization when people are trying to protect delta from flood. Also, floods can happen more easily in dense network. Since dense river network can strengthen the effect of changing flood flow and sea-level rise, people may need to worry more about the disastrous result brought about by climate change within dense river network.



---

*Source: 3928920-2015-12-28.xml*