Flooding in Delta Areas under Changing Climate: Response of Design Flood Level to Non-Stationarity in Both Inﬂow Floods and High Tides in South China

: Climate change has led to non-stationarity in recorded floods all over the world. Although previous studies have widely discussed the design error caused by non-stationarity, most of them explored basins with closed catchment areas. The response of flood level to nonstationary inflow floods and high tidal levels in deltas with a dense river network has hardly been mentioned. Delta areas are extremely vulnerable to floods. To establish reliable standards for flood protection in delta areas, it is crucial to investigate the response of flood level to nonstationary inflow floods and high tidal levels. Pearl River Delta (PRD), the largest delta in South China, was selected as the study area. A theoretical framework was developed to quantify the response of flood level to nonstationary inflow floods and the tidal level. When the non-stationarity was ignored, error up to 18% was found in 100-year design inflow floods and up to 14% in 100-year design tidal level. Meanwhile, flood level in areas that were ≤ 22 km away from the outlets mainly responded to the nonstationary tidal level, and that ≥ 45 km to the nonstationary inflow floods. This study will support research on the non-stationarity of floods in delta areas.


Introduction
Since the early twentieth century, climate change has altered hydrological cycles all over the world. Flood records are no longer stationary [1][2][3][4][5]. In particular, delta areas have experienced intensified flooding due to the combined impacts of the inflow floods from the upstream watershed and the rising high tidal level induced by sea-level rise (SLR). Meanwhile, livable climate and convenient sea transportation have contributed to increasing population and economy in delta areas around the world [6,7]. Thus, delta zones have become extremely vulnerable to floods. It is important to understand the response of flood level to nonstationary inflow floods and high tidal levels in delta areas, especially ones with a dense river network.
There were some previous studies discussing the variation in design tidal levels under the impact of SLR. Trends in recording series, as well as the correlation between these trends and other climate factors, were explored using sophisticated extreme value analysis methodology [8]. To calculate reliable at-site frequency of the increasing tidal level, many previous researchers generated a stationary series before carrying out frequency analysis. Some restored the stationary time series by removing the increment caused by SLR [9][10][11], while some tried to generate stationary samples by duplicating the characteristics of historical tidal processes, which were assumed to be stationary [12]. Other researchers assumed that the rate of high tidal levels was time-dependent, and conducted frequency analysis on the recorded series with a time-correlated variable [8][9][10][11][12][13]. [12]. Other researchers assumed that the rate of high tidal levels was time-dependent, and conducted frequency analysis on the recorded series with a time-correlated variable [8][9][10][11][12][13].
Design flood levels in delta areas were also analyzed by statistical methods in past studies. Joint distribution was the most widely used approach to statistically estimate flood levels in delta areas, which was simultaneously affected by multiple forcing factors like rainstorms, storm surges, upstream floods, and high tides. To take non-stationarity caused by environmental changes into concern, joint distribution incorporated the nonstationary characteristics of changing factors into marginal distribution functions. However, in this way, both the difficulty of parameter estimation and the uncertainty of results were significantly increased [14]. Another commonly used method for delta flood protection is to first calculate two design flood levels, under the impact of either the inflow floods or the high tidal level, and choose the higher one as the design value [15]. This method is simple and convenient, but it neglects the interconnection between inflow floods and high tides.
This research aims to estimate the response of flood level to the impact of inflow floods and downstream tidal level, when the time series of these two factors are both nonstationary. To achieve this purpose, Pearl River Delta (PRD) in South China was selected as the study area. This paper is arranged as follows. In Section 2, a brief introduction of the study area, dataset and methodology is provided. In Section 3, the response of the flood level to the nonstationary inflow floods and the high tidal level is estimated, and its spatial pattern is discussed. Conclusions are made in Section 4. The investigation in PRD can be used as reference for estimations of the flood level response in other delta areas with similarly complicated environmental conditions.

Study Area and Dataset
Pearl River Delta, which covers an area of 39,380 km 2 , is the largest delta in South China ( Figure 1). It has a humid subtropical climate. The monsoon period in this area lasts from April to September. Due to the terrain altitude, the dense river network develops from the northwest to the southeast. Water flow in PRD comes from West River and North River, and all the water in the delta flows into Pearl River Estuary through multiple outlets. Based on its source of inflow floods water, PRD area can be further divided into a West River sub-delta, and a North River sub-delta. Floods in two upstream river basins mainly occur between June and August. Meanwhile, extreme tides in the estuary, induced by typhoons or storm surges, mainly occur around September [16].  PRD is the megaregion where Macao is located [17]. In 2014, the population of cities in PRD surpassed that of Chicago and it was ranked in the top 30 in the world [18]. In 2015, the total income Water 2017, 9, 471 3 of 16 of only two cities in PRD, i.e., Guangzhou and Foshan, reached $20.6 million [19]. With its fast growth in economy and population, PRD is extremely vulnerable to flood risk.
In this research, recorded peak flood volumes in the flood season, i.e., from June to August, in Wuzhou and Shijiao were used. So were the high tidal levels in seven outlet stations ( Figure 1, Table 1). All the data were collected during the period from 1958 to 2011.

