Introduction

With rapid global climate change and the development of urbanization, urban flood and water logging problems are becoming increasingly prominent worldwide. To reduce the impact of flood disasters on urban areas, it is urgent to propose a set of efficient and stable urban flood disaster simulation methods to simulate rainstorms and flood processes in urban areas and provide a decision-making basis for urban flood control and drainage, rescue and relief (Wang et al. 2017; Xu et al. 2020).

The urban rainstorm and flood simulation methods can be classified as hydrology, hydrodynamics and hydrodynamic coupling methods (Dams et al. 2015; Wijeratne et al. 2023). The hydrology method was the earliest method to be proposed (Bell et al. 2016); this method involves rainfall-runoff generation and flow convergence calculations based on the hydrology principle in a basin with simple structure and high efficiency. Unfortunately, only the flow process at the outlet of the basin can be determined, and the hydraulic characteristics cannot be determined at a specific location. The hydrodynamic method is based on basin grid cells to carry out surface runoff generation calculations (Tansar et al. 2022). This kind of method has high computational accuracy but low computational efficiency. The hydrology and hydrodynamic approach combines the advantages of the above methods and takes the sub-basin as the hydrological response unit to calculate surface runoff generation and flow convergence. The hydrodynamic method is used to calculate the flow after water enters the underground pipe network. The overflow from underground to the surface is calculated by using a two-dimensional shallow water model. This method has both a good hydrology base and high computational efficiency (Zhao et al. 2018; Hou et al. 2017).

In view of the coupling methods for one- or two-dimensional models, several scholars have studied and constructed coupling models based on coupling theory (Yu et al. 2014; Hsu et al. 2000). Amirhossein Nazari1 et al. (2023) proposed the integrated SUSTAIN-SWMM-MCDM approach for the optimal selection of LID practices in urban stormwater systems. The stormwater management model (SWMM) was used for rainfall-runoff and hydraulic modeling. Six scenarios of different combinations of LID practices were developed. In 2018, Ali et al. (2023) studied parallelization of AMALGAM algorithm for a multiobjective optimization of a hydrological model. In 2021, Ghodsi et al. (2021) modeled the effectiveness of rain barrels, cisterns, and downspout disconnections for reducing combined sewer overflows in a city-scale watershed. In 2022, Hashemi and Mahjouri (2022) carried out a global sensitivity analysis-based design of low-impact development practices for urban runoff management under uncertainty. In 2019, Macro et al. (2019) proposed a new multiobjective optimization tool for green infrastructure planning with SWMM.

Many urban rainfall flood simulation methods based on traditional hydrology, which involve fewer model parameters and a fast calculation speed, are applied during the urban planning and design stage. Unfortunately, characterizing the process and mechanism of surface water logging is difficult with such models, and the simulation accuracy is limited (Jenkinson 2010). Hydrodynamic methods directly characterize the physical process of surface flow generation and confluence through solving two-dimensional shallow water equations; these methods can provide more intuitive and powerful support for research on the generation process and mechanism of urban floodwater logging and have greater scientific value and practical value (Hefzul Bari et al. 2016). However, due to urban complexity, it is easy for high-precision urban rainfall flood models to cause numerical algorithm instability, and the calculation time is too long, which reduces their practical use value. Usually, simplified models are needed for this kind of model, but simplified models cannot perfectly characterize small-scale topographic features such as buildings, roads and other structures (Pugliese et al. 2022). In addition, in most studies that use two-dimensional surface hydrodynamic models to simulate urban waterlogging, the influence of underground pipe networks on waterlogging may often be neglected, and the water deduction method may be simply used to depict the confluence and drainage functions of surface floods. In recent years, one- and two-dimensional coupled computing models, such as MIKE-SWMM, have the ability to calculate and couple two-dimensional slope flow generation and non-constant flow in pipe networks (Schubert et al. 2017; Shariat et al. 2019; Sheikh et al. 2021). However, because of its accuracy and stability, complex urban surface conditions are difficult to meet, and the long calculation time makes it difficult to meet the needs of complex applications, such as urban water logging early warning systems and real-time displays (Dumitrache, et al. 2023).

