Nonlinear Flood Responses to Tide Level and Land Cover Changes in Small Watersheds

: Regarding global warming, the threat of ﬂooding is projected to increase due to the change in intensity and frequency of single drivers and ampliﬁcation caused by multi-driver interactions. This interaction becomes more complicated in developing regions with rapidly changing land cover. As a result, demands on ﬂood risk management are rising especially in small watersheds, which are more vulnerable to driver disturbances compared with large watersheds. Existing studies focused on large watersheds rather than small watersheds. However, the ﬁndings derived from large-scale analysis cannot be transferred to small watersheds directly. This research investigated the ﬂood responses in the Yonghe River Watershed (YRW) (63.8 km 2 ) in Guangzhou, China, considering the impact of land cover change. The YRW experienced a disastrous compound ﬂood on 22 May 2020. A hydrodynamic model integrating the Hydrologic Engineering Center’s Hydrologic Modeling System and River Analysis System (HEC-HMS and HEC-RAS, respectively) was established and calibrated using the inundation depths observed during the ﬂood. Model analysis using multiple scenarios showed that the watershed is river-dominated, and ﬂood responses to the three factors are nonlinear but with different increasing rates. The response curves for tide levels and land cover changes increase faster at high values, whereas the rainfall intensity curves vary slightly. These ﬁndings highlight the importance of integrating tidal impacts into ﬂood risk management, even in river-dominated coastal watersheds. The study further recommends that in small watersheds, 50% imperviousness is an indicator of the urgent demand for ﬂood risk management measures.


Introduction
Globally, floods are among the most hazardous disasters, accounting for 44% of total events during 2000-2019, which affected 1.6 billion people [1].Moreover, due to global warming, the frequency and intensity of heavy precipitation events have increased since the 1950s over most land areas, and sea levels have also shown an evident rising trend [2].Consequently, the adverse impacts of the interactions between inland and oceanic processes are greatly enhanced.Compound flooding is driven by multiple mechanisms, including intense rainfall, riverine flow, astronomic tides, storm surges, and waves, either simultaneously or in close succession.Their occurrence probability increases in low-lying coastal regions, where the population and assets are crowded.
Compound-flood analyses are typically conducted using statistical or hydrodynamic models.The abundance of observational data enables statistical models to reveal the dependence among extreme rainfall, riverine flow, and storm surges and estimate the probability of co-occurrence [3][4][5][6][7][8][9].A typical example could be a tropical cyclone event that delivers storm surges and intense precipitation together [10,11].Nasr et al. [12] found that along the United States coastline, the highest dependence exists for surge waves, followed by surge precipitation, surge discharge, wave precipitation, and wave discharge.Ghanbari et al. [13] revealed that the northeast Atlantic and western Gulf Coasts experience the Land 2023, 12, 1743 2 of 17 highest compound flood probability under current conditions of climate and anthropogenic activities.In Europe, the Mediterranean coasts currently experience the highest compound flood probability from precipitation and storm surges, and a high probability is projected along the northern parts of the European coast [14,15].In China, the dependence on storm surges and heavy precipitation is higher at southeastern tide gauges (latitude < 30 • N) than at northern tide gauges, while several sites have significant dependence on the tropical cyclone season [16].The dependence of different compound flood drivers can provide macrotrends at continental or national scales; however, they are still insufficient for flood risk management in a specific estuary or watershed, which relies on high-resolution water level data.The GIS-based bathtub model simply overlays the water level from distinct sources and has been widely applied in large-scale studies owing to its easy-to-use features.However, this tends to cause a bias in flood risk because the interaction between different mechanisms is ignored [17].Compared to deterministic bathtub models, hydrodynamic models can produce more accurate results, such as water depth and flood extent [18][19][20], and assist in the identification and mapping of the dominant forces, which is critical for policymaking in flood risk management for coastal regions.
Globally, the transition of flood mechanisms from ocean-dominated to river-dominated varies, depending not only on the intensity of individual driving forces but also on other controlled variables, including riverbed elevation, bottom roughness, estuary geomorphology, wind conditions, and phase lag [21][22][23][24][25]. Furthermore, a specific estuary region can be delineated into three types of flooding zones: hydrological zone (river-dominated), tidal zone (ocean-dominated), and transition zone [19].Estuaries act as frequency filters, where low-frequency waves decay at a lower rate and propagate more inland than high-frequency waves [26].For instance, astronomical tidal waves from the Yangtze River Estuary propagate 700 km inland [26], whereas backwater zones can extend up to 500 km in the Mississippi River [22].Under a nonstationary climate, the return periods of storm tides reduce significantly [25], and thus the inland propagation of low-frequency tides is likely to be more frequent, imposing new challenges for river-dominated regions.
Unlike large watersheds, the impact of compound flooding on small watersheds with a spatial coverage of less than hundreds of square kilometers has received limited attention.Since small watersheds are typically independent hydrological units in estuarine regions, they are influenced by ocean processes.However, this interaction has not been adequately represented in large-scale analyses owing to their coarse resolution.Compared with the main river watersheds, small watersheds have a shorter response time at the hourly scale and much smaller peak discharges, suggesting their high vulnerability to backwater effects from ocean processes.Therefore, knowledge learned from modeling at the estuary scale should be transferred to small watersheds with caution.In addition, urbanization is still a rapid process in developing regions.Changes in urban hydrology, such as an increase in peak discharge and an advance in peak time, can influence the intensity of compound flooding.Therefore, the impact of land cover changes on compound flooding and its interaction with other driving forces should be considered.
This study aimed to explore flood responses in small developing coastal watersheds where demands on flood risk management are rising rapidly.A watershed in Guangzhou China that experienced an extreme flood event on 22 May 2020 was analyzed.A hydrodynamic model was first established and then calibrated using flood-stage investigation data.The model was further applied to investigating flood responses to rainfall, tide, and land cover changes.Finally, the model uncertainties and the possible transfer of learned knowledge in this case to other small watersheds were discussed.