Methodology
In this research, a methodology was developed to estimate the impact of the nonstationary inflow floods and the high tidal level on flood level. First, a time-varying moments (TVM) model was applied to evaluate the design inflow floods from the upstream river basins and the high tidal levels around the delta estuaries. Then, the influence of the nonstationary inflow floods and the high tidal level on flood level was estimated through a one-dimensional hydrodynamic model using the TVM-estimated design inflow floods and the high tidal level as its inputs. Four scenarios were created by combining the nonstationary inflow floods and the high tidal level, and then the nonstationary water level in PRD was estimated under these scenarios. The whole schedule of this combined methodology is shown in the flowchart presented in Figure 2. All the parts that are considered nonstationary are highlighted in blue.
Water 2017, 9,471 3 of 16 PRD is the megaregion where Macao is located [17]. In 2014, the population of cities in PRD surpassed that of Chicago and it was ranked in the top 30 in the world [18]. In 2015, the total income of only two cities in PRD, i.e., Guangzhou and Foshan, reached $20.6 million [19]. With its fast growth in economy and population, PRD is extremely vulnerable to flood risk.
In this research, recorded peak flood volumes in the flood season, i.e., from June to August, in Wuzhou and Shijiao were used. So were the high tidal levels in seven outlet stations ( Figure 1, Table 1). All the data were collected during the period from 1958 to 2011.

Methodology
In this research, a methodology was developed to estimate the impact of the nonstationary inflow floods and the high tidal level on flood level. First, a time-varying moments (TVM) model was applied to evaluate the design inflow floods from the upstream river basins and the high tidal levels around the delta estuaries. Then, the influence of the nonstationary inflow floods and the high tidal level on flood level was estimated through a one-dimensional hydrodynamic model using the TVM-estimated design inflow floods and the high tidal level as its inputs. Four scenarios were created by combining the nonstationary inflow floods and the high tidal level, and then the nonstationary water level in PRD was estimated under these scenarios. The whole schedule of this combined methodology is shown in the flowchart presented in Figure 2. All the parts that are considered nonstationary are highlighted in blue.

Time-Varying Moments Estimating Nonstationary Inflow Floods and High Tidal Levels
The time-varying moments method, which has been widely applied in the statistical estimation of nonstationary peak flood volume, was used in this investigation to calculate the nonstationary design value of inflow floods and high tidal level [20]. In this method, to incorporate the non-stationarity, the first two moments of the time series were initially assumed to have trends that are either linear (L) or parabolic (P) [21,22].
Four conditions were generated by assembling the L or P trend in the first two moments, labeled A to D [22]: (1) when only the first moment, i.e., the mean of the time series, is time-dependent, whether the trend is L or P, the condition is named A; (2) when only the second moment, i.e., the standard deviation of the time series, is time-dependent, whether the trend is L or P, the condition is named B; (3) when both the mean and the standard deviation have time-dependent trends in, and the standard deviation is a product of, the mean and a constant value, i.e., the coefficient of variation (Cv), the condition is named C; (4) when both the first two moments have trends, but they are not directly correlated, the condition is named D. By combining Scenarios A-D and two types of trends, i.e., L and P, eight time-trend models were generated. All models are listed in Table 2, where m represents the mean of the nonstationary time series that is defined in relation to both the mean of a stationary time series, m, and time, t. Likewise, σ denotes the nonstationary standard deviation and is considered to be associated with both the stationary standard deviation σ and time, t. In all functions, a m , a σ , b m , and b σ denote the constants.
Note: P3 denotes the Pearson type III.
According to studies carried out in the 1990s, Pearson type III performs the best among all commonly used probability distribution functions for either the flood flow or the high tidal level in PRD [23][24][25][26][27][28][29][30]. Therefore, Pearson type III was selected as the probability density function (PDF) used in design value estimation in this research. Based on the eight time trend models, parameters of Pearson type III were all defined to be time-dependent variables. All these parameters were estimated by the maximum likelihood method, where t was used as a variable. The goodness of fit of PDFs, i.e., Pearson type III with parameters of different time trend models, was tested by the Akaike Information Criterion (AIC).

One-Dimensional River Network Hydrodynamic Model
The one-dimensional hydrodynamic model has been widely used in flood level estimation. It is comparatively more accurate in representing the cross-sectional area of channels and requires only a small amount of field data to set up the model. Also, it can be quickly set up, and the computation is fast [31]. Therefore, the one-dimensional hydrodynamic model was chosen as the basic tool to estimate the flood level in this research.

