An urban pluvial flood simulation model based on diffusive wave approximation of shallow water equations

Extreme natural disasters have been more frequent in recent years due to many factors, such as global climate change. Urban floods caused by sudden heavy rainstorms have become serious problems in many cities, where floods have resulted in heavy economic loss and have also led to casualties. Risk analysis is needed to support disaster prevention and mitigation. Numerical simulation is an effective tool for quantitative risk analysis, and for this reason studying urban rainstorm flood simulation will be helpful for disaster response.

Many kinds of numerical simulation methods have been used in flood simulation (Yin et al. 2011; Vacondio et al. 2012; Hillebrand et al. 2016; Rai et al. 2016). Among them, the methods based on shallow water equations (SWEs) are widely used because they have a clear physical meaning and achieve a balance between simulation speed and result accuracy. Some of these methods solve fully dynamic SWEs (Caleffi et al. 2003; Liang et al. 2007). Smith & Liang (2013) utilize GPU-based technology for fast urban flood modeling for the fully dynamic SWEs. From these studies, it can be seen that the methods for solving fully dynamic SWEs are complex, particularly when applied to complex topographies found in the real world. To overcome this problem, some simplified methods have been established. Diffusion wave approximation is one of the most important of the simplified solving methods (Fernández-Pato & García-Navarro 2016). It ignores the acceleration effect and assumes that the gravitational force and the resistance force are always in equilibrium, thus greatly simplifying the solving process. There are many flood simulation models based on the diffusion wave approximation, such as LISFLOOD-FP (Bates & De Roo 2000) and JFLOW (Bradbrook et al. 2004). Also, some recent studies deal with model reduction to achieve faster urban flood simulation (Chen et al. 2012a, 2012b).

However, most of the current flood simulation models are designed for fluvial floods. Many important aspects of urban pluvial floods, such as building considerations and rainfall considerations, are not given enough attention. Improper handling of these details may introduce errors or even result in simulation failure. Therefore, a simulation model designed for urban pluvial floods is needed to make numerical simulation and risk analysis more accurate.

In this paper, an urban pluvial flood simulation model, UPFLOOD, is established based on diffusion wave approximation of SWEs. The model takes into account all the characteristics of urban pluvial floods, and includes many improvements in simulation details. The model is validated through case studies and experimental data.

For open channel flow the Manning Formula is much more widely used than the Darcy-Weisbach Equation. Therefore, in this paper we use Equation (3) rather than Equation (4).

The simulation region is divided into a grid. For simplicity, all the cells are squares with side length of w (m). Figure 1 shows a grid cell (Cell C) and its 8 neighboring cells. Of the 8, only Cells N, E, S, and W are directly connected to Cell C, through Edges CN, CE, CS, and CW respectively. Water can only flow into or out of Cell C via these edges.

Figure 2 shows a flowchart of the simulation. For each time step, the flowrate at each cell edge is calculated using the water depth in the cells. Then, the time step is determined from the water depth and flow velocity. Next the water depth for each cell is updated using the flowrates at the cell edges. If the simulation end condition is satisfied (commonly this is when the simulation time exceeds a given value), the simulation is complete; otherwise, a new time step is entered and the simulation continues. Details of these steps are described below.

The water surface gradient at cell edges is calculated in a local coordinate system. The cell stencil for a cell edge is shown in Figure 3.

When calculating using Equation (13), if certain cells (Cells 1, 2, 5, and 6) are unavailable (out of boundary or occupied by buildings), a problem will occur since the water surface elevation in these cells is unavailable. In this case, the method for calculating needs to be modified. Equation (13) consists of two terms: the first term is the finite difference of water elevation between Cells 1, 3, and 5; the second term is that for Cells 2, 4, and 6. Take the first term as an example. Cell 3 should always be available (otherwise there will be no need to calculate the flowrate at that edge, as it will always be zero). Figure 4 shows all 4 situations for Cells 1, 3, and 5.

For these 4 cases, the first term should be calculated using ⁠, ⁠, and 0 respectively. The second term can be obtained using a similar technique.

The flowrate is positive when water flows from Cell 3 to Cell 4, and negative when water flows from Cell 4 to Cell 3. The water depth d in this equation is calculated using Equation (14). Similarly, flowrates at other cell edges can be obtained. Flowrates are positive if water flows are in +x or +y direction; negative if water flows are in −x or −y direction.

To ensure accuracy, a maximum time step (s) is set. If the time step obtained by Equation (17) is greater than ⁠, the time step will be set to ⁠.

