Assessment of the effects of natural and anthropogenic drivers on extreme flood events in coastal regions

With climate change leading to sea level rise and frequent occurrence of extreme weather events, extreme flood disasters occur frequently in coastal cities, especially under the combined influence of heavy rainfall, fluvial floods, and storm surges. Compound floods from multiple sources are prone to superimposed effects, greatly aggravating flooding. In addition, the influence of anthropogenic activities such as man-made flood control projects in coastal regions also need to be explored. Using the Ling River Basin on the east coast of China as a case study, this study investigated the impacts of natural and human drivers on the spatiotemporal distribution of extreme floods based on a 1D–2D coupled hydrodynamic model. The hydrodynamic model was calibrated and validated based on historical flood events, and the validated model was further used to reproduce the devastating flood caused by Typhoon Lekima in 2019. The results show that the model has great capability and accuracy in simulating both the flood occurrence process and its inundation extent under natural and anthropogenic influences. The results also indicate that the high amount of discharge from the upper reaches of the Ling River was the main driver that caused the catastrophic consequences, and 83% of the inundation is caused by this factor. Comparatively, the failure of flood protection measures, which was a human driver, aggravated the influence of floods in local regions. Assessing the individual and superimposed effects of natural and human drivers on flooding may help enhance the resilience of flooded coastal areas.


Introduction
Coastal flooding is one of the most frequently occurring and destructive natural disasters, causing heavy casualties and losses and impacting many low-lying coastal areas around the world. Climate change has caused sea level rise and increased the frequency of hurricanes/typhoons and extreme flood events in coastal areas, such as Hurricane Katrina in the United States in 2005, Typhoon Haiyan in the Philippines in 2013 and Typhoon Lekima in China in 2019.
On a global scale, climate warming will very likely cause more flood inundations in coastal regions in the future. Coastal urban areas are floodplains with low-lying terrains and face the influence of multiple natural drivers ) such as sea level rise, storm surge, rainfall, river discharge, and astronomically high tide level, any of which can directly lead to flooding in the region. In reality, heavy rains and high tide levels often occur simultaneously (Zhou and Jie 2019). Furthermore, the combination of heavy rains, upstream river floods, and storm surges will greatly aggravate the severity of the floods (Xu et al. 2019). Therefore, costal compound flood hazards have increasingly attracted the interest of researchers (Bushra et al. 2019;Kew et al. 2013;Hurk et al. 2015;Wahl et al. 2015), and their threat emphasizes the importance of taking a holistic approach to assessing flood risk in coastal areas (Pasquier et al. 2019).
In addition to the natural environment and climate change, human intervention and socioeconomic development also have a negligible impact on spatiotemporal changes in flood inundation in coastal regions. With the accelerating process of urbanization, anthropogenic drivers constantly change the underlying surface conditions, and the proportion of impervious area has increased year by year. This alters the natural hydrological process and results in large flood peaks in urban areas (Zhang et al. 2016a, b).
Meanwhile, coastal cities are densely populated and have highly developed economies. Once floods occur, the losses are magnificent. Humans have carried out a series of flood protection measures (FPM) to reduce flood risk, such as the construction of reservoirs, embankments, roads, railways, bridges and other infrastructure in urban areas. However, these measures may sometimes lead to additional hazards in the case of extreme floods, such as reservoir bursts, dike bursts, and culvert blockages (Gan et al. 2018). Therefore, assessing the positive and negative contributions of human influence has become an important perpetual concern.
Researchers have studied the dependence and occurrence probability between the drivers of coastal floods from a statistical point of view, but they cannot explain the impact of these drivers on the risk of disaster-causing flooding nor identify specific areas in coastal areas that are vulnerable to flooding. The development of hydrodynamic models provides an opportunity for this (Ciullo et al. 2019;Gilles et al. 2012;Lee et al. 2018;Li et al. 2016;Schmid-Breton et al. 2018;Sharma et al. 2018). Hydrodynamic modeling has been widely used in flood simulation and risk analysis (Jung and Jasinski 2015;Koue et al. 2018;Li et al. 2017;Mardani et al. 2020;Ramos-Fuertes et al. 2014;Zhang et al. 2015;Zhou et al. 2018). The most common model is one-dimensional (1D) hydrodynamic simulation. The 1-D model is convenient for modeling and fast calculation and is suitable for describing channel processes. However, it does not behave well in simulating dynamics over floodplains, which have 2-D characteristics (Carrivick 2006;Moya et al. 2016;Pathirana et al. 2011). Considering the advantages of the 1D model in the channel and the 2D model over the floodplain, 1D/2D coupled models (Patel et al. 2017) were developed. These coupled models can dynamically represent the interactions among coastal regions, urban areas, rivers and floodplains. In recent years, they have been successfully used to simulate the impact of river discharge, extreme rainfall, sea level rise and marine storm surge on floods in coastal areas (Bermúdez et al. 2019;Gallien et al. 2014;Xu et al. 2019).
Although previous studies have investigated the causes of extreme compound flood events, they usually mainly considered the natural drivers alone. However, studies considering the impacts of both anthropogenic and natural aspects in highly developed coastal cities are rare. Thus, using the Ling River Basin as the study area, a 1D/2D coupled hydrodynamic model with real flood disaster investigation data is used in this study to investigate the compound impacts of multiple natural drivers, including storm surge, inland fluvial flooding and heavy rainfall, as well as the anthropogenic local FPM, on flood risk. Additionally, we also aimed to determine the relative influence of anthropogenic and natural drivers, which will provide useful information for the government and decision-makers.