Therefore, considering the problems that exist in previous studies, the main objectives in this paper are to (1) improve SWMM as a foundation for realizing urban flood hydrology and hydrodynamic model coupling; (2) propose coupling methods for one- and two-dimensional models and provide a robust coupling strategy for urban rainstorm and flood multidimensional and multiprocess simulation; (3) construct an urban flood hydrology and hydrodynamic coupling model (HHDCM); (4) construct SWMM simulation results as the input source of the underground pipe network to form a continuous urban flood calculation system; (5) select July.20 extreme rainstorm in Zhengzhou, Henan Province, China, as the sample to verify and calibrate the constructed model; and (6) simulate the urban flood inundation situation caused by July.20 extreme rainstorm; and finally, draw the urban inundation distribution map.

Methodology

Improvement of the SWMM

Model structure and parameters

Figure 1 shows the improved basic structure of the SWMM. The first step for SWMM is to divide the study area into several sub-catchment areas. In each sub-catchment, the surface runoff generation and confluence are calculated, and each sub-catchment area can be divided into permeable and impervious surfaces. Surface runoff enters the pipeline network or river channel through slope convergence and then flows to the outlet of the drainage area after confluence (Baek et al. 2020; Balali et al. 2014; Barbosa et al.2012).

Fig. 1
figure 1

Improved basic structure of the SWMM

Full size image

Table 1 shows the 12 basic parameters of the SWMM. These parameters can be classified as geometric parameters or rating parameters according to different determination methods, in which the geometric parameters are directly obtained through measurement, and the rating parameters need to be determined and verified by the measured urban hydrology process.

Table 1 Basic parameters of the SWMM
Full size table

Among the geometric parameters, the characteristic width W is the characteristic width of the sub-catchment area and is an important and sensitive parameter. According to the geometric relationship, W includes three parts, i.e., W1, W2 and W3, which can be determined according to the following formula: 1, 2 and 3

$${W}_{1}=W$$
(1)
$${W}_{2}=\frac{{A}_{2}}{{A}_{2}+{A}_{3}}W$$
(2)
$${W}_{3}=\frac{{A}_{3}}{{A}_{2}+{A}_{3}}W$$
(3)

where \({W}_{1}\) is the width of impervious area A1 (m); \({W}_{2}\) is the width of permeable area A2 (m), where there is no water storage in the low surface area (m); and \({W}_{3}\) is the width of permeable area A3(m), where there is a certain amount of water storage in the low surface area (m). A1, A2 and A3 can be determined by the area and percentage of the sub-catchment area. The key problem is to determine the total characteristic width W in a sub-catchment area. Generally, W can be determined by formula (4).

$$W=C\sqrt{\left({A}_{1}+{A}_{2}+{A}_{3}\right)}$$
(4)

where C is the shape coefficient, which is generally 0. 2–0. 5.

Surface flow generation process calculation (Platz et al. 2020; Rezazadeh et al. 2019)

Based on the concept of infiltration capability, Horton proposed an infiltration pattern. When the rainfall intensity is less than or equal to the ground infiltration capacity, all the rainfall infiltrates underground, and no surface runoff is generated. When the rainfall intensity is greater than the ground infiltration capacity, the infiltration rate is equal to the infiltration capacity. The rainfall above the infiltration capacity will form surface runoff. Therefore, the runoff generation pattern can be calculated by formula (5):

$${R}_{s}=\underset{i>{f}_{p}}{\overset{T}{\int }}(i-{f}_{p}){\text{d}}t$$
(5)

where Rs is the surface runoff (mm); i is the rainfall intensity (mm/h); fp is the infiltration capacity (mm/h); t is the time (h); and T is the total rainstorm duration(min).

Surface runoff confluence calculation (Saeid Eslamian et al. 2023; Du et al. 2012)

In the SW, each sub-catchment area is generalized as a nonlinear reservoir for confluence calculations. In confluence calculations, it is assumed that each sub-catchment area is a rectangle with a width W (called the characteristic width). The nonlinear reservoir water balance equation can be expressed by the partial differential Eq. (6):

$$\frac{\partial Q}{\partial t}=i-e-f-q$$
(6)

where \(\frac{\partial Q}{\partial t}\) is the water storage variation at time t (m3); i is the rainfall intensity (mm/s); e is the surface evaporation rate (mm/s); f is the infiltration rate (mm/s); and q is the surface runoff rate (mm/s).

Assuming that the surface runoff in each sub-catchment area is regarded as uniform flow passing through a rectangular channel with width W, depth d − ds, and slope S, the surface runoff Q can be calculated by using Manning’s formula (7):

$$ Q = \frac{1.49}{n}S^{{{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}}} R_{x}^{{{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}}} A_{x} $$
(7)