Model Setup and Input Data
The one-dimensional hydrodynamic model used in the present study was developed from the one created by the Center for Water Resources and Environment (2015), and the whole model was set up in Fortran [32]. This model adopted the Saint-Venant continuity equation and the longitudinal momentum balance equation as the governing equations. The implicitly weighted four-point Preissmann scheme [33] was used to solve these governing equations.
The model set the time-step as 10 min, and a space lag from 500 to 2500 m. In all, the model included 196 channels and 144 nodes ( Figure 3). The upstream inflow from Gaoyao (West River) and Shijiao (North River), as well as the tidal levels in seven outlets ( Figure 1, Table 1), were set as its boundaries, and the recorded data in these stations were used as the input to run the model. momentum balance equation as the governing equations. The implicitly weighted four-point Preissmann scheme [33] was used to solve these governing equations. The model set the time-step as 10 min, and a space lag from 500 to 2500 m. In all, the model included 196 channels and 144 nodes ( Figure 3). The upstream inflow from Gaoyao (West River) and Shijiao (North River), as well as the tidal levels in seven outlets ( Figure 1, Table 1), were set as its boundaries, and the recorded data in these stations were used as the input to run the model. Conducting simultaneous field survey over a large basin like PRD region is extremely time-consuming and expensive. Meanwhile, although large-scale field measurement has been carried out by the Guangdong Hydraulic Bureau every five years, the data are beyond our access. Owing to the limitations of survey data, this research applied the topographic data collected in the flood season of 1998 to develop a 1D model. All topographic maps were obtained from the Guangdong Hydraulic Bureau, in a measuring scale of 1:5000. The hourly flood flow and water level were simultaneously detected at 40 stations during a typical flood in PRD. Data collected from 12:00 p.m. on 25 June 1998 to 12:00 p.m. on 28 June 1998 at hydrological stations in Gaoyao and Shijiao were used to drive the model and calibrate the parameters. The fluvial water level at each sampling frequency was linearly interpolated to the model time-step.  Hydraulic Bureau every five years, the data are beyond our access. Owing to the limitations of survey data, this research applied the topographic data collected in the flood season of 1998 to develop a 1D model. All topographic maps were obtained from the Guangdong Hydraulic Bureau, in a measuring scale of 1:5000. The hourly flood flow and water level were simultaneously detected at 40 stations during a typical flood in PRD. Data collected from 12:00 p.m. on 25 June 1998 to 12:00 p.m. on 28 June 1998 at hydrological stations in Gaoyao and Shijiao were used to drive the model and calibrate the parameters. The fluvial water level at each sampling frequency was linearly interpolated to the model time-step.

Model Validation
The hydrodynamic model was validated before estimating the flood level at PRD stations. The coefficient of riverbed roughness was first set between 0.016 and 0.035 [33]. Next, after the model had been running for 30 days, the parameters were updated by the results obtained on day 30. The aforementioned procedures were repeated more than 10 times before the simulation results stabilized. The model parameters were finally set when the difference between the simulation results of day 1 and day 30 remained within 1%.
To assess the performance of the model, measured and model-estimated amplitudes for eight constituents were compared in six stations, i.e., Makou, Sanshui, Tianhe, Baijiao, Xiaheng, and Shaluowei, and the result was shown in Figure 4. Difference bands were both plotted at ±0.025 and ±0.05 m, with an overall RMSE of 6.87 cm. For each of these six stations, recorded and estimated flood levels were also compared. The average absolute errors in the six stations were below 0.05 m, with a minimum of 0.001 m at Baijiao. To further measure the model efficiency, an index (α) was applied, was calculated by the following equations: where Z is the estimated flood level; Z 0 , the recorded flood level; and i, the number of hours. When α ≥ 75%, the simulation of the model was predicted to be valid. The value of α is 95% in Makou, 83% in Sanshui, 90% in Tianhe, 92% in Baijiao, 81% in Xiaheng, and 92% in Shaluowei. All αs are above 80%, and the model has been proved to be valid.

Model Validation
The hydrodynamic model was validated before estimating the flood level at PRD stations. The coefficient of riverbed roughness was first set between 0.016 and 0.035 [33]. Next, after the model had been running for 30 days, the parameters were updated by the results obtained on day 30. The aforementioned procedures were repeated more than 10 times before the simulation results stabilized. The model parameters were finally set when the difference between the simulation results of day 1 and day 30 remained within 1%.
To assess the performance of the model, measured and model-estimated amplitudes for eight constituents were compared in six stations, i.e., Makou, Sanshui, Tianhe, Baijiao, Xiaheng, and Shaluowei, and the result was shown in Figure 4. Difference bands were both plotted at ±0.025 and ±0.05 m, with an overall RMSE of 6.87 cm. For each of these six stations, recorded and estimated flood levels were also compared. The average absolute errors in the six stations were below 0.05 m, with a minimum of 0.001 m at Baijiao. To further measure the model efficiency, an index (α) was applied, was calculated by the following equations: where is the estimated flood level; , the recorded flood level; and i, the number of hours. When α 75%, the simulation of the model was predicted to be valid. The value of α is 95% in Makou, 83% in Sanshui, 90% in Tianhe, 92% in Baijiao, 81% in Xiaheng, and 92% in Shaluowei. All αs are above 80%, and the model has been proved to be valid.