Study area
The Ling River Basin, located in Zhejiang Province, China, was studied. The geographical location is between 28°40 0 -29°04 0 N and 120°49 0 -121°41 0 E. The landform of the study area is dominated by hills, and the terrain slopes from northwest to southeast. The basin is surrounded by mountains on three sides and is adjacent to the East China Sea on one side. The drainage of the study area mainly relies on the city's natural water system, the Ling River (as shown in Fig. 1). The Ling River originates from the Shifeng and Yong'an tributaries in the upper reaches. Both are mountainous rivers that are characterized by meandering channels, steep bathometry slopes, rapid currents and high banks. The middle reaches are the flat Datian Plain and Yicheng Plain. The Datian Plain is located in the urban area of the city of Linhai, which has a dense population and a developed economy. The Yicheng Plain is hilly and lowlying and is often waterlogged when it rains. The downstream is affected by the tides of the East China Sea.
Protection levels vary among river embankments in the study area. The embankment near downtown was built to protect the area against a 50-year flooding standard, but the other embankments meet the standards of only 20-year or 10-year floods. Some areas along the Yong'an and Shifeng tributaries have natural banks without any river embankment. Downtown Linhai contains an ancient town with a history of more than 1600 years. The ancient town wall is not only a historical and cultural relic but also has flood protection functions today. The ancient town wall is an important part of the current FPM system. The Ling River Basin was recently hit by the landing of Typhoon Lekima in 2019, which has been the strongest typhoon since 1949. Typhoon Lekima made landfall on the coast of Zhejiang Province on the early morning of August 10th, 2019. Rain began in Linhai from 8:00 on August 8th and stopped at 15:00 on the 10th. The maximum three-day areal rainfall in Linhai was 371.2 mm. Heavy rainfall led to a significant increase in the water level in the Ling River, and severe flash floods were triggered in some towns. Among these towns, the ancient town in Linhai suffered the most, as it was flooded for 3 days.
In summary, the upper reaches of the study area are threatened by mountain floods, and the lower reaches are affected by East China Sea tides. The economically developed and densely populated downtown Linhai is located in the middle reaches. Once flooded, casualties and economic losses are severe. The study area often suffers compound flood disasters and is a high-risk area of flood disasters.

