Abstract

Urban floods caused by sudden heavy rainstorms are becoming more frequent and causing serious problems in many cities. Developing methods to simulate urban rainstorm floods is helpful in disaster prevention and mitigation. In this paper, we establish an urban pluvial flood simulation model based on diffusion wave approximation of shallow-water equations. The model takes full account of the characteristics in urban pluvial floods, and includes many improvements in simulation details. These details include building a consideration method, a rainfall consideration method, and so on. A new calculation method of water surface gradient is established, which 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 are verified through simple cases with analytical solution and experiments with measured data. The results show that the new model is more accurate than common diffusion wave approximation models. A new treatment to avoid ‘checkerboard oscillation’ is established. In comparison with existing methods, the new method proved to be the best. A proof of concept shows that the new model can deal with complex situations and is helpful for urban drainage system planning.

INTRODUCTION

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.

METHODS

Basic principle

The one-dimensional fully dynamic SWEs are: 
formula
(1)
 
formula
(2)
where t is time (s), x is a spatial coordinate (m), d is water depth (m), q is specific discharge (m2/s), s is the source term (m/s), including rainfall, infiltration and drainage, g is gravitational acceleration (=9.8 m/s2), is the slope of the ground (dimensionless), is the friction slope (dimensionless). Two common formulas for calculating are based on the Manning Formula (in Equation (3)) and the Darcy-Weisbach Equation (in Equation (4)): 
formula
(3)
 
formula
(4)
where k is a unit conversion factor (=1 m1/3 s−1), n is the Manning coefficient (dimensionless), V is depth averaged flow velocity (m/s), f is the Darcy friction factor (dimensionless).

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).

In diffusion wave approximation of SWEs, it is assumed that the gravitational force and frictional resistance force are always in equilibrium. Thus the momentum equation (Equation (2)) becomes: 
formula
(5)
Extending Equations (1) and (5) to the two-dimensional form: 
formula
(6)
 
formula
(7)
where is the vector of specific discharge (m2/s), is the gradient of the water surface (dimensionless), and is a vector of the friction slope (dimensionless), the friction slope vector and the velocity vector are in the same direction, because the resistance force and water flow velocity are in opposite directions. Therefore, , where is the vector length of , is the vector length of . can be calculated using the one-dimensional formula (Equation (3) or Equation (4)), as long as flow velocity V is replaced by . According to Equation (7), 
formula
(8)
the left and the right side should be in the same direction, so and are in opposite directions: 
formula
(9)
The vector lengths of the left and the right side in Equation (8) should be equal: 
formula
(10)
is a function of , while is not directly correlated with . Therefore, can be solved using Equation (10). When using the Manning Formula, the analytic solution is: 
formula
(11)
can be obtained by substituting into Equation (9), and can be calculated. Then, can be calculated using Equation (6). Therefore, the water depth can be updated for the next time step, and the new flow velocity can be obtained. The simulation is done by repeating this process.

Algorithm overview

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 1

A grid cell and its neighboring cells.

Figure 1

A grid cell and its neighboring cells.

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.

Figure 2

Flowchart of simulation.

Figure 2

Flowchart of simulation.

Calculate flowrates at cell edges

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.

Figure 3

Cell stencil for a cell edge.

Figure 3

Cell stencil for a cell edge.

The perpendicular and tangential components of the water surface gradient are obtained using the following approximations: 
formula
(12)
 
formula
(13)
where , , , , , and are water surface elevations in Cells 1, 2, 3, 4, 5, and 6, respectively.

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.

Figure 4

All four situations of three cells.

Figure 4

All four situations of three cells.

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

Water flow velocity can be calculated using the methods that are introduced in the ‘Basic principle’ Section (e.g. using Equation (11)). During the calculation, the water depth d is: 
formula
(14)
where h is the water surface elevation (m), and z is the ground elevation (m). In any cell, z + d = h.
After is calculated, the x component of the flow velocity can be obtained using: 
formula
(15)
The flowrate at that edge (m3/s) is: 
formula
(16)

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.

Determine the time step

Time step, , is determined using: 
formula
(17)
where C is the Courant number (dimensionless), and and are maximum flow velocity and maximum water depth across the whole simulation area. In order to ensure stability, C should be less than 1. When C is no more 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 the beginning is not negative. This is important for the simulation, as negative water depth will cause wrong simulation results (e.g., the flow velocity calculated by Equation (11) will no longer be a real number due to negative water depth). We show the proof of this later. In practical applications, it is suggested that C should be no more than 0.25, so that non-negative water depth can always be guaranteed and the simulation will never crash.

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 .

