Effect of Land Use and Drainage System Changes on Urban Flood Spatial Distribution in Handan City: A Case Study

: This study simulated urban ﬂooding under various land use and drainage system conditions and described the process of historical ground–underground construction and its inﬂuence on spatial variations in waterlogging, taking Handan City as an example. The obtained results can provide support for urban water security and sustainable urban water resource management. The land use change, represented by the expansion of sealed surfaces, has a positive impact on the distribution and the volume of ﬂood in Handan City, while the drainage system has the opposite effect. The ﬂooding distribution changes over decades reveal that ﬂooding risk is reduced in most areas by improved drainage conditions but exacerbated in impervious areas and riversides due to increasing impermeable areas, the rapid draining of pipes, and poor outlet conditions. This study demonstrates how the dual changes in land use and drainage pipeline networks affect urban ﬂooding distribution; we suggest considering land use and the extension of drainage pipelines in future construction.


Introduction
Against the background of rapid economic growth in China, the rate of urbanization in China increased dramatically from 26% in 1990 to 64.72% in 2021 [1]. Rapid urbanization causes many environmental problems, among which urban water security is becoming more and more notable [2]. Taking the 2021 Henan floods as an example, China's Henan province suffered from a catastrophic flooding disaster from 17 to 23 July 2021 as a result of heavy rainfall. Especially on 20 July, Zhengzhou, the provincial capital of Henan, experienced a severe river flood and concurrent urban waterlogging. A record amount of rain fell within one hour, reaching 201.9 mm, breaking the previous record in China, which had stood since 1951. The flood disaster caused 380 deaths and led to a direct economic loss of CNY 41 billion (USD 5.67 billion) after investigation and evaluation in Zhengzhou [3]. Urbanization has transformed the processes of the urban hydrological cycle, which, in turn, has intensified the frequency and intensity of extreme hydrological events and exerted a far-reaching impact on urban water security [4,5]. Therefore, it is essential to evaluate the influence of urbanization on urban floods and identify the impact factors.
Urbanization has affected processes of the urban hydrological cycle in the following ways: (a) expanding impermeable surfaces, such as sidewalks and buildings, may increase surface runoff and reduce groundwater recharge [6][7][8][9]; (b) lowering surface roughness by reducing the vegetation distribution and flattening the surface may accelerate the surface runoff process [10,11]; (c) hardening river bottoms may accelerate the river confluence process [12,13]; (d) building drainage pipelines may accelerate pipeline confluence [14,15]. Hence, the dual changes in ground (such as expanding impermeable surfaces and lowering surface roughness) and underground construction (such as hardening river bottoms and building drainage pipeline networks) must be considered in urban stormwater management.

Study Area
Handan City is located in southern Hebei Province, China (Figure 1a). At the end of 2021, the city had a total population of 9.36 million, ranking second in Hebei; GDP totaled CNY 411.4 billion (USD 56.94 billion), ranking fourth in Hebei [37]. Handan lies in the North China Plain, in a continental semi-arid climate zone, with strong monsoonal influence. The mean daily temperatures range from −0.9 • C in January to 27.3 • C in July, while the annual mean temperature is 14.3 • C. The average annual precipitation is 502 mm, and the rainy season generally occurs in July and August.
Surrounded by the ring road, the study area was the central area of Handan City (Figure 1b). The area of this region is 91.7 km 2 , and the geographical coordinates are 36 • 20 -36 • 44 N and 114 • 03 -114 • 40 E. The topography of Handan City is generally high in the south and west and low in the north and east (Figure 1c). With its typical northern climate characteristics and medium-sized area, the downtown Handan area represents a typical medium-sized, rapidly developing city in northern China. typical medium-sized, rapidly developing city in northern China.
The study area is a high-risk area for urban waterlogging problems. A typical example is the "7·18 Handan Heavy Rainstorm event" in 2016, which caused the collapse of 36,000 houses, affecting 2.04 million people, with 108 reported dead and missing, and a direct economic loss of CNY 18.9 billion (USD 2.62 billion). Urban flooding management is of great significance for Handan City in the prevention and mitigation of natural disasters. However, only a few researchers have performed studies related to urban stormwater management in Handan [38][39][40][41][42][43].

Data
In order to simulate the rainfall runoff of the study area, the input data included the following: (1) weather data providing the input rainfall data; (2) land use data providing the infiltration parameters for the runoff generation process and the roughness parameters of different land covers for the overland flow concentration process; and (3) drainage network data providing the concentration paths and their physical parameters for the concentration process. Additionally, to validate the simulation, social media data were used to represent the observed data.

Weather Data
The design storm curve in this study was provided by the water conservancy bureau of Handan, ascertained from the Chicago hyetograph method developed by Tongji University [44]. The study area is a high-risk area for urban waterlogging problems. A typical example is the "7·18 Handan Heavy Rainstorm event" in 2016, which caused the collapse of 36,000 houses, affecting 2.04 million people, with 108 reported dead and missing, and a direct economic loss of CNY 18.9 billion (USD 2.62 billion). Urban flooding management is of great significance for Handan City in the prevention and mitigation of natural disasters. However, only a few researchers have performed studies related to urban stormwater management in Handan [38][39][40][41][42][43].