where n is the surface roughness coefficient, S is the average slope of the sub-catchment area (%), \({A}_{x}\) is the crossing section area of the nonlinear reservoir (m2), and \({R}_{x}\) is the hydraulic radius (m). The water section area \({A}_{x}\) and hydraulic radius \({R}_{x}\) are calculated separately by the following formulas (see formula (8) and formula (9), respectively):

$${A}_{x}=W(d-{d}_{x})$$
(8)
$${R}_{x}=d-{d}_{x}$$
(9)

By substituting formulas (8) and (9) into formula (7), the surface runoff equation can be obtained as formula (10).

$$ Q = \frac{1.49}{n}WS^{{{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}}} \left( {d - d_{x} } \right)^{5/3} $$
(10)

By dividing Q by the area of the sub-catchment area, the surface runoff flow q can be calculated via formula (11) and formula (12).

$$ q = \frac{{1.49WS^{1/2} }}{{A_{n} }}\left( {d - d_{x} } \right)^{5/3} $$
(11)
$$ \frac{\partial Q}{{\partial t}} = i - e - f - \alpha \left( {d - d_{x} } \right)^{5/3} $$
(12)

in which, \(\alpha = \frac{{1.49WS^{1/2} }}{{A_{n} }}\)

where α is the confluence coefficient on the slope surface.

In the above formula, let i * = i − e − f; then, i * is the net rain, and the above formula can be rewritten as follows:

$$ \frac{\partial Q}{{\partial t}} = i \times \left( { - \alpha \left( {d - d_{x} } \right)^{5/3} } \right) $$
(13)

This is a nonlinear differential equation. Usually, an approximate solution can be obtained by using the Newton–Raphson iteration method (Nazeri Tahroudi et al. 2023).

Pipeline network flow confluence model (Lu et al. 2019)

The pipeline flow confluence process can be approximately considered flood wave movement without lateral flow input. The SWMM can provide three methods (constant flow method, flood wave movement and flow dynamic wave method) for pipeline network flow simulation. The constant flow and flood wave movement methods are simplified methods that have high computational efficiency but low application potential. The flow dynamic wave method is based on the one-dimensional (St.Vennat) flow equation, which has the best application effect. Therefore, the principle and calculation method of the flow dynamic wave method are emphasized here. The governing equation for the flow dynamic wave method is the Saint-nan equation, as shown in the following equations:

$$\frac{\partial Q}{\partial t}+\frac{\partial A}{\partial t}=0$$
(14)
$$ gA\frac{\partial H}{{\partial x}} + \frac{{\partial \left( {Q^{2} /A} \right)}}{\partial x} + \frac{\partial Q}{{\partial t}}gAS_{f} = 0 $$
(15)

where Q is the flow rate (m3/s), A is the crossing section area (m2), H is the water depth (m), g is the gravitational acceleration (m/s2), and Sf is the friction resistance ratio, which can be calculated by using the Manning formula(16).

$$ S_{f} = \frac{J}{{gAR^{4/3} }}Q\left| V \right| $$
(16)

where J = gn2, n is the integral roughness of the pipeline network, R is the hydraulic radius (m), and v is the section flow speed (m/s). The absolute value means that the friction resistance force direction is opposite to the flow direction.

From \(\frac{{Q}^{2}}{A}={v}^{2}A\), Eq. (17) can be obtained:

$${\text{gA}}\frac{\partial H}{\partial x}+2Av\frac{\partial v}{\partial x}+{v}^{2}\frac{\partial A}{\partial x}+\frac{\partial Q}{\partial t}+gA{S}_{f}=0$$
(17)

By substituting Q = Av into formula (17) and multiplying by Q = Av, we can obtain the following equation:

$${\text{Av}}\frac{\partial v}{\partial x}=-v\frac{\partial A}{\partial t}-{v}^{2}\frac{\partial A}{\partial x}$$
(18)

By combining formula (17) and formula (18), we can obtain formula (19).

$${\text{gA}}\frac{\partial H}{\partial x}-2Av\frac{\partial A}{\partial t}-{v}^{2}\frac{\partial A}{\partial x}+\frac{\partial Q}{\partial t}+gA{S}_{f}=0$$
(19)

The solution of Eq. (19) needs to be combined with the following continuity Eq. (20):

$$\frac{\partial H}{\partial t}=\frac{\sum {Q}_{i}}{\omega }$$
(20)

where H is the water level or head at node (m), Qi is the flow rate through the node (m3/s), and ω is the free water surface area at node (m2).

By combining formula (19) and formula (20), the pipeline flow and water head at the node in period ∆t can be determined. Because it is a partial differential equation, it can be solved by the finite differential method.

