Abstract
Hydrodynamic models were commonly used for flood risk management in urban area. This paper presents initial efforts in developing an urban flood inundation model by coupling a one-dimensional (1D) model with a two-dimensional (2D) model to overcome the drawbacks of each individual modelling approach, and an additional module is used to simulate the rainfall-runoff process in study areas. For the 1D model, the finite difference method is used to discretize the Saint-Venant equations. An implicit dual time-stepping method (DTS) is then applied to a 2D finite volume model for an inundation simulation to improve computational efficiency. A total of four test cases are applied to validate the proposed model; its performance is demonstrated by a comparison with an explicit scheme and previously published results (an extensive physical experiment benchmark case, a vertical linking example, and two real drainage cases with actual topography). Results demonstrate that the proposed model is accurate and efficient in simulating urban floods for practical applications.
1. Introduction
Over the past decades, there have been an increasing number of extreme rainfall and flood events occurring globally in relation to climate change, which paralyze cities and result in serious social and economic consequences. In addition, small-scale urban flooding occurs frequently throughout the world in relation to poor drainage. Consequently, there is a growing need for storm water management strategies to reduce urban flood risk.
Hydrodynamic models are solid tools in urban storm water management. Numerical simulations of urban floods are important when scientifically planning and designing urban drainage systems and providing efficient urban flood disaster control and management strategies. Although many models have been developed for river and coastal flooding [1–3], urban flooding models have not yet been adequately developed; this partly attributes to the complex flow processes in urban areas when inundation occurs. As there is a one-dimensional (1D) flow property in sewer systems, the highly efficient one-dimensional (1D) model is the most commonly used tool to simulate the hydraulic performance of urban drainage systems. A number of sewer models are used widely, such as MOUSE [4], Infoworks CS [5], and Storm Water Management Model (SWMM) [6]. However, 1D models can only predict the surcharge volume from a drainage system in terms of the ground’s surface. Therefore, to accurately simulate detailed urban flood propagation and inundation of the ground’s surface, coupling of the 1D and 2D hydrodynamic model is necessary to consider flow interactions between underground pipes and the ground’s surface, and the 1D model and 2D model in the coupled model are used to simulate flows in the pipe/river drainage system and surface inundation, respectively. In recent years, several coupled 1D-2D models for urban floods have been developed and applied [7–10].
2D models are widely used to simulate urban flood inundation [11–14]. However, 2D calculation can be time consuming, and large variations in cell sizes may occur because of the complex urban underlying surface or local refinement, which leads to lower efficiency because the time step is limited by numerical stability and is largely determined by the smallest cells in the mesh. Although the time steps of 1D and 2D models are both limited by numerical stability, the time steps adopted in 1D models can usually be much larger than those of 2D models. However, most existing 1D-2D coupled models require the same time step in 1D and 2D models, which may result in lower efficiency of the 1D models. To improve the drawback of the inefficient 2D model, large-sized grids and simplified models are commonly adopted, which can then lead to a loss of accuracy.
The aim of this paper is to develop a coupled model for urban flood inundation by using an implicit dual time-stepping method (DTS) [15]. DTS is adopted to a 2D model to improve run time efficiency. The method is fully implicit, matrix-free, and easily implemented, meanwhile providing the possibility of parallel computing. With the application of the DTS, arbitrarily large time steps of the 2D model are permitted and computational efficiency of both 1D and 2D model is significantly improved. In addition, the coupled model has greater stability due to the use of an implicit scheme and is therefore particularly suitable for urban flood simulations where there is a complex underlying surface.
In this study, the 1D and 2D governing equations for free surface and pressurised flows are presented together with details of numerical schemes. The proposed model is then validated against four numerical test problems, and in the final section the study is summarized and concluded.
2. Numerical Schemes
2.1. Rainfall-Runoff Model
To simulate the rainfall-runoff process in urban areas, the study area is firstly divided into subcatchment areas; these represent land areas that collect rainwater and allow infiltration and drainage to a specific node. These subcatchments are then divided into impervious and pervious areas to enable the separate calculation of runoff, and the Horton infiltration method [16] is used to simulate the infiltration process in pervious area. Rainwater collected by the subcatchments finally drains to the drainage network, and the 1D model is used to calculate the flow routing in the drainage network.
2.2. 1D Model
The 1D model is used to simulate flows in the urban drainage system where free surface and pressurised flows coexist. Free surface flow can be expressed mathematically by the 1D shallow-water equation (SWE), and if a fictitious Preissmann slot is introduced, the same set of governing equations can be used for full pipe flow without a free surface, known as pressurised flows. The continuity and momentum equations of the 1D model are given bywhere denotes time; is the Cartesian coordinates; is discharge; is the wet cross-section area; is the water level and water head for free surface flow and pressurised flow, respectively; is gravity acceleration with a value of 9.81 m/s2; is the momentum correction factor; and is the friction slope.
Firstly, the governing equations are solved by the finite difference method using Preissmann’s four-point implicit difference scheme in individual channels or pipes. The discretization of Preissmann’s scheme is as follows: where and θ are weighting coefficients, and are the time step and space step, respectively, is the index of time level, is the index of locations, and is the general function.
Then, the water level and flow discharge at each section of a channel or pipe can be denoted by using water level at the junctions (as shown in Figure 1) if a recursive method is adopted. Flow discharge and the water level at junctions can be solved by adopting an iteration procedure. Finally, the water level and flow discharge of each section are updated by applying the obtained state value of junctions.
However, Preissmann’s scheme is not able to simulate supercritical flows, and to overcome this limitation the supercritical flow is treated as a subcritical flow by reducing the influence of the convective momentum term in the momentum equation. The convective momentum term can be expanded to two parts,where is the speed of a section.
A reduction coefficient based on the Froude number is then introduced to gradually reduce the first part of the convective momentum term.
2.3. 2D Model
The 2D model for urban surface flood inundation simulation is based on two-dimensional shallow-water equations in a conservative form. If we neglect wind effects and the Coriolis term, the governing equations can be written aswhere denotes time; and are the Cartesian coordinates; , , and denote vectors of the conserved flow variables and fluxes in the - and -directions, respectively; and is the source vector containing the bed slope source and friction source . These vectors are expressed aswhere , , , and are the water depth, depth-averaged velocity in the - and -directions, and bottom elevation, respectively; and are bed slopes in the - and -directions, respectively; and and are friction slopes in the - and -directions, respectively.
A Godunov-type finite volume method with an HLLC flux function is adopted and the variables are located at the geometric centres of cells. The 2D SWE is integrated over the control volume and expressed as
Using the second-order backward difference scheme at each cell in time discretization and applying Green’s theorem, (6) can be written aswhere is the area of cell i; is the unit outward vector normal to the boundary; and are the index and length of the faces of cell i, respectively; is the numerical flux through the edge; is the physical time level; is the source term of cell i; and are the solutions from the current and previous time steps; and is the unknown solution of the next time step to be solved.
Equation (7) can be solved by introducing a pseudo-time variable ,
The process of iteratively solving (8) is known as inner iteration. As the inner steady problem converges, the residual term converges to zero and the solution of (7) advances from to . The inner steady problem can be solved using a variety of schemes, and an effective implicit matrix-free LU-SGS algorithm [17] is adopted in this study. To obtain high efficiency and prevent possible endless loops, a convergence tolerance and a maximum number of iterations are used to control the inner iteration scheme. The iterative process is terminated only if the residual term is less than the convergence tolerance or the number of iterations exceeds the maximum number of iterations.
3. Interaction between the Sewer and Overland Flow
In this study, the manhole (junction) function points of the flow exchange between the surface and subsurface systems (as shown in Figure 2). Surcharge occurs when water flows from the sewer drainage system to the ground surface system (DTG) at a time when the water level in the manhole reaches ground elevation. In addition, water flows from the ground’s surface into the sewer drainage systems (GTD) when the water level, or the water head in a manhole, is lower than the ground water level.
There is currently no specific theory or widely accepted methods used to accurately describe the overflow from manhole to surface. Therefore, in this study the discharge of DTG is computed during the iteration process of solving the manhole water levels and is given bywhere is the area of the manhole; and are the water level at the current and previous time steps, respectively; is the inflow or outflow of the branch around the manhole; and is the number of the branches around the manhole.
The discharge of GTD through manholes can be described using either a weir equation or an orifice equation [3, 18]. In this paper, the weir equation, which handles both free and submerged flows, is adopted and is given bywhere is the discharge coefficient; is the perimeter of the node; and are the water depth in the node and ground surface above the top of the node, respectively; and is the exchange discharge.
4. Application and Discussion
In this section, the capability of the proposed model is validated using four test case applications: one extensive physical experiment benchmark case, one vertical linking example and two real drainage cases with actual topography.
4.1. Toce River Experiment with Blocks for 2D Model
In urban areas, densely spaced buildings have a significant effect on flood propagation and inundation. To provide validation data for urban flood models, a series of experiments were conducted for the IMPACT (Investigation of extreMe flood Processes And unCertainTy) project by CESI (Centro Elettrotecnico Sperimentale Italiano) [19]. These experiments were extensively applied to test the applicability of models in high density urban areas [11, 20, 21]. In the experiments, a cluster of blocks was placed in an instrumented physical model of a short reach of the Toce River valley in Italy on a 1 : 100 scale to simulate the blockage effects of an urban district.
One set of experiments was used to verify the 2D model presented in this paper, and the terrain of the domain and the distribution of blocks are illustrated in Figure 3. Although the possibility of a very thin film of water or moisture wetting the concrete bed due to previous runs during the experiment exists, the water film has little effect on the results, so the whole domain was assumed initially dry and in the given inflow boundary (illustrated in Figure 4) the model was run for 60 s. Experimental data from four observation stations were used for validation. Two stations (P3 and P4) located at the upper reaches of the blocks and two stations (P9 and P10) located between the building and the downstream of the blocks, respectively. Figure 5 shows the model’s predictions and experimental measurements at these four observation stations. The numerical results of the water depth at these locations are in good agreement with the experimental data and the predictions of a Godunov-type 2D flood model developed by Kim et al. [21].
(a)
(b)
(c)
(d)
The explicit MUSCL-Hancock scheme (MHS) is more efficient than the commonly used Runge-Kutta scheme and is thus chosen here to be representative of explicit schemes. The essential idea of the MUSCL-Hancock scheme is to use cell averages to predict the values of the conserved quantities at cell edges at and then use these predicted values to update the solution to [22, 23]. A comparison was made between the DTS and the MHS, where the adaptive step method was used to obtain the most efficient performance. Courant number was used to determine the time steps dynamically for both MHS and DTS. Table 1 shows the relative run time of various cases computed by DTS and MHS with different Courant numbers, and a comparison of numerical results is presented with different Courant number in Figure 6. It can be seen from Figure 6 and Table 1 that the run time of DTS are only 37% of that determined with MHS with no other significant differences observed. Courant number must be less than 1.0 for an explicit scheme and thus the time steps of the MHS are largely limited. Because DTS is unconditionally stable, large time steps can be adopted. The results also suggested that the efficiency of DTS improved as the Courant number increased; however, the effect was not obvious over a certain range.
4.2. Real Drainage Case for 1D Model
To verify the ability, reliability, and precision of the 1D model presented in this paper in relation to dealing with an actual urban drainage problem, a real drainage system of a small community in Changsha City, China, was simulated. Results of the present model were then compared with those of the SWMM model [6], which is one of the most complete and widely used models globally for planning, analysis, and design related to drainage systems within urban areas. The community was built in the 1990s and covers a total area of 11.7 hectares, 56% of which are impermeable, and the average gradient is 0.3%. As illustrated in Figure 7, the entire area was divided into 23 subcatchments, and 23 trunk pipes were identified for modelling.
Two designed storms (labelled R1 and R2) were calculated by the model, and Figure 8 illustrated the discharge hydrographs of the outfall. The total rainfall of R1 and R2 was 58 mm and 75 mm, respectively. It can be seen from Figure 8 that the results of present 1D model in this paper are almost identical with the results of SWMM model, which indicated that the proposed model has reliable numerical accuracy. It should be noted that two models have some nuances in the calculation of the flow peak, which can possibly be attributed to the different numerical methods used in the two models.
(a)
(b)
4.3. Vertical Linking Case
A vertical linking case was designed to assess the rationality and feasibility of the coupling strategies of the 1D and 2D models presented in this paper. As shown in Figure 9, the pipe network was formed by six pipes and six nodes, and all nodes except Nodes 1 and 6 could exchange water with the surface. The pipes and nodes were numbered (their attributes are listed in Tables 2 and 3, resp.). A square closed floodplain with a length of 200 m was situated at the top of Nodes 2, 3, 4, and 5; the elevation of the floodplain was 0 m. The Manning coefficient for the floodplain was set to 0.025; the floodplain and pipe network were initially dry and the inlet flow at Node 1 was increased from 0 to 1 m3/s in the first 10 minutes and then maintained at a constant rate, while a free boundary was given at Node 6.
The coupled model presented in this paper was used to simulate this case, while Infoworks ICM, which is widely used commercial drainage software that includes a 1D model and a 2D model, was also used to simulate this case for a comparison.
Figure 10 illustrates the distribution of the flow field and water level throughout the floodplain; Tables 4 and 5 show the discharge in each pipe and the water depth at each node after a steady state is reached, respectively. From these figures and tables it can be seen that the water flows from the pipe to the surface through Node 2, as the pipe diameters of Pipes 2 and 3 are relatively small and the drainage capacity is insufficient, but that water flows from the floodplain to the pipe network through Nodes 3, 4, and 5. The given inlet flow at Node 1 maintained at a constant rate at last, so once a balance of the two kinds of flow is reached, which means the discharge of water flowing in and out is equal, the flow in the floodplain was eventually stabilized. From a comparison of the results in these tables, it is evident that the pipe discharge and node water depth simulated by the two models are basically in agreement, which indicates that the coupled model presented in this paper gives reliable results. The differences in the results simulated by the models may be due to different coupling algorithms and meshing. However, it is considered that overall the model results have certain rationality.
(a)
(b)
With respect to the water volume balance, the accumulated inflow volume should be equal to the sum of the accumulated outflow volume and the water volume flowing onto the floodplain, if the change of water volume inside the channel is ignored. Similarly, the volume in the floodplain should be equal to the net exchange water through four nodes. Figure 11(a) shows the history of the change process in total volume of GTD, total volume of DTG, and the water volume within the floodplain. Figure 11(b) shows the history of the change process in total inflow volume, the sum volume of the outflow, and the water in the floodplain. From Figures 10(a) and 10(b), it can be seen that the sum of the GTD and water volume within the floodplain is equal to the total volume of the DTG and that the inflow volume is equal to the sum of the outflow volume and the water volume within the floodplain; this indicates good water conservation properties.
(a)
(b)
4.4. Real Drainage Case for Coupled 1D and 2D Model
The coupled 1D and 2D model proposed and established in this paper was further tested by application to a real urban rainstorm simulation case in the Xinhepu Community, Guangzhou City. Complete data of the terrain and pipe network are available for this community and it was thus selected as a benchmark case to assist the verification of the proposed model. The total area of the study area is 12.27 hectares and 85% of the area is impervious. Subcatchment discretization was based on the layout of streets, buildings, and terrain features and the entire study area was divided into 94 subcatchments. In all, 381 conduits were identified for modelling through streamlining and consolidation; these are shown in Figure 12(a). To make an accurate calculation of the flood depth caused by a rainstorm, the surface 2D model considered only the street area; the housing area was assumed to be watertight and thus the water flowed only onto the street. The assumption is reasonable because the flood water will not flow into the building when the water depth is not very deep. Figure 12(b) shows the distribution and terrain of the streets in the calculation range and it can be seen that the elevation of the whole area became gradually lower from north to south.
(a)
(b)
Four rainstorms with different return periods were simulated to assess the drainage capacity and waterlogging risks of the community. A return period is a statistical measurement typically based on historic data denoting the average recurrence interval of the rainfall or flood over an extended period of time. Details of pipes reached full drainage capacity (FDC) and surface flooding are presented in Tables 6 and 7, respectively, where the following is evident. () With an increase in the return period of the rainstorm, the rate of the FDC pipes improved and there was an increase in the maximum velocity and flooded area of the streets. However, there was little change in the length of time that pipes remained in FDC state; this could be because that the drainage capacity of partial pipes was poor and that the FDC state occurred during light rainfall, thereby demonstrating that increases in rainfall intensity have little effect on the length of time that pipes remain in a FDC state. () When rains with return periods of 1 and 2 years occurred, a smaller number of pipes reached FDC and the water ran at a low velocity in the streets with a depth below 0.10 m. However, with a return period of 10 years, almost one-third pipes were at a FDC state and there was greater surface flooding in the streets, thereby indicating an increased risk of flooding. However, overall, rainstorm runoff was eliminated quickly and did not cause a large area of flooding, which indicates the community has good drainability.
Figure 13 shows the maximum water depth and FDC time of pipes of a rainstorm with a return period of 10 years. As can be seen from Figure 13, the FDC time of pipes is mainly concentrated within 15 min and the water depth within the streets is within 0.2 m, which is completely consistent with the results of Table 6. The drainage network and facilities of Xinhepu Community have been upgraded in recent years and the overall drainage drainability is good; therefore, the results of the model are in good agreement with the actual situation of the region.
(a)
(b)
5. Discussion and Conclusion
In this paper, a coupled 1D and 2D model for urban flood inundation using an implicit DTS method was presented and analysed. The following conclusions are made.
An implicit dual time-stepping method was applied to the 2D model to improve run-time efficiency. Results show that the run time can be reduced by more than 60% with no significant loss of accuracy, compared to the explicit scheme. The effectiveness of the present model in simulating floods in a dense urban area was verified through application to an experimental case. The model shows good performance, which indicates its capability of accurately and efficiently handling complex flows involving urban buildings.
A vertical linking case was designed to assess the rationality and feasibility of coupling strategies and the water volume balance. The results computed by the present model were compared to those obtained by Infoworks ICM; only subtle differences were noted, which indicates that the coupling strategies of the present model have certain rationality. A water quantity analysis also shows that the model has good water conservation properties.
The effectiveness and robustness of the present model were verified through application to an experimental case of an urban flood, a real drainage case for the 1D model, a vertical linking case, and an actual case with irregular topography and complex pipe network. The model once again showed a good performance, which indicates that it is capable of accurately and efficiently simulating complex water flows in urban areas.
Weir formula and orifice formula are the most common methods to calculate flow interaction including GTD and DTG, but they do not always give accurate results in some specific situations, which is a problem faced by almost all urban flood models. The presented 1D-2D coupled model adopted a weir formula to compute the discharge of GTD and a method to compute the discharge of DTG during the iteration process. This combination of the two methods is slightly more stable, but there is no evidence that this combination is more applicable than weir formula and orifice formula. Therefore, more studies about flow interaction mechanism and calculation method are very necessary and should be considered in future work.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the National Natural Science Foundation of China (Project no. 50979062), the International S&T Cooperation Projects, the Ministry of Science and Technology of China (no. 2012DFG21780), and IWHR Research & Development Support Program (JZ0145B222017).