Data
In order to simulate the rainfall runoff of the study area, the input data included the following: (1) weather data providing the input rainfall data; (2) land use data providing the infiltration parameters for the runoff generation process and the roughness parameters of different land covers for the overland flow concentration process; and (3) drainage network data providing the concentration paths and their physical parameters for the concentration process. Additionally, to validate the simulation, social media data were used to represent the observed data.

Weather Data
The design storm curve in this study was provided by the water conservancy bureau of Handan, ascertained from the Chicago hyetograph method developed by Tongji University [44].
Equation (1) is a common form of the design rainstorm intensity formula: where i is the mean rainstorm intensity, mm/min; T is the design rainfall return period, a; t is the rainfall duration, min; and A 1 , b, B, and c are parameters related to local rainstorm characteristics and need to be solved. Among them, A 1 is a parameter representing the rainfall intensity; B is a dimensionless parameter representing the variation in the rainfall intensity; b and c are constants independent of the return period. The design drainage frequency is 50%, which can provide reasonable safety for traffic and pedestrians at a reasonable cost. The formula only applies to rainfall continuing for less than 3 h, and the recurrence interval for storm sewer design is 2 years, according to the local design standard for wastewater engineering; therefore, a typical design rainfall process lasting for 1 h with a return period of 2 years was selected. Equation (2) is the design rainstorm intensity formula from the Handan Hydrology Bureau with experimental parameter values. The experienced value of A 1 is 6.025 mm/min, B is 0.86, b is 7.22, and c is 0.528. Thus, Equation (2) is obtained.

Land Use Data
The land use data of Handan in 1987Handan in , 1997Handan in , 2007, and 2017 were interpreted using remote sensing data from Google Maps. The spatial resolution of the downloaded remote sensing images of the study area was up to 1 m. For these remote sensing data, eCognition software was used to perform supervised classification, and the corresponding classification results and remote sensing images were corrected and coordinate-transformed. Based on the land surface properties of Handan City, the land use types were classified into 6 classes, i.e., farmland, grassland, water, woodland, bare land, and impervious surface ( Figure 2). Equation (1) is a common form of the design rainstorm intensity formula: where is the mean rainstorm intensity, mm/min; is the design rainfall return period, a; is the rainfall duration, min; and , b, B, and c are parameters related to local rainstorm characteristics and need to be solved. Among them, is a parameter representing the rainfall intensity; is a dimensionless parameter representing the variation in the rainfall intensity; and are constants independent of the return period.
The design drainage frequency is 50%, which can provide reasonable safety for traffic and pedestrians at a reasonable cost. The formula only applies to rainfall continuing for less than 3 h, and the recurrence interval for storm sewer design is 2 years, according to the local design standard for wastewater engineering; therefore, a typical design rainfall process lasting for 1 h with a return period of 2 years was selected. Equation (2)  (2)

Land Use Data
The land use data of Handan in 1987Handan in , 1997Handan in , 2007, and 2017 were interpreted using remote sensing data from Google Maps. The spatial resolution of the downloaded remote sensing images of the study area was up to 1 m. For these remote sensing data, eCognition software was used to perform supervised classification, and the corresponding classification results and remote sensing images were corrected and coordinate-transformed. Based on the land surface properties of Handan City, the land use types were classified into 6 classes, i.e., farmland, grassland, water, woodland, bare land, and impervious surface ( Figure 2).  Land use change is a major issue in rapid globalization and urbanization, particularly in rapidly developing countries and regions throughout Asia. After experiencing  Table 1, with the matrix depicted in Table 2.  Table 1 summarizes the changes in different types of land use over the past 30 years. The overall trend shows that the decreases in farmland and grassland are 24.46% and 13.85%, respectively, and the increase in impermeable area is 44.13%. Over the past three decades, the impermeable surface has increased significantly, from 32.44% to 76.57%.
From the land use transfer matrix (Table 2), the types of land use with the greatest changes are farmland, grassland, and impermeable surface. The increasing impermeable areas are mainly farmland (20.02 km 2 ), bare land (11.52 km 2 ), and grassland (9.24 km 2 ). Moreover, of all the converted land use types, the conversion rates to impermeable surface are as follows: 72.99% of farmland, 69.19% of grassland, 70.06% of water, 70.16% of woodland, and 78.53% of bare land in 1987 converted to imperious surface by 2017. It can be inferred that buildings and pavements with rapid urbanization are occupying all other types of land use.