Two-dimensional hydrodynamic model

The two-dimensional shallow water equation cannot maintain the balance between the internal force term and the source term under complex terrain conditions. To precisely maintain balance and ensure simulation accuracy in complex terrain, this paper proposed the following improved control Eq. (21):

$$\frac{\partial u}{\partial t}+\frac{\partial f}{\partial x}+\frac{\partial g}{\partial y}=S$$
(21)

in which \(u = \left( {\begin{array}{*{20}l} \eta \hfill \\ {uh} \hfill \\ {vh} \hfill \\ \end{array} } \right)\) \(f = \left( {u^{2} h + \frac{1}{2}g(\eta^{2} - 2\eta z_{b} )} \right)\)

$$ g = \left( {v^{2} h + \frac{1}{1}g(\eta^{2} - 2\eta Z_{b} )} \right)\;s = \left( {\begin{array}{*{20}l} { - C_{f} tb\sqrt {u^{2} + v^{2} } - g\eta \frac{{\partial Z_{b} }}{\partial x}} \hfill \\ { - C_{f} vt\sqrt {u^{2} + v^{2} } - g\eta \frac{{\partial Z_{b} }}{\partial y}} \hfill \\ \end{array} } \right) $$

where η is the water surface height (m); η = h + zb and h are the water depth (m); zb is the riverbed elevation (m); t represents the time (s); x and y represent the coordinates along the two directions of the plane (m); u, f, g, and S represent the fluid variable, x-direction flux, y-direction flux, and the vector of the source term, respectively; q is a point source term relevant to factors such as rainfall; v represents the velocity (m/s); \(\frac{{\partial Z_{{\text{b}}} }}{\partial x}\) and \(\frac{{\partial Z_{{\text{b}}} }}{\partial y}\) represent the bed slopes in the x and y directions, respectively; and Cf represents the bed roughness coefficient. The above governing models can be solved by using the Godunov-type finite volume method.

Construction and implementation methods of the urban HHDCM coupling model

HHDCM coupling structure

The influencing factors of urban floods are complex and include surface and underground drainage networks and river boundary conditions, which involve complex water flow exchange mechanisms among the surface, underground drainage network and river channel. A one-dimensional hydrodynamic model has the advantages of a relatively simple solution, high computational efficiency and limited modeling data. Although two-dimensional hydrodynamic models have the advantages of detailed results and are suitable for dealing with uncertainty flows, these models still have deficiencies, such as complex calculation processes, low calculation efficiency, high requirements for modeling data and difficulty in generalizing hydraulic structures. The water flow process in urban areas is complex and includes not only flows with obvious one-dimensional properties, such as pipe networks, rivers and streets, but also flows with obvious two-dimensional properties, such as surface and street intersections. It has been proven in practice that any kind of model could unavoidably involve some shortcomings in urban flow simulation. Therefore, it is necessary to couple one dimension with a two-dimensional model to construct an urban flood hydrology hydrodynamic coupling model (HHDCM). The basic structure of the HHDCM is shown in Fig. 2.

Fig. 2
figure 2

Structure of the hydrology and hydrodynamic coupling flood model (HHDCM)

Full size image

HHDCM coupling and implementation

Model coupling in the vertical direction

One- and two-dimensional model coupling in the vertical direction, namely model coupling from the surface to the subsurface, could mainly focus on flow exchange between the land surface and urban underground drainage pipe networks. In the model simulation, the node is the only channel for carrying out flow exchange between the underground drainage network and the land surface, and the purpose of vertical coupling is to calculate flow exchange at the nodes. Assuming that the node water head is H1D and the surface grid water level corresponding to the node is H2D, the flow exchange in the vertical direction can be divided into three cases according to the relationship between the status of H1D and H2D.

  1. 1.

    H1D > H2D, the flow in the pipe network system overflows up to the surface through the node;

  2. 2.

    H1D < H2D, water flows from the surface to the underground drainage pipe network;

  3. 3.

    H1D = H2D, if the surface has no accumulated water or the height of the node head is below the surface height, then there is no flow exchange between the surface and underground drainage pipe network. Therefore, only the first two cases need to be calculated.

Due to the shortage of vertical flow exchange mechanisms and calculation methods, the calculation methods are limited. Generally, the weir or orifice flow formula is used to calculate the exchange flow at a vertical node. Therefore, this paper adopted both the weir and orifice flow formulas to calculate the vertical flow exchange.

(1) Node overflow (H1D > H2D).

