The Coupling Effect of Flood Discharge and Storm Surge on Extreme Flood Stages: A Case Study in the Pearl River Delta, South China

The low-lying Pearl River Delta in South China is subject to severe flood threats due to watershed floods, sea level rise, and storm surges. It is still unknown to what extent and how far inland storm surges and sea level rise impact the extreme flood stages. This study investigated the coupling effect of flood discharge and storm surge on the extreme flood stages in the Pearl River Delta by using on site observations and simulations generated by the Hydrologic Engineering Center-River Analysis System model. The results show that flood discharges dominated the flood stages in the middle and upper Pearl River Delta, while the storm surges had maximum impact near the river mouth. The storm surges and flood stages showed a significant increase after 2002 in the Hengmen waterway. The design flood stages for the post-2002 period were 0.23–0.89 m higher than the pre-2002 ones at Hengmen at the six return periods from 5 to 200 years examined in this study. Their difference declined toward the upper waterway and reduced to zero about 23 km away from the Hengmen outlet. The coincidence of extreme flood discharges and storm surges further escalates the extreme flood stages in the lower 30 km of estuarine waterways. Our results quantify the severe threats due to sea level rise and intensified storm surges in the lower Pearl River Delta, and are significant for urban planning and designing and managing flood control facilities in the Pearl River Delta and in other coastal fluvial deltas.