Study Area
The average annual precipitation in Guangzhou is approximately 1800 mm, which is unevenly distributed over time.More than 80% of the annual precipitation occurs during the rainy season from April to September [27].Guangzhou is also influenced by Land 2023, 12,1743 semidiurnal mixed tides in Lingding Bay.Regarding land cover changes, the built-up area in Guangzhou has increased rapidly, from 297.5 km 2 in 2000 to 1380.6 km 2 in 2022.
The Yonghe River Watershed (YRW) is located to the east of the central urban area of Guangzhou, China (Figure 1).The YRW covers two metro stations, Guanhu and Xinsha, that were completely flooded on 22 May 2020 (Figure 2).Moreover, some communities and villages in the downstream area were severely flooded, to depths of >1.5 m.
unevenly distributed over time.More than 80% of the annual precipitation occurs during the rainy season from April to September [27].Guangzhou is also influenced by semidiurnal mixed tides in Lingding Bay.Regarding land cover changes, the built-up area in Guangzhou has increased rapidly, from 297.5 km 2 in 2000 to 1380.6 km 2 in 2022.
The Yonghe River Watershed (YRW) is located to the east of the central urban area of Guangzhou, China (Figure 1).The YRW covers two metro stations, Guanhu and Xinsha, that were completely flooded on 22 May 2020 (Figure 2).Moreover, some communities and villages in the downstream area were severely flooded, to depths of >1.5 m.
The watershed covers a relatively small area of 63.8 km 2 .The Yonghe River is 21.1 km long and originates from the hilly area in Huangpu District, flows through Zengcheng District, and ends at the East River, one of the main tributaries of the Pearl River.The percentage of impervious surfaces in the YRW increases from 13.9% in 2000 to 54.4% in 2020.The average slope of YRW is 5.71 degree and quartiles are 0.46, 2.06, 7.56, and 58.46 degrees.

Topography
Accurate topography is critical in flood modeling because low-quality DEM data may lead to large uncertainties in the results, such as flooding area and inundation depth.The topographic data were acquired using airplane-based LiDAR, and the spatial resolution of the resulting DEM data was 2 m.The watershed covers a relatively small area of 63.8 km 2 .The Yonghe River is 21.1 km long and originates from the hilly area in Huangpu District, flows through Zengcheng District, and ends at the East River, one of the main tributaries of the Pearl River.The percentage of impervious surfaces in the YRW increases from 13.9% in 2000 to 54.4% in 2020.The average slope of YRW is 5.71 degree and quartiles are 0.46, 2.06, 7.56, and 58.46 degrees.

Model Data 2.2.1. Topography
Accurate topography is critical in flood modeling because low-quality DEM data may lead to large uncertainties in the results, such as flooding area and inundation depth.The topographic data were acquired using airplane-based LiDAR, and the spatial resolution of the resulting DEM data was 2 m.

Land Cover
Land cover in 2020 was extracted from the Gaofen 1 Satellite image (multi-spectrum, 8 m spatial resolution) using an object-oriented classification approach.The classification results included water, vegetation, bare earth, and construction land.The overall accuracy was 97.7%, and the Kappa coefficient was 0.96.

