Variation of uncertainty of drainage density in flood hazard mapping assessment with coupled 1D–2D hydrodynamics model

Coupled 1D–2D hydrodynamic models are widely utilized in flood hazard mapping. Previous studies adopted conceptual hydrological models or 1D hydrodynamic models to evaluate the impact of drainage density on river flow. However, the drainage density affects not only river flow, but also the flooded area and location. Therefore, this work adopts the 1D–2D model SOBEK to investigate the impact of drainage density on river flow. The uncertainty of drainage density in flood hazard mapping is assessed by a designed case and a real case, Yanshuixi Drainage in Tainan, Taiwan. Analytical results indicate that under the same return period rainfall, reduction in tributary drainages in a model (indicating a lower drainage density) results in an underestimate of the flooded area in tributary drainages. This underestimate causes higher peak discharges and total volume of discharges in the drainages, leading to flooding in certain downstream reaches, thereby overestimating the flooded area. The uncertainty of drainage density decreases with increased rainfall. We suggest that modeling flood hazard mapping with low return period rainfalls requires tributary drainages. For extreme rainfall events, a lower drainage density could be selected, but the drainage density of local key areas should be raised.


Introduction
In coupled one-dimensional and two-dimensional hydrodynamics models (1D-2D model), the channel flows are simulated in a 1D model, and the floodplain flows are simultaneously simulated in 2D. Two-dimensional hydrodynamics schemes include simple-volume conservative storage-filling algorithms, diffusive wave scheme and fully dynamic shallow water modelling. Since 1D-2D models have high accuracy and short calculation time, it is widely used in Meso (regional) scale and Micro (local) scale flood hazard assessment (Apel et al. 2009;de Moel et al. 2015;Werner 2004;Werner et al. 2005). Wilson et al. (2007) adopted 1D-2D model LISFLOOD-FP and topographic data from the Shuttle Radar Topography Mission to derive the inundation of Amazonian seasonally flooded wetlands. Price and Vojinovic (2008) proposed the use of the concept of digital city to manage urban flood disasters, based on the tropical island of St Maarten, one of five land area of the Netherlands Antilles. The flood hazard maps were simulated with the flooding model MIKE 11 package developed by DHI Water and Environment. Timbadiya et al. (2015) used a 1D-2D model (MIKE 11 and MIKE 21) to forecast the water level of lower Tapi River in India. The flood hazard maps in Tainan, Taiwan were simulated using the 1D-2D model SOBEK for flood risk management (Doong et al. 2016). Wu et al. (2017) applied a 1D-2D model, coupling the SWMM and LISFLOOD-FP models, to predict the future flooding scenarios in Dongguan, China. Yang et al. (2018aYang et al. ( , 2018b) used the SOBEK model to evaluate the flood risk transfer effects due to land development in lowlands.
Several sources of uncertainties in 1D-2D models have been widely studied, analyzed and summarized. The main sources of uncertainties for the models include model selection and parameter settings, as well as hydrological and geological data of the models (Bales and Wagner 2009;de Moel et al. 2015;Merwade et al. 2008;Teng et al. 2017). Werner (2004) compared 1D (SOBEK-RIVER), 2D (DELFIT-FLS) and 1D-2D (SOBEK-OVER-LAND FLOW) models to predict flood stages of the River Saar in Germany. The research results demonstrated 1D and 1D-2D models generated more reliable simulation results than 2D model. Werner et al. (2005) compared SOBEK 1D model and coupled SOBEK 1D-2D model to predict floodplain inundation in the towns of Usti and Orlici in Czech Republic. Their research results show that the simulation results of 1D model and 1D-2D model in the flooded area were similar. However, the 1D-2D model showed better performance in comparison of distributed water level observations. Apel et al. (2009) compared linear interpolation methodology, 1D-2D model (LISFLOOD-FP) and 2D models for the risk analysis in the municipality of Eilenburg in Saxony, Germany. The simulation results of the 1D-2D and 2D models were highly consistent with the observed flooding extents and depths. However, the long calculation time required by the 2D model was a major obstacle in the calibration process. Dimitriadis et al. (2016) employed a Monte-Carlo approach to analyze the uncertainty of input discharge, longitudinal and lateral gradients, roughness coefficient and the grid size in three 1D and quasi-2D models (HEC-RAS, LIS-FLOOD-FP, and FLO-2d). Their research result revealed that the channel and floodplain friction, and inflow discharge were the main uncertainties in flood propagation.
Another major source of model uncertainty is the model input, which include both hydrological and geological data. Many researchers have examined the design rainfall, inflow hydrograph, lateral discharge, channel and floodplain geometry, as well as the initial condition of soil moisture. Wu et al. (2010) developed a risk analysis model, combined with multivariate Monte Carlo simulation and Advance First-Order Second-Moment method. The model was applied to assess the risk of underestimating flood peak flows due to the uncertainties in rainfall information (rainfall depth, duration and storm pattern), and the uncertainties in the parameters of the rainfall-runoff model (Sacramento Soil Moisture Accounting model, SAC-SMA). Wu et al. (2011) presented a risk analysis model to consider the risk of damage to the flood protection structure of the Keelung River due to uncertainties in hydrological and hydraulic analysis, including design rainfall event, freeboard of level and the flood-diversion channel. Wong et al. (2015) used a 1D-2D model (LISFLOOD-FP) to study the effect of river channel cross-section and longitudinal-section geometry on flood dynamics. Their results indicated that uncertainty on channel longitudinal-section variability only affected the local flood dynamics, but did not significantly affect the friction sensitivity or inundation mapping. These effects were negligible, as they were less important than other uncertainties, such as boundary conditions. Parameter values and initial condition impact a hydrological model results. Pradhan et al. (2020) built a distributed hydrological model of the sub-basin of Reynolds Creek Experimental Watershed using the GSSHA (Gridded Surface Sub-surface Hydrological Analysis) model, and evaluated the impact of different resolutions of satellite imagery-based soil moisture on the discharge simulation.
The grid resolution and the post-processing quality significantly influence the flood simulation results. Adeogun et al. (2015) proposed a 1D-2D model, coupling 1D sewer model (SWMM) and 2D inundation model, and analyzed the sensitivity of the mesh resolution and roughness. The research results revealed that a higher-resolution mesh enhanced the flooding simulation results, but lengthened the calculation time. Noh et al. (2018) developed a hybrid parallel code (H12) of the 1D-2D model to accelerate the calculation speed of the flood simulations with a hyper-resolution grid. Their research results showed that the hyper-resolution grid modeling accurately described the flooding situation in urban areas, while the coarser resolution grid modeling led to local isolation and distortion of the flooded area, due to the lack of topographical details. Meesuk et al. (2015) presented Multidimensional Fusion of Views-Digital Terrain Model (MFV-DTM) to integrate groundview Structure from Motion (SfM) observations with top-view LiDAR data. A 1D-2D model was applied to simulate an extreme urban flood event that occurred on June 10, 2003 in Kuala Lumpur, Malaysia. The research results demonstrated that MFV-DTM modeling could represent a true flood situation better than the standard LiDAR-DTM modeling and Filtered LiDAR-DTM modeling.
The drainage density in the model influences the flood simulation results for flood hazard mapping. Yildiz (2004) used a physically based spatially distributed hydrological model to study the sensitivity of flow simulation to the drainage density, taking the Monongahela River Basin in the United States as an example. Pallard et al. (2009) took the Po River in Italy as an example, and applied a spatially distributed hydrological model to investigate the relationship between the drainage density and the mean, standard deviation coefficient of variation and coefficient of skewness of the annual maximum peak flow. Ogden et al. (2011) and used the GSSHA model to evaluate the impact of impervious area, drainage density, width function and subsurface storm drainage on peak flood flows using the Dead Run watershed in Lanzhou State, USA as a case. Those studies, which were based on the concept hydrological model or 1D model, found that a larger drainage density leads to an increase in peak discharge. However, adjusting the drainage density in a flood simulation alters not only discharge, but also the flooded area and location. This issue has not been discussed in previous studies. Therefore, this investigation designed two experiments using the SOBEK model to measure the effect of drainage density on two-dimensional flooding simulation.

SOBEK model
The hydrological and hydraulic model was developed with the SOBEK model (version 2.13) developed by Deltares in the Netherlands. The SOBEK model has several modules, including Rainfall-Runoff, 1D FLOW-Rural and Overland Flow-2D (Deltares 2017). The flooding simulation process in this study has two phases. The first phase is 1 3 the hydrological phase, in which the Rainfall-Runoff module is run to convert rainfall into runoff, and then the sub-catchment discharge is computed. The second phase is the hydraulic phase. After the runoff is introduced into the channel, the 1D-2D module is run to perform the calculation with complete Saint Venant equations. For one-dimensional and two-dimensional flow, continuity and momentum equations are calculated numerically using the Delft scheme. When the discharge exceeds the flood capacity of the channel, the water overflows the channel, and then the dynamic flow on the surface is simulated through the 2D module.
This study employed the US Soil Conservation Service Curve Number method to compute abstractions. In the rainfall-runoff model, the US Soil Conservation Service's dimensionless unit hydrograph was employed to transform rainfall into sub-catchment discharge. The SCS CN Method is described as follows.
where P e denotes rainfall excess (mm); P represents depth of precipitation (mm); S indicates potential maximum retention (mm), and CN is the curve number. Lag time is computed according to the SCS Lag Equation, shown as follows: where T 1ag denotes lag time (h); L represents length of the longest drainage path (m); H indicates average catchment slope (%); S is the potential maximum retention (mm), and CN denotes the curve number. The SCS Dimensionless Unit Hydrograph Model is adapted as the Rainfall-Runoff model to estimate the discharge.
where T lag represents the lag time of peak flow (h); T c indicates the time of concentration (h); T p is the time to peak flow (h); T r denotes the recession limb time (h); T b represents the runoff time base (h); Q p indicates the quantity of peak flow (cms); A is the area of catchment (km 2 ), and R denotes the quantity of effective rainfall (mm). (1)

Experiment I
Five different models were designed (as shown in the Table 1). Model A included Drainages 2-years, 5-years, 10-years, 50-years, and 100-years. Model B included Drainages 5-years, 10-years, 50-years, and 100-years, and so on (as shown in Fig. 1). A 2-years return period flood could be drained by Drainage 2-years. A 5-years return period flood could be drained by Drainage 5-years, and so on. Cross-sections of the designed channel are shown in Table 2.
Each model was simulated with 2-years, 5-years, 10-years, 50-years, and 100-years return period rainfalls. In each scenario, the flooded area, as well as the peak discharge (Q A ) and the total volume of discharge at Point A (V A ) were calculated and compared.

Study area
This study used the catchment of Yanshuixi Drainage in Tainan, Taiwan as a case study. The Taiwan Strait, Zengwun River and Yanshuei River are situated in the west, north and south of the catchment, respectively. The catchment area is 108.454 km 2 , and the terrain is inclined from northeast to southwest. The elevation of the catchment area is between Elevation Level (E.L.) − 0.2 m and 9.5 m. The drainage length is 19.20 km, and the average slope is about 1/7,000. The river export is about 315 m wide, and most of the river bank is concrete revetment (see Fig. 2). The land usages in the catchment of Yanshuixi Drainage are agriculture (42.4%), forest (1.8%), traffic (9.9%), water (5.9%), build-up (19.3%), public (2.0%), recreation (3.0%), mining (0.1%) and other (15.6%). Due to the flat and low terrain, the drainage outlets are affected by the tidal level, and typhoons and heavy rains often cause floods downstream.

Model parameters
The DEM with an elevation accuracy of 10-20 cm was output from the Airborne LiDAR point cloud run by the Ministry of Interior, Taiwan. The grid resolution of DEM was set to 20 m in the experiments. The parameters of the flow path and the slope in the catchment were estimated from DEM, and the CN was estimated from land use. The CN of each sub-catchment area was based on weighted calculation of different land use areas. The surface friction parameters (White-Colebrook K n ) were derived according to land use. The SOBEK model was configured with regional drainages, farmland drainages, and related hydraulic structures, such as gates, pumping stations, detention basins and bridges. The channel cross-sections were set according to the measurement data. The channel friction parameters (Manning n) were set according to the channel type.
The exit boundary condition of the model is based on the tide level at the Sicao tide station. Observed rainfall data from the Quantitative Precipitation Estimation and Segregation Using Multiple Sensor (QPESUMS) are obtained from 6 Dopper radars of the Central Weather Bureau. The interval of the weather information was 10 min, and the grid resolution was 1.3 km (Chiou et al. 2005;Wang et al. 2016).

Model calibration and verification
The water level data in the Yanshueixi Drainage were accumulated from three Tainan City Government water level stations (see Fig. 2), which measure the water level using ultrasonic sensors with an accuracy of 2.5%. The degree of bias between the observed and simulated water depths was measured by coefficient of efficiency (CE). The equation of CE is shown below: where H sim (t i ) denotes the simulated water depth (m) at time t; H obs (t i ) represents the observed water depth (m) at time t; H obs t i indicates the average observed water depth (m). Storm 0611 in 2016 and Tropical Storm Haitang in 2017 were utilized for model calibration and verification, respectively. The CN value and overland roughness k n of model calibration were the same as those of model verification. Channel roughness n was adjusted for model calibration. Figure 3 depicts the comparison between the observed and simulated water levels of three water level stations. The CE values of Storm 0611 were 0.847, 0.994, and 0.960 at Station 1, 2, and 3, respectively. The CE values of Tropical Storm Haitang were 0.916, 0.996, and 0.983 at Station 1, 2, and 3, respectively. It indicated that the model accuracy was high. Figure 4 illustrates the calibrated CN of sub-catchment, overland roughness kn and channel roughness n.

Drainage density
In this study, drainage density (D d ) is the total length of all the drainages in a watershed divided by the total area of the watershed.
where L is total length of all the drainages in a watershed (km); A is total area of the watershed (km 2 ). Two types of drainages, regional drainage and farmland drainage, were set in models. The total length of regional drainages in Model I was 75.370 km. In Model II, in addition to regional drainages, 50.143 km of farmland drainages was added. Table 3 lists the channel densities of Model I and II, which are 0.695 km −1 , and 1.157 km −1 , respectively.

Flooded areas comparison
The differences of simulated flooded areas between Model I, and II were evaluated using three indices, defined as follows. Figure 5 illustrates the schematic diagram. where

Experiment I
Tables 4, 5, and 6 depict the results of Experiment I, and Fig. 6 display the simulated flooded areas under 2-100-years return period rainfall. Under the 2-year return period rainfall, Drainage 2-years to Drainage 100-years had greater flood capacities than the 2-year return period. Therefore, no flooding occurred in Model A to Model E, and the Q A and V A of all models were 856 cms and 13,795,000 m 3 . Under the 5-year return period rainfall, only Drainage 2-years had a lower flood capacity than the 5-year return period flood. Since Model A contained Drainage 2-years, 25 ha of flooding occurred. The other models did not include Drainage 2-years, so no flooding occurred. At the same time, the Q A and V A in Model A were 2111 cms and 38,164,000 m 3 , which were smaller than those of other models (2114 cms and 38,199,000 m 3 ). Since Model A contained 25 ha of flooding, the flood water was retained on the surface, resulting in a decrease in drainage.
Under the 10-year return period rainfall, the 10-year return period flood was greater than the flood capacities of Drainage 2-years and Drainage 5-years. In Model A, 300 ha of flooding occurred in Drainage 2-years and Drainage 5-years. Since Model B contains Drainage 5-years, 212 ha of flooding occurred. Conversely, Drainage 2-years and Drainage 5-years were not set in other models, so no flooding occurred.  Because SOBEK is a 1D-2D model, the discharges from Rain-Runoff model are input into the drainages first. If the discharge is less than the flood capacity of the drainage, then flooding does not occur. If the discharge is greater than the flood capacity of the drainage, then flooding occurs. Therefore, the reduction in tributary drainages results in two effects on the simulation results of the 1D-2D model.
Effect 1: Under the same return period rainfall, fewer tributary drainages in a model (lower drainage density) indicates a smaller flooded area.
Effect 2: Under the same return period rainfall, fewer tributary drainages in a model (lower drainage density) indicates an increased peak discharge and total volume of discharge in the drainages.
Previous works have noted that a larger drainage density leads to an increase in flow rate (Ogden et al. 2011;Pallard et al. 2009;Yildiz 2004). However, our research results indicate that an increase in drainage density does not necessarily lead to an increase in flow. Previous investigations adopted the spatially distributed hydrological model or the GSSHA model, and therefore determined surface and underground discharges are drained to the outlet of the catchment. The increase in channels shortens of the concentration time, causing the peak flow rate to rise with the increased drainage density. However, not all runoff is drained to the catchment outlet in the 1D-2D model. A simulation scenario with a high drainage density is modeled with dense and narrow upstream channels. Floodwater is retained in the upstream when the upstream channels are flooded, reducing discharge from outlets.

Experiment II
The analysis results are shown in Table 7, Figs. 7, 8, and 9. Table 7 shows that in all simulated rainfall scenarios, Model II have larger flooding areas than Model I, and A(I) is approximately 83.9-96.9% of A(II). However, P(I ∩ II) is only 49.7-66.1%. Thus Models I and II have similar total flooded areas, but significantly different flooded locations, and this difference increases with an increase in rainfall.
The A(II − I) of 2-years, 5-years, 10-years, 50-years, and 100-years return period rainfalls were 86.32 ha, 162.16 ha, 329.56 ha, 499.80 ha and 672.32 ha, respectively. However, P(II − I) of 2-years, 5-years, 10-years, 50-years, and 100-years return period rainfalls were 31.2%, 23.7%, 27.4%, 18.2% and 19.8%, respectively (as shown in Table 7). A(II − I) occurred because farmland drainages were set up in Model II, but not in Model I, so some flooding that occurred in farmland drainages could not be shown in Model I. This caused an underestimate of the flooded area in Model I (Effect 1 had been proved in experiment I). The flood capacity of farmland drainage is almost 2-5-years return period, and the flood capacity of regional drainage is 10-years return period. Under 2 return period rainfall, the majority of flooding occurred in farmland drainages, and P(II − I) = 31.2%. However, the P(II − I) values during the 50-years or 100years return period rainfalls were only 18.2% and 19.8%, respectively, because most of the discharges far exceeded the drainage capacities of regional drainages and farmland drainages. In general, A(II-I) increased as the total rainfall rose, and P(I − II) decreased significantly as the total rainfall rose (see Figs. 8 and 9).
The A(I − II) of 2-years, 5-years, 10-years, 50-years, and 100-years return period rainfalls were 53.12 ha, 128.68 ha, 162.52 ha, 428.72 ha, and 512.88 ha, respectively. However, P(I − II) of 2-years, 5-years, 10-years, 50-years, and 100-years return period rainfalls were 19.2%, 18.8%, 13.5%, 15.6% and 15.1%, respectively (see Table 7). A(I − II) occurred because farmland drainages were set up in Model II, but not in Model I, so some flooding in farmland drainages could not be shown in Model I. Therefore, the peak discharges and volumes of discharges in drainages in Model I were greater than those in Model II (Effect 2 was shown in Experiment I). In Model I, the higher discharges in drainages caused the overflow occurrence at certain downstream reaches, resulting in an overestimate of the flooded area. In general, A(I − II) increased as the total rainfall rose, and P(I − II) declined slightly with increasing total rainfall (see Figs. 8 and 9). Overall, P(I − II) and P(II − I) decreased with increasing rainfall, and P(I ∩ II) rose with rising rainfall (see Fig. 9), indicating that the uncertainty of drainage density declined with rising rainfall.
These experimental results indicate that the drainage density in the 1D-2D model had a direct impact on the accuracy of the flood hazard mapping. The analytical results reveal a model with higher drainage density was more similar to the real terrain, and generated more accurate simulation results. However, obtaining detail topographical data of tributary drainages is difficult, expensive and time-consuming. Therefore, in the flooding hazard mapping, we recommend evaluating the drainage density with reference to two factors, the simulated rainfall scenario and the local key areas.
In the model, the determination of drainage density is directly related to the simulated rainfall situation. A higher drainage density should be selected for a lower simulated rainfall level. Conversely, a lower drainage density could be selected given a higher simulated rainfall level. Flood hazard mapping is generally simulated with extreme historical or designed rainfall events, and the models are typically calibrated and verified by the hydrological stations located in rivers (Timbadiya et al. 2015;Wilson et al. 2007).
In this case, the models are normally set up in with only vertical and horizontal sections of the mainstreams. When the extreme flow is greater than the flood discharge capacity of the river, large-scale flooding occurs on both sides of the mainstream, causing flooding in low-lying areas. Conversely, a model for the flood hazard mapping with low return period rainfall, such as 2-years, 5-years, and 10-years should establish tributary drainage, such as regional drainage or farmland drainage, otherwise it might underestimate the simulated flooded area, and skew the location of flooded area. We recommend raising the drainage density of local key areas in the model, such as important preservation objects and settlements, in the process of flood hazard mapping. Since the flood hazard map is an important basis for the subsequent decision of disaster prevention measures, the accuracy of flooding simulation in some key areas is very important. Therefore, if the drainage density of the entire survey area cannot be increased owing to limitations of funding and data acquisition, then the drainage density can be increased in key regions.
Another method of reducing the uncertainty of drainage density is to perform hydrodynamic calculations on high-resolution two-dimensional DEM data using the 2D model. This method is computationally intensive, so it should be combined with other high-speed computing technologies, such as GPU (Morales-Hernández et al. 2021).
The latest technology can be employed to improve the accuracy of the flood hazard mapping. A large number of tributary drainage data need to be collected, saved and managed in order to improve the accuracy of flooding simulation. The management platform of the basic data for the model is very important. Performing surveys manually in the field is quite time-consuming and expensive. The digital terrain data can be obtained using drones, and processed through appropriate post-processing technology (Meesuk et al. 2015). As well as the calibration and verification of the model, distributed flooding observation data are also required to improve the accuracy of the model, and especially to confirm the accuracy of the flooded area. Smart water level gauges and Internet of Things technology can be efficiently used to obtain a large number of continuous and accurate flooding data (Chang et al. 2018). The application of drones for large-scale identification of range of flooding is also a significant direction for future research (Yang et al. 2020).

Conclusion
Drainage density has a significant impact on the accuracy of the 1D-2D model simulation, but has not been discussed in previous studies. Simplifying the tributary drainages has two main effects under the same rainfall scenario. Under the same return period rainfall, fewer tributary drainages in a model indicates a smaller flooded area (Effect 1), as well as an increased peak discharge and total volume of discharge in the drainages (Effect 2). In the case study of Yanshuixi Drainage, Effect 1 caused an underestimate of the flooded area, and Effect 2 caused an overestimate of the flooded area. The impact of these two effects on the simulation results gradually decreased as the rainfall increased. These characteristics can be utilized to select the appropriate drainage density. If the model is to simulate rainfall with a low return period, then the tributary drainage data cannot be simplified. However, for flood hazard mapping with extreme events, if only mainstreams are set in the model, then the tributary drainages at important locations can be raised to improve the accuracy of the simulation. Manual survey of tributary drainages is costly and time-consuming. For large-scale flood hazard mapping, airborne radar or drones can quickly obtain a large amount of tributary drainage data, making an important future research direction. Future studies will set the drainage density as a variable, so that the appropriate drainage density can be selected according to different rainfall scenarios. The function of flooding depth and flooded area can be established with artificial intelligence models to accelerate the estimation of flooding range and depth.