Now we prove that when C is less than 0.25 and the source term s is ignored, the water depth will never be negative after a time step, as long as the water depth at beginning is not negative.

In this paper, we re-implemented two existing diffusion wave approximation models for comparison.

When water flow is not perpendicular to the cell edges, errors will be introduced (it can be proved that the flowrate obtained using Equation (26) will be greater than that obtained using Equation (16)). Therefore, the UPFLOOD model that is established in this paper should be more accurate than the LISFLOOD-FP model.

JFLOW (Bradbrook et al. 2004) is a two-dimensional overland flooding simulation model that is based on a diffusion wave approximation of SWEs. Unlike our UPFLOOD model, flowrates in the JFLOW model are calculated at the center of the cells, rather than at the cell edges. For each cell, the steepest water surface slope vector (equivalent to water surface gradient in this paper) is obtained by comparing the water elevation in this cell with those in the four neighboring cells. Once the flowrate vector has been calculated, the flowrates to the neighboring cells in the x and y directions are set to be the x and y components of the vector, respectively.

Because each cell has only one flowrate vector, water can only flow out of a cell in two orthogonal directions. In some cases, this will introduce errors. When simulating floods using the JFLOW model, if the water surface elevation in a cell is higher than that in both east and west cells (or north and south cells), water will only flow out from the cell to one of the neighboring cells (the one with the lower water surface elevation). However, in the real world, water will flow in two opposite directions in that situation. This means that the JFLOW model will introduce extra simulation errors under certain conditions. A detailed simulation case will be analyzed later. In a plain flood simulation, the topography is simple and such problems seldom occur. However, in urban pluvial flood simulations with complex topography, this problem may be significant.

In urban pluvial flood simulations, many details should be treated properly, otherwise significant errors may be introduced and the accuracy of the simulation results cannot be ensured.

In this paper, buildings are treated as solid blocks. Grid cells that are occupied by buildings are excluded from the simulation area. Since water cannot flow in or out of a building cell, the flowrate between a building cell and a normal cell is always zero.

In most of the cells, the rainfall source term is equal to the rainfall intensity I (m/s, converted from common units such as mm/h). However, in cells that are adjacent to building cells, the rainfall source term is slightly greater than the rainfall intensity. The difference is equal to the total number of cells occupied by the building divided by the total number of cells surrounding the building. In this way, the total rainfall in the simulation area will not decrease due to buildings in this area. This consideration is also consistent with actual situations that most buildings directly discharge rain water from the roof to the adjacent ground.

For simplicity, the extra resistance caused by rainfall is ignored in this paper unless specified otherwise.

The ground infiltration source term can be calculated using empirical formulas such as Horton's equation, and the Green-Ampt method (Mays 2010). Although infiltration parameters for different kinds of soil are available, data for common urban ground surfaces, such as concrete and asphalt, are still lacking. For this reason, in this paper we consider urban ground surfaces to be impervious, and infiltration is ignored.

As negative source terms (ground infiltration source term and drainage source term) exist, the total source term may be negative. In the real world, a negative source term means the cell is drying, it cannot cause negative water depth. In order to make sure that the water depth in cells is always non-negative and thus ensure stability of the simulation, additional treatment is needed. As previously mentioned, when the Courant number is less than 0.25 and the source term is ignored, the water depth will never be negative at the end of a time step, as long as the water depth is not negative at the beginning of the time step. Therefore, if a negative water depth occurs, it must be caused by a negative source term (drainage or infiltration). In this situation, a simple but effective method is adopted: if the water depth in a cell is negative at the end of a time step, adjust drainage and/or infiltration rate so that the water depth of this cell is exactly zero.

In fluvial flood simulations, the process of wetting and drying is important because the flooded region is usually an important result. Many useful methods have been established (Bradbrook et al. 2004; Yu & Lane 2006). However, in an urban pluvial flood simulation, almost no useful data can be obtained through counting wet grid cells, because all the cells are wet during the rain. Also, the model established in this paper can work properly no matter whether the cells have positive water depth or zero water depth. Therefore, no special treatment of drying and wetting is needed. In this way, artificial parameters are not needed and the simulation becomes more ‘natural’. In addition, the simulation procedure can be simplified and thus the simulation speed is slightly increased.