Rainfall and Tide Levels
Rainfall data were collected from the official website of the Water Resources Department of Guangdong Province (http://210.76.80.76:9001/Map/Map.aspx?id=2, accessed on 6 October 2020).In 2020, the storm event occurred from 8 p.m. on 21 May to 7 a.m. on 22 May and peaked at midnight (Figure 3).The entire watershed was divided into subareas based on the locations of the weather stations using the Thiessen polygon method.Rainfall depths at the five stations are listed in Table 1.

Topography
Accurate topography is critical in flood modeling because low-quality DEM data may lead to large uncertainties in the results, such as flooding area and inundation depth.The topographic data were acquired using airplane-based LiDAR, and the spatial resolution of the resulting DEM data was 2 m.

2.2.2.. Land Cover
Land cover in 2020 was extracted from the Gaofen 1 Satellite image (multi-spectrum, 8 m spatial resolution) using an object-oriented classification approach.The classification results included water, vegetation, bare earth, and construction land.The overall accuracy was 97.7%, and the Kappa coefficient was 0.96.

Rainfall and Tide Levels
Rainfall data were collected from the official website of the Water Resources Department of Guangdong Province (http://210.76.80.76:9001/Map/Map.aspx?id=2, accessed on 6 October 2020).In 2020, the storm event occurred from 8 p.m. on 21 May to 7 a.m. on 22 May and peaked at midnight (Figure 3).The entire watershed was divided into subareas based on the locations of the weather stations using the Thiessen polygon method.Rainfall depths at the five stations are listed in Table 1.At approximately the same time, the tidal level peaked (Figure 4).Huangpu was the nearest tide station to the study area, and tidal data from Huangpu during the storm event were collected from the National Marine Data and Information Service's website (https://www.cnss.com.cn/tide/,accessed on 6 October 2020).At approximately the same time, the tidal level peaked (Figure 4).Huangpu was the nearest tide station to the study area, and tidal data from Huangpu during the storm event were collected from the National Marine Data and Information Service's website (https://www.cnss.com.cn/tide/,accessed on 6 October 2020).

Model Setup
The flood process was modeled by coupling the Hydrologic Engineering Center's Hydrologic Modeling System (HEC-HMS) and the River Analysis System (HEC-RAS).The study area was divided into multiple catchments, each of which had an outlet for the main channel to the Yonghe River.The HEC-HMS simulated the hydrological processes in each catchment and provided the flow rate at the outlets, which served as the boundary conditions for the HEC-RAS.To interact with ocean processes, the boundary conditions also included the tidal level.Finally, HEC-RAS simulated the overflow of floodwater along the main channel and calculated the flooded area and inundation depth using the simulated flood stage and high-quality terrain data.

HEC-HMS
Runoff generation was calculated using the Soil Conservation Service Curve Number (SCS-CN) method, in which only one parameter, CN, is required.For one catchment, the CN value equals the mean CN value of different land cover types weighted by area (Equation (1)).

Model Setup
The flood process was modeled by coupling the Hydrologic Engineering Center's Hydrologic Modeling System (HEC-HMS) and the River Analysis System (HEC-RAS).The study area was divided into multiple catchments, each of which had an outlet for the main channel to the Yonghe River.The HEC-HMS simulated the hydrological processes in each catchment and provided the flow rate at the outlets, which served as the boundary conditions for the HEC-RAS.To interact with ocean processes, the boundary conditions also included the tidal level.Finally, HEC-RAS simulated the overflow of floodwater along the main channel and calculated the flooded area and inundation depth using the simulated flood stage and high-quality terrain data.

HEC-HMS
Runoff generation was calculated using the Soil Conservation Service Curve Number (SCS-CN) method, in which only one parameter, CN, is required.For one catchment, the CN value equals the mean CN value of different land cover types weighted by area (Equation ( 1)).
where CN is curve number of the catchment, n is the number of land cover types in the catchment, and A i and CN i represent the area and curve, respectively, for the ith land cover type.Runoff volumes generated in the catchments were routed to the catchment outlet using the SCS Unit Hydrograph approach, which has a principal parameter, lag time.The lag time can be defined as follows: Lag = l 0.8 (S+1) 0.7 517.217Sl p 0.5 (2) where Lag is the time (h), l is the flow length (m) of the longest flow path, Sl p indicates the average slope (%), and S represents the maximum potential retention (mm), equal to (1000/CN) − 10, where CN represents the number of SCS curves.The flow length and average slope can be obtained using the HEC-GeoHMS 10.7 software, which is a geospatial hydrology toolkit in ArcMap.In the HEC-HMS model, the base flow was set as zero.The base flow in such a small watershed is limited and its impact on rainfall event simulations was represented by model calibration implicitly.

HEC-RAS
The geometric model of the Yonghe River was established based on a 2 m LiDAR DEM using HEC-GeoRAS, which is also a toolkit in ArcMap.Cross sections were defined at 100 m intervals along the main channel and cross-section lines were extended outwardly until watershed boundary.Floodplains lower than overflow water level are considered to be flooded.The energy loss comprised contraction or expansion and friction losses.The cross sections of the Yonghe River changed gradually, especially in the downstream areas.Thus, contraction and expansion coefficients were set to 0.1 and 0.3, respectively.The evaluation of friction losses required Manning's roughness coefficients for the main channel and floodplains.Roughness coefficients were defined separately for the main channel, left floodplain, and right floodplain.The floodplains were empirically divided into upstream and downstream parts based on significant differences in land cover.The roughness values were determined during model calibration using the 22 May flood event.

Model Calibration 2.4.1. Inundation Depths
A field investigation was conducted on 25 May 2020 using surveying equipment, including the global navigation satellite system's real-time kinematic (GNSS RTK) and total station equipment.GNSS RTK can determine the three-dimensional coordinates of a single point with centimeter precision.The key to accurate positioning is the real-time corrections that the equipment requires from a continuously operating reference station (CORS) system.When satellites are inadequately available due to blockages by buildings and canopies, the GNSS RTK cannot function well; then, the alternative is "GNSS RTK + total station".The GNSS's receiver was first moved to a position where accurate positioning was to be achieved.A total station was then used to survey the height difference between the flood mark and the GNSS receiver.Finally, the elevation of the floodmark above the 1985 National Elevation Datum was calculated using the height difference and GNSS RTK observations.In this manner, three sites were investigated, namely Guanhu Station, Xinjie Village (close to Xinsha Station), and Lixin Road (Figure 1b), and their inundation depths were observed to be 0.87, 1.65, and 1.0 m, respectively.

Roughness Coefficients
The roughness coefficients used for the model calibration are listed in Table 2.The upstream floodplains exhibited a wider range of roughness coefficients; thus, more samples were collected from there.In total, 432 combinations of roughness coefficients were obtained.the observed inundation depths at the three sites considered in the field investigation, as follows: where h sim (i) and h obs (i) represent the simulated and observed inundation depths, respec- tively.The roughness coefficients with the lowest MD among all 432 combinations were applied for the subsequent scenario analyses.The calibrated model was then applied to examine the responses of rainfall intensity, tidal level, and land cover change, and to identify the dominant flood mechanism in coastal watersheds.

Scenarios Design 2.5.1. Rainfall Scenarios
Since most of the precipitation from the 22 May event was concentrated over 3 h, rainfall scenarios were designed for a duration of 3 h.According to the local intensity-durationfrequency relationships released by the Guangzhou Bureau of Water Authority [28], the rainfall depths for typical return periods are listed in Table 3.The rainfall depth for each return period was distributed over time by the Chicago curve [29], where the rising and falling limbs can be represented as: where q(t b ) and q(t a ) are the time series of rainfall intensity in the rising and falling limbs, respectively; t b and t a are the times before and after the peak rainfall intensity, respectively; r is the ratio of the peak rainfall intensity time to the total duration and was set to 0.33; A 1 is the rainfall depth (mm) for the one-year return period; F is the parameter of rainfall depth variation; T is the return period (year); b is used to adjust the rainfall duration; and n is related to rainfall attenuation.Parameters including F, b, and n were determined by calibration using historical rainfall events.The values of F, b, and n were 0.438, 11.259, and 0.750, respectively [28].

Tide Level Scenarios
The tide levels at the Huangpu Station for the five typical return periods [30] are listed in Table 4.The tidal gauge station, marked in Figure 1a, is only 20 km away from the outlet of Yonghe River.As a result, the tidal processes at the outlet of Yonghe River can be represented by processes at the tidal gauge station.The tidal levels in Table 4 were directly applied as boundary conditions to the Yonghe River.2020), which provides high-resolution (30 m) land cover products with quality control.It has an overall accuracy greater than 85% for construction land and residential areas.Land cover for 2025 and 2035 (30 m spatial resolution) was predicted using the Future Land Use Simulation Model Software (FLUS 2.4) [31].The FLUS model is an integrated model for multitype land-use scenario simulations that couples human and natural effects [32].

Perturbation Analysis
As rainfall intensity, tide level, and land cover imperviousness have different dimensions, the contribution of one influential factor can be measured by the relative change in flooding indicators, flooded area or inundation depth, divided by the relative change in the factor.First, a baseline scenario was set, i.e., 5-year design rainfall, 2-year design tide, and land cover in 2020.The model was then run by varying one factor and controlling for the other two.The relative change in the flood indicators can be defined as where β result is the relative change, and R base and R i are the model results for flooded area, and inundation depth under the baseline and perturbation scenarios, respectively.Similarly, the relative change in influential factors can be defined as where β f actor is the relative change, and F base and F i are influential factors under the baseline and perturbation scenarios, respectively.
The ratio β result / β f actor for one factor represented the measure of its impact on flooding.The larger the ratio, the more powerful the factor.

Model Performance
In the model calibration, the best combination of roughness coefficients produced the smallest error of the mean inundation depth, 0.06 m (Table 5).The maximum inundation depth simulated using the best combination is shown in Figure 5. Compared with the upstream area, the downstream area was flooded more severely, and the inundation depth was more than 1.5 m in some places.

Model Performance
In the model calibration, the best combination of roughness coefficients produced the smallest error of the mean inundation depth, 0.06 m (Table 5).The maximum inundation depth simulated using the best combination is shown in Figure 5. Compared with the upstream area, the downstream area was flooded more severely, and the inundation depth was more than 1.5 m in some places.

Responses to Rainfall Intensity Change
The model was run by varying the rainfall scenario based on the baseline scenario (5year design rainfall, 2-year design tide, and land cover in 2020), and the two indicators, average inundation depth and flooded area, were computed for all scenarios (Figure 6).A significant monotonic increasing trend with return period was observed for both indicators.Average inundation depth increased significantly by 63.8%, from 0.58 m at 5-year scenarios to 0.95 m at 100-year scenarios.Flooded areas rose from 7.43 km 2 at 5-year scenarios to 12.97 km 2 at 100-year scenarios, showing a spatial expansion of 74.6%.
The increases in both the average inundation depth and flooded area exhibited nonlinear trends among the five return periods.The increase in inundation depth between the LX GH XJ

Responses to Rainfall Intensity Change
The model was run by varying the rainfall scenario based on the baseline scenario (5-year design rainfall, 2-year design tide, and land cover in 2020), and the two indicators, average inundation depth and flooded area, were computed for all scenarios (Figure 6).A significant monotonic increasing trend with return period was observed for both indicators.Average inundation depth increased significantly by 63.8%, from 0.58 m at 5-year scenarios to 0.95 m at 100-year scenarios.Flooded areas rose from 7.43 km 2 at 5-year scenarios to 12.97 km 2 at 100-year scenarios, showing a spatial expansion of 74.6%.5-and 20-year return periods (0.191 m) was larger than that between the 20-and 100-year return periods (0.176 m); similarly, increases in flooded areas were 2.986 km 2 between the 5-and 20-year return periods and 2.546 km 2 between the 20-and 100-year return periods, indicating a decrease in the rising rate.

Responses to Tide Level Change
The responses to tide levels are shown in Figure 7.Although the changes in the two indicators caused by variations in tide levels were much smaller than those in the rainfall The increases in both the average inundation depth and flooded area exhibited nonlinear trends among the five return periods.The increase in inundation depth between the 5-and 20-year return periods (0.191 m) was larger than that between the 20-and 100-year return periods (0.176 m); similarly, increases in flooded areas were 2.986 km 2 between the 5-and 20-year return periods and 2.546 km 2 between the 20-and 100-year return periods, indicating a decrease in the rising rate.

Responses to Tide Level Change
The responses to tide levels are shown in Figure 7.Although the changes in the two indicators caused by variations in tide levels were much smaller than those in the rainfall scenarios, they also showed an increasing trend.Average inundation depth increased by 6.9%, from 0.58 m at 5-year scenarios to 0.62 m at 100-year scenarios.Flooded areas rose from 7.43 km 2 at 5-year scenarios to 7.82 km 2 at 100-year scenarios, with a spatial expansion of 5.2%.

Responses to Tide Level Change
The responses to tide levels are shown in Figure 7.Although the changes in the two indicators caused by variations in tide levels were much smaller than those in the rainfall scenarios, they also showed an increasing trend.Average inundation depth increased by 6.9%, from 0.58 m at 5-year scenarios to 0.62 m at 100-year scenarios.Flooded areas rose from 7.43 km 2 at 5-year scenarios to 7.82 km 2 at 100-year scenarios, with a spatial expansion of 5.2%.
Compared with rainfall intensity, the increase in both indicators due to tide level changes exhibited an opposite nonlinear trend.The increase in inundation depth between the 2-and 20-year return periods (0.008 m) was smaller than that between the 20-and 100year return periods (0.036 m).Similarly, the increase in flooded areas was 0.153 km 2 between the 2-and 20-year return periods and 0.232 km 2 between the 20-and 100-year return periods, indicating an increase in the rate of increase.Compared with rainfall intensity, the increase in both indicators due to tide level changes exhibited an opposite nonlinear trend.The increase in inundation depth between the 2-and 20-year return periods (0.008 m) was smaller than that between the 20-and 100-year return periods (0.036 m).Similarly, the increase in flooded areas was 0.153 km 2 between the 2-and 20-year return periods and 0.232 km 2 between the 20-and 100-year return periods, indicating an increase in the rate of increase.

Responses to Land Cover Change
Due to rapid urbanization in the study area, the percentage of impervious land increased from 13.9% in 2000 to 54.4% in 2020 (Table 6); however, this number is predicted to further increase in the future.The land cover product of the FLUS model predicted imperviousness to increase by 61.6% in 2025 and 70.4% in 2035, which would be slower than that in the first decade of the 21st century.Since the urbanization process is spatially heterogeneous in watersheds, the downstream area would be occupied at a faster pace than the upstream area.
Flood responses to land cover changes are shown in Figure 8.The magnitude of change in flood indicators was between those in the rainfall-and tide-level scenarios.Both indicators showed an overall increasing trend.Average inundation depth increased by 17.2%, from 0.57 m in 2000 to 0.67 m in 2035.Flooded areas rose from 7.15 km 2 in 2000 to 8.53 km 2 in 2035, with a spatial expansion of 19.3%.Responses to land cover changes exhibited nonlinear trends like the responses to the tide levels.The rate of increase during 2020-2035 was found to be larger than that during 2000-2020, while the increase in inundation depth between 2000 and 2020 (0.009 m) was smaller than that between 2020 and 2035 (0.088 m).Increases in flooded area were 0.28 km 2 and 1.1 km 2 for 2000-2020 and 2020-2035, respectively, suggesting a stronger turn in the rising rate than the nonlinear trend in the responses to tide level.8.The magnitude o change in flood indicators was between those in the rainfall-and tide-level scenarios.Both indicators showed an overall increasing trend.Average inundation depth increased by 17.2%, from 0.57 m in 2000 to 0.67 m in 2035.Flooded areas rose from 7.15 km 2 in 2000 to 8.53 km 2 in 2035, with a spatial expansion of 19.3%.Responses to land cover changes ex hibited nonlinear trends like the responses to the tide levels.The rate of increase during 2020-2035 was found to be larger than that during 2000-2020, while the increase in inun dation depth between 2000 and 2020 (0.009 m) was smaller than that between 2020 and 2035 (0.088 m).Increases in flooded area were 0.28 km 2 and 1.1 km 2 for 2000-2020 and 2020-2035, respectively, suggesting a stronger turn in the rising rate than the nonlinea trend in the responses to tide level.

Comparison of Nonlinear Responses
Since changes in a single factor in the given scenarios did not show an equal rate o change in other factors, the analysis of a single factor cannot reveal the responses pre cisely.Instead, a perturbation analysis based on the baseline scenarios, 5-year design rain fall, 2-year design tide, and land cover in 2020 was conducted to compare the influences of the three factors on the flooding indicators.The results are shown in Figure 9.

Comparison of Nonlinear Responses
Since changes in a single factor in the given scenarios did not show an equal rate of change in other factors, the analysis of a single factor cannot reveal the responses precisely.Instead, a perturbation analysis based on the baseline scenarios, 5-year design rainfall, 2-year design tide, and land cover in 2020 was conducted to compare the influences of the three factors on the flooding indicators.The results are shown in Figure 9.
The impact of a factor can be represented by the slope, which is defined as the relative change in the flood response divided by the relative change in the factor (Table 7).Generally, both flooding indicators were most influenced by rainfall intensity, followed by land cover and tide levels.Considering the curve segments that connect the designed scenarios, considerable variations can be observed in the slope of the segments in each curve.The rainfall intensity curve varied the least among segments, whereas the other two curves showed stronger nonlinearity, increasing faster for larger values of factor.The curve corresponding to land cover could be divided into two phases according to the scenario of 2020, suggesting that the considered watershed is currently at a critical point regarding land cover imperviousness, after which the flood risk will increase rapidly.The curves of tide level showed strong rising trends in the last phase, S4-S5, with slope of 0.70 (flood depth) and 0.33 (flood area).The impact of the tide level at extreme tides was greater than, or at least comparable to, land cover change.Considering the continuously rising sea level and the accelerating global water cycle, the negative impact of tides may increase.
If a diagonal line starting from the original point is drawn in the upper right portion, most curves are under the line, suggesting that relative changes in influential factors do not induce the same amount of relative changes in flood depth and flood area.A probable reason for these outcomes is depth-area-volume relationship of detention areas such as floodplains, which could be simplified as circular cones.If water is continuously added to a circular cone, increases of water area and water depth triggered by a constant volume show a declining trend.Moreover, curve shapes are subject to floodplain topography.Table 7. Slopes for responses to rainfall intensity, land cover, and tide level.S1-S5 denote the five scenarios designed for a given factor.For example, rainfall intensity scenarios S1, S2, S3, S4, and S5 correspond to 5-, 10-, 20-, 50-, and 100-year return periods, respectively.Phase S1-S2 denotes response change from 5-to 10-year return periods.The depth and area represent the average inundation depth and flooded area, respectively.Table 7. Slopes for responses to rainfall intensity, land cover, and tide level.S1-S5 denote the five scenarios designed for a given factor.For example, rainfall intensity scenarios S1, S2, S3, S4, and S5 correspond to 5-, 10-, 20-, 50-, and 100-year return periods, respectively.Phase S1-S2 denotes response change from 5-to 10-year return periods.The depth and area represent the average inundation depth and flooded area, respectively.If a diagonal line starting from the original point is drawn in the upper right portion, most curves are under the line, suggesting that relative changes in influential factors do not induce the same amount of relative changes in flood depth and flood area.A probable reason for these outcomes is depth-area-volume relationship of detention areas such as floodplains, which could be simplified as circular cones.If water is continuously added to a circular cone, increases of water area and water depth triggered by a constant volume show a declining trend.Moreover, curve shapes are subject to floodplain topography.

Impact of Tides on River-Dominated Watersheds
According to the model analysis, compound flooding in the Yonghe River was primarily dominated by precipitation, and this conclusion is consistent with previous studies on the Pearl River Delta (PRD), where flood discharges dominated the flood stages in the middle and upper PRD, while storm surges had the maximum impact near the river mouth [33].Since the YRW is located in the upper part of the PRD, the impact of tides in this region is smaller than that near the river mouth.Nevertheless, the increase in tidal impacts should receive more attention in long-term policymaking.Qiu et al. [34] reported that extreme sea level has been the dominant driver of compound coastal-fluvial flooding in the PRD over the past 60 years, and its impact is projected to rise with global warming in the future.
In the past two decades, multiple typhoons have caused record-breaking storm surges at the mouth of the Pearl River, leading to severe damages to coastal cities in Guangdong, such as Zhuhai, Shenzhen, and Guangzhou.For example, at almost all hydrological stations near the Pearl River mouth, water levels recorded by Mangkhut (No. 1822) and Hato (No. 1713) ranked first and second in the history since hydrological records became available, respectively.In Zhongda Station, the three largest water levels were 8. Similar studies have been conducted in other regions.An analysis in Haikou, a city in southern China, showed that the return period of storm tides was reduced by 10.3% under the nonstationary scenario [25].Unusual long-distance inland propagation can be observed in the Yangtze River [26] and the Mississippi River [22].Extreme tidal levels under the nonstationary scenario will cause small-and medium-sized watersheds to face greater flood hazards than ever before.
In summary, flood risk management in coastal watersheds that are currently riverdominated should also consider the impact of tides, particularly for extreme water levels.As suggested in this study, the flood risk increases faster at high tide levels owing to the nonlinear feature of the flood response to the tide.

Timing for Flood Risk Management in Developing Regions
Flood risk management measures, particularly engineering measures, require substantial financial support, which is even more challenging for developing regions that struggle to balance multiple criteria, including economic growth, pollution control, disaster reduction, etc.The demand and supply of flood risk management are always imbalanced, which further amplifies flood risk by simultaneously increasing hazards and exposure.Rapid changes in land cover are the main reason for the increase in flood risk in developing regions.First, the transformation from a natural landscape to construction or residential land not only leads to losses in natural infiltration and detention capacity but also reduces surface roughness.Consequently, it increases the magnitude of flood discharge and leads to earlier peak timing, exacerbating the flood hazards.Second, the continuous growth of population and industry significantly increases exposure to floods, creating a new demand for flood protection.
The Environmental Kuznets Curve (EKC) provides a reference guide for the timing of flood risk management.The EKC was initially developed to describe the relationship between economic progress and environmental degradation over time as an economy progresses.The curve is an inverted U-shape, in which pollution increases with initial development but declines with further economic progress.Under this framework, the YRW was found to be at a turning point currently.As of 2020, more than half of the area was covered by impermeable surfaces, and metro systems have stretched to this area, suggesting that its socioeconomic condition has grown adequately.In contrast, the flood control capability of levees (less than a 10-year return period) and the design standard of the drainage system (approximately a 1-year return period) in this area are insufficient.Furthermore, the nonlinear flood responses to land cover changes (Table 7) confirm that increases in construction land in the future will increase the risk of flood hazards.
Based on the analysis conducted in this study, 50% imperviousness is recommended as an indicator of the urgent demand for flood risk management measures in small watersheds.A previous study in the Yangtze River Delta [35] explored the impact of urbanization on floods in a basin with a 1371 km 2 area and found that urbanization impacts can be better understood at the sub-basin level.Suburban areas may experience a high flood risk with urbanization expansion, even though the impact at the basin level may not be remarkable.Taking Guangzhou as an example, its built-up area in 2022 was 1380.6 km 2 , which accounted for 18.56% of the total area of 7434.4 km 2 .However, highly developed small watersheds such as the YRW are expected to face flood hazards more severe than ever before.The findings of the study on the PRD are consistent with conclusions at the sub-basin level in the Yangtze River Delta.The results for these two world-class delta regions are informative for recognizing the optimal timing for flood risk management in other small coastal watersheds.

Uncertainties and Limitations
First, the hydrodynamic model did not include the drainage process.During urbanization, the agricultural landscape was transformed into a mixture of manufacturing plants, newly built residential communities, and original villages and farmlands.The design standard for drainage in this suburban area is lower than that in the downtown areas.Moreover, this is quite common in villages where drainage depends mainly on overland flow.The exclusion of the drainage process is likely to impact simulations with small return periods but has limited influence on scenarios with large return periods.
Second, the model was calibrated using only the flood marks without timestamps.Without other events for validation, the model may be overfitted and, thus, its ability to generalize may be questioned.This concern can be largely addressed by the fine-scale data used in this study, including topography and land cover.In addition, the absence of flood peak time may produce uncertainties in model calibration.It should be noted that flood depth is time-dependent.If a hydrograph is used for model calibration, deviation of flood peak time would reduce the model performance significantly.Nevertheless, in a small watershed like YRW, flood processes are quick and, thus, it is acceptable to calibrate a model only with maximum flood depth.
Third, the roughness parameters remained unchanged among the five land cover scenarios, which is likely to underestimate the nonlinearity of flood responses.In the first (2000) and second (2010) land cover scenarios, using roughness parameters from the 2020 scenario likely overestimated flood responses.In contrast, underestimation may have occurred in future land cover scenarios (2025 and 2035).Consequently, the nonlinearity of flood responses to land cover change was underestimated to some extent.However, only one flood event was available for model calibration; thus, the roughness parameters for the other four scenarios could not be precisely determined.
Finally, this study does not include the impact of climate changes, which are influencing rainfall intensity and tide level across the globe.As mentioned in Section 4.1, record-breaking tide levels have already reduced the return periods.The impact is the same for rainfall intensity.Based on model results in this study and climate change trends, it can be anticipated that tide level will have a continuously rising impact on flood hazard from a long-term perspective.

Benefits and Future Work
Population and economic growth in coastal regions create an urgent demand for flood risk regulation measures.According to China's statistical yearbooks, the shares of coastal provinces in the total population were 41.33%, 43.25%, and 45.05% in 2000, 2010, and 2020, respectively, showing a clear move towards coastal areas.Furthermore, coastal provinces have accounted for 76.22% of the net population increase over the past two decades.Regarding the economy, the gross domestic product (GDP) of coastal provinces increased from USD 7894 trillion in 2000 to USD 76,703 trillion in 2020, presenting nearly a tenfold expansion.
The conclusions on the nonlinear flood responses are an important reminder for small watersheds that have not experienced extreme flood events like the one that occurred on 22 May 2020.According to incomplete statistics, out of more than 50,000 small-and medium-sized river watersheds in China, 85% are located along the coasts of river basins [36].Compared with the significantly enhanced flood prevention in the main rivers, small-and medium-sized river watersheds are weak points in flood risk management systems.Approximately 70-80% of flood-related losses in China are caused by floods in small-and medium-sized watersheds [37].
In developing regions, the lack of hydrological data is common, making it difficult to understand flood mechanisms and assess the intensity of future flood hazards.Although extreme flood events, such as the 22 May 2020 flood in Guangzhou, opened a window to understanding flood mechanisms, the rare occurrence of extreme events makes them infeasible for most watersheds.In future research, small watersheds should be classified into limited groups using machine learning techniques based on environmental factors and socioeconomic conditions.Thus, knowledge learned from limited extreme events can be transferred to other data-scarce small watersheds to enhance future flood risk management.In addition, the impact of climate changes should be incorporated to produce possible future flood hazard scenarios.

Conclusions
Compound flooding caused by multiple drivers simultaneously or in close succession poses a challenge to flood risk management in coastal regions owing to nonlinear interactions among drivers.Rapid land cover change leads to an increase in peak discharge and an advance in peak time; therefore, it plays a significant role in flood risk management in developing coastal regions.This study investigated the nonlinear responses of compound floods to rainfall intensity, tide levels, and land cover changes in a developing coastal watershed in Guangzhou, China, which is currently undergoing rapid urbanization.Regarding the impact on inundation depth and flood area, the most influential factor was rainfall intensity, followed by land cover change and tide level.Flood responses to all three factors showed a nonlinear increasing pattern but with different rates of increase.The response curves of the tidal level and land cover change increased faster at high values, whereas the rainfall intensity curve varied slightly.These findings highlight the importance of integrating tidal impacts into flood risk management even in river-dominated coastal watersheds.
In this study, the transfer of learned knowledge was discussed.Currently, the extreme water levels and long-distance inland propagation of storm surges are creating larger flood hazards in small watersheds, which are further amplified by rapid urbanization.Based on the findings of this study, 50% imperviousness is recommended to be used as an indicator of the urgent demand for flood risk management measures in small watersheds.
Additionally, this study had some limitations.First, the model did not incorporate the drainage process, which may have led to uncertainties in the simulation of small return periods.Second, the model was not validated after calibration because of a lack of historical flood data.Third, the analysis of the simulation results did not consider the probability of the co-occurrence of different drivers.Therefore, in future research, nonlinear flood responses should be tested across multiple spatial scales by considering additional combinations of drivers and their co-occurrence probabilities.Furthermore, knowledge learned from extreme flood events regarding environmental factors and socioeconomic conditions should be transferred to similar watersheds selected by machine learning techniques.

Figure 1 .
Figure 1.Study area location (a) and landscape (b).Huangpu (HP) and Zengcheng (ZC) are the two districts in eastern Guangzhou.The field survey sites A, B, and C represent Lixin Road, Guanhu Station, and Xinjie Village, respectively.Rain gauges in and around the study area include Jiancun (JC), Jiuyu (JY), Yeling (YL), Guangzhou Toyota (GT), and Xintang (XT).The watershed is empirically divided into upstream and downstream areas by the red solid line.The landscape image was taken from the Gaofen 1 Satellite in 2020.

Figure 1 .
Figure 1.Study area location (a) and landscape (b).Huangpu (HP) and Zengcheng (ZC) are the two districts in eastern Guangzhou.The field survey sites A, B, and C represent Lixin Road, Guanhu Station, and Xinjie Village, respectively.Rain gauges in and around the study area include Jiancun (JC), Jiuyu (JY), Yeling (YL), Guangzhou Toyota (GT), and Xintang (XT).The watershed is empirically divided into upstream and downstream areas by the red solid line.The landscape image was taken from the Gaofen 1 Satellite in 2020.Land 2023, 12, x FOR PEER REVIEW 4 of 17

Figure 2 .
Figure 2. Guanhu Metro Station during (a) and post (b) the flood.The photo (a) was taken on 22 May 2020, and sourced from Weibo, a public social media platform in China.The photo (b) was captured by the research team on 25 May 2020, when a rescue team was pumping flood water out from the metro station.

Figure 2 .
Figure 2. Guanhu Metro Station during (a) and post (b) the flood.The photo (a) was taken on 22 May 2020, and sourced from Weibo, a public social media platform in China.The photo (b) was captured by the research team on 25 May 2020, when a rescue team was pumping flood water out from the metro station.

Figure 2 .
Figure 2. Guanhu Metro Station during (a) and post (b) the flood.The photo (a) was taken on 22 May 2020, and sourced from Weibo, a public social media platform in China.The photo (b) was captured by the research team on 25 May 2020, when a rescue team was pumping flood water out from the metro station.

Figure 3 .
Figure 3. Temporal distribution of rainfall depth at rain gauges in and around the study area.

Figure 4 .
Figure 4. Water levels at Huangpu Tide Station under the National Elevation datum of 1985 from 21 May to 22 May.

Figure 4 .
Figure 4. Water levels at Huangpu Tide Station under the National Elevation datum of 1985 from 21 May to 22 May.

Figure 6 .
Figure 6.Average inundation depth and flooded area under rainfall scenarios in typical return periods.

Figure 6 .
Figure 6.Average inundation depth and flooded area under rainfall scenarios in typical return periods.

Figure 6 .
Figure 6.Average inundation depth and flooded area under rainfall scenarios in typical return periods.

Figure 7 .Figure 7 .
Figure 7. Average inundation depth and flooded area under tide level scenarios in typical return periods.

Figure 8 .
Figure 8.Average inundated depth and inundated area under different land cover scenarios.

Figure 8 .
Figure 8.Average inundated depth and inundated area under different land cover scenarios.

Figure 9 .
Figure 9. Influence of three influential factors (rainfall intensity, tide level, and land cover change) on average inundation depth and flooded area.

Figure 9 .
Figure 9. Influence of three influential factors (rainfall intensity, tide level, and land cover change) on average inundation depth and flooded area.
28 m in 2018 (Mangkhut, 1822), 7.79 m in 2017 (Hato, 1713), and 7.73 m in 2008 (Hagupit, 0814).At Dashi Station, corresponding water levels were 8.19 m in 2018, 7.95 m in 2017, and 7.79 m in 2008.These new extreme records have significantly changed the distribution of the water levels, thus reducing the return period of storm tides.

Author Contributions:
Conceptualization, H.H. and X.W.; methodology, H.H., Y.P. and C.W.; formal analysis and investigation, H.H., Y.P. and C.W.; writing-original draft preparation: H.H. and Y.P.; writing-review and editing, H.H., Y.P., C.W. and X.W.All authors have read and agreed to the published version of the manuscript.Funding: This research was supported by Guangdong Basic and Applied Basic Research Foundation (grant #2022A1515011494), and the Science and Technology Program of Guangzhou (grant #202201011726).

Table 1 .
Total rainfall depth and maximum three-hourly rainfall depth of adjacent rain gauges.
2.4.3.Parameter OptimizationEach combination was input to simulate the 22 May flood event, and its performance was measured by the mean deviation (MD) of the simulated inundation depths from

Table 3 .
Average intensity and rainfall depth for 3 h designed storms with typical return periods in Guangzhou, China.

Table 4 .
Tide levels for the five typical return periods at Huangpu Station based on the 1985 National Elevation Datum.Besides the 2020 land cover data, data from four additional years (2000, 2010, 2025, and 2035) were used.Land cover data for 2000 and 2010 were obtained from the Geographical Information Monitoring Cloud Platform (http://www.dsac.cn/DataProduct/Detail/200804, accessed on 6 October

Table 5 .
Model's performance in the calibration.

Table 5 .
Model's performance in the calibration.

Table 6 .
Percentage of impervious surface from 2000 to 2035.