Introduction
The low-lying Pearl River Delta (PRD) in South China is subject to severe flood threats due to flood discharge and the backwater effect of spring high tides and storm surges caused by tropical cyclones (Zhang et al. 2017;Chen 2019). Flood risk in the PRD has been exacerbated due to rapid urbanization and anthropogenic alterations ). Complex levees have been built to resist the flood water and protect the low-lying parts of the delta (Yin et al. 2020). Behind these levee systems, more and more towns, commercial centers, and manufacturing enterprises have been developed during the past four decades (Wang et al. 2020). The levee systems shrank the natural river floodplain and reduced the floodwater storage and discharge rate, leading to higher flood stages in the waterways. Meanwhile, dam construction in the upstream and sand dredging to improve navigation in the main waterways of the river often resulted in riverbed downcutting and a larger backwater effect during high tide and storm surges. Although both riverbed downcutting and the levee system help to mitigate the impact of more frequent floods and storm surges in the low-lying delta, it is still unknown how efficiently they can work during extreme cases such as the coincidence of severe but infrequent floods and storm surges.
The river networks in the Pearl River Delta have extremely complex hydrodynamic characteristics that are influenced by nonlinear interactions between river and tide dynamics (Cai et al. 2018). In the flood season when river flow dominates over the tide, the average water levels are mainly controlled by the freshwater discharge, especially in the middle and upper parts of the PRD. The large river discharge dramatically dampens down the tidal amplitude and flood tide duration by the risen water level slope, resulting in a larger tidal asymmetry, the periodic difference between the falling and rising tidal periods (Cai et al. 2014a(Cai et al. , 2014bJi et al. 2020). The duration of the rising tide is usually shorter than the falling tide, and such difference (asymmetry) becomes larger in the middle and upper delta due to the nonlinear friction effect of freshwater discharge . The tidal asymmetry is one of the critical factors affecting sediment transport and the formation of shoals and tidal flats (Brown and Davies 2010).
Natural variations and cumulative, intense human activities have engendered significant changes in the hydrodynamics of the PRD (Zhang et al. 2010). There was no evident change in annual water discharge in the Pearl River from 1957 to 2004, but the sediment load showed a significant decreasing trend due to dam construction in the middle and upper reaches of the Pearl Rivers (Zhang et al. 2008). Inconsistent with freshwater discharge variations, the annual maximum water level exhibited a decreasing trend in the upper PRD but an increasing trend in the middle and lower river reaches from the 1950s to 2005 (Zhang, Yan et al. 2009). Meanwhile, both tidal range and flood tidal duration revealed an increasing trend in the upper delta from the 1950s to 2005. The turning point occurred in the 1990s, coinciding with the downcutting of the riverbed due to the long-term and large-scale sand excavation prior to 1990 (Luo et al. 2007). The lowered riverbed and the rising sea level altered the hydrodynamics in the Pearl River Delta by increasing the water depth and reducing friction, leading to a lower flood stage and larger tidal prism, range, and duration in the upper delta (Zhang et al. 2010). In contrast, the larger tidal prism and upstream movement of the tidal limit obstructed flood discharge and escalated the flood stage in the middle and lower PRD (Zhang et al. 2017;Liu et al. 2019). Besides tide and flood discharges, storm surges caused by tropical cyclones often delivered an abrupt upswing of flood stages and posed a severe threat to the lower PRD (Li and Li 2013;Chen 2019).
Accurate determination of the flood stage profile along the main waterways is crucial for designing and supervising flood control infrastructures in the Pearl River Delta. In 2002, the Department of Water Resources (DWR) in Guangdong Province issued a suite of design flood stages at eight frequencies (return periods) of 20% (5a), 10% (10a), 5% (20a), 3.3% (30a), 2% (50a), 1% (100a), 0.5% (200a), and 0.33% (300a) along the river networks in the PRD (DWR 2002). These design flood stages are the legal guide for local and regional urban development and have played an important role in practice. Compared to the design flood stages issued by DWR in 1982, the flood stages at Makou Station in the upper PRD decreased by 0.7-1.0 m for different frequencies, while the flood stages at Ronggui and Xiaolan Stations in the middle PRD rose about 0.5 m. Zhang et al. (2017) identified similar alteration patterns of extreme water levels in the PRD by flood frequency analysis. Nonetheless, the DWR's design flood stages of 2002 in the lower PRD had been easily exceeded recently by the annual maximum water levels during storm surges. For example, the design flood stage was only 3.59 m for the largest return period of 300 years, while the recorded maximum water levels at the Hengmen outlet were 3.77 m at 12:00 on 23 August 2017 for typhoon Hato and 4.03 m at 17:00 on 16 September 2018 for typhoon Mangkhut (Chen 2019). It is urgent to diagnose the reason causing the large difference between the current peak water levels and the 2002 design flood stages (DWR 2002).
Taking the Pearl River Estuary, South China, as an example, the primary objective of this study is to estimate the coupling effect of the intensified storm surges on the extreme flood stages. This offers a more practical and conservative guidance for the planning, design, and management of flood control facilities and aids urban development planning in the Great Bay Area of Guangdong-Hongkong-Macao and in other similar riverine deltas. In situ water level observations for historical representative flood discharges and storm surges were used to analyze the spatiotemporal variations of flood stages along the West River and Hengmen waterways in the Pearl River Delta. A one-dimension hydrodynamic model, the Hydrologic Engineering Center-River Analysis System (HEC-RAS), was established to model the coupling effect of storm surges and flood discharges and to estimate the design and extreme flood stages along the Xiaolan-Hengmen waterways in the lower PRD as a case study.

Study Area and Data
There are complex and connected river channels in the Pearl River Delta. The hydrodynamics of these channels are affected by both fresh water discharge and tidal forcing. This study chose one of the simple waterways to demonstrate the coupling effect of flood and storm surge on the extreme flood stages.

Study Area
The Pearl River Delta is located in South China, and is mainly fed by the West River, North River, and East River (Fig. 1), which contribute an annual stream flow of 7,124 m 3 /s, 1,465 m 3 /s, and 719 m 3 /s respectively . The West River and North River interchange water via the Sixianjiao waterway, which is the first bifurcation of both rivers in the PRD. Makou and Sanshui are their control stations after both rivers enter the Pearl River Delta. Nanhua is the second bifurcation of the West River, whose discharge is near evenly diverted between the Xihai and Donghai waterways in normal flow conditions. The flow discharge in the Donghai waterway is further diverted into the Xiaolan, Jiya, and Ronggui waterways and finally pours into the bay of Lingding Yang via the Hengmen and Hongqili outlets (Fig. 1). The Xiaolan and Jiya waterways converge into the Hengmen waterway upstream of Donghe Station.
The target waterways in this study include Donghai, Jiya, Xiaolan, and Hengmen waterways (Fig. 1). The Donghai, Jiya, and Xiaolan waterways are in the middle and lower parts of the PRD, and are severely affected by the coupling effect of floods, tide, and storm surges . The average spring tide range in the waterways is 1.0-1.7 m and is controlled by the Pacific oceanic tide propagation and the freshwater discharge (Zhang et al. 2010;Cai et al. 2018). There are on average 2.8 landfalls of tropical cyclones annually that can generate severe storm surges in this area . The total length for the Donghai-Xiaolan-Hengmen waterways is 55 km from the upstream Nanhua Station to the Hengmen outlet. The width of Xiaolan waterway is about 190 m on average, whereas the widths of the upstream Donghai and the downstream Hengmen waterways are much wider, about 1000 m and 600 m on average, respectively. The mean depth of the waterway is about -10 m, with four deep pits (-20--30 m) in the Donghai waterway and one pit (-15 m) near the confluence of the Xiaolan and Jiya waterways in the Hengmen waterway. The hydrological stations of Nanhua and Rongqi are used for setting the upstream boundary conditions for the hydraulic modeling, while Hengmen is used to determine the downstream boundary condition. Other stations (unfilled circles) are used for model calibration and validation (Fig. 1).

Data Description
The data used in this study are summarized in Table 1. The bathymetry data of waterways were surveyed in 1999 by the Department of Water Resources (DWR) of Guangdong Province and converted into a 5 m grid. The water level data with 5-minute interval are mainly used for setting the boundary condition of model calibration and validation. The flow velocity data were measured near Hengmen using a ship-borne Acoustic Doppler Current Profiler on 18-20 June 2019. The hourly water level and flow data and the peak storm surge data are used to compute the extreme water levels for different frequencies, for statistical analysis of historical events, and for the base water levels of scenario simulations.

Design Flood Stages
The design flood stages of main waterways in the PRD were officially issued by the DWR of Guangdong Province in 2002 and were used to set the upstream boundary condition at Nanhua and Rongqi Stations, and provide a reference for the flood stages simulated in this study. They were simulated by a one dimensional numerical model using the 1999 bathymetry data in the delta. These data include 99 main waterways, 144 bifurcations/confluences, 2,425 cross sections, and a total length of 1,510 km (DWR 2002). The flood stages were simulated separately at two boundary conditions, the flood-dominated and tide-dominated, at eight frequencies (return periods) of 0.33% (300a), 0.5% (200a), 1% (100a), 2% (50a), 3.3% (30a), 5% (20a), 10% (10a), and 20% (5a). Makou and Hengmen are the upstream and downstream boundary stations. For the flood-dominated simulation, the upstream boundary condition was computed from annual maximum flow and then converted into water levels using the stage-flow curve at Makou; the downstream water levels were computed from the in situ observations of tidal levels at a near-shore station of Sanzao and then converted into the water levels at Hengmen and other outlets by the correlation method (DWR 2002).
For the tide-dominated simulation, the upstream boundary condition was set as the multiyear mean of annual maximum flow (27,697 m 3 /s) at Makou. The downstream water levels were computed from the in situ observations of tidal levels at Hengmen and other outlets. The combination of maximum water levels from both conditions was taken as the design flood stage for the waterways in the Pearl River Delta. The design flood stages at six stations (Fuzhouhe, Xiaolan, Shakou, Hengjingzha, Donghe, and Hengmen) were used to compare with those obtained in this study (Fig. 1). All elevation and water level data were adjusted to the 1985 Yellow Sea Geodatum employed in this study.

Historical Flood and Storm Surge Analysis
Water level and flow recorded for historical events are the best sources to understand the hydrodynamics along the waterways during floods and storm surges. Figure 2 illustrates the flood stages and astronomic tide from the Hengmen outlet, through Xiaolan and Nanhua, to Makou ( Fig. 1), the control station of West River, during the periods of two typical flood discharges and storm surges.
The duration of flood discharge increases with intensity. For instance, the water levels beyond 6 m at Makou lasted about 100 hours (4 days) for the 1-year (return period) flood that occurred on 6-15 June 2020 (Fig. 2a), while the water levels beyond 8 m lasted over 200 hours (9 days) for the 5-year flood event that occurred on 1-9 September 1988 (Fig. 2b). This indicates that the intense flood discharge has high probability to encounter the spring high tide, for example, the spring tide at the beginning (Lunar: 16 April 2020) and end (Lunar: 29 July 1988) of both flood events. Even with a much smaller flow (water level) at Makou, the peak water level at the Hengmen outlet during the 2020 flood was higher than that during the 1988 flood, which was caused by the backwater effect of the spring high tide. Moreover, both flood events dampened the tidal range and greatly escalated the low tide at the Hengmen outlet ( Table 2). The tide signal at Nanhua and Xiaolan almost disappeared when their water levels were higher than 4 m ( Fig. 2a, b).
For the two storm surges at return periods of 10 years on 9 September 1983 ( Fig. 2c) and 50 years on 16 September   (Fig. 2d), the peak water levels caused by the storm surges lasted only a few hours and were largely determined by the landfall time of each typhoon and the astronomical tide phases. For example, the increased water levels at Hengmen due to storm surges were 1.71 m for the 1983 typhoon Ellen and 2.86 m for the 2018 typhoon Mungkhut, while their peak water levels were 3.31 m and 3.94 m, respectively (Table 2). During the 1983 Ellen storm surge, there was a small flood discharge at Makou (water level = 3.78 m, Q = 12,300 m 3 /s), which obstructed the upstream propagation of the storm surge, with no storm surge signal showing at Makou (Fig. 2c). In contrast, there was strong storm surge of 1.75 m at Makou during the 2018 Mungkhut typhoon when the freshwater discharge was low (daily mean water level = 1.53 m, Q = 4700 m 3 /s). In addition, the water levels at Nanhua and Xiaolan were quite similar and lower than that at Hengmen during low flow and the large storm surge of the 2018 Mungkhut, while the water levels at Nanhua and Xiaolan had evident differences and were higher than water levels at Hengmen, which was due to the coupling effect of floods, spring tide, and storm surges during the 1983 Ellen storm surge. Such coupling effects on the flood stages were largely underestimated in the design flood stages issued by the DWR in 2002, and is explored in the succeeding analysis of this study.

Model Setting
The one-dimensional HEC-RAS hydrodynamic model was established to simulate the coupling effect between flood discharge and storm surge. Various boundary conditions (water levels) were designed to simulate the extreme water level profile along the study waterways based on the observed water levels during historical floods and storm surges.

Model Calibration and Validation
The HEC-RAS model was developed by the Hydrologic Engineering Center (HEC) of the U.S. Army Corps of Engineers, and has been widely used in channel hydrologic simulation, including flow and water level simulation, dam breaks, inundation simulation, and risk management (USACE 2010;Parhi 2018;Islam et al. 2019). The HEC-RAS model is based on the continuity equation in Eq. 1 and momentum equation in Eq. 2, which describe mass conservation and momentum variation respectively (USACE 2010).
where V is the mean velocity of a cross section, qz/qx is the water surface slope, Q and A T are the flow and total flow area at the midpoint of the control volume, q l is the lateral inflow per unit length, x is distance along the channel, t is time, z is water surface elevation, and S f is friction slope. In Eq. 2, three main forces are consideredpressure, gravity, and friction.
The HEC-RAS model was established along the Donghai, Jiya, Xiaolan, and Hengmen waterways using the bathymetry data, which were combined with the land elevation and levee height data (Fig. 1). The left and right banks were constrained by the levees. The cross sections were extracted with a mean interval of 200 m using the HEC-GeoRAS plug-in in ArcMap, and a smaller interval of 20 m was used in the bend and at the bifurcation and confluences. Expansion and contraction coefficients of cross sections near bridges use the default values suggested in the model user manual, that is, 0.3 and 0.1 on general cross sections, 0.5 and 0.3 on the upper and lower cross sections of a bridge (USACE 2010).
The HEC-RAS model was calibrated for the 72-hour unsteady flow during the spring tide on 18-20 August 2019 and driven by the 5-minute water levels recorded at the  (USACE 2010). The Nash-Sutcliffe efficiency coefficient (NSE) and the Mean Absolute Difference (MAD) between the simulated and observed water levels at the other 8 stations were used to optimize the model parameters (Fig. 1). The NSE and MAD were computed using Eqs. 3 and 4, respectively (Nash and Sutcliffe 1970;Pinos and Timbe 2019).
where N is number of water levels in sequence, h obs is the water level in observation sequence, h sim is the water level in simulation results, h obs is the mean water level in observation sequence, w t is the weight.
The MAD is also used as the calibration parameter for the velocity.
After calibration, the mean flow velocity from the HEC-RAS model is in good agreement with the in situ measured values at Hengmen,for example, NSE = 0.85, MAD = 0.11 m/s (Fig. 3). The HEC-RAS model was then validated using the hourly water levels recorded at Nanhua, Rongqi (upstream), and Hengmen (downstream) as boundary conditions for the unsteady flow during the storm surge on 15-18 September 2018 and the flood discharge on 8-12 June 2020. The NSE and MAD between the simulated and observed water levels at the other 8 stations were calculated to evaluate the model performance (Fig. 1).
The Mangkhut super typhoon first landed on the Philippines on 15 September 2018, then moved northeast through the South China Sea and landed again on Taishan in Guangdong Province on 16 September. It was the largest one of 2018, with a maximum wind speed of 45 m/s and air pressure of 955 hPa at landfall. The highest storm surge was 3.39 m at Sanzao near the landfall and 2.86 m at Hengmen (Table 2, Fig. 2d), which broke the tidal record at Hengmen and was equivalent to a 50-year storm surge event (COIN 2018). Fortunately, the storm surge occurred during the neap tide and did not meet spring high tide or flood although it caused severe damage and significant loss of life and property.
The subtropical high was stagnating over South China and brought heavy rains in the upper basins of the North and West Pearl River in early June 2020, forming the number one flood peak on 7-15 June during the spring tide (Lunar 16-24 April) (Fig. 2a). The peak water levels at Makou and Nanhua rose to 7.4 m and 5.8 m, respectively. Both the 2020 flood and 2018 storm surge events offered us good opportunities to validate our model performance and to set the base water levels for modeling extreme floods and storm surges.

Scenario Simulations
After the HEC-RAS model had been calibrated and validated, it was applied to simulate the extreme flood stages under various scenarios: floods, storm surges, and the coupling of floods and storm surges. Since the flood duration of extreme flood events with frequency smaller than 0.2 is usually longer than 10 days (Fig. 2b), the discharge of flood water usually encounters the spring tide. Thus, the boundary conditions of flood simulations were set during the spring tide. The in situ water levels during the storm surge on 7-11 September 1983 and the flood on 8-12 June 2020 were used as the base water levels to derive the hourly water levels after the extreme water levels were set for different frequencies at the three boundary stations (Table 3).
For the flood simulations, the upstream boundary water levels at Nanhua and Rongqi were set as the same design flood stages of DWR (2002) at six frequencies (return periods) of 20% (5a), 10% (10a), 5% (20a), 2% (50a), 1% (100a), and 0.5% (200a) ( Table 3). As demonstrated in Fig. 2b, the hourly water levels at Nanhua and Rongqi were set as stable during the 120-hour simulations. The high tide water levels at Hengmen (downstream) were set based on the high tide water level (2.53 m) of the 1-year flood on 9 June 2020, and adding the increase of water levels at the  (Table 2). In other words, the flood discharge has a much larger magnitude on increasing the water levels for low tide than for high tide. After setting the water levels of high tide and low tide, the hourly water levels during the 120-hour simulation period were generated by a two-step conversion proportional to the in situ recorded water levels on 8-12 June 2020 (Fig. 2a).
For the storm surge simulations, the extreme water levels at Hengmen were first calculated from the annual maximum water levels from 1954 to 2019 using the Gumbel distributions at six frequencies (Table 3). The peak water levels recorded at Nanhua, Rongqi, and Hengmen during the 10-year storm surge on 9 September 9 1983 were used as the base water levels, and the peak water levels at the other return periods were generated by adding their difference of extreme water levels in different frequencies at Hengmen (Table 3). The hourly water levels at the three boundary stations during the 120-hour simulation period were derived from the peak water levels proportionally based on the in situ recorded water levels on 7-11 September 1983, which also represented a random combination of landfall storm surge and tide phases (Fig. 2c).
For the coupling simulations of floods and storm surges, the upstream boundary conditions at Nanhua and Rongqi were the same as the flood simulations (Table 3). At the downstream Hengmen boundary station, the high tide levels were the same as the storm surge simulation, while the low tide levels were same as the flood simulations. The hourly water levels at Hengmen during the 120-hour simulations were derived from the high tide and low tide levels proportionally against the in situ recorded water levels on 7-11 September 1983, which presents a general encounter probability between storm surges and tidal phases (Fig. 2c).
The design water levels along the waterways were the combinations of the simulated maximum water levels for the flood-dominated simulations and the storm surge dominated simulations. The former were set as the flood discharges in six frequencies at the upstream stations of Nanhua and Rongqi that encounter the 5-year storm surge at Hengmen, and the latter were set as 5-year flood discharge that encounters the storm surges in six frequencies ( Table 3). The extreme flood stages were simulated when the same frequency floods and storm surges occur coincidentally.

Results
In this section, we present the results for model calibration and validation during a spring tide, flood and storm surge, and the scenario simulations.   (0.06 m) to 0.94 (0.14 m) at stations along the Jiya waterway (Fig. 4). Figures 5 and 6 demonstrate the validation results for the flood on 8-12 June 2020 and the storm surge on 15-18 September 2018. The HEC-RAS model effectively simulates the water level variations during the flood discharge and storm surge. For the 2020 flood, the simulated water levels were generally lower than the observations in the Xiaolan waterway, but agreed well with the observations in the Jiya and Hengmen waterway. The NSEs (MAD) vary from 0.78 (0.10 m) to 0.81 (0.22 m) at stations in the Xiaolan waterway, while they vary from 0.6 (0.06 m) to 0.96 (0.28 m) at stations in the Jiya and Hengmen waterway. For the 2018 Mungkhut storm surge, the simulated water levels overall agreed well with the observations in phase and magnitude for all stations on both waterways, except for lower peak water levels. The NSEs (MAD) vary from 0.88 (0.06 m) to 0.98 (0.18 m). This indicates that the model is well calibrated and can be used to simulate the flood discharges and storm surges in different scenarios.

Scenario Simulations
The scenario simulations include flood only, storm surge only, flood dominated, storm surge dominated, and the coincident of flood and storm surge with same frequency. The design flood stages of DWR (2002) at several stations along the Xiaolan-Hengmen waterways were illustrated as a reference. Figure 7 presents the simulated maximum water levels during floods and storm surges and the coupling (coincidental) effects of floods and storm surges in six frequencies (return periods). The flood discharges dominate the maximum water levels along the upper 40 km beginning at Nanhua, and the simulated maximum water levels were in good agreement with the design flood stages of DWR (2002). In contrast, the storm surges dominate the maximum water levels near the river mouth and affect a range of up to 20 km upstream starting at Hengmen (Fig. 7). The simulated maximum water levels during storm surges were much higher than the design flood stages of DWR (2002) at the stations of Donghe and Hengmen outlet. If the flood and storm surge occur simultaneously, the coupling effect of flood and storm surge further escalates the extreme water levels and the affecting ranges increase with the

Design Flood Stages
Although the coupling effect of flood and storm surge generates much higher peak water levels near the coastal areas, they seldom occur at the same time, especially the low frequency events. The low-frequency floods usually occur in the first flood season ranging from April to June when the cold front system and the South China Sea Monsoon could produce basin-wide storms and severe flood (Wang et al. 2019), while large-scale typhoons normally form from July to September . Following the practice of the design flood stages issued by DWR in 2002, the flood-dominated scenarios are set as the floods encountering a 5-year storm surge, and the storm surge dominated scenarios are set as the storm surges encountering a 5-year flood. The maximum water levels from both scenarios form the design flood stages for the same return period. For the flood-dominated simulations (Fig. 8a), the maximum water levels from the coupled simulation of 5-year flood against the 5-year storm surge agreed well with the design flood stages of DWR. For the other frequency floods, the maximum water levels agreed with the design flood stages of DWR within the upper 15 km waterways, but are slightly lower near the Hengmen outlet. For the storm surge-dominated simulations (Fig. 8b), the storm surges control the maximum water levels, which are much higher than the design flood stages of DWR within the lower 20 km of the waterway near the Hengmen outlet.
The combinations of maximum water levels from both flood-dominated and storm surge-dominated simulations form the design flood stages at the same frequency along the Dohai-Xiaolan-Hengmen waterways. The new design flood stages generally agree with those of DWR in 2002 within the upper 20 km waterways near Nanhua, but are much higher along the lower 30 km waterways, especially  (Fig. 8c). Moreover, the extreme water levels were generated when a flood discharge encounters a storm surge at the same frequency during spring tide (Fig. 8d), although such cases seldom occur and the combined frequency is much smaller.

Discussion
Our results reveal that the flood stages along the waterways near the Hengmen outlet were up to 0.89 m higher than the design flood stages issued by DWR in 2002 (Fig. 8c, Table 4). The increasing flood stages are consistent with the results of Zhang et al. (2017) that all stations in the lower PRD had significant increments of flood stages by 0.35, 0.59, and 1.06 m for return periods of 10, 20, and 50 years respectively in a comparison between pre-1980 and post-1980 periods. Compared to the design flood stages issued by DWR in 1982, the flood stages of 2002 at Xiaolan in the middle PRD rose about 0.5 m (DWR 2002). This suggests that the extreme flood stages along the estuary waterways in the lower PRD have been rising during the past half century and the area faces severe flood threat (Zhang et al. 2017).
The extreme flood stages in the estuary waterways are subject to the forcing of river discharges and tides, especially the storm surges, although their time scales differ widely as tidal waves vary hourly, while high flood discharges can last days or weeks (Leonardi et al. 2015;Ji and Zhang 2019). The continuous increments of flood stages are likely caused by the intensified human activities in the delta and the global sea level rise.

Sea Level Rise
Sea level rise is one of the most direct factors that escalate the extreme flood stages in the estuary waterways of the Pearl River Delta (Silva et al. 2020). From 1906 to 2014, the mean high tide in the Lingding Yang Bay had increased by 0.40 m from 1. 58 m in 190658 m in , to 1.62 m in 193658 m in , 1.66 m in 195458 m in , 1.82 m in 197158 m in , and 1.98 m in 201458 m in (Gong et al. 2019. The mean sea level rise rate is 3.6 mm/year in the coast of Guangdong Province during 1980-2019, plus the monthly anomaly of 156 mm, for example, in May 2019 above the annual mean (DNRG 2020). This indicates that sea level rise could explain up to 0.3 m (about 30%) of the increment of extreme flood stages in the estuary waterways since 1980. If the current rising rate continues, the mean sea level in the estuary of the Pearl River is projected to rise with an extra 108 mm by 2050 and 288 mm by 2100. Moreover, the rising sea level reduces the riverbed/seafloor and lateral friction and increases the propagation of the tidal and storm surge waves along the estuary waterways, resulting in a nonlinear or higher rise of extreme flood stages greater than suggested by the recorded sea level rise (Zhang et al. 2010).

Downcutting of Riverbed and Channel Dredge
From 1986-2003, over 8.7910 8 m 3 of sand was excavated from the PRD, resulting in an average downcutting of the depths of 0.59-1.73 m in the main channels of the West River (Luo et al. 2007), which in turn altered the variations of the spatiotemporal distributions of extreme flood stages in the PRD (Zhang et al. 2010). For instance, the flow ratio at Sansui was 21.3% pre-1991 and 25.6% post-1991 (Zhang, Yan et al. 2009). There was a dramatic downshift of the flow-stage rating curves in the upper delta between the 2008 June flood and the 1998 June flood . Riverbed dowcutting caused longer flood tide duration, larger tidal prism, and greater range especially in the upper PRD (Zhang et al. 2010), which resulted in more total flood discharge and anomalous high water levels in the middle and lower delta ). In addition to dredging activity in the river channels, largescale channel dredging occurred in the Lingding Yang Bay, where water areas deeper than 10 m increased by 36 km 2 during 1988-2008 (Wu et al. 2014). Between 2012 and 2013, the local deepening of channel depth was about 5 m to accommodate ever larger ships in Guangzhou Harbor, far beyond the magnitude of natural evolution (Wu et al. 2016). The abnormal high storm surges and water levels in the lower PRD during typhoon Hato in 2017 and 2018's typhoon Mangkhut were an integrated effect of these human alterations and sea level rise (Chen 2019). Nevertheless, it is still challenging to quantify the specific amount of water level increase due to the downcutting of riverbed and channel deepening in the Lingding Yang Bay by sand excavation, navigation dredging, or reduced sediment load resulting from reservoir construction.

Land Reclamation
The land area of the delta grew by 1,160 km 2 from 1850 to 2015, while the water area and volume decreased by 35% and 39%, respectively (Wu et al. 2018). Further dramatic land use changes occurred during 1988-2008, when nearly 200 km 2 of the land reclaimed from the intertidal and shallow depths was mostly used for industrial development (Wu et al. 2014). At the same time, the mean water depth increased by 15% due to reduction in shallow water and the large-scale navigation channel dredging in the Lingding Yang Bay (Wu et al. 2014). The eastern outlets of the Pearl River, for example, extended seaward 8-10 km during the 1980s-2005 (Zhang et al. 2010). The national tidal station of Hengmen outlet was located on the coastline of the Lingding Yang Bay, but now is about 8 km away from the bay coastalline due to the land reclamation of Cuiheng District in Zhongshan and the Wanqingshan in Guangzhou, and became an inner part of the Hengmen waterways (Fig. 1). The land reclamation extended the flood discharge path and in turn elevated the flood stages in the middle and lower part of waterways of the Pearl River Delta.

Variations of Stream Flow and Storm Surges
There were no significant changes in annual precipitation from 1951 to 2005 throughout the Pearl River Basin (Zhang, Xu, Becker et al. 2009), but the annual maximum stream flow had a slightly upward trend at the Gaoyao Station, which is a reference station of flow and flood for the West River (Zhang, Yan et al. 2009). The increased maximum stream flow did not raise the maximum water level at Makou in the upper PRD due to riverbed downcutting, but did escalate the flood stages in the waterways of the middle and lower delta (Zhang, Yan et al. 2009). However, the contribution of the increased maximum flood discharge to the increment of the extreme water levels at Hengmen outlet was near negligible. Moreover, the maximum stream flow is expected to decrease in the near future after several large reservoirs are in operation in the upper West River.
Generally, flood discharge lasts days and weeks and dominates the extreme water level profile along the middle and upper waterways in the PRD, while the storm surges dominate the extreme water levels in the estuary waterways near the outlet, although these events last only a few hours (Wu et al. 2020). The maximum wind speed of tropical cyclones showed an increased trend in the West North Pacific over the period 1981-2006, although the occurrence frequency showed a decreasing trend (Elsner et al. 2008). Table 5 demonstrates that the extreme storm surges at Hengmen rose by 0.28 to 0.87 m, and the extreme water levels rose by 0.17 to 0.63 m between pre-2002 and post-2002 periods in the six return periods from 5 years to 200 years. Compared to the boundary conditions of tide-dominated simulations of DWR (2002) at Hengmen, the extreme water levels applied in this study were 0.23 to 0.89 m higher in the six return periods (Table 4), which are the primary cause that contributes to the differences of design flood stages between this study and the 2002 design flood stages (DWR 2002). Their difference declined toward the upper waterway and was close to zero around the Hengjingzha Station, which is 22.8 km away from the Hengmen outlet (Table 4).
The extreme flood stages and the design flood stages are assumed to be the same at Hengmen. Their difference increased toward inner waterway and reached a peak near the Donghe Station, which is 4.5 km away from the Hengmen outlet, and then reduced toward the upper waterways beyond Hengjingzha (Fig. 7, Table 4). Nevertheless, extreme floods in the entire river watershed and the storm surges from the off-coast area are separate events and seldom occur at the same time. In practice, the design flood stages were the maximum combinations of water levels simulated by the extreme flood discharge and the coastal storm surges (DWR 2002). The extreme water levels are illustrated as comparison and guidance to handle the uncertainty of design flood stages and grace heights in designing flood control facilities.

Conclusion
This study analyzed the spatiotemporal variations of historical flood discharges and coastal storm surges from in situ water level observations along the West River and Hengmen waterways in the Pearl River Delta. The flood discharge events last over 10 days and dominate the extreme flood stages of waterways in the middle and upper PRD, and dramatically escalate the low tide levels and reduce the tidal range at the Hengmen outlet. In contrast, the storm surges only last a few hours but dominate the extreme flood stages near the estuary waterway.
Considering the abnormal high water levels caused by storm surges in the lower PRD in recent years, we established the HEC-RAS hydraulic model to simulate the coupling effect of extreme flood discharges and storm surges along the Donghai, Xiaolan/Jiya, and Hengmen waterways during spring tide. Generally, the simulated water levels agreed well with the in situ observations at the eight stations for calibration during the spring tide in August 2019 and for validation during the 1-year flood in June 2020 and the 50-year storm surge in September 2018.
As for the scenario simulations, the storm surges and flood stages showed a significant increase between pre-2002 and post-2002 periods in the Hengmen waterway. The coupling effect of flood discharges and storm surges further escalates the extreme flood stages in the lower 30 kms of waterways although they are separate events and seldom occur at the same time. Our results quantify the severe threats due to sea level rise and intensified storm surges in the lower delta, and are significant for urban planning and for managing and designing flood control facilities in the PRD and in other similar fluvial deltas.