The inertia terms are ignored in the diffusion wave approximation, which leads to the ‘checkerboard oscillation’ phenomenon when cells are small and the time step is large (Hunter et al. 2005). In one time step, the volume of water flow between two adjacent cells is so great that the water surface elevations of these two cells reverses in the next time step and oscillation occurs. This phenomenon is bad for simulation accuracy. In fluvial flood simulations, this phenomenon seldom occurs because the water is coming from only one side (the river side) and therefore the water depth gradient has one dominating direction. However, in urban pluvial floods, still water surface is common, where the water surface gradient is very small and is easy to be disturbed. Also, in urban rainstorm flood simulations, cells are often very small to represent the complex landforms caused by buildings and infrastructure. Therefore, the ‘checkerboard oscillation’ phenomenon is common in urban rainstorm flood simulations. Hunter et al. (2005) pointed out that there are two simple solutions to prevent the ‘checkerboard oscillation’ phenomenon. The first is to add a flow limiter, the second is to use adaptive time steps.

It should be noted that in the original text of Hunter et al. (2005), the limitation is twice as large as that of Equation (35). However, that will be not enough for small ‘checkerboard oscillations’ caused by numerical errors to decay down to a flat free surface. So in this paper, Equation (35) is used to solve this problem.

After time step limitations of all edge cells are taken into consideration, the maximum possible time step value can be obtained. The time step is determined according to the Courant number condition and time step limitations.

In this way, there will be a lower bound for Equation (35), and the problem of extremely small time steps will be solved. However, an artificial parameter and formula will be introduced, which will introduce extra errors and reduces the ability of the simulation to reflect the nature of water flow.

In this paper, a new method, the advanced flow limiter method, is established to prevent the ‘checkerboard oscillation’ phenomenon.

This is approximately equivalent to stopping water flow after time ⁠, when the water surface elevations in these two cells are equal. In this way, water surface elevations in two cells will not reverse and the ‘checkerboard oscillation’ can be avoided.

This method can be seen as a generalized flow limiter. It uses not only water surface elevation but also the change rates of water surface elevation. The artificial influence introduced by this method is much lower than for traditional flow limiters. Taking stable flow on a gentle slope as an example, the change rates in all grid cells will be very small (may be even zero). If a traditional flow limiter method is used, the flowrate will probably be limited because the gradients of the water surface elevations are small, thus additional errors are introduced. However, if the advanced flow limiter method is used, such problems will not exist because the limit condition will never be triggered, since the ‘chasing time’ will be very large.

Three different methods for the prevention of the ‘checkerboard oscillation’ phenomenon (flow limiter method, adaptive time step method, and advanced flow limiter method that is established in this paper) will be tested and compared later in this paper.

Three different flood simulation models: LISFLOOD-FP model, JFLOW model, and the UPFLOOD model established in this paper, are tested with simple cases. The results obtained by the different models are compared to theoretical solutions, and the effectiveness of these models is evaluated.

In all three cases, the size of the grid cells is set to 5 m × 5 m. Ground infiltration and drainage are ignored. The Manning coefficient is set to 0.01. For convenience, variables X and Y are used to denote the grid cell number in x and y directions.

Figure 5 shows the sketch map of Case 1. The simulation area covers 50 m × 50 m, and is divided into a 10 × 10 grid. The ground slope is 0.01 radians (=0.573 degrees). The ground is high in the west and low in the east. Stable slope flow with a water depth of 0.1 m is simulated. The initial condition is: water depth across the whole simulation area is 0.1 m when t = 0. The boundary condition is: water depth is always 0.1 m in cells where X = 1 or X = 10. The simulation output result is the total flowrate across the thick red line in Figure 5.

Theoretically, according to the Manning formula, the slope flow velocity is 2.154 m/s. Considering water depth and total length of the red line, the total flowrate should be 10.77 m³/s. Table 1 lists the simulation results from the three models.

In this simple case, all three models gave the correct result, showing the effectiveness of these models.

Figure 6 shows the sketch map of Case 2. The simulation area covers 50 m × 50 m, and is divided into a 10 × 10 grid. The ground slope is 0.01 radians (=0.573 degrees). The ground is high in the northwest and low in the southeast. Stable slope flow with a water depth of 0.1 m is simulated. The initial condition is: water depth across the whole simulation area is 0.1 m when t = 0. The boundary condition is: water depth is always 0.1 m in cells where X = 1, X = 10, Y = 1 or Y = 10. The simulation output result is the total flowrate across the thick red line in Figure 6.

Theoretically, according to the Manning formula, the slope flow velocity is 2.154 m/s. Considering water depth and total length of the red line, the total flowrate should be 10.66 m³/s. Table 2 gives the simulation results for the three models.