Update water depth at cells

Using Cell C in Figure 1 as an example, the change rate in water depth is: 
formula
(18)
The water depth at the end of a time step will be: 
formula
(19)
This completes our introduction to the simulation implementation method.

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.

Take Cell C as an example: only 4 cell edges are directly connected to the cell. Since the flow velocity is never greater than , the maximum possible value of flowrates at these cell edges is: 
formula
(20)
where d is calculated using Equation (14). When water flows out from Cell C, the water surface elevation in Cell C must be higher than that in the neighboring cell, thus 
formula
(21)
Ignoring the source term, the minimum rate of change in water depth can be achieved when water flows out of the cell in four directions at maximum flowrate, 
formula
(22)
Time step 
formula
(23)
Considering that , , therefore: 
formula
(24)
so the water depth at the end of the time step is 
formula
(25)
This means that the water depth will never be negative after a time step. This completes the proof.

Comparison with similar simulation models

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

LISFLOOD-FP model

LISFLOOD-FP (Bates & De Roo 2000) is a floodplain simulation model that is based on a diffusion wave approximation of SWEs. The structures of the LISFLOOD-FP model and the UPFLOOD model are similar. Both models calculate flowrates at cell edges and then update the water depth in the cells for each time step. In the LISFLOOD-FP model the water flow is always assumed to be perpendicular to the cell edges. Therefore, the flowrate can be calculated by a simpler equation: 
formula
(26)
Since the original literature of the LISFLOOD-FP model (Bates & De Roo 2000) did not give the expression for water depth d at cell edges explicitly, Equation (14) is used in this paper.

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 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.

Yu & Lane (2006) provided an equation for calculating equivalent water depth: 
formula
(27)
where and are water depth calculated using Equation (14) in the x and y directions (m), and and are x and y components of the steepest water surface slope vector (dimensionless).

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.

Algorithm details

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.

Building consideration method

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.

Source term consideration method

The source term mainly includes terms for rainfall, ground infiltration, and drainage. The total source term can be expressed as: 
formula
(28)

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.

In some studies it has been observed that rainfall may introduce extra resistance to overland flow (Emmett 1970; Katz et al. 1995). Considering this phenomenon, Equation (7) should be modified to: 
formula
(29)
where is the vector of friction slope caused by rainfall (dimensionless). As with , is also in the opposite direction to the flow velocity. Therefore, Equation (10) will become: 
formula
(30)
If rainfall is seen as a branch inflow with a horizontal initial velocity of zero, an expression of can be obtained through the momentum theorem: 
formula
(31)
Hsu et al. (2000) have used this expression. Using this expression and the Manning formula, Equation (30) can be converted to: 
formula
(32)
This is a quadratic equation of . Since the quadratic coefficient is positive and the constant term is negative, there should be a positive root and a negative root. As , the positive root is the correct solution. That is to say, the new model established in this paper is valid even when the extra resistance caused by rainfall is considered. This demonstrates the extensibility of this new model.

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.

The drainage source term can be calculated using an orifice discharge formula, as long as the drainage pipe system works properly. 
formula
(33)
where C is the discharge coefficient (dimensionless), and A is the total effective cross-sectional area of drainage openings in a cell (m2).

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.

Wetting and drying process method

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.

Prevention of the ‘checkerboard oscillation’ phenomenon

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.

If the flow limiter method is used, the flowrate at cell edges will be limited. Taking Edge CE as an example, when a flow limiter is used, the absolute value of the flowrate will become: 
formula
(34)
If the adaptive time step method is used, the time step is limited according to the water surface gradient. Taking Edge CE as an example using the adaptive time step method, a limit condition for the time step will be: 
formula
(35)

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.

Both of these methods can eliminate ‘checkerboard oscillation’, but they have drawbacks. Using a flow limiter will probably underestimate flowrates and introduce errors. Using the adaptive time step method will result in a very small time step when the water surface is flat (as can be seen from Equation (35), the time step limitation will tend to zero if surface gradient tends to zero), resulting in a sharp drop in simulation speed. A workaround is to use a modified formula for flow velocity when the difference in water surface elevation between two cells is less than a given threshold : 
formula
(36)

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.