Considering the flow state in the pipe network, the overflow can be calculated by using the orifice flow formula. The calculation steps are as follows:

  • Assuming that the node allows for overflow, the inspection well area can be considered to be an overflow storage area; then, the water head H1D can be calculated. The surface grid water level H2D at the overflow location can be obtained by using two-dimensional model simulation results. The exchange flow can be calculated by using formula (22).

    $${Q}_{n,s}={C}_{0}{A}_{mh}\sqrt{2g({H}_{1D}-{H}_{2D})}$$
    (22)

    where \({Q}_{n,s}\) is the overflow at the current time step, m3/s; \({C}_{0}\) is the orifice flow coefficient, the value range is [0,1]; \({A}_{mh}\) is the node water storage area, m2; and g is the gravitational acceleration, m/s2.

  • According to the node water storage calculated by using the SWMM model, the overflow in the next time step should not exceed the node total storage. Moreover, the maximum allowable flow Qem should be set. Therefore, to ensure model stability, Formula (23) is used to limit the overflow.

    $$ Q_{n,s} = {\text{min}}\left( {V_{mh} /t_{n + 1} ,{ }Q_{{em,{ }}} Q_{n,s} } \right) $$
    (23)

    where \({V}_{mh}\) is the node storage, m3; and \({t}_{n+1}\) is the next time step.

  • The overflow volume can be considered the external outflow of the SWMM model node and the source term of the two-dimensional grid cell. Then, the overflow at the next time step can be calculated.


(2) Node return flow.