Increasing Inflow Floods
Trends in the inflow floods of PRD were estimated by Mann-Kendall method [34,35]. As is shown in Table 1, inflow floods from West River had significantly increased since 1958. Since 1950, global warming and atmosphere circulation abnormities have led to fewer but more intense tropical cyclones in West River [36][37][38]. Meanwhile, the extreme precipitations (over 99 percentile) in West River have increased by more than 45 mm since 1959 [39]. As a result, the flood flow in lower West River became significantly higher in the flood season [40].
Based on the trend test results, design flood in West River was estimated by TVM, and CL was detected to be the best-fit time-trend model for the increasing inflow floods ( Table 2). Thirty-year simple moving mean and standard deviation of annual peak flood flow series in West River was analyzed and displayed in Figure 5. Since the first two moments of the flood series both had linear trends, the results of AIC test were proven to be reliable.

Increasing Inflow Floods
Trends in the inflow floods of PRD were estimated by Mann-Kendall method [34,35]. As is shown in Table 1, inflow floods from West River had significantly increased since 1958. Since 1950, global warming and atmosphere circulation abnormities have led to fewer but more intense tropical cyclones in West River [36][37][38]. Meanwhile, the extreme precipitations (over 99 percentile) in West River have increased by more than 45 mm since 1959 [39]. As a result, the flood flow in lower West River became significantly higher in the flood season [40].
Based on the trend test results, design flood in West River was estimated by TVM, and CL was detected to be the best-fit time-trend model for the increasing inflow floods ( Table 2). Thirty-year simple moving mean and standard deviation of annual peak flood flow series in West River was analyzed and displayed in Figure 5. Since the first two moments of the flood series both had linear trends, the results of AIC test were proven to be reliable.  The 100-, 50-, and 20-year design flood flow in West River in 2050 were then estimated. To compare the design volumes with and without taking non-stationarity into consideration, the relative difference E was calculated as follows: where X T denotes the nonstationary design value, and X 0 denotes the stationary one. The gap between a 100-year and a 20-year X T is ∆X T , and the difference between a 100-year and a 20-year X 0 is ∆X 0 . Both ∆X T and ∆X 0 were calculated and compared. The greater the value of ∆X T or ∆X 0 , the steeper the slope of the upper tail (and the greater the pressure in flood protection. (The "upper tail" here specifically refers to the part of the flood frequency curve with an exceedance probability lying between 1% and 5%.) Results showed that by ignoring the increasing tendency, the 100-year peak flood flow was underestimated by 17.51% (9825 m 3 /s) (Table 3). Meanwhile, the relative underestimation (E) was augmented when the exceeding probability reached 5%, i.e., by 17.91% (9417 m 3 /s) for the 50-year flood, and by 18.52% (8809 m 3 /s) for the 20-year flood. In addition, ∆X T (9553 m 3 /s) was larger than ∆X 0 (8537 m 3 /s), indicating that non-stationarity in the flood had contributed to a steeper upper tail. Table 3. Comparison of design flood.

Varied High Tidal Level
High tidal levels in seven outlets were tested, and tidal levels in all outlets but Nansha had upward tendencies. Somehow, only the tendencies in Hengmen, Denglongshan, Huangjin, and Xipaotai reached a significance of 90% (Table 1). Rising sea-level rise at the Pearl River Estuary pushed up the high tidal levels, and severe sand excavation and channel regulation from the 1980s to the early 21st century around the upstream of Nansha neutralized the rising trend [41].
Nonstationary estimations were carried out in the six outlets that had significant trends. In each station, best-fit models was detected, and the nonstationary design value was compared with the stationary one (Tables 2 and 3). Results showed that, whether the return period was 100, 50, or 20 years, design levels were all underestimated in Hengmen, Denglongshan, and Huangjin, and the increase in the return period had led to a greater underestimation, which was displayed in the value of E. Meanwhile, the 100-, 50-, and 20-year Es ranged from 9.33% (0.25 m) to 11.16% (0.27 m) in Hengmen, 14.07% (0.37 m) to 16.88% (0.39 m) in Denglongshan, and 14.06% (0.45 m) to 20.15% (0.55 m) in Huangjin. Both X T and Es decreased along the coastline from the southwest to the northeast, i.e., from Huangjin and Denglongshan to Hengmen. Furthermore, since ∆X T reached 0.24 m in Hengmen, 0.30 m in Denglongshan, and 0.37 m in Huangjin, it is easy to conclude that the slope of upper tail was steeper in the southwest (Huangjin) than in the northeast (Hengmen).
The direction of the ocean current affected the impact of SLR on the tidal level. Based on this criteria, Wu [42] divided the PRD outlets into three groups. In the first group, to which Huangjin belonged, high tides in outlet were mainly influenced by the current from the Indian Ocean. In the second group, to which Denglongshan belonged, high tidal levels were affected by a substantial impact of ocean currents from both the Pacific Ocean and the Indian Ocean. In the third group, to which Hengmen belonged, high tidal levels were solely influenced by the Pacific Ocean current.
The surrounding environment of all these outlet stations affected the ocean current impact as well. Take Huangjin as an example: since it was located in a bay surrounded by tiny islands, the tidal vibration affected by the current impact was amplified. In Hengmen, conversely, the tidal vibration induced by the current effect was reduced, due to the land lying to the east of Pearl River Estuary [42,43].

Response of the Design Flood Level to Nonstationary Inflow Floods and the High Tidal Level
Design flood levels in 22 stations were estimated under six scenarios, combining nonstationary inflow floods and nonstationary high tidal level (Table 4, Figure 6). The design flood level estimated with these nonstationary inputs was called the nonstationary flood level, whereas the one without considering non-stationarity was named the stationary flood level. For all 22 stations, nonstationary and the stationary flood levels were compared under Scenarios 1 and 4.
Water 2017, 9, 471 9 of 16 tidal vibration affected by the current impact was amplified. In Hengmen, conversely, the tidal vibration induced by the current effect was reduced, due to the land lying to the east of Pearl River Estuary [42,43].

Response of the Design Flood Level to Nonstationary Inflow Floods and the High Tidal Level
Design flood levels in 22 stations were estimated under six scenarios, combining nonstationary inflow floods and nonstationary high tidal level (Table 4, Figure 6). The design flood level estimated with these nonstationary inputs was called the nonstationary flood level, whereas the one without considering non-stationarity was named the stationary flood level. For all 22 stations, nonstationary and the stationary flood levels were compared under Scenarios 1 and 4.  In 21 out of 22 stations, the nonstationary flood level was higher than the stationary. Only in Hengshan was the stationary flood level higher. This is because Hengshan was located just 18 km away from the outlet at Xipaotai, and its flood level was highly affected by the nonstationary design tidal level in Xipaotai. Due to the result in Table 1, the nonstationary design tidal level in Xipaotai was lower than the stationary one. Therefore, the stationary flood level in Hengshan was higher than the nonstationary one.  In 21 out of 22 stations, the nonstationary flood level was higher than the stationary. Only in Hengshan was the stationary flood level higher. This is because Hengshan was located just 18 km away from the outlet at Xipaotai, and its flood level was highly affected by the nonstationary design tidal level in Xipaotai. Due to the result in Table 1, the nonstationary design tidal level in Xipaotai was lower than the stationary one. Therefore, the stationary flood level in Hengshan was higher than the nonstationary one.
To further assess the response of flood level to the nonstationary inflow floods and high tidal level in PRD, the 22 stations were divided into three groups by the K-means clustering method, based on the absolute difference between their nonstationary and stationary flood levels (Table 5). Among the Under both scenarios, five stations, i.e., Baiqing, Da'ao, Lianyao, Zhuyin, and Zhuzhou, belonged to the first group. These five stations were located between 25 and 42 km away from the delta outlets, and their flood levels were similarly affected by the high tidal levels. Meanwhile, the responses in Banshawei, Sanshanjiao, and Fengmamiao stations were the lowest under both scenarios. The distance between these stations and their outlets was within the area of 15-33 km. Outlets of these stations, i.e., Wanqingshaxi and Sanshanjiao, exhibited no significant trend in the high tidal level.
When the scenario changed from 1 to 4, stations located in the West River sub-delta, i.e., Makou, Baijiao, Jiangmen, Nanhua, and Tianhe, fell from groups with higher mean values to ones with lower mean values. Meanwile, the stations in the North River sub-delta, i.e., Zidong and Sanduo stations, jumped from the group with lower mean (Group 2) to the group with a higher one (Group 1). Alteration in geometries in the upstream river channel increased the flood flow distribution ratio in the North River sub-delta. Under Scenario 1, the flow distribution ratio in Sanshui, i.e., the income of the North River sub-delta, was 0.7% (300 m 3 /s), higher than under Scenario 4 ( Table 3). In other words, the increase in the upstream flood flow in the North River sub-delta was greater under Scenario 4, while the increment of income flood flow in the West River sub-delta dropped by 300 m 3 /s.

Response of the Flood Level along the River Channel
The flood level in each station along the same river channel was connected into a flood line, and the 100-year nonstationary flood line was compared with the stationary one in two representative river channels (Figures 6 and 7). Makou-Zhuyin channel was located in the West River sub-delta, where there were five stations along the way, i.e., Makou, Tianhe, Jiangmen, Da'ao, and Zhuyin. Sanshui-Sanshanjiao channel was located in the North River sub-delta, and Sanshui, Sanduo, and Sanshanjiao were located along the river channel from the delta income to the outlet. To further assess the response of flood level to the nonstationary inflow floods and high tidal level in PRD, the 22 stations were divided into three groups by the K-means clustering method, based on the absolute difference between their nonstationary and stationary flood levels (Table 5). Among the three groups, Group 1 had the highest mean value, and Group 3 the lowest. If a station belonged to the same group under both scenarios, it was marked in blue (Table 5). Under both scenarios, five stations, i.e., Baiqing, Da'ao, Lianyao, Zhuyin, and Zhuzhou, belonged to the first group. These five stations were located between 25 and 42 km away from the delta outlets, and their flood levels were similarly affected by the high tidal levels. Meanwhile, the responses in Banshawei, Sanshanjiao, and Fengmamiao stations were the lowest under both scenarios. The distance between these stations and their outlets was within the area of 15-33 km. Outlets of these stations, i.e., Wanqingshaxi and Sanshanjiao, exhibited no significant trend in the high tidal level.
When the scenario changed from 1 to 4, stations located in the West River sub-delta, i.e., Makou, Baijiao, Jiangmen, Nanhua, and Tianhe, fell from groups with higher mean values to ones with lower mean values. Meanwile, the stations in the North River sub-delta, i.e., Zidong and Sanduo stations, jumped from the group with lower mean (Group 2) to the group with a higher one (Group 1). Alteration in geometries in the upstream river channel increased the flood flow distribution ratio in the North River sub-delta. Under Scenario 1, the flow distribution ratio in Sanshui, i.e., the income of the North River sub-delta, was 0.7% (300 m 3 /s), higher than under Scenario 4 (Table 3). In other words, the increase in the upstream flood flow in the North River sub-delta was greater under Scenario 4, while the increment of income flood flow in the West River sub-delta dropped by 300 m 3 /s.

Response of the Flood Level along the River Channel
The flood level in each station along the same river channel was connected into a flood line, and the 100-year nonstationary flood line was compared with the stationary one in two representative river channels (Figures 6 and 7). Makou-Zhuyin channel was located in the West River sub-delta, where there were five stations along the way, i.e., Makou, Tianhe, Jiangmen, Da'ao, and Zhuyin. Sanshui-Sanshanjiao channel was located in the North River sub-delta, and Sanshui, Sanduo, and Sanshanjiao were located along the river channel from the delta income to the outlet.  In either channel, the nonstationary flood line was higher than the stationary one, and the absolute gap between two flood lines gradually narrowed from the delta incomes to the outlets (Figure 7a,b). The relative difference E in the two channels varied dramatically in the spot 60 km away from the outlet (Figure 7c,d). When it came around the area that was 60 km away from the outlet, E kept dropping along the Sanshui-Sanshanjiao channel, while it increased sharply in Makou-Zhuyin channel. In Sanshui-Sanshanjiao channel, since the high tidal level in the outlet displayed no significant upward trend (Table 1), the increase in the design flood line mainly responded to the upstream floods' increment. Thus, the impact of the nonstationary inflow floods decreased smoothly along the channel solely due to the reduction in the elevation. On the contrary, however, the increment of the 100-year tidal level reached 14.07% in the outlet of Makou-Zhuyin channel (Table 3). Affected by the interactive impact of increasing inflow floods and high tidal level, the relative difference between flood lines increased in the area that was 60 km away from the outlet.
Both the absolute and the relative difference between the nonstationary and stationary flood line declined sharply in Jiangmen, caused by the flow distribution ratio variation in its upstream, i.e., the node where Tianhe and Nanhua parted. A large portion of the flood flow (45.53%) from Makou went through Tianhe, located 13 km away from Jiangmen. However, the width-depth ratio in Tianhe was only half that in Makou [44]. Therefore, the relative gap between nonstationary and stationary flood level in Tianhe increased. Nevertheless, the width-depth ratio in Jiangmen is 2.4 times that in Tianhe [44]. Thus, the relative gap between the nonstationary and the stationary design level in Jiangmen was narrowed abruptly from Tianhe. In either channel, the nonstationary flood line was higher than the stationary one, and the absolute gap between two flood lines gradually narrowed from the delta incomes to the outlets (Figure 7a,b). The relative difference E in the two channels varied dramatically in the spot 60 km away from the outlet (Figure 7c,d). When it came around the area that was 60 km away from the outlet, E kept dropping along the Sanshui-Sanshanjiao channel, while it increased sharply in Makou-Zhuyin channel. In Sanshui-Sanshanjiao channel, since the high tidal level in the outlet displayed no significant upward trend (Table 1), the increase in the design flood line mainly responded to the upstream floods' increment. Thus, the impact of the nonstationary inflow floods decreased smoothly along the channel solely due to the reduction in the elevation. On the contrary, however, the increment of the 100-year tidal level reached 14.07% in the outlet of Makou-Zhuyin channel (Table 3). Affected by the interactive impact of increasing inflow floods and high tidal level, the relative difference between flood lines increased in the area that was 60 km away from the outlet.
Both the absolute and the relative difference between the nonstationary and stationary flood line declined sharply in Jiangmen, caused by the flow distribution ratio variation in its upstream, i.e., the node where Tianhe and Nanhua parted. A large portion of the flood flow (45.53%) from Makou went through Tianhe, located 13 km away from Jiangmen. However, the width-depth ratio in Tianhe was only half that in Makou [44]. Therefore, the relative gap between nonstationary and stationary flood level in Tianhe increased. Nevertheless, the width-depth ratio in Jiangmen is 2.4 times that in Tianhe [44]. Thus, the relative gap between the nonstationary and the stationary design level in Jiangmen was narrowed abruptly from Tianhe.

Leading Factors of Non-Stationarity and the Flood Level Response
Nonstationary flood levels under different scenarios were further calculated and compared. To compare the response of the flood level to nonstationary inflow floods and that to the nonstationary tidal level, index X was developed and calculated by the following equation: where Z i (i = 1, 2, 3) is the nonstationary design flood level simulated under scenario i. When X ≤ 0.5, the flood level in PRD mainly responded to the nonstationary high tide, whereas when X ≥ 5, the flood level mainly responded to the nonstationary inflow floods.
In general, X in each station was positively correlated to the distance between this station and the delta outlet ( Figure 8). Xs of stations that were less than 22 km away from the PRD outlet were all below 0.5. Meanwhile, when the distances between stations and the delta outlet were more than 45 km, all Xs were above 5. It can be deducted that in PRD, in stations that are more than 45 km away from the outlets, the flood level was mainly influenced by the nonstationary inflow floods. Furthermore, flood levels in stations that are located within 22 km of the delta outlets were strongly affected by the nonstationary high tidal level.

Leading Factors of Non-Stationarity and the Flood Level Response
Nonstationary flood levels under different scenarios were further calculated and compared. To compare the response of the flood level to nonstationary inflow floods and that to the nonstationary tidal level, index X was developed and calculated by the following equation: where (i = 1, 2, 3) is the nonstationary design flood level simulated under scenario i. When 0.5, the flood level in PRD mainly responded to the nonstationary high tide, whereas when 5, the flood level mainly responded to the nonstationary inflow floods.
In general, X in each station was positively correlated to the distance between this station and the delta outlet ( Figure 8). Xs of stations that were less than 22 km away from the PRD outlet were all below 0.5. Meanwhile, when the distances between stations and the delta outlet were more than 45 km, all Xs were above 5. It can be deducted that in PRD, in stations that are more than 45 km away from the outlets, the flood level was mainly influenced by the nonstationary inflow floods. Furthermore, flood levels in stations that are located within 22 km of the delta outlets were strongly affected by the nonstationary high tidal level.  As for the flood level in stations lying between 22-45 km from delta outlets, their response to either nonstationary inflow floods did not increase when they became farther away from the delta outlets. Take Zhuzhou, for example; though it is located 39 km away from the outlet, its X was the same as that of Banshawei, which was 33 km away from the outlet. Since the impact of extreme tides decreased when the distance between the station and delta outlets increased, X in Zhuzhou was supposed to be higher than that in Banshawei. On the other hand, however, Zhuzhou was much farther away from the entrance of PRD than Banshawei ( Figure 8). Therefore, the impact of upstream change decreased to a greater extent in Zhuzhou than in Banshawei. Moreover, the outlet of Zhuzhou had a more severe upward tendency than Banshawei (Table 1). Thus, the non-stationarity of the high tidal level in the outlet of Zhuzhou exerted a more severe impact on the flood level in Zhuzhou.
Although Rongqi was located 2 km farther from the delta outlet than Xiaolan [44], X in Rongqi was the same as in Xiaolan (Figure 8). Simply according to the distance between stations and the delta outlets, X in Rongqi was supposed to be higher than in Xiaolan. However, the estimated X was inconsistent with our assumption. This was due to the impact of a dense river network and the directions of the development of river channels [37,44]. Rongqi was located on a river way almost parallel to the latitude. Since the channel was in accordance with the direction of the reduction of the elevation, flow speed in this channel slowed down. Meanwhile, the region where Rongqi was located had a higher density of river networks than other spots in PRD, and the complexity of river channels further weakened the impact of nonstationary high tidal levels, whose effect decreased from the outlets to the incomes. Based on Xs in Zhuyin and Sanshanjiao, it was concluded that the estimated water level in Sanshanjiao was more affected by non-stationarity in the upstream flood flow. Somehow, their distance to the delta outlets was almost the same, i.e., 25 km for Zhuyin and 24 km for Sanshanjiao. As was mentioned in Section 3.1.1, there was an increase in the upstream flow in the North River sub-delta with the rising input flood flow in the West River. Meanwhile, Sanshanjiao was closer to the income station than Zhuyin. The impact of nonstationary inflow floods was thus maintained. Likewise, although the distances from Lanshi and Ma'an to their outlets were both 27 km, Lanshi had an X much higher than that of Ma'an. Since Ma'an was located at the downstream of a node where the flood from both Makou and Sanshui joined, the variance in the upstream flood change exerted little impact. The distance between the station and the income station was the main cause. As for the flood level in stations lying between 22-45 km from delta outlets, their response to either nonstationary inflow floods did not increase when they became farther away from the delta outlets. Take Zhuzhou, for example; though it is located 39 km away from the outlet, its X was the same as that of Banshawei, which was 33 km away from the outlet. Since the impact of extreme tides decreased when the distance between the station and delta outlets increased, X in Zhuzhou was supposed to be higher than that in Banshawei. On the other hand, however, Zhuzhou was much farther away from the entrance of PRD than Banshawei ( Figure 8). Therefore, the impact of upstream change decreased to a greater extent in Zhuzhou than in Banshawei. Moreover, the outlet of Zhuzhou had a more severe upward tendency than Banshawei (Table 1). Thus, the non-stationarity of the high tidal level in the outlet of Zhuzhou exerted a more severe impact on the flood level in Zhuzhou.
Although Rongqi was located 2 km farther from the delta outlet than Xiaolan [44], X in Rongqi was the same as in Xiaolan (Figure 8). Simply according to the distance between stations and the delta outlets, X in Rongqi was supposed to be higher than in Xiaolan. However, the estimated X was inconsistent with our assumption. This was due to the impact of a dense river network and the directions of the development of river channels [37,44]. Rongqi was located on a river way almost parallel to the latitude. Since the channel was in accordance with the direction of the reduction of the elevation, flow speed in this channel slowed down. Meanwhile, the region where Rongqi was located had a higher density of river networks than other spots in PRD, and the complexity of river channels further weakened the impact of nonstationary high tidal levels, whose effect decreased from the outlets to the incomes. Based on Xs in Zhuyin and Sanshanjiao, it was concluded that the estimated water level in Sanshanjiao was more affected by non-stationarity in the upstream flood flow. Somehow, their distance to the delta outlets was almost the same, i.e., 25 km for Zhuyin and 24 km for Sanshanjiao. As was mentioned in Section 3.1.1, there was an increase in the upstream flow in the North River sub-delta with the rising input flood flow in the West River. Meanwhile, Sanshanjiao was closer to the income station than Zhuyin. The impact of nonstationary inflow floods was thus maintained. Likewise, although the distances from Lanshi and Ma'an to their outlets were both 27 km, Lanshi had an X much higher than that of Ma'an. Since Ma'an was located at the downstream of a node where the flood from both Makou and Sanshui joined, the variance in the upstream flood change exerted little impact. The distance between the station and the income station was the main cause.

Conclusions
Here, the response of flood level in a dense river network to nonstationary inflow floods and high tidal level was studied. Pearl River Delta (PRD) in South China was selected as the study area. Non-stationarity in inflow floods and tidal level was first evaluated by statistical methods such as that of Mann-Kendall, and the design value was further calculated by the time-varying moment method. The flood level estimated with nonstationary input was calculated under six scenarios in a one-dimensional hydrodynamic network model. The flood level estimated with nonstationary inputs was compared with that with stationary inputs, and the response of the flood level to non-stationarity was studied. The methodology utilized in this research for estimating the flood level response to non-stationarity can be used as a reference to investigate the design flood level change in other deltas with a dense river network. It was found that the flood level response was highly affected by nonstationary inflow floods and the high tidal level induced by environmental changes.
In PRD, not considering the non-stationarity in incoming flood flow would have resulted in an 18% underestimation of the 100-year flood level, and an underestimation up to 14% (0.37 m) in the 100-year tidal levels among the multiple outlets. The underestimation of incoming flood flow or downstream high tidal level would have further caused underestimation of flood level along river channels. This situation will not only challenge the existing standard for the local flood protection but also lead to a higher degree of loss in both economy and life. For deltas that have endured intense environmental changes just like PRD, the design flood level needs to be updated by taking the non-stationarity of hydrological inputs into account.
In this paper, the spatial pattern of the flood level response was discovered. The flood level in places located 25-42 km away from the outlets in the West River sub-delta responded the most to the non-stationarity in both inflow floods and high tidal level. Meanwhile, in the North River sub-delta, places lying 15-33 km away from delta outlets had the minimum response to the nonstationary inflow floods and the high tidal level. Due to the change in the geometries of upstream river channel, the North River sub-delta had a greater distribution ratio of inflow floods. This resulted in a greater increase in the flood level in the North River sub-delta, with a rise in the nonstationary flood and tidal level of a smaller return period. The site in the West River sub-delta situated 60 km away from the outlet was an anomalous spot since the at-site flood level merely responded to the non-stationarity.
Based on the leading factor that the flood level responded to, PRD was divided into three regions. In general, regions within 22 km of delta outlets responded the most to the increasing tidal level induced by sea level rise, and spots that were located more than 45 km away from outlets were more affected by increasing inflow floods than by increasing high tidal levels. Stations that were located 22-45 km away from the delta outlets were influenced by both the sea level rise and the inflow increase. As for these stations, distance to the income stations, river network density, and the directions of river channels could all affect the extent of flood level response to the mutual impact of the nonstationary inflow floods and high tidal level. How these factors affect the water levels and their extents at different stations remains to be studied further. Meanwhile, the boundary of the incoming flood impact or that of the tidal level in other deltas can be detected in the same way.