Data
The spatial data (i.e., elevation, land use map, and soil distribution map), which were provided by the Geographic Information Center of Zhejiang Province, China, were prepared as inputs for the hydrodynamic model. High- quality topographic data are required to ensure the accuracy of the flood analysis model results. Digital Elevation Model (DEM) data of the study area with a resolution of 30 m were obtained. The extent, water depth and flow velocity of floods are very sensitive to roughness parameters. To express the spatial surface roughness, a grid file of land use/coverage with a resolution of 30 m in the study area was collected. The land use/cover of each category (such as bare land, residential area, industrial area, and water area) in the grid cell was identified. In addition, other relevant geographic information data in the area were also collected, such as rivers, roads, dikes, and basic water conservancy engineering facilities. High-precision channel geometry data are very important for the determination of the Manning coefficient of the riverbed roughness. In this study, a total of 279 river cross-sections in the basin were collected by the Linhai Water Resources Bureau. The hydrological boundary is an important input condition for the simulation. Because the study area is affected by floods from multiple sources, the boundary conditions included not only regional rainfall data but also the river discharge upstream and tide gauge data of downstream estuaries. The discharge data were measured at the Baizhi'ao hydrological station (along the Yong'an tributary) and the Shaduan hydrological station (along the Shifeng tributary) in the upper reaches of the Ling River. The rainfall data of four rainfall stations and tide level data from the Haimen station at the entrance of the sea were collected from the Linhai Water Resources Bureau.
Historical flood disaster data were used to validate the hydrodynamic model. In this study, four typhoon events that seriously impacted the east coast of China from 2007 to 2020 were selected. Through the Linhai Water Conservancy Bureau, flood disaster investigation data from Typhoon Wipha in 2007, Typhoon Sudiro in 2015, Typhoon Lekima in 2019 and Typhoon Hagupit in 2020 were collected and used to test the precision and accuracy of the model.

Methodology
The process of this study can be divided into three parts as shown in Fig. 2: hydrodynamic model setup and evaluation, driver scenarios design and assessment of the drivers' effect. The methods associated with the first two are described subsequently in Sect. 3. The results of the third part are presented in Sect. 4.

Model principle
The original model used in this study is a one-and twodimensional (1D-2D) coupled hydrodynamic model, which has been applied and validated in many areas along the coast of China (Zhang et al. 2020;Zhang and Tang. 2021). The details of the model can be seen in Appendix A.
The 1D river model is used to quickly simulate the hydrodynamic process of the river network in the basin. The 1D Saint-Venant equations were used for simulating flows in the river channel, and the four-point implicit Preissmann scheme was used for discretization. When the water flow overflows the channel and spreads to the land, the water flow shows obvious 2D characteristics, and the flow pattern is more complicated (Bermúdez et al. 2019). The floodplain is divided into 2D unstructured triangular grids to realize spatial dispersion, and a two-dimensional shallow water model is constructed to simulate the flood evolution outside the river (Shen et al. 2019). The unstructured finite-volume method is used for the numerical resolution of 2D shallow water equations (Yoon and Kang. 2004). The numerical fluxes are computed by the Roe algorithm (Fraccarrllo et al. 2003). The method of limited water depth for wetting/drying treatment is used (Brufau et al. 2002;Xu et al. 2012).
The coupling of the 1D river model and the 2D floodplain model is implemented through the hydraulic connection conditions at the dike breach (Dutta et al. 2007). At the location of the breach, weir flow equations are used to calculate the flow and water level at the breach, and 1D and 2D models are coupled through the water balance relationship (Zhang et al. 2010). The key to model coupling is the interaction of hydraulic information between river channels and the land. The hydraulic parameters of the cross-section are calculated by the 1D river network model, and the average hydraulic parameters of the 2D grid are calculated by the 2D floodplain model. The grid boundary at the breach is connected to the cross-section of the river channel. The 1D river network model provides the discharge Q as a boundary condition for the 2D model, and the discharge Q is distributed at the water inlet of the 2D grid computing unit. The water levels in the 2D model grids varies, and the average water level Z of the grid at the 2D coupling boundary is returned to the 1D model to calculate the discharge at the next time step. In this way, the coupling of 1D and 2D models is realized. The coupling principle is shown in Fig. 3.