Both the weir and orifice flow equations are used to calculate the return flow. The specific calculation steps are as follows:

  • The node water head H1D and the grid water level H2D at the corresponding position are obtained by using SWMM and two-dimensional model simulation results. Based on the water level difference between the surface and the node, the return flow can be calculated by using Formula (24).

    $${Q}_{s,n}=\left\{\begin{array}{c}{c}_{w}\omega {h}_{2D}\sqrt{2g{H}_{2D}} {H}_{1D}\le {Z}_{2D}<{H}_{2D}\\ {c}_{0}{A}_{mh}\sqrt{2g({H}_{2D}-{H}_{1D})} {Z}_{2D} \le {H}_{1D} <{H}_{2D}\end{array}\right.$$
    (24)

    where \({c}_{w}\) is the weir flow coefficient, for which the value range is [0,1]; Qs, n, is the node return flow (m3/s); h2D is the surface water depth (m); \(\omega \) is the node circumference or rainwater outlet width (m); and Z2D is the ground elevation (m).

  • To ensure model stability, Formula (25) is used to limit the surface return flow.

    $$ Q_{s,n} = {\text{min}}\left( {{ }Q_{{s,n,{ }}} ,Q_{e,m} ,{ }V/t_{n + 1} } \right) $$
    (25)

    where V is the total water amount in the cell connected to the inspection well (m3).

The calculated return flow can be input to the SWMM model as the node external inflow; at the same time, the return flow can be considered the source item and substituted into the two-dimensional model to carry out the following step.

Model coupling in the horizontal direction

To solve the problem of flow exchange between the channel and ground, one- and two-dimensional model couplings are used to address flow exchange in the horizontal direction. According to the difference between the coupling position and flow at the coupling point, the surface one- and two-dimensional model couplings are divided into forward and lateral connections.

A forward connection means that the channel is connected to a two-dimensional area passing through the upstream and downstream regions. The coupling point is located on two sides of the river channel. The flow could be exchanged in the two-dimensional calculation area through two ends of the river channel. The flow direction at the junction point is consistent with the flow direction in the river channel, which involves the boundary conditions of the channel upstream and downstream.

The lateral connection refers to the link point between the river channel and the two-dimensional area passing through both sides of the river channel. The coupling position is located on both sides of the river channel. The water flows from two banks of the river into a two-dimensional area or into the river channel from a two-dimensional model calculation area. If the flow direction at the connection point is not consistent with the channel flow direction, a certain angle would usually form, which would not involve the upstream and downstream boundary conditions of the river channel. Therefore, two different connection modes usually require different connection strategies and computational methods. The mutual boundary method and weir flow formula method can be used to calculate the forward connection and lateral connection flow exchanges, respectively.


(1) The forward connection.

In the forward direction, the flow exchanges within a two-dimensional area through both ends of the channel. The forward connection can be calculated by using one- and two-dimensional models that mutually provide boundary conditions. The key problem for this method is determining the boundary conditions of two models at the junction point. The specific ideas and steps are as follows:

  • The two-dimensional model takes the downstream outflow of SWMM as the flow boundary condition, namely:

    $$ Q_{1D,n} = \mathop \sum \limits_{k = 1}^{M} q_{k,n + 1} l_{k} $$
    (26)

    where Q1D and n are the flow rates of the connecting section between the channel and two-dimensional area (m3/s), M is the number of grid edges connecting two dimensions and the river channel, lk is the side length of the kth grid edge (m), and qk is the unit width flow of the kth grid edge (m2/s). Based on the flow boundary calculation methods in the two-dimensional model, the Manning formula was adopted to distribute the flow to each grid edge at the junction point; see formula (27).

    $$ q_{i,n + 1} = \frac{{Q_{1D,N} \left( {h_{2D}^{5/3} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{M} \left( {h_{2D}^{5/3} l} \right)}} $$
    (27)

  • The weighted average water level of the junction unit can be used as the water level boundary condition of the SWMM.

    $$ Z_{1D,n + 1} = \mathop \sum \limits_{k = 1}^{M} \frac{{z_{k,n} l_{k} }}{L} $$
    (28)

    where Z1D, n + 1 is the water level boundary condition in the channel in the next step (m); z is the grid water level (m); and L is the total length of the forward connection boundary (m).


(2) The lateral connection.

The purpose of lateral connections is to simulate flow exchange at the junction point of two banks and in the computational area of a two-dimensional model. When the computational unit is a period of river reach, the channel flow is simulated with a one-dimensional model. When the calculation unit is a grid, the surface flow is simulated with a two-dimensional model. Therefore, the fixed boundary condition can be applied at the junction point between the two-dimensional model and the river channel. The specific steps for the lateral connection calculation are as follows.

  • When calculating water exchange, the connecting units of one- and two-dimensional models should be set separately in advance. Usually, the grid that connects a two-dimensional area and a river channel can be considered a unit. The key problem is to calculate the grid water level Hc and the channel water level Hr. The grid water level Hc can be obtained by using a two-dimensional model. The river water level needs to be interpolated in the simulation results of the SWMM. To obtain the river channel water level Hr, the upstream and downstream river channel water levels of the nodes in the calculated river reach must be extracted from the SWMM calculation results; then, the relevant river water level can be calculated via linear interpolation.

  • When the lateral exchange flow is calculated according to the surface grid water level Hc and the corresponding river channel level Hr, four conditions may exist:

    1. 1.

      When both Hc and Hr are less than the top elevation of the channel embankment Ze, no lateral flow exchange occurs, i.e., Q = 0;

    2. 2.

      When Hc > Hr and max (Hc, Hr) > Ze, lateral flow exchange occurs, and the flow direction is from the two-dimensional area to the channel;

    3. 3.

      When Hc < Hr and max (Hc, Hr) > Ze, lateral flow exchange occurs, and the flow direction is from the channel to the two-dimensional area;

    4. 4.

      When Hc = Hr > Ze, in practice, the flow direction should be determined according to the surface and river channel flow directions. When the weir flow formula is used to calculate the flow, the initial flow velocity is not considered, thus, Q = 0. The flow exchange Q can be calculated by using the weir flow formula (29).

    $$ Q = \left\{ {\begin{array}{*{20}l} {0.35b_{e} h_{{{\text{max}}}} \sqrt {2gh_{{{\text{max}}}} } } \hfill & {\frac{{h_{{{\text{min}}}} }}{{h_{{{\text{max}}}} }} \le \frac{2}{3}} \hfill \\ {0.91b_{e} h_{{{\text{min}}}} \sqrt {2g(h_{{{\text{max}}}} - { }h_{{{\text{min}}}} )_{{}} } } \hfill & {\frac{2}{3}{ } < { }\frac{{h_{{{\text{min}}}} }}{{h_{{{\text{max}}}} }} \le 1} \hfill \\ \end{array} } \right. $$
    (29)

    where \(h_{max} = max\left( {H_{r} ,H_{c} } \right) - Z_{e} ;\;h_{min} = min(H_{r} ,{ }H_{c} ) - Z_{e}\)

    where Hr and Hc are the water levels upstream and downstream of the weir, respectively (m), which are the water level of the river channel and the two-dimensional grid, respectively; Ze is the elevation of the weir, which can be set according to the actual situation (m); and be the width of the weir, which is generally the length of the junction between the grid and the channel (m).

HHDCM coupling verification-taking the Zhengzhou July.20 extreme rainstorm as an example

Background

Zhengzhou is the provincial capital city of Henan Province, China. The topography around Zhengzhou is greater in the southwest and less in the northeast. The natural height difference from east to west is more than 30 m. The surface slope in the southwest area is large, and the surface slope in the central area is relatively gentle. The average annual precipitation in Zhengzhou is 635.6 mm, and the rainfall is mainly concentrated from June to September of each year. Flood season rainfall accounts for more than 60% of the annual total rainfall. There are more than 30 rivers in the city, one each in the Yellow River and the Huaihe River basin. The Yellow River water system mainly includes the main stream of the Yellow River, the Luo River, the Si Shui River and the Kushui River. The Huaihe River system mainly includes the Jinshui River, Suoxu River, Shuangji River and Jialu River (see Fig. 3).

Fig. 3
figure 3

River system in Zhengzhou

Full size image

By the end of 2022, the urban area of Zhengzhou was projected to expand to 1055.27 km2, and the urbanization rate reached 73.4%. In 2022, the total population of Zhengzhou reached 10.136 million, and the total population in the central area exceeded 7 million. According to statistics, since 2006, the average annual extreme rainstorm and water logging disasters have occurred 15 times/year, causing an average annual economic loss of more than 200 million yuan (approximately 30.7mllion USD). On July 20, 2021, Zhengzhou experienced a historic extreme rainstorm (i.e., July.20 rainstorm) (Xinhua News Agency. 2022). From 16:00 to 17:00 on July 20, the hourly rainfall was 201.9 mm. From 20:00 on July 19 to 20:00, the 1-day (24 h) rainfall amounted to 552.5 mm. From 20:00 on July 17 to 20:00 on July 20, 3 days (72 h) of rainfall amounted to 624.1 mm. The rainfall depth per unit area is 617.1 mm, which causes severe waterlogging in most of the urban area (The State Council 7.20 Disaster Investigation Team. 2022).

Basic materials

According to statistics, during recent years, most rainstorm centers occurred in the Jinshui district of Zhengzhou, where severe flood disasters occurred. Most waterlogging points in the Zhengzhou urban area are normally located in Jinshui District. Fortunately, sufficient rainfall observation stations are distributed in the Jinshui area. Therefore, we selected a typical area in Jinshui District for this study. The selected study area is 64.42 km2. Based on the topography, buildings and drainage pipeline distribution, the entire study area was divided into 7700 sub-catchments. The one-dimensional model contains 3928 pipes and open channels, 3471 inspection wells, 12 drainage pump stations, and 24 flow outlets. The model calculation grid is an unstructured grid. The buildings within the study area were considered to be non-flooded areas. To ensure that the building area was not flooded in the simulation, the buildings were excluded from the grid area when dividing the grids. To obtain the building outline, a high-precision remote sensing map was used in this paper. The river embankment is set up as a lateral coupling boundary. Water exchange between the river channel and land surface is permitted. The computational step in the one- and two-dimensional models was 0.5 s. According to the July.20 extreme rainstorm distribution features, 17 maximum short duration rainfall sequences were selected (Table 2). The data and criteria used in this study are derived from the officially published hydrology records.

Table 2 Statistics on July.20 extreme rainstorm with different short duration
Full size table

Twelve electronic water meters were installed in the study area to monitor the waterlogging situation. Twelve liquid level meters were installed in the pipelines to monitor the water depth in the pipeline inspection wells (Fig. 4).

Fig. 4
figure 4

A sketch map of the flood inundation distribution

Full size image

Model verification results and analysis

According to the July.20 rainstorm process (Table 2), the water depths of the inspection wells in the main pipeline were selected for model verification. The Nash efficiency coefficient (ENS) (formula 30) was adopted to evaluate the model verification results.

$${\text{ENS}}=1-\frac{\sum_{{\text{t}}=1}^{{\text{t}}}{({{\text{Q}}}_{0}^{{\text{t}}}-{{\text{Q}}}_{{\text{m}}}^{{\text{t}}})}^{2}}{\sum_{{\text{t}}=1}^{{\text{T}}}{({{\text{Q}}}_{0}^{{\text{t}}}-\overline{{{\text{Q}} }_{0}})}^{2}}$$
(30)

where ENS is the Nash coefficient, Qo refers to the observed value at time t, Qm refers to the simulated value at time t, Qt (superscript) represents the value at time t, and \(\overline{{Q }_{0}}\) represents the average observed value.

The ENS value ranges from negative infinity to 1. An ENS value close to 1 indicates good simulation performance, and the model confidence is high. In addition, if the ENS is far less than 0, the model is unreliable.

Generally, for urban areas with very complex underlying surface features, if ENS > 0. 6, then the model can be considered to have good precision. The calculation results reveal that the ENS of the flood depth process in the inspection wells is 0. 914 and 0. 602, respectively (Fig. 5), which indicates that the model simulation accuracy is desirable.

Fig. 5
figure 5

Comparison of measured and simulated flood hydrographs in the inspection wells

Full size image

The flood location is also an important criterion for model rating. Table 3 shows the comparative analysis of the simulated and measured flood locations, which shows that the measured flood locations are perfectly consistent with the simulated water depth except at site S4.

Table 3 Comparison between simulated and measured flood locations in July.20
Full size table

Figure 4 shows the urban inundation depth distribution caused by July.20 rainstorm, which shows that the simulation results basically agree with the measured results.

Table 4 shows the differences between the simulated and measured floodwater depths, which shows that according to a comparison between the measured and simulated results, the maximum errors were discovered at points S3 and S10. The relative error at point S3 is − 12.5%, and the relative error at the S10 point is 12.5%. It can also be seen that the simulation errors at most points are less than 10%. The simulation relative errors are less than 5% for more than 30% of the measurement points.

Table 4 Error analysis of the model simulation results
Full size table

Therefore, the urban flood model studied in this paper can reflect the real situation of rainstorms flood disasters and has good reliability and accuracy.

Discussions

  • The two-dimensional shallow water equation cannot maintain the balance between the internal force term and the source term under complex terrain conditions. To precisely maintain balance and ensure simulation accuracy in complex terrain, this paper proposed improved control equations. The control equations can be solved by using the Godunov-type finite volume method, which can save calculation time and improve accuracy.

  • Although two-dimensional hydrodynamic models have the advantages of detailed results and are suitable for dealing with uncertain flow, these models still have deficiencies, such as complex calculation processes, low calculation efficiency, high requirements for modeling data and difficulty in generalizing hydraulic structures. In this paper, we coupled a one-dimensional model with a two-dimensional model to construct an urban flood hydrology hydrodynamic coupling model (HHDCM).

  • HHDCM coupling and implementation were carried out both from the vertical and horizontal directions, which solved the problem of flow exchange between the channel and surface. The flood simulation results were improved.

  • According to the difference between the coupling position and flow at the coupling point, the surface one- and two-dimensional model couplings are further divided into forward and lateral connections. The forward connection can be calculated by using one- and two-dimensional models that mutually provide boundary conditions. The lateral connection is used to simulate flow exchange at the junction point of two banks and in the computational area of the two-dimensional model. Therefore, the fixed boundary condition can be applied at the junction point between the two-dimensional model and the river channel.

  • Compared with previous methods (Lu et al. 2019), the simulation accuracy clearly improved. Table 4 shows that the maximum relative error for flood depth is 12.5%.

Conclusions

In this paper, based on the improved SWMM and hydrodynamic models, one- and two-dimensional model coupling was studied from the lateral, forward and vertical direction. The urban flood hydrology hydrodynamic coupling model (HHDCM) was constructed by combining hydrology simulation with a pipe network and node hydrodynamic model, a river channel hydrodynamic model and two-dimensional urban flood disaster models.

  1. 1.

    The SWMM structure, model parameters and surface confluence process simulation were improved, and relevant SWMM parameters were proposed.

  2. 2.

    The original pipe network confluence model was improved, and a new two-dimensional hydrodynamic model was constructed, which improved the simulation accuracy of complex topography.

  3. 3.

    The coupling structure and implementation method of the urban SWMM and two-dimensional hydrodynamic model were studied, and a detailed coupling form of the SWMM and two-dimensional hydrodynamic model was proposed.

  4. 4.

    Model coupling from different directions (vertical, horizontal, forward, lateral, and node flow exchange) was carried out to realize desirable model coupling results between the SWMM and two-dimensional hydrodynamic model.

  5. 5.

    A total of July.20 rainstorms in Zhengzhou were used as a study case, and different rainfall durations during the rainstorm process were selected for model verification and calibration. The verification results reveal that the simulated flood depth and area are consistent with the measured results. The maximum relative error is 12.5%. The verification results verify that the model has good accuracy and reliability and can be applied to simulate urban rainstorms and flood disasters.

  6. 6.

    The studied results would play an important reference role in urban extreme flood simulation and flood control. Some parameters could be further verified in practice in the future.