Take Cell C as an example: the change rate of the water depth can be obtained using Equation (18). Similarly, the change rates of the water depths in other cells can also be obtained, for example in Cell E. These change rates are also change rates of water surface elevation, as the ground elevations are constants. For Cells C and E, we introduce ‘chasing time’: 
formula
(37)
If , the equation will represent the time needed for the water depth in these two cells to become equal. If , the water depth in two cells will inverse in one step, and the ‘checkerboard oscillation’ phenomenon will occur. At that time, the flowrate between the two cells should be adjusted to avoid this oscillation. In this paper, the adjusted flowrate is: 
formula
(38)

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.

RESULTS AND DISCUSSION

Verifying the different simulation models

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.

Case 1: Orthogonal slope

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.

Figure 5

Sketch map of Case 1.

Figure 5

Sketch map of Case 1.

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 m3/s. Table 1 lists the simulation results from the three models.

Table 1

Simulation results of Case 1

Model Flowrate (m3/s) Result 
LISFLOOD-FP 10.77 Pass 
JFLOW 10.77 Pass 
UPFLOOD 10.77 Pass 
Model Flowrate (m3/s) Result 
LISFLOOD-FP 10.77 Pass 
JFLOW 10.77 Pass 
UPFLOOD 10.77 Pass 

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

Case 2: Oblique slope

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.

Figure 6

Sketch map of Case 2.

Figure 6

Sketch map of Case 2.

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 m3/s. Table 2 gives the simulation results for the three models.

Table 2

Simulation results of Case 2

Model Flowrate (m3/s) Result 
LISFLOOD-FP 12.68 Failed 
JFLOW 10.66 Pass 
UPFLOOD 10.66 Pass 
Model Flowrate (m3/s) Result 
LISFLOOD-FP 12.68 Failed 
JFLOW 10.66 Pass 
UPFLOOD 10.66 Pass 

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.

Case 3: Watershed under rainfall

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).

Figure 7

Sketch map of Case 3.

Figure 7

Sketch map of Case 3.

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.

Table 3

Simulation results of Case 3