Model setup
3.1.2.1 One-dimensional (1D) river network model setup Based on the data of the Ling River Basin, a 1D river network model is constructed to simulate the streamflow process (Fig. 4). Modeling is performed for the main rivers, including the Ling River and its four main tributaries (i.e., Shifeng, Yong'an, Yicheng, and Datian). Due to the lack of cross-section data, the centralized inflow mode is used for the other smaller tributaries. The centralized inflow mode adds the tributary flow to the corresponding nodes of the main river. The generalized section interval of the 1D river is approximately 200-2000 m, with a total of 279 cross-sections, and the total length of the generalized river channel is approximately 176 km.  The boundary conditions input into the model include two categories: The upper boundary of the 1D river network is selected from the Shaduan hydrological station and the Baizhi'ao hydrological station, and the lower boundary is set at the Haimen tide gauge station.
3.1.2.2 Two-dimensional (2D) floodplain model setup In this study, the maximum possible inundation extent was defined as the entire study area of the 2D floodplain, and the boundary of the study area is delineated based on DEM. An unstructured triangular grid (Papaioannou et al. 2016) that is easy to adapt to complex terrain is selected to generalize the ground surface, and the total number of grids divided in the study area is 13931. In this study, the water blocking and diversion effects (Zhang et al. 2016a, b) of linear features (such as highways, railways, provincial roads, dikes, etc.) are considered and treated as internal grid division constraints to reflect their influences on water flow and flood routing process. In the highly populated and developed urban center area, to accurately simulate the inundation extent and depth, the grid is densified with an average grid area of 0.0004 km 2 . The grid is assigned elevation values according to the DEM, as shown in Fig. 5.

Model calibration
Calibration process is very important to the accuracy of model. Previous studies (Bharali and Misra 2021;Barati et al. 2012;Akbari and Barati 2012) have shown that the sensitivity of Manning's roughness coefficient is critical to the simulation results. In this study, the roughness coefficient is taken as the main calibration parameter in the process of model calibration.
Our strategy is to calibrate Manning's roughness coefficient separately according to the characteristics of each river in the study area. Firstly, we set the initial values (0.2-0.65 in this study) by referring to the hydraulic design manual and other references (Zhao and He 2010). Secondly, using try and error approach, the roughness values that best fit the observed water level process were determined.
Using the measured discharge data during the Sudiro flood (from 0:00 on August 9, 2015, to 20:00 on August 11, 2015) and the Hagupit flood (from 13:00 on August 3, 2020, to 22:00 on August 5, 2020) as the model calibration input, the simulated hydrological process is compared with the measured water level process, as shown in Fig. 6a, b. It  can be seen that the parameters of the model can accurately reflect the flood routing process. Both the Yong'an tributary and Shifeng tributary are mountainous rivers, whose riverbeds mostly consist of sand, gravel and pebbles. The Shifeng tributary runs through the canyon shoal in Linhai, with a meandering channel course and alternating deep pools and shoals. The upstream terrain is steep and undulating. The channel is rich in sand and gravel, and the roughness is considered to be 0.065. In the upper reach of the Yong'an tributary, the riverbed is steep and the flow is rapid. Since the riverbed is stony and uneven, the roughness is considered to be 0.054. Downstream of the Yong'an tributary and Shifeng tributary to the confluence of the Ling River, the terrain is relatively gentle, the river channel is relatively flat and straight, and the roughness of the river is considered to be 0.030. The river channel from the confluence to the middle reaches of the Ling River is clean and straight without beach and pool, and the roughness is considered to be 0.028. The Datian tributary, Yicheng tributary and the lower reaches of the Ling River tributary are all plain river networks, and the roughness of the river channel is considered to be 0.030 (Fig. 6c).
The roughness values for the floodplain grids are calculated based on the land cover types. The roughness values for different types of land use are obtained from Liu et al. (2017). The influence of buildings on roughness was not considered in this analysis since house occupation ratios are unknown and considered to be insignificant in the study area. According to the relevant regulations of land use classification, six different land use types are extracted, and the roughness value corresponding to each land use type is shown in Table 1.
In addition to the key calculation parameter of roughness coefficient, the model structure parameters need to be set, including time step, initial conditions and internal boundary conditions. Considering the stability and time requirements of the model calculation, the time step of the model calculation is considered to be 30 s. The initial condition parameter, that is, the depth of ponding on the surface at the beginning of the model, is set to 0 in this study, assuming that there is no ponding on the surface. The inner boundary condition which is the judged water depth (m) of the dry and wet interface is considered to be 0.1.

Model verification
The model results for river flow are validated during the the Typhoon Wipha flood event in 2007. Compared with Typhoon Soudelor in 2015, this typhoon is more intense. Using the measured discharge data during the Wipha flood as the model's upper boundary input and the tide level data at Haimen as the lower boundary input, the simulated hydrological process is compared with the measured water level process, as shown in Table 2. Table 2 shows that the difference between the simulated water level and the measured water level at the validation points (C1-C6 in Fig. 4) along the Yong'an tributary, Shifeng tributary and Ling River is within 0.10 m, indicating that the model simulation is good, and the parameter selection and cross-section generalization are all reasonable.

Scenarios of natural drivers
When the inundation depth depends only on the river flow, a 100-year discharge will produce a 100-year depth. However, this is not the case when multiple drivers are involved in the flooding process, which is common in coastal urban areas. During the Typhoon Lekima, the measured peak discharge of the Shaduan hydrological station is 3680 m 3 /s, which is the largest discharge in the record of the Shifeng tributary. The measured peak discharge of the Baizhi'ao hydrological station is 6900 m 3 /s, which is the second-largest discharge in the record of the Yongan tributary. The time difference between the peak discharges of the two tributaries is only 1 h. Due to Typhoon Lekima, the maximum three-day areal rainfall in Linhai reaches 371.2 mm, breaking the local short-term rainstorm record. Meanwhile, the lower reaches of the Ling River encounters the astronomically high spring tide, and the measured highest tide level at the Ximen tide gauge reaches 10.98 m. According to the actual situation above, the extreme flooding caused by Typhoon Lekima is considered to be a comprehensive disaster event affected by multiple sources, including heavy rainfall, fluvial floods, and high tide downstream.
To understand the independent or superimposed impact of a single factor and combined factors of natural drivers on flood inundation, four scenarios are designed (Table 3). For each of these scenarios, the boundary conditions are adjusted to apply the three single drivers mentioned above and the combined drivers.

Scenarios of human drivers
Besides natural drivers on flooding, the impact of human drivers (such as failure of local FPM) is also worth exploring. The ancient town wall is part of the Ling River embankment. The field investigation after the Typhoon Lekima flood show that as the water level of the Ling River continued to rise, the flooding outside the gates became too deep and greatly exceeded their capacity. Consequently, the Chaotian Gate, Wangjiang Gate, Zhenning Gate, Xingshan Gate, and Jingyue Gate were broken successively. The flood entered the ancient town through these five breaches, resulting in an average inundation depth of more than 1.8 m. To evaluate the impact of the embankment breach on the flooding, two scenarios of FPM were designed as follows.
(1) In the scenario of a dike-break flood, according to the gates that were actually broken, each gate serves as a breach (the location of the breach is shown in Fig. 7) to simulate the flooding process under the actual breach condition.
(2) In the scenario of the normal function of the FPM, we assume that all functions of the existing FPM are normally operated.
By simulating and comparing the flood consequences under these two scenarios, the change in flood inundation driven by the failure of FPM is evaluated.

Evaluation of model performance in reproducing the Lekima flood event
To evaluate the model's ability to reproduce the Lekima flood, the accuracy of the model results is evaluated by using the measured data during the flooding period and the field survey data after the flood.

Assessment of inundation extent
The maximum inundation extent simulated by the model is compared with the extent observed by satellite images. Figure 8a shows the satellite remote sensing map of the flood inundation provided by the Academy of Aerospace Information, Chinese Academy of Sciences and the National Satellite Meteorological Center. The blue area indicates the flooded area after the Lekima, and the light yellow area indicates the mud flooded area. Figure 8b shows the model simulation results, in which the blue area is the maximum flood inundation extent. Figure 8c shows a comparison between the model and observations. The comparison shows that the model-simulated flood extent and the satellite-observed flood extent tend to be consistent in their spatial distribution, and the inundation extent simulated by the model is slightly larger than that shown in the satellite image. The reason for the difference is that the model result expresses the real-time maximum flooded area of the study area under Typhoon Lekima, and all flooded areas are reflected. The remote sensing image is taken after the flood event. According to the simulation results, the flood inundation area is distributed in some villages on both sides of the Yong'an tributary and the Shifeng tributary, as well as the area of the Yicheng Plain and the Datian Plain, which are located on both sides of the Ling River. More than 70% of the main body of the inundated area is concentrated in the urban area near the Linghai urban areas (the red wire frame area in Fig. 8b). Analyzing the topographic map of the study area, it can be found that the Linhai urban area is located in a long and narrow area of flat land between the mountains and that this special terrain is one of the major causes of floods in this area.

Assessment of simulated inundation depth
By comparison with typical flood mark points, the accuracy of the simulated inundation depth in the floodplain area is verified. A total of 18 typical flood mark points with measured water level data are selected (Fig. 8c). The comparison results are shown in Table 4. The differences between the measured water depth and simulated water

Evaluation of the impact of natural drivers
The simulated inundation areas corresponding to different water depth classifications are shown in Table 5. A comparison showed that in the results of case 4, the average flood depth reaches 2.53 m, and the area with a flood depth greater than 3 m reaches 46.05 km 2 , which causes serious disaster to the study area. Moreover, this comprehensive effect is not a simple summation of the single effects of each driver of upstream fluvial floods, heavy local rainfall and downstream high tide. Under different water depth ranks, this superimposed effect makes the difference between an increasing and decreasing flood depth. Under the influence of the single driver of rainfall (case 2) or high tide (case 3), the inundation area is mainly concentrated in the shallow water depth range (0.5-1 m). As the water depth increases, the inundation extent decreases gradually, and the proportion of inundation area corresponding to deep water ([ 2 m) is the smallest. Especially under the influence of high tide drivers, the submerged surface is very small and concentrates only in the area with the lowest elevation. Comparatively drivers from various sources, and its inundation law is roughly the same as that driven by upstream fluvial floods. The area submerged by deep water ([ 2 m) occupies more than half of the total submerged area. However, due to the complexity of the interaction mechanism among various drivers, this study cannot quantify the specific contribution of individual natural drivers. Figure 10 shows the difference in the flood inundation extent under different natural driver scenarios. Due to the differences in drivers, the submerged depth and extent disproportionately vary. The driver of the downstream high tide level has the smallest contribution to the regional flood inundation. The impact caused by the rainfall driver is greater than that at high tide, but the inundation area is much smaller than the impact caused by the combined conditions on Typhoon Lekima. Comparing the results driven by upstream fluvial floods, it can be seen that during   Typhoon Lekima, fluvial floods were the main driver that led to serious disasters in the Linhai urban area. Compared with other studies on the spatial pattern of storm surge, rainfall and river-driven coastal composite floods, it can be found that although the dominant drivers might be storm surges (Pasquier et al. 2019) or rainfall (Zellou and Rahali 2019) in different regions, the general consensus is that in coastal areas with rivers, river discharge contributes to the incurred additional flood damage of the inundation process (Kumbier et al. 2018;Shin et al. 2019), and from the estuary to the upper reaches of the river, the sensitivity of river discharge increases significantly. Especially in the southeastern coastal areas of China, single flow estuary areas that are similar to the Ling River Basin are the most typical; when quantifying the impact of flood drivers, the river discharge during typhoons is of particular concern.

Assessment of the impact of failure of FPM
A comparison of the two FPM scenarios is shown in Table 6, which shows that the difference in submerged area and depth under any water depth classification in the two scenarios is close. Both the flood extent and depth caused by gate breaches do not significantly aggravate inundation in the whole study area. There have been several breaches in river embankments, but they have not further aggravated the catastrophic impact of fluvial floods. To discover why this is the case, we further analyze the flooding process under the two scenarios (Fig. 11).
Fourteen hours after the start of the simulation, thefloods continuously advance from the upstream Shifeng tributary and Yong'an tributary to the Ling River, and some areas on both sides are submerged because the upper reaches of the Ling River can defend against only 20-year floods, whereas the flood caused by Lekima has a return period of 80 years (Chinese official report). The water level in the river is higher than that in river dikes, and the overflowing water spreads to the floodplain. Meanwhile, Fig. 11a, e show that under the two scenarios, the ancient town is not submerged due to the protection of the ancient city wall. In addition, floods overflow from the Daqing tributary (a small tributary connected to the Ling River). After the 18th hour, the flood overflowing from the Daqing tributary in scenario 2 (Fig. 11f) gradually spread to the southwest and submerge part of the densely built area. In addition, in scenario 1 (Fig. 11b), it can be seen that the water inflow from the breach of the gate continues to flood into the town, and nearly half of the ancient town is submerged. After the 23rd hour, the overflow flood in scenario 2 (Fig. 11g) already submerges more than half of the streets outside the ancient town in Linhai, but the ancient town is not yet submerged. In the results of scenario 1 (Fig. 11c), it can be seen that the flood that enters the ancient town from the breaches further expands the inundation area and tends to converge with the overflow flood from the northeast. At the 45th hour, the flooding stabilizes, and the final flooding range under the two schemes is very close.
Through the comparison of the inundation process under these two scenarios, it is concluded that even if there are no breaches of the city gates, because the existing embankment standard is lower than the high water level caused by Typhoon Lekima, a large number of floods overflow the embankment top from the upper reaches of the Ling River and its tributaries (such as the Daqing tributary) and finally cover the whole ancient town.
In addition to the inundation process, the changes in water depth at the five gates under the two scenarios are also compared (Fig. 12). The submerged water depths of all gate locations in scenario 2 are all lower than the water depth in scenario 1, and the difference in the maximum water depth at Chaotian Gate reaches 1.8 m. In terms of flood arrival time, the time at which the flood water begins to appear at the gate location in scenario 2 is 8-15 h later than in scenario 1. By analyzing the inundation water depth and flood arrival time, it can be found that due to the wall gate breach, the human driver significantly exacerbates the inundation of the ancient town, which indirectly makes it the most serious disaster area in the typhoon flood event.
Although in the ultrahistorical extreme flood event, the partial failure of FPM has little effect on the overall flooding, it has a significant impact on the depth of the flood water and flood diffusion path in parts of the area. Especially when this affected area has an extremely high economic or historical and humanistic value, the disastercausing contribution of by human drivers deserves special attention.

Conclusion
The overall goal of this study is to deepen the understanding of coastal flood drivers and flood inundations using extreme flood events in the Ling River Basin on the eastern coast of China as a case study. The following conclusions are drawn: 1. Combining the 1D-2D hydrodynamic coupling model with the postflood survey data, the calibrated model can satisfactorily reproduce the flood inundation process and results under Typhoon Lekima in 2019, which indicates the applicability and rationality of this method for flood simulation under the combined effect of multisource drivers in complex underlying surface areas. 2. By comparing the simulation results under different natural driver scenarios in the Lekima flood event, it can be seen that although rainfall in the urban area increases the surface runoff and the waterlogging and river flood cannot be discharged smoothly from the estuary due to the jacking of the high tide level, the high discharge of the upstream is the key natural driving force of the destructive disaster. 3. The impact of human drivers on flooding needs to be given more attention. The lessons from the Typhoon Lekima flood event show that even if there are few extreme events in the historical records, coastal cities vulnerable to floods may need to strengthen at least local FPM (such as embankment reinforcement) to protect their important densely populated areas. 4. Although the model simulation results are for a specific case, analyzing the roles of different drivers in historical extreme flooding may be helpful for understanding the relationship between the impact of multiple drivers, such as sea level rise, population growth and urban development in coastal areas, and flood inundation in the future. 5. One limitation of this study was the shortage of observation data, this paper only takes the extreme flood event of Typhoon Lekima as an example to carry out analysis. In future research, we will collect more flood events in the study area and conduct a more comprehensive evaluation and analysis of the driving force of multi-source floods in the area.

Appendix A
Numerical model

The 1D hydrodynamic model
The Saint-Vinant equations describing the one-dimensional flow of the river channel are: where q is the lateral discharge, g is the gravity acceleration, Q,A and Z are the flow discharge through A, the cross sectional area and the water level, V x is the component of the side inflow velocity in the direction of the water flow, which can generally be approximated to zero,K is the flow modulus, and a is the momentum correction coefficient.
The Eq. (1) are numerically discretized using the fourpoint linear implicit difference format, and the difference equations of any river reach are obtained as: Taking the water level of the first node and the water level of the last node as free variables, the three-coefficient chasing method is used to eliminate the water level and flow of the middle section, and finally two equations of the relationship between the flow of the first and last section and the water level of the first and last nodes are derived.
The two equations forms are as follows: in which: Z ðIÞ is the water level of the first node. Z ðJÞ is the water level of the last node. Q L1 and Q L2 are the first and last cross-section flow. From the back to the front, the flow of this section is expressed as a linear function of the water level of the cost section and the water level of the end node, and the recursion equation is as follows: Similarly, starting from the first reach, try to express the section flow as a linear function of the cost section water level and the first node water level.
Therefore, the formula can be obtained from the above recursive formula. When calculating the recurrence formula, a, b, n, h, g and c six catch-up coefficients need to be saved. Once the water levels of the first and last nodes are obtained, the flow rate of the same section can be calculated by using the Eqs. (4) and (5): Solve simultaneously to get: After obtaining Z i substitute it into the Eq.
(3) to get Q i .

The 2D hydrodynamic model
The 2D shallow water equations are used as the governing equations for flood routine in the land. 2D shallow water equations in conservation form are as follows, where h is the water depth; u and v are the depth averaged components of the velocity U along the x and y coordinates respectively, g is the gravity acceleration; t is a time variable, s ox ¼ Àoz b =ox s oy ¼ Àoz b =oy s ox and s oy are bed slopes in the x-and y-directions, z b is the river bottom elevation, s fx and s fy are the friction slopes in the x-and y-directions, n is Manning's roughness coefficient. Neglecting diffusion of momentum due to wind effects and the Coriolis term.
In the numerical discretization of the above equations, a cell-centered finite volume method is formulated for Eq. (8) over a triangular control volume is adopted. The variables exist in the center of the element and the boundary of the grid is the control. where the dependent variables of the system are stored at the center of the cell and represented as piecewise constants. Equation (8) can be written as: where F ¼ ðE; HÞ, and the integral form of the Eq. (9) over a fixed volume V can be introduced, Z U i is the average value of the conservative variables over the volume V i at a given time, Using Green's formula, the surface integral in Eq. (10) is transformed into line integral, as follows, in which, A i is the area of the control volume, oV i is the boundary of the grid, b S is the integral value of the source term and n is the unit outward normal vector.
The second term on the left in Eq. (12) can be described as follow.
where m is the number of cells edges, F i;j is the interface flux passing through the boundary j of the grid i,, l ij is the length of the boundary j of the grid i. F i;j can be solved by the approximate Riemann solution in Roe format. 3. The 1D-2D hydrodynamic model coupling A weir-type connection is used at the coupling section of 1D model and 2D model to realize the flow exchange between surface and river. Actual overflow per unit width (q) is calculated using Eq. (14) below. Then, the total discharge can be obtained multiplying overflow per unit width (q) by the width of the levee failure.
In which, q is overflow per unit width, h 1 = max (Zu, Zd)-Zb, h 2 = min (Zu, Zd)-Zb, Zu, Zd is the water levels inside and outside the breach, respectively; Zb is the elevation of the top of the breach, g is gravitational acceleration constant.

Wetting and drying treatment
In this paper, the method of limited water depth is used to deal with moving boundary problem. When dealing with moving boundaries, the limited depth method divides the grid into dry grid, wet grid and semi-dry grid.
Dry grid: at time n, if the grid water depth h \ htol1 and the water depth of adjacent units is also less than htol1, no flow and momentum will pass through the common boundary. If the water depth of all adjacent units is less than htol1, the unit is dry.
Semi-dry grid: at time n, if the grid water depth htol1 \ h \ htol2 and the water depth of adjacent units is less than htol2, there is only flow flux but no momentum flux on the adjacent interface.
Wet grid: at time n, the grid water depth H [ htol2.
Note: htol1 is the dry-wet threshold in the two-dimensional control parameters, htol2 defaults 5 times htol1.