In this case, a significant error is found in the LISFLOOD-FP model results. The fundamental reason for this is that flow velocity components that are parallel to cell edges are ignored.

Figure 7 shows the sketch map of Case 3. The simulation area covers 45 m × 50 m, and is divided into a 9 × 10 grid. The western part of simulation area has a slope of 0.02 radians (=1.146 degrees), low in the west and high in the east. The eastern part of simulation area has a slope of 0.01 radians (=0.573 degrees), low in the east and high in the west. A rainstorm, with duration of 15 min and rainfall intensity of 100 mm/h, is simulated. The initial condition is: water depth across the whole simulation area is 0 when t = 0. The simulation output results are total water volume in west region (cells where X > 5) and east region (cells where X < 5).

As the water depth is shallow and the slope is steep, it can be assumed that the water surface gradient is equal to the ground slope. According to Manning's formula, the flow velocity (proportional to flowrate) is proportional to the square root of the slope. Therefore, the ratio between the water flow from the center cells (X = 5) to the west cells (X < 5) and water flow from the center cells to the east cells (X > 5) is ⁠. The total volume of rain water that falls into the center cells (10 grid cells) is ⁠. Finally, all the water flows to the east cells or to the west cells. The total volume of water that flows to the west cells and east cells is and ⁠, respectively. The rain water that falls into the east cells and west cells stays in the east cells or the west cells, respectively. Finally, the total volume of water in the west region and east region will be and ⁠, respectively. Table 3 lists the simulation results for the three models.

In this case, the results obtained by the JFLOW model are significantly different from the theoretical solutions. The reason for this is discussed above in the ‘JFLOW model’ Section.

The following conclusions can be drawn from these simple case tests: (1) all three models are valid for overland flow simulation; (2) the LISFLOOD-FP model may introduce significant errors if the direction of water flow and cell edges are not orthogonal; (3) the JFLOW model may introduce significant errors at local high points in the terrain; (4) the UPFLOOD model developed in this paper produced good results in all cases, indicating its effectiveness.

Three simulation models (UPFLOOD, LISFLOOD-FP, and JFLOW) are tested with the artificial rainfall experiment carried out by Cea et al. (2010). The case with Geometry S20 and Hyetograph Q25T40 is selected as the test case. Figure 8 compares the simulation results obtained by three models and the experimental data. All the simulation results fit with experimental data. There are observable differences among simulation results obtained by different models, especially at the beginning of the experiment (as shown in Figure 9).

Figure 10 compares the root mean square difference between simulation results and experimental data during the whole experiment (150 seconds). It can be seen that the UPFLOOD model that is established in this paper obtained better results than existing diffusion wave approximation models (LISFLOOD-FP and JFLOW). Figure 11 shows the water depth distribution obtained using the UPFLOOD model. The results have a good fit with results obtained using fully dynamic SWE model (Cea et al. 2010).

A simple case in which the ‘checkerboard oscillation’ phenomenon will easily occur was simulated. Simulations having a large time step (Courant number = 0.1) were carried out using different ‘checkerboard oscillation’ prevention methods: no treatment; flow limiter method; adaptive time step method (h_lin = 1 mm); and the advanced flow limiter method established in this paper. A simulation with a small time step (Courant number = 0.005) was also carried out and the result is used as a reference.

The simulation area covers 40 m × 40 m, and is divided into a 40 × 40 grid. The size of each grid cell is 1 m × 1 m. The Manning coefficient is set to 0.01. There are two small drainage points in the simulation area, with a discharge coefficient C = 0.65, and cross-sectional area A = 0.05. One of the drains lies in the northeast corner, the other in the southeast corner. Initially, there is no water in the simulation area. The rainfall intensity is 200 mm/h and total simulation time is 60 min. Figure 12 shows the simulation results obtained with the different methods.

It can be seen that serious ‘checkerboard oscillation’ will occur in the case where the time step is large and nothing else is done. The flow limiter method can eliminate ‘checkerboard oscillation’, but by limiting the flow, resulted in the simulated water depth being much higher than it should be. The adaptive time step method got very good results, these being nearly the same as the results obtained for a small time step. The performance of the advanced flow limiter method is also good. Figure 13 shows the average water depth errors (compared to simulation results for a small time step) of the different treatment methods.