Drainage Network Data
The drainage pipe network data provide the pipeline concentration parameters (length, shape, depth, and Manning's roughness coefficient of the conduit). When considering the ground and underground concentration progress, rivers are taken as open channels incorporated into the drainage network for calculation because solidification of the river bottom was completed. The pipeline concentration parameters of the river were set based on parameters such as the shape of the river cross-section, the feathering factor at the bottom of the river, and the slope of the river. The outlet node of the drainage to the river was set as the end node of each pipe. The water flow exchange between pipe and river followed the one-dimensional Saint-Venant equation.
Five main rivers run through the study area: Fuyang River, Qin River, Zhu River, Zhizhang River, and Shuyuan River. Qin River flows from west to east; Fuyang River flows from south to north; Zhizhang River bypasses the urban area from southwest to northeast around the central area. From the southernmost part of the river system upstream of the city, there is a confluence of the three rivers at Zhangzhuang Bridge, where there is a sluice gate; when heavy rains come, the sluice gate is shut to prevent incoming floods, and the upstream water flows into Zhizhang River to the southeast to protect the urban area from upstream flood. Thus, the upstream water parameter was set to 0 during the severe rainstorms in the SWMM.
To perform hydraulic calculations of the pipe network, according to the principle of area equivalence, multiple pipes were merged into one pipe, and those with an inner diameter of less than 300 mm were omitted to simplify the original structure. The simplified drainage network of the study area is shown in Figure 3. river was set as the end node of each pipe. The water flow exchange between pipe and river followed the one-dimensional Saint-Venant equation. Five main rivers run through the study area: Fuyang River, Qin River, Zhu River, Zhizhang River, and Shuyuan River. Qin River flows from west to east; Fuyang River flows from south to north; Zhizhang River bypasses the urban area from southwest to northeast around the central area. From the southernmost part of the river system upstream of the city, there is a confluence of the three rivers at Zhangzhuang Bridge, where there is a sluice gate; when heavy rains come, the sluice gate is shut to prevent incoming floods, and the upstream water flows into Zhizhang River to the southeast to protect the urban area from upstream flood. Thus, the upstream water parameter was set to 0 during the severe rainstorms in the SWMM.
To perform hydraulic calculations of the pipe network, according to the principle of area equivalence, multiple pipes were merged into one pipe, and those with an inner diameter of less than 300 mm were omitted to simplify the original structure. The simplified drainage network of the study area is shown in Figure 3. Drainage networks include the natural or manual removal pathways of ground and underground water from an area. The canal system was designed for 100-year frequency floods; there are six storage facilities with a total volume of 14 million m 3 . Figure 4 shows the drainage network of the study area in 1987, 1997, 2007, and 2017, while Table 3 lists the sub-catchments, nodes and links, and total pipe length of different years. The drainage network was initially constructed using GIS data provided by the East Sewage Treatment Plant of Handan. Drainage networks include the natural or manual removal pathways of ground and underground water from an area. The canal system was designed for 100-year frequency floods; there are six storage facilities with a total volume of 14 million m 3 . Figure 4 shows the drainage network of the study area in 1987, 1997, 2007, and 2017, while Table 3 lists the sub-catchments, nodes and links, and total pipe length of different years. The drainage network was initially constructed using GIS data provided by the East Sewage Treatment Plant of Handan.   The construction of the drainage system has rapidly expanded during the last 30 years. As shown in Table 3, the total pipe lengths grew from 99.77 km in 1987 to 251.05 km in 2017. In this study, the development of the drainage system was applied in urban flood simulations, and the relationships between land cover changes and spatial variations in waterlogging were analyzed.

Stormwater Management Model
The stormwater management model (SWMM), which has the advantages of being open-source and fast to compute, developed by the United States Environmental Protection Agency (US-EPA) in 1971, is a dynamic rainfall-runoff-subsurface runoff simulation model used for single-event to long-term (continuous) simulation of the surface/subsurface hydrology quantity and quality from primarily urban/suburban areas. For a large area such as Handan City, the SWMM can calculate the accumulated waterlogging amount of the node in each catchment. In order to visualize the results, the obtained data were equally distributed to each sub-catchment, and then the water surface was smoothed to determine the spatial distribution of waterlogging.  The construction of the drainage system has rapidly expanded during the last 30 years. As shown in Table 3, the total pipe lengths grew from 99.77 km in 1987 to 251.05 km in 2017. In this study, the development of the drainage system was applied in urban flood simulations, and the relationships between land cover changes and spatial variations in waterlogging were analyzed.

Stormwater Management Model
The stormwater management model (SWMM), which has the advantages of being open-source and fast to compute, developed by the United States Environmental Protection Agency (US-EPA) in 1971, is a dynamic rainfall-runoff-subsurface runoff simulation model used for single-event to long-term (continuous) simulation of the surface/subsurface hydrology quantity and quality from primarily urban/suburban areas. For a large area such as Handan City, the SWMM can calculate the accumulated waterlogging amount of the node in each catchment. In order to visualize the results, the obtained data were equally distributed to each sub-catchment, and then the water surface was smoothed to determine the spatial distribution of waterlogging.
Surface runoff can be generated by rainfall, snowfall, or the melting of snow or glaciers. According to the land use and surface drainage trend, a watershed can be divided into several sub-catchment areas, and the runoff process is calculated according to the characteristics of each sub-catchment area. Runoff in the SWMM is simulated for single or continuous rainfall and snowmelt [45]. The dynamic wave method is used to consider flow characteristics during surface flooding and the pressure effects caused by the inverse slope in some parts of a catchment [46]: where R is the runoff of a sub-catchment; S 1 is the area of the pervious surface; S 2 is the area of the impervious surface with depression storage; S 3 is the area of the impervious surface with no depression storage; i (t) is the function of rainfall intensity with time; e (t) is the function of evaporation intensity with time; f (t) is the function of infiltration intensity with time; and D 1 and D 2 are the depression storage of the pervious surface and impervious surface with depression storage, respectively. The kinematic wave method cannot calculate backwater, backflow, or pressurized flow; thus, it is only applicable to tree drainage networks. The dynamic wave method can describe the storage, sink, inlet and outlet losses, backflow, pressurized flow, and backwater of the pipe; thus, it is suitable for complex annular drainage networks. The elevation difference between the outfall bottom of the drainage network and the river bottom in this study was generally between 0.5 m and 1.0 m in Handan City; backflow and full pipe flow phenomena can occur during heavy storms. Therefore, this study used the dynamic wave method to simulate the drainage confluence in the main urban area of Handan.

Determination of Model Parameters
The SWMM can simulate the runoff process generated from sub-catchments, flow depth in open channels and closed pipes, and the ponding quantity for inlet nodes. The model requires three types of parameters: runoff-producing parameters, overland flow concentration parameters, and pipeline concentration parameters.
Sub-catchments are hydrologic units of land, the topography and drainage systems of which dispose of direct surface runoff to a single discharge point. The properties associated with a sub-catchment in the SWMM are as follows: area, width, slope (%), imperviousness (%), soil properties, Manning's n for impervious and pervious surfaces for overland flow, depth of depression storage in pervious and impervious areas, and the proportion of impervious area with no depression storage (%).
The 3D Analyst tools in ArcGIS 10.2 (ArcGIS, ESRI, Redlands, CA, USA) were used to generate the slope of each point from the digital elevation model (DEM) data.
Width is defined as sub-catchment divided by the length of stream flow. The length of overland flow is difficult to determine accurately; therefore, sub-catchment width parameters are subject to calibration to obtain a close match between observation and calculation. The sub-catchments are classified into two types by observing their shapes: right-angled trapezoids and others. Most smaller and numerous sub-catchments, located in the central zone of the study area, are right-angled trapezoids, whereas the larger and fewer subcatchments, located in the edge zone, exhibit various shapes. Hence, the length of the smaller sub-catchments can be calculated by computing geometry, and the length of larger sub-catchments can be measured in ArcGIS 10.2.
The sub-catchment impervious roughness is taken from the user's manual, and pervious roughness is calculated by area-weighted averaging for varying surfaces. The Manning's roughness values of each surface are listed in Table 4. The rate of infiltration is a function of soil properties in the sub-catchment area. The maximum infiltration rate, minimum infiltration rate, and decay constant α of the Horton infiltration model were 76.2 (mm/h), 3.81 (mm/h), and 0.0006, respectively, according to research performed in Tianjin [47], a municipality near the study area.
Inverse elevation and maximum depth are two pipeline parameters in the SWMM. The inverse elevations of each node were provided by the local sewage treatment plant, and the maximum depth was calculated by subtracting the invert elevation from the ground elevation of DEM data.

Validation
The rapid nature of urban flooding from storms makes the accurate surveying of water depths difficult to observe, thus hindering the validation of urban flood models [48]. Nevertheless, social media is now an open platform for collecting emergency information [49,50]. Scraping data through a Web Spider bot is an important way to harvest verifiable data. A heavy rainstorm in Handan triggered floods on 19 July 2016, leaving many validation-supporting footprints on the Internet. Weibo, a widely used social media platform in China, was chosen as the scraping dataset. The crawling targets were "rain" and "waterlogging," and the time range was defined as 19 July 2016, 0:00 to 20 July 2016, 12:00. The 646 related posts contained valid information, such as photographs, text, upload locations, and time stamps. In order to protect the privacy information of uploaders, the exchangeable image files were deleted while posting the message, resulting in the absence of the original shooting location of the photo, which meant that the uploading location could have been different from the waterlogging location. In order to confirm the location, it was necessary to identify landmarks in the image. Geographic locations are essential to validate the model. The geographic position could be located based on unique landmarks in the photograph through a mobile mapping application, such as Google Maps or Gaode Maps.
The input data of the SWMM contained the rainfall process of 19 July 2016, the land use map in 2016, and the drainage pipeline in 2016. The Web Spider harvested several thousand related photographs, but only 11 included precise location information. The model simulation results showed that all 11 points had accumulated water ( Figure 5). This indicates that the model inversion results are generally consistent with the actual situation.

Parameter Sensitivity Analysis
Parameter sensitivity analysis involves identifying the input parameters that exert considerable impact on the outcomes [51]. The Morris method, a widely used approach for parameter screening, is applied to determine the influence of land use changes and drainage pipeline changes on flood risk. The modified Morris method calculates the average value of the Morse coefficient to obtain the sensitivity factor S, as shown in Equation (4): where S is the sensitivity index; n is the number of times the SWMM was run; is the outcome of the SWMM run for the i t time; is the input for the i run; is the initial value of the outcome; and is the initial value of the input. According to previous research [52], the sensitivity indices are ranked into four classes: when | | ≥ 1.00, the parameter has very high sensitivity; when 0.20 ≤ | | < 1.00, the parameter has high sensitivity; when 0.05 ≤ | | < 0.02 , the parameter has medium sensitivity; and when 0.00 ≤ | | < 0.05, the parameter has low sensitivity [53].

Urban Waterlogging Risk Level
The impacts of water depth on different traffic subjects, such as pedestrians, trams, motor vehicles, electric vehicles, and automobiles, were used as the basis for classifying the risk. In total, five levels of risk are classified ( Table 5). The thickness of an athletic shoe sole is generally about 3 cm; thus, it is believed that a water depth of less than 3 cm does not affect pedestrians in walking. The height from the ankle to the ground is generally 10 cm; therefore, it is believed that a water depth between 3 cm and 10 cm will affect pedestrian walking. The critical value of water wading depth for electric vehicles is 30 cm [53], and the heights of the air outlets and inlets of automobiles are around 30 cm and 70 cm,

Parameter Sensitivity Analysis
Parameter sensitivity analysis involves identifying the input parameters that exert considerable impact on the outcomes [51]. The Morris method, a widely used approach for parameter screening, is applied to determine the influence of land use changes and drainage pipeline changes on flood risk. The modified Morris method calculates the average value of the Morse coefficient to obtain the sensitivity factor S, as shown in Equation (4): where S is the sensitivity index; n is the number of times the SWMM was run; Y i is the outcome of the SWMM run for the i t time; X i is the input for the i run; Y 0 is the initial value of the outcome; and X 0 is the initial value of the input. According to previous research [52], the sensitivity indices are ranked into four classes: when |S| ≥ 1.00, the parameter has very high sensitivity; when 0.20 ≤ |S| < 1.00, the parameter has high sensitivity; when 0.05 ≤ |S| < 0.02, the parameter has medium sensitivity; and when 0.00 ≤ |S| < 0.05, the parameter has low sensitivity [53].

Urban Waterlogging Risk Level
The impacts of water depth on different traffic subjects, such as pedestrians, trams, motor vehicles, electric vehicles, and automobiles, were used as the basis for classifying the risk. In total, five levels of risk are classified ( Table 5). The thickness of an athletic shoe sole is generally about 3 cm; thus, it is believed that a water depth of less than 3 cm does not affect pedestrians in walking. The height from the ankle to the ground is generally 10 cm; therefore, it is believed that a water depth between 3 cm and 10 cm will affect pedestrian walking. The critical value of water wading depth for electric vehicles is 30 cm [53], and the heights of the air outlets and inlets of automobiles are around 30 cm and 70 cm, respectively.
When the inlet port is submerged, traffic is basically paralyzed. Thus, 30 cm and 70 cm were set to be the critical heights of different flooding risk levels. Table 5. Urban flooding risk level classification.

Risk Level Water Depth (cm) Extent of the Threats to Different Means of Transportation
There is little impact on traffic. II 3-10 The water depth is higher than the insteps, thus slowing down pedestrians.

10-30
When the rail surface water depth exceeds 30 cm, electric vehicles must be shut down. IV  The water depth is higher than the outlet, thus impeding motor vehicles. V >70 The water depth is higher than the inlet port, basically paralyzing traffic.

Effect of Land Use Change
The distributions of flooding areas simulated using different land use in 1987, 1997, 2007, and 2017 are shown in Figure 6. The results were obtained directly from SWMM simulations, which used a two-year rainfall event and the drainage system in 2017. respectively. When the inlet port is submerged, traffic is basically paralyzed. Thus, 30 cm and 70 cm were set to be the critical heights of different flooding risk levels.

Risk Level Water Depth (cm) Extent of the Threats to Different Means of Transportation
There is little impact on traffic. Ⅱ 3-10 The water depth is higher than the insteps, thus slowing down pedestrians. Ⅲ 10-30 When the rail surface water depth exceeds 30 cm, electric vehicles must be shut down. Ⅳ

30-70
The water depth is higher than the outlet, thus impeding motor vehicles. Ⅴ >70 The water depth is higher than the inlet port, basically paralyzing traffic.

Effect of Land Use Change
The distributions of flooding areas simulated using different land use in 1987, 1997, 2007, and 2017 are shown in Figure 6. The results were obtained directly from SWMM simulations, which used a two-year rainfall event and the drainage system in 2017. The distributions of flooding water depths are also shown in Figure 6. The high-risk areas corresponding to level V were mainly located in the center of Handan City in 1987, then grew to the north and south in 1997, and spread to the whole study area in the subsequent 20 years. Increasing areas of waterlogging shared the same growing region of the impervious surface ( Figure 2).
As shown in Table 6, the areas with flood depths of more than 3 cm were positively correlated with the percentage of the impervious surface. The total flood volume exhibited the same trend as the water depth. The flooding volume is an indicator of the destruction of municipal infrastructure [54], and the locations of flood nodes represent high-risk areas. With the proportion of impervious surface increasing by 136.11% every 10 years, the total flood volume increased by 387.87%, from 7.89 × 10 4 m 3 to 3.85 × 10 5 m 3 . The land use The distributions of flooding water depths are also shown in Figure 6. The high-risk areas corresponding to level V were mainly located in the center of Handan City in 1987, then grew to the north and south in 1997, and spread to the whole study area in the subsequent 20 years. Increasing areas of waterlogging shared the same growing region of the impervious surface ( Figure 2).
As shown in Table 6, the areas with flood depths of more than 3 cm were positively correlated with the percentage of the impervious surface. The total flood volume exhibited the same trend as the water depth. The flooding volume is an indicator of the destruction of municipal infrastructure [54], and the locations of flood nodes represent high-risk areas. With the proportion of impervious surface increasing by 136.11% every 10 years, the total flood volume increased by 387.87%, from 7.89 × 10 4 m 3 to 3.85 × 10 5 m 3 . The land use change, represented by the expansion of the impervious surface, exerts a considerable impact on the spatial distribution of waterlogging and the volume of flooding.  Figure 7 shows the spatial variations in waterlogging using the drainage system. As detailed in Section 3.1, the results were obtained from SWMM simulations with biennial rainfall events and land use situations in 2017. The deep waterlogging areas occupied most of Handan City in 1987 but shrunk to the central zone in 2017, corresponding to the gradual reduction in areas with high-risk levels. The results reveal that a more developed drainage system may improve flood conditions in most areas, especially those located far from rivers. As more and more new pipelines transport water to rivers, the distribution of floods in the study area has become more concentrated near the outlets.   Figure 7 shows the spatial variations in waterlogging using the drainage system. As detailed in Section 3.1, the results were obtained from SWMM simulations with biennial rainfall events and land use situations in 2017. The deep waterlogging areas occupied most of Handan City in 1987 but shrunk to the central zone in 2017, corresponding to the gradual reduction in areas with high-risk levels. The results reveal that a more developed drainage system may improve flood conditions in most areas, especially those located far from rivers. As more and more new pipelines transport water to rivers, the distribution of floods in the study area has become more concentrated near the outlets. The changes in different water depth areas in Table 7 show that, with the improvement in drainage conditions, the shallowest and deepest areas increased rapidly, from 44.24% to 75.24% and from 0.12% to 1.08%, respectively. Additionally, the areas with the deepest water depths are distributed at riverside areas or pipe outlets. As shown in Table  7, the total flood volume is negatively correlated with the total pipe length. With the total pipe length increasing by 151.63% to 251.05 km, the total flood volume decreased by 27.56% from 5.31 × 10 5 m 3 . The results reveal that a more effective drainage system may The changes in different water depth areas in Table 7 show that, with the improvement in drainage conditions, the shallowest and deepest areas increased rapidly, from 44.24% to 75.24% and from 0.12% to 1.08%, respectively. Additionally, the areas with the deepest water depths are distributed at riverside areas or pipe outlets. As shown in Table 7, the total flood volume is negatively correlated with the total pipe length. With the total pipe length increasing by 151.63% to 251.05 km, the total flood volume decreased by 27.56% from 5.31 × 10 5 m 3 . The results reveal that a more effective drainage system may reduce the total flood volume, although it exacerbates the flooding risk at riverside areas and the pipe outlet.  Figure 8 shows the changes in flood areas of different land use structures and drainage systems over the past three decades. The SWMM uses biennial rainfall events to simulate the amount of flooding at each node. When combining the effects of land use change with drainage network change, the study results show that the flood conditions have deteriorated, especially in the central zones, and most of the flooding nodes are located in areas with impervious surfaces. This result is partially similar to a previous study [55], which demonstrated that the flood volume is more sensitive to urban development and drainage construction in central urban zones. The outlet in the SWMM is a node through which water can flow away from a pipeline to a river. Due to poor outlet conditions, flooding near the Fuyang River is growing worse in the eastern region. The topography of Handan is relatively gentle and could not meet the requirements for a drainage slope; therefore, the outlet of the drainage pipelines may be submerged by river floods during heavy rainfall. Under such conditions, the drainage system cannot discharge into the river; thus, the water is accumulated around the riverside where the outlet cannot transfer water to the river.   Figure 8 shows the changes in flood areas of different land use structures and drainage systems over the past three decades. The SWMM uses biennial rainfall events to simulate the amount of flooding at each node. When combining the effects of land use change with drainage network change, the study results show that the flood conditions have deteriorated, especially in the central zones, and most of the flooding nodes are located in areas with impervious surfaces. This result is partially similar to a previous study [55], which demonstrated that the flood volume is more sensitive to urban development and drainage construction in central urban zones. The outlet in the SWMM is a node through which water can flow away from a pipeline to a river. Due to poor outlet conditions, flooding near the Fuyang River is growing worse in the eastern region. The topography of Handan is relatively gentle and could not meet the requirements for a drainage slope; therefore, the outlet of the drainage pipelines may be submerged by river floods during heavy rainfall. Under such conditions, the drainage system cannot discharge into the river; thus, the water is accumulated around the riverside where the outlet cannot transfer water to the river.  As shown in Table 8, the flooding volume increases slightly year by year, from 3.31 × 10 5 m 3 to 3.85 × 10 5 m 3 . The proportion of area with a water depth above 70 cm increased from 0.11% to 1.08%, and that between 30 cm and 70 cm increased from 0.41% to 0.93%. The changes in the total flooding volume over the years indicate that the construction of drainage systems is not consistent with the land use changes. Meanwhile, the spatial distributions of waterlogging indicate that the improved drainage system may trigger more severe flooding over impervious surfaces and the riverside due to rapidly draining pipes and poor outlet conditions.

Urban Flooding Volume Distribution Pattern Analysis
The percentage of impervious surface and length of drainage pipe represents the land use and the drainage system, with sensitivity indices of 1.09 and −1.08, respectively. The sensitivity of land use change parameters is slightly higher than the drainage pipe change parameters, but both are highly sensitive. This result is consistent with the results of Li et al. [56], who observed that the N-Imperv (percent impervious %) and Manning's roughness coefficient of pervious and impervious areas had a greater effect on the total runoff volume. The development of drainage systems contributes similarly to the land use pattern, which means that more effort is needed for the overall drainage system development.
In order to describe the flooding spatial distribution patterns shown in Figures 6-8, the correlations of imperviousness, total pipe lengths, and flooding volume are plotted in Figure 9. It can be seen in Figure 9a that as the imperviousness of the city increases from 32.43% to 76.57%, the flood volume in the city increases from 7.8 × 10 4 m 3 to 3.85 × 10 5 m 3 , and there is a positive correlation between impervious surface and total flood volume; the linear correlation coefficients reached 0.9892. Changes in land use type and imperviousness have direct impacts on surface runoff [57]. As shown in Table 8, the flooding volume increases slightly year by year, from 3.31 × 10 5 m 3 to 3.85 × 10 5 m 3 . The proportion of area with a water depth above 70 cm increased from 0.11% to 1.08%, and that between 30 cm and 70 cm increased from 0.41% to 0.93%. The changes in the total flooding volume over the years indicate that the construction of drainage systems is not consistent with the land use changes. Meanwhile, the spatial distributions of waterlogging indicate that the improved drainage system may trigger more severe flooding over impervious surfaces and the riverside due to rapidly draining pipes and poor outlet conditions.

Urban Flooding Volume Distribution Pattern Analysis
The percentage of impervious surface and length of drainage pipe represents the land use and the drainage system, with sensitivity indices of 1.09 and −1.08, respectively. The sensitivity of land use change parameters is slightly higher than the drainage pipe change parameters, but both are highly sensitive. This result is consistent with the results of Li et al. [56], who observed that the N-Imperv (percent impervious %) and Manning's roughness coefficient of pervious and impervious areas had a greater effect on the total runoff volume. The development of drainage systems contributes similarly to the land use pattern, which means that more effort is needed for the overall drainage system development.
In order to describe the flooding spatial distribution patterns shown in Figures 6-8, the correlations of imperviousness, total pipe lengths, and flooding volume are plotted in Figure 9. It can be seen in Figure 9a that as the imperviousness of the city increases from 32.43% to 76.57%, the flood volume in the city increases from 7.8 × 10 4 m 3 to 3.85 × 10 5 m 3 , and there is a positive correlation between impervious surface and total flood volume; the linear correlation coefficients reached 0.9892. Changes in land use type and imperviousness have direct impacts on surface runoff [57]. In Figure 9b, it can be seen that as the total pipe length increases from 99.1 km to 251.05 km, the total flood volume decreases from 5.31 × 10 5 m 3 to 3.85 × 10 5 m 3 , which shows a negative correlation; the linear correlation coefficients reached 0.8654. This result is different from previous studies [55]; however, we believe that this finding is reasonable because, with the construction of urban pipe networks, stagnant urban water will drain In Figure 9b, it can be seen that as the total pipe length increases from 99.1 km to 251.05 km, the total flood volume decreases from 5.31 × 10 5 m 3 to 3.85 × 10 5 m 3 , which shows a negative correlation; the linear correlation coefficients reached 0.8654. This result is different from previous studies [55]; however, we believe that this finding is reasonable because, with the construction of urban pipe networks, stagnant urban water will drain more smoothly, resulting in a decrease in the flood volume. Another possible reason is that the effects of land use changes and urban pipe networks on urban stormwater are discussed separately in this paper using the control variables approach.
In Figure 9c, the effects of urban land use change and urban pipe network construction on urban stormwater are not linear, and the combined effect of both increases the amount of urban waterlogging year by year, although the magnitude of change is significantly reduced compared with the single-factor case, especially between 2007 and 2017. The increase is very slight, which indicates that the aggressive expansion period of urban construction has ended, and urban construction from the perspective of water security has become more scientific and orderly.
For a better understanding of the information in Figure 8, all statistics of each subcatchment for each year and the flooding volume were plotted. Figure 10 shows that the percentage of impervious area of each sub-catchment, the area of each sub-catchment, and the drainage network condition exert a considerable impact on the flood volume. With the land use development, the impervious area of each sub-catchment appears to be rising, and with the drainage pipeline network construction, the area of each sub-catchment reveals a generally decreasing trend. Waterlogged sub-catchments of more than 1.5 × 10 4 m 3 appeared in 1987 and 1997 due to the larger area of each catchment. The flood amount for each waterlogged sub-catchment has fallen during the last 30 years due to the improved drainage network conditions. Nonetheless, this does not signify safer water security conditions because the total amount increases each year. more smoothly, resulting in a decrease in the flood volume. Another possible reason is that the effects of land use changes and urban pipe networks on urban stormwater are discussed separately in this paper using the control variables approach. In Figure 9c, the effects of urban land use change and urban pipe network construction on urban stormwater are not linear, and the combined effect of both increases the amount of urban waterlogging year by year, although the magnitude of change is significantly reduced compared with the single-factor case, especially between 2007 and 2017. The increase is very slight, which indicates that the aggressive expansion period of urban construction has ended, and urban construction from the perspective of water security has become more scientific and orderly.
For a better understanding of the information in Figure 8, all statistics of each subcatchment for each year and the flooding volume were plotted. Figure 10 shows that the percentage of impervious area of each sub-catchment, the area of each sub-catchment, and the drainage network condition exert a considerable impact on the flood volume. With the land use development, the impervious area of each sub-catchment appears to be rising, and with the drainage pipeline network construction, the area of each sub-catchment reveals a generally decreasing trend. Waterlogged sub-catchments of more than 1.5 × 10 4 m 3 appeared in 1987 and 1997 due to the larger area of each catchment. The flood amount for each waterlogged sub-catchment has fallen during the last 30 years due to the improved drainage network conditions. Nonetheless, this does not signify safer water security conditions because the total amount increases each year. As Figures 9 and 10 show, the flood volume increased greatly along with urbanization from 1987 to 2007. However, the flood volume increased with the same urbanization speed from 2007 to 2017. This indicates that the strategy of sustainable development in Handan has worked well in the field of urban planning, urban construction, and urban flood security.
Several limitations are acknowledged in this study. Although the land use development and drainage system construction data were obtained from multiple sources (e.g., remote sensing data and local drainage conditions) to ensure the reliability of the results, the simulation still has limitations due to a lack of input data (e.g., spatial distribution and hydrological conditions of the soil, high-resolution DEM data, hydrological measurements of ground and underground drainage systems, and operation conditions of pumps). The availability of soil spatial distribution is an important factor in choosing a As Figures 9 and 10 show, the flood volume increased greatly along with urbanization from 1987 to 2007. However, the flood volume increased with the same urbanization speed from 2007 to 2017. This indicates that the strategy of sustainable development in Handan has worked well in the field of urban planning, urban construction, and urban flood security.
Several limitations are acknowledged in this study. Although the land use development and drainage system construction data were obtained from multiple sources (e.g., remote sensing data and local drainage conditions) to ensure the reliability of the results, the simulation still has limitations due to a lack of input data (e.g., spatial distribution and hydrological conditions of the soil, high-resolution DEM data, hydrological measurements of ground and underground drainage systems, and operation conditions of pumps). The availability of soil spatial distribution is an important factor in choosing a suitable infil-tration model, which is crucial for calculating flood volumes. With higher resolution, the DEM data and land use data provided detailed information for two-dimensional hydraulic models for a more realistic simulation of flood distributions. Additionally, the social media information was a faulty substitute for historical water depth observations and hydrological records to validate the simulation results. The lost information on accurate water depth, shooting time, and location could undermine the reliability of validation. Thus, the lack of measuring stations constrains the building and verifying model of urban flood simulation in rapidly developing cities in recent years in China.

Conclusions
In response to the question of how the dual changes in land use and drainage systems affect the distribution of flooding, this study interpreted remote sensing images of Handan City from the past three decades. Then, the dual changes in hydrological conditions caused by ground and underground development were investigated, a stormwater simulation model was developed, and the results were verified using social media data in order to identify the combined effects of spatial variations in waterlogging and volume.
In this case study of Handan City, against the backdrop of rapid economic growth, there has been a boom in both ground and underground construction. In addition, land use change, represented by the expansion of sealed surfaces, has exerted a positive impact on the distribution and amount of flooding. While the impermeable area has increased from 29.74 km 2 to 70.22 km 2 , the proportion of area susceptible to submergence under water depths over 70 cm has increased from 0.55% to 1.08%, and the total volume of flooding has increased from 7.89 × 10 4 m 3 to 3.85 × 10 5 m 3 . Furthermore, the construction of the drainage network has increased from 99.77 km to 251.05 km, and the flooding volume has reduced from 5.31 × 10 5 m 3 to 3.85 × 10 5 m 3 ; however, the proportion of area susceptible to submergence under water depths over 70 cm has increased from 0.12% to 1.08%. The dual changes in land use and drainage pipeline network have increased the total flooding amount from 3.31 × 10 5 m 3 to 3.85 × 10 5 m 3 . Finally, the proportion of areas susceptible to submergence under water depths of more than 70 cm has increased from 0.11% to 1.08%, while those susceptible to submergence between 30 cm and 70 cm of water grew from 0.41% to 0.93%. These joint effects show that the drainage system can mitigate the effect of land use change, although at the same time also aggravate the flooding conditions in impervious and riverside areas, especially in central zones.
The main findings of the research are as follows: (1) both surface and underground construction are highly sensitive to the distribution of flooding; (2) the land use change is a positive factor influencing the total flood volume, whereas the drainage system is negative, and the waterlogging risk increases in impervious areas and riverside areas of central zones; (3) as discussed previously, the simulation results still come with uncertainties due to condition simplification and data limitations in this study. The results of this study suggest that sustainable planning for drainage systems in rapidly urbanizing cities requires greater consideration of the long-term developments of city-wide land use, which necessitates cross-sectoral and cross-region cooperation.