Model Volume in west region (m3Volume in east region (m3Result 
LISFLOOD-FP 28.67 27.58 Pass 
JFLOW 31.25 25.00 Failed 
UPFLOOD 28.67 27.58 Pass 
Model Volume in west region (m3Volume in east region (m3Result 
LISFLOOD-FP 28.67 27.58 Pass 
JFLOW 31.25 25.00 Failed 
UPFLOOD 28.67 27.58 Pass 

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.

Brief summary

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.

Test of different models with experimental data

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 8

Comparison of simulation results and experimental data.

Figure 8

Comparison of simulation results and experimental data.

Figure 9

Comparison of simulation results and experimental data (first 30 seconds).

Figure 9

Comparison of simulation results and experimental data (first 30 seconds).

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).

Figure 10

Comparison of error among three simulation models.

Figure 10

Comparison of error among three simulation models.

Figure 11

Water depth distribution obtained by UPFLOOD model (t = 40 s).

Figure 11

Water depth distribution obtained by UPFLOOD model (t = 40 s).

Testing ‘checkerboard oscillation’ prevention methods

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 (hlin = 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.

Figure 12

Simulation results obtained with different methods.

Figure 12

Simulation results obtained with 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.

Figure 13

Average water depth errors for the different treatment methods.

Figure 13

Average water depth errors for 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.

Figure 14

Time cost of different treatment methods.

Figure 14

Time cost of different treatment methods.

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.

Proof of concept

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.

Figure 15

Simulation area.

Figure 15

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 m2.

Figure 16

Elevation map of simulation area.

Figure 16

Elevation map of simulation area.

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.

Figure 17

Distribution of water depth at the end of rainstorms.

Figure 17

Distribution of water depth at the end of rainstorms.

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.

Figure 18

Existing drains and additional drains in proposed plan.

Figure 18

Existing drains and additional drains in proposed plan.

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.

Figure 19

Comparison of inundation area for the different plans.

Figure 19

Comparison of inundation area for the different plans.

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.

CONCLUSIONS

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.

ACKNOWLEDGEMENTS

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).

REFERENCES

REFERENCES
Bates
P. D.
&
De Roo
A. P. J.
2000
A simple raster-based model for flood inundation simulation
.
Journal of Hydrology
236
,
54
77
.
doi: 10.1016/S0022-1694(00)00278-X
.
Bradbrook
K. F.
,
Lane
S. N.
,
Waller
S. G.
&
Bates
P. D.
2004
Two dimensional diffusion wave modelling of flood inundation using a simplified channel representation
.
International Journal of River Basin Management
2
,
211
223
.
doi: 10.1080/15715124.2004.9635233
.
Caleffi
V.
,
Valiani
A.
&
Zanni
A.
2003
Finite volume method for simulating extreme flood events in natural channels
.
Journal of Hydraulic Research
41
,
167
177
.
doi: 10.1080/00221680309499959
.
Cea
L.
,
Garrido
M.
&
Puertas
J.
2010
Experimental validation of two-dimensional depth-averaged models for forecasting rainfall-runoff from precipitation data in urban areas
.
Journal of Hydrology
382
,
88
102
.
doi: 10.1016/j.jhydrol.2009.12.020
.
Chen
A. S.
,
Evans
B.
,
Djordjević
S.
&
Savić
D. A.
2012a
A coarse-grid approach to representing building blockage effects in 2D urban flood modelling
.
Journal of Hydrology
426–427
,
1
16
.
doi: 10.1016/j.jhydrol.2012.01.007
.
Chen
A. S.
,
Evans
B.
,
Djordjević
S.
&
Savić
D. A.
2012b
Multi-layered coarse grid modelling in 2d urban flood simulations
.
Journal of Hydrology
470–471
,
1
11
.
doi: 10.1016/j.jhydrol.2012.06.022
.
Emmett
W. W.
1970
The Hydraulics of Overland Flow on Hillslopes (Report No. 662-A)
. United States Government Printing Office
,
Washington DC
.
Fernández-Pato
J.
&
Garcı́a-Navarro
P.
2016
2D Zero-Inertia model for solution of overland flow problems in flexible meshes
.
Journal of Hydrologic Engineering
21
(
11
).
doi: 10.1061/(ASCE)HE.1943-5584.0001428
.
Hillebrand
G.
,
Klassen
I.
&
Olsen
N. R. B.
2016
3D CFD modelling of velocities and sediment transport in the Iffezheim hydropower reservoir
.
Hydrology Research
48
.
doi: 10.2166/nh.2016.197
.
Hsu
M. H.
,
Chen
S. H.
&
Chang
T. J.
2000
Inundation simulation for urban drainage basin with storm sewer system
.
Journal of Hydrology
234
,
21
37
.
doi: 10.1016/S0022-1694(00)00237-7
.
Hunter
N. M.
,
Horritt
M. S.
,
Bates
P. D.
,
Wilson
M. D.
&
Werner
M. G.
2005
An adaptive time step solution for raster-based storage cell modelling of floodplain inundation
.
Advances in Water Resources
28
,
975
991
.
doi: 10.1016/ j.advwatres.2005.03.007
.
Katz
D. M.
,
Watts
F. J.
&
Burroughs
E. R.
1995
Effects of surface roughness and rainfall impact on overland flow
.
Journal of Hydraulic Engineering
121
,
546
553
.
doi: 10.1061/(ASCE)0733-9429(1995)121:7(546)
.
Liang
D.
,
Lin
B.
&
Falconer
R. A.
2007
Simulation of rapidly varying flow using an efficient TVD-MacCormack scheme
.
International Journal for Numerical Methods in Fluids
53
,
811
826
.
doi: 10.1002/fld.1305
.
Mays
L. W.
2010
Water Resources Engineering
.
John Wiley & Sons
,
London
,
UK
.
Rai
P. K.
,
Chahar
B. R.
&
Dhanya
C. T.
2016
GIS-based SWMM model for simulating the catchment response to flood events
.
Hydrology Research
48
,
384
394
.
doi: 10.2166/nh.2016.260
.
Smith
L. S.
&
Liang
Q.
2013
Towards a generalised GPU/CPU shallow-flow modelling tool
.
Computers & Fluids
88
,
334
343
.
doi: 10.1016/j.compfluid.2013.09.018
.
Vacondio
R.
,
Rogers
B.
,
Stansby
P.
&
Mignosa
P.
2012
SPH modeling of shallow flow with open boundaries for practical flood simulation
.
Journal of Hydraulic Engineering
138
,
530
541
.
doi: 10.1061/(ASCE)HY.1943-7900.0000543
.
Yin
Z.
,
Jie
Y.
,
Xu
S.
&
Wen
J.
2011
Community-based scenario modelling and disaster risk assessment of urban rainstorm waterlogging
.
Journal of Geographical Sciences
21
,
274
284
.
doi: 10.1007/s11442-011-0844-7
.