The adaptive time step method performs the best, followed by the advanced flow limiter method. The flow limiter method and the no treatment method do not perform well, which intuition would suggest. For numerical simulations, both result accuracy and simulation speed should be taken into consideration. Otherwise, shortening the time step would solve all problems and no extra treatment method would be needed. Figure 14 shows the time that simulations using different treatment methods take. All simulations were carried out on the same personal computer.

Simulations using a large time step are much faster than simulations using a small time step. Of the large time step simulation methods, the adaptive time step method, which has the lowest simulation error, is the slowest. The advanced flow limiter method has relatively low result error (about 20% of the flow limiter method and about twice the adaptive time step method) and a short simulation time (about 20% of the adaptive time step method, and similar to the flow limiter method). It should therefore be a good method for preventing the ‘checkerboard oscillation’ phenomenon. Later in this paper, all the simulations will use the advanced flow limiter method.

This section is a proof of the concept that the UPFLOOD model can deal with complex situations and is helpful for urban drainage system planning.

A region in Beijing is selected as the study area. It stretches 1.1 km from north to south and 1.5 km from east to west covering an area of 1.65 square kilometers. To avoid any influence caused by boundaries, the simulation area was extended by 50 meters in all directions. Figure 15 shows the region of study and the extended simulation area.

The topography is a very important factor in floods. In this paper, Digital Elevation Model (DEM) data are provided by the Beijing Institute of Surveying and Mapping. Figure 16 shows the elevation map of the simulation area. The position of every drain was carefully mapped on site. The investigation found that there are about 620 drains in this region and most of them have a cross sectional area of about 0.17 m².

In practice, rainstorm scenarios should be selected according to meteorological data. In this paper, for simplicity, three simple rainstorm scenarios are selected: 200 mm/h rainfall of 15 min duration, 100 mm/h rainfall of 30 min duration, and 50 mm/h rainfall of 60 min duration. In all these simulations, the most serious inundation occurs just at the end of the rainstorm. Figure 17 shows the distribution of water depth at the end of the storm.

These rainstorms have the same total amount of rain, but the shorter, heavier storm caused more serious inundation. That is because the amount of water draining out via the drainage system is limited in short, strong storms.

It is obvious that inundation points for different rainstorm scenarios are similar. Basically, they are all local low-lying areas. Setting up additional drains in these areas may be an ideal way to mitigate rainstorm floods. Figure 18 shows a proposed plan for a drainage system upgrade. The number of drains in the simulation area will increase from 620 to 648.

After the proposed upgrade, the total drainage capacity will increase by 4.5% (=(648–620)/620 × 100%). As a comparison, another proposal is to improve each existing drain to increase the total drainage capacity by 10%. Simulations are carried out for these different plans. Figure 19 shows the comparison of inundation area, which in this paper is defined as the total area where water depth is greater than 0.15 m.

It can be seen that adding more drains in the seriously inundated areas is much more effective than evenly improving each existing drain. This observation underlines the fact that plans made with the help of numerical simulations are much better than those produced without this work. This result also shows that the UPFLOOD model is useful for urban pluvial flood prevention and mitigation.

In this paper, an urban pluvial flood simulation model, UPFLOOD, was established based on diffusion wave approximation of SWEs. The model takes full account of characteristics in urban pluvial floods, and includes many improvements in simulation details. These details include a building consideration method, a rainfall consideration method, and so on. A new calculation method of water surface gradient is established. The tangential gradient at cell edge, which is often ignored in existing diffusion wave approximation models (e.g. LISFLOOD-FP), is calculated using a stencil. This method is suitable for complex topology in urban pluvial flood simulation and can reduce unnecessary simulation error introduced by calculation methods. The accuracy and stability of the model were verified through simple cases with analytical solutions and experiments with measured data. The results show that the UPFLOOD model is more accurate than common diffusion wave approximation models (LISFLOOD-FP and JFLOW). A new treatment to avoid ‘checkerboard oscillation’, the advanced flow limiter method, is established. This method can be seen as a generalized flow limiter method. It uses not only water surface elevation but also the change rates of water surface elevation. The artificial influence introduced by this method is much lower than for traditional flow limiters. In comparison with existing treatment methods (flow limiter method and adaptive time step method), the new method proved to be the best. A proof of concept shows that the UPFLOOD model can deal with complex situations and is helpful for urban drainage system planning.

This work was supported by National Key R&D Program of China (Grant No. 2018YFC0809900), the National Natural Science Foundation of China (Grant No. 71774093) and China Clean Development Mechanism Foundation (Grant No. 2013049).