Hydrodynamic Simulation of Urban Waterlogging Based on an Improved Vertical Flow Exchange Method

Abstract

In the 1D–2D coupled simulation of urban waterlogging, the calculation process of vertical flow exchange is independent from the 1D hydraulic calculation, resulting in a failure to consider the node head and pipe flow during the exchange flow calculation, which may lead to irrational results and further affect the stability of the model calculation. However, setting an upper limit for the exchange flow may introduce excessive subjective factors into the simulation process. In this study, a vertical flow exchange method based on the water balance of nodes is proposed. When a node is in an overloaded state, the calculation of vertical flow exchange at the node is integrated into the 1D hydraulic simulation process, thus taking into consideration the influence of the node head and pipe flow when calculating vertical flow exchange. Additionally, the iterative solution method used in the 1D hydraulic model ensures numerical harmony between the vertical flow exchange, node head and pipe flow, thus ensuring the stability of the coupled calculation. For the non-overloaded nodes, the calculation of the vertical flow exchange was conducted using a variable-head orifice discharge formula, enabling the consideration of changes in the surface water depth during the calculation of the node backflow. Using the InfoWorks ICM model as a benchmark, a comparative analysis of case simulation results demonstrated that the improved vertical flow exchange method was able to accurately and stably simulate the process of vertical flow exchange. When used with the improved vertical exchange method, the coupled model gave simulation results that closely matched those of the benchmark model.

1. Introduction

In recent years, with the increasing frequency of extreme rainfall disasters, more cities have implemented planning and construction for drainage and flood control. In this process, urban hydrodynamic models have been widely applied [1,2,3,4]. Early hydrodynamic models focused on simulating surface runoff processes and hydrodynamic processes in underground drainage systems, but they were unable to simulate surface waterlogging [5,6]. In order to address this deficiency, Hsu M. H. (2000) used the SWMM (Storm Water Management Model) as a one-dimensional (1D) model coupled with a two-dimensional (2D) diffusive overland flow model, in order to achieve a coupled urban waterlogging simulation [7]. However, in this method, the flow exchange between models was unidirectional. When the drainage network was overloaded, water flowed from the drainage system to the ground surface; when the drainage system had sufficient discharge capacity, surface water on the ground could not flow back to the drainage system. In further research, a series of coupled models capable of bidirectional exchange between an underground drainage system and surface water were developed [8,9,10]. Although coupled models have the advantages of detailed results and adeptness in handling simulations of the bidirectional exchange of water flows, their computational processes are complex and their efficiency is low. In order to achieve a balance between simulation accuracy and computational efficiency, some researchers explored waterlogging simulation methods that could avoid solving hydrodynamic equations yet still maintain good accuracy, such as with flood analysis algorithms based on GIS [11], an enhanced inundation model for flood calculations [12], and some simple models of cellular automata [13]. However, these methods can only determine the extent of waterlogging and water depth, but are unable to simulate the hydrodynamic processes of waterlogging, so their application scope is limited.

Currently, the coupled 1D–2D model is the main method for simulating urban waterlogging. Vertical flow exchange is a key segment in the coupled process, and it significantly influences the accuracy of the simulation results and the water balance between the 1D and 2D models [14,15,16]. The current calculation method for vertical flow exchange mainly uses the head difference between the water levels of nodes in the 1D model (H_1D) and the surface water levels of the 2D model (H_2D) as the driving head, and orifice or weir flow formulas are utilized to calculate the flow exchange between the 1D and 2D models [17].

The typical calculation formula for node overflow is as follows [2]:

$$Q_{of}=c_o{A}_{mh}\sqrt{2g(H_{1D}-H_{2D})} \mathrm{i}\mathrm{f}{ H}_{1D}>H_{2D }$$

(1)

where Q_of is the calculated overflow value (m³/s); c_o is the orifice discharge coefficient; A_mh is the manhole area of nodes (m²); g is the gravitational acceleration (m/s²); H_1D is the water level of the 1D model’s nodes (m); and H_2D is the water level of the 2D model (m).

The typical calculation formula for the node backflow is as follows [2]:

$$Q_{bf}=\left\{\begin{array}{c}{c}_{w}W{H}_{2D}\sqrt{2g{H}_{2D}} \mathrm{if}{ H}_{1D}{\le Z}_{2D}<H_{2D} \mathrm{and} Q_{bf}<Q_{bm} \\ c_{\mathrm{o}}{A}_{mh}\sqrt{2g\left(H_{1D}-H_{2D}\right)} \mathrm{if}{ Z}_{2D}\le H_{1D}<H_{2D} \mathrm{and} Q_{bf}<Q_{bm}\\ Q_{bm} \mathrm{if}{ Q}_{bf}\ge Q_{bm} \end{array}\right.$$

(2)

where Q_bf is the calculated value of the backflow (m³/s), and in order to achieve a stable simulation, the value of Q_bf cannot exceed the maximum backflow rate of Q_bm; c_w is the weir discharge coefficient; W is the perimeter of the node or the width of the rainwater outlet (m); and Z_2D is the ground elevation where the node is located (m).

When calculating node backflow or intercepted discharge by inlets, in order to consider the water interactions between the surface and drainage system more accurately, water energy upstream of the node or inlet (sum of the head difference and velocity head) is used as the driving head for Formula (2) [18]. Furthermore, some studies have discussed the relationship between the discharge coefficient and the surface flow pattern, and have concluded that the orifice and weir discharge coefficient have a significant correlation with the Froude number of surface flow [18,19]. The vertical water exchange method takes into consideration that velocity head and flow pattern are often used for detailed 1D models of drainage systems, including the full drainage structure composed of conduits, manholes, gullies and inlets. These are coupled with 2D models in a fully distributed way [20,21], where the runoff volumes are estimated and applied directly to the elements of the 2D model of the overland surface, and exchange with the 1D model through inlets. However, due to limitations in data availability and computational power, the Fully Distributed models are generally suitable for small-scale cases. For medium to large-scale models, the Semi-Distributed Urban Stormwater models remain the preferred choice. In Semi-Distributed models, conceptual, empirical or physical-based methods transform runoff routing into inflow hydrographs, which are then applied to the selected computational nodes of the drainage system. Not every inlet is modeled, but they are clustered to computational ones. In this case, the discharge coefficient has been decoupled from specific inlets and has become a subjectively determined parameter that needs to be calibrated using actual measurement data.

The determination of discharge coefficients for flow exchange equations over a range of hydraulic conditions and inlet types has a significant impact on the simulation results. Recent studies related to water exchange through manholes or inlets using the weir and orifice-type equations have been implemented using physical model experiments [22,23,24] and numerical simulations [25]. The values of the discharge coefficients are influenced by the inlet geometry and the approaching flow characteristics (flow depth and Froude number). Rubinato et al. [19] studied water exchange of scaled circular manholes under subcritical flow with a Froude number of 0.151–0.691; the results revealed that under subcritical flow conditions, water exchange is sensitive to hydraulic heads within interaction nodes and surface water. Hence, the uncertainty related to head loss in hydraulic structures, surface roughness and other parameters in urban hydrology may have significant implications to water exchange and should be considered in the selection of weir/orifice formulas and the determination of discharge coefficients. Under supercritical flow conditions, the flow velocity parameter has to be considered while determining the discharge coefficients. Cosco et al.‘s study [18] highlighted the strong relationship between discharge coefficient ranges and the upstream Froude number under supercritical flow, and formulated power functions in order to express the weir/orifice discharge coefficients as a function of Froude number. Ali Zaiter et al. (2024) gave a comprehensive review about “equations and methodologies of inlet drainage system discharge coefficients” [26]. Ali’s review summarized orifice and weir discharge coefficient ranges of rectangular inlets, grated circular inlets and circular non-grated inlets under subcritical flow and supercritical flow, respectively, and concluded that there is no uniform pattern or trend to determine the discharge coefficient of each inlet as it may vary significantly with the change in Froude number, inlet geometry, grate geometry, grate orientation, flow velocity and water depth. Therefore, weir and orifice discharge coefficients have to be calibrated and adjusted experimentally, taking into consideration the crucial factors that affect the inlet discharge capacity [27]. Russo et al. (2015) analyzed extreme floods in Barcelona using the 1D–2D coupled model created by InfoWorks ICM, and gave a detailed description of how to calculate and calibrate the discharge coefficients of manholes and storm inlets [28].

The calculation of flow exchange through orifice or weir formulas is not only complex, but also involves significant uncertainty in the selection of discharge coefficients. Moreover, the introduction of an artificially set value of Q_bm adds subjectivity and uncertainty to the simulation of coupled models. The fundamental reason for the instability in the coupling process when not using a limitation for Q_bm is that the calculation of the backflow is independent of the 1D model. Thus, interactions among variables such as the vertical flow exchange, node head and conduit flow are neglected, which may lead to irrational backflow values. When an irrational value is involved in the coupling process, computational instability may occur. Fan Y. (2017) suggested determining the Q_bm based on the simulation stability of the model and the actual situation of the node [29]. According to this method, the appropriate value of Q_bm can only be obtained through trial and error during the simulation process. In order to implement a stable coupled simulation without artificially set values, Peng G. (2022) calculated the backflow while considering the free space of a manhole and connected conduits [30]; however, this was a compromise solution and the interactions among the vertical flow exchange, node head and conduit flow were still neglected.

This study introduces the principle of node inflow and outflow balance into the calculation process for vertical flow exchange, and the calculation of vertical flow exchange is integrated into the 1D model in order to consider interactions among the vertical flow exchange, node head and conduit flow so that harmonious results can be obtained. This improved method simplifies the calculation process for vertical flow exchange and enhances the stability of the coupled process.

2. Coupled 1D–2D Hydrodynamic Model

2.1. Hydrodynamic Models

The SWMM is a hydrologic–hydraulic water quality simulation model that utilizes dynamic modeling techniques [31]. As an open-source model, it is widely used in the field of urban stormwater simulation. The SWMM is suitable for the 1D hydrodynamic simulation of drainage systems in urban areas. In this study, the SWMM is used as a 1D model for the hydrological simulation of sub-catchments and the hydrodynamic calculation of drainage systems.

Conservative shallow water equations are used as the governing equations in the 2D model, and they are solved using a self-developed simulation module with the Godunov scheme for the finite volume method [32]. Conservative shallow water equations are adopted as the governing equations in the 2D model. The governing equations can be written in matrix form, as shown in Formulas (3) and (4).

$$\frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial y}=S_f+S_b$$

(3)

where U, F and G represent the conservation vector, the flux vectors in the x and y directions, respectively, and S_b and S_f are the source terms representing the bed slope and friction effects, respectively. In addition, t, x, and y represent time and the coordinate axes of a Cartesian coordinate system. The specific forms of U, F, G, S_b and S_f are shown in Formula (4).

$$\begin{array}{c}U=\left[\begin{array}{c}h\\ hu\\ hv\end{array}\right],F=\left[\begin{array}{c}hu\\ {hu}^2+\frac{1}{2}g{h}^2\\ huv\end{array}\right]\end{array},G=\left[\begin{array}{c}hv\\ huv\\ {hv}^2+\frac{1}{2}g{h}^2\end{array}\right],S_b=\left[\begin{array}{c}0\\ -gh\frac{\partial z_b}{\partial x}\\ -gh\frac{\partial z_b}{\partial y}\end{array}\right],S_f=\left[\begin{array}{c}0\\ -C_{f}u\sqrt{u^2+v^2}\\ -C_{f}v\sqrt{u^2+v^2}\end{array}\right]$$

(4)

where h represents water depth (m); u and v represent the flow velocities in the x and y directions, respectively (m/s); z_b represents the absolute elevation of the ground surface (m); and $C_f=\frac{{gn}^2}{h^{1/3}}$, in which n is the Manning coefficient. Additionally, water level η is used in the second order spatial reconstruction and the non-negative water depth reconstruction—η = h + z_b.

Formulas (3) and (4) are solved herein by a finite volume Godunov-type numerical scheme proposed by Liang (2010), which has an emphasis of constructing a well-balanced scheme and evaluating the friction source term in order to ensure non-negativity of water depth [33,34]. The monotonic upstream-centered scheme for conservation laws (MUSCL) linear reconstruction method [35] with a MINMOD slope limiter [36] is used along with a Harten, Lax and van Leer approximate Riemann solver, with the contact wave restored (HLLC) for the calculation of interface fluxes [37]. The second-order temporal accuracy is achieved by using a Runge–Kutta time integration [38].

During the coupling process of the 1D and 2D models, two main aspects are involved: time synchronization and flow exchange. For time synchronization, this study used the 1D model time step as the benchmark and achieved time synchronization by adjusting the 2D model time step. For flow exchange, this study only considered vertical flow exchange, while lateral and frontal flow exchanges were neglected. Closed boundary conditions were adopted at the boundaries of the surface grid for the 2D model.

2.2. Time Synchronization in the Coupled 1D and 2D Models

Since the time step of the 1D model (Δt_1D) is often larger than that of the 2D model (Δt_2D), Δt_1D is taken as a benchmark. By adjusting Δt_2D, the synchronization of the time steps in the 1D and 2D models is achieved. During the simulation process, when the 1D model advances by one time step, the 2D model runs n steps. The first n − 1 steps of the 2D model are calculated using the Courant–Friedrichs–Lewy (CFL) criterion, and the value of Δt_2D for the nth step is calculated as Δt_1D − $\sum _1^{n-1}{\Delta t}_{2D}$. The synchronization method for the 1D and 2D models is presented in Figure 1.

3. Improved Vertical Flow Exchange Method

The three typical cases of vertical flow exchange at nodes are shown in Figure 2.

In Figure 2, cases (a) and (b) represent situations of overloaded nodes, where the water head inside the node (same as the water level of the surface cell) is greater than the maximum at the upper edge of the connected conduit section, which is denoted as H_P. Since the 1D model generally ignores the volume inside node wells, this implies that the node water head H_J exceeds the ground elevation H₀. In such cases, the calculation of conduit flow in the 1D model should consider the surface water level of the 2D model at the node location. The inflow and outflow of the node are calculated according to the conduit flow values. The difference between the inflow and outflow is considered as the vertical flow exchange between the 1D and 2D models. If the value of the vertical flow exchange is positive, this indicates that the node is an overflow node. If the value is negative, this indicates that the node is a backflow node. If the value is equal to 0, this indicates that there is no vertical flow exchange in this node. Case (c) represents the situation of a non-overloaded node, where the node water head H_J is less than the maximum at the upper edge H_P of the section of connected conduits. In this case, the calculation of pipe flows in the 1D model was independent of the surface water level, and the backflow of the node was calculated using the orifice flow formula and was added as the inflow to the node.

3.1. Calculation of the Vertical Flow Exchange

The nodes of the 1D model are associated with the mesh cells of the 2D model in which they are located, with the mesh cell area serving as the ponded area for each node.

Figure 3 illustrates the relationship between the ponded area of the node and the mesh cell. The dark blue cell represents the mesh cell corresponding to the location of the node. Before the coupling simulation begins, the ponded area of the corresponding node is updated using the area of this cell. This cell not only participates in the calculations of the 2D model, but also contributes to the computation of the vertical flow exchange. The remaining cells are ordinary cells that only participate in the calculations of the 2D model.

3.1.1. Calculation of the Vertical Flow Exchange in Overloaded Nodes

When a node is in an overloaded state, the calculation of the flow in the connected conduits needs to consider the influence of the surface water level. In 1D hydraulic simulations, before the calculation of the dynamic wave in each time step begins, the water level of the mesh cell in which the node is located is used to update the node’s water head in the 1D model. Then, the 1D hydraulic calculation of the dynamic wave is performed. The SWMM module uses an iterative method to calculate node water heads and conduit flow rates, with the convergence of node water heads as the stopping condition for iteration, and in each iteration, the vertical flow exchange is updated by applying the principle of water conservation. This ensures numerical harmony among the node water heads, conduit flows and vertical exchange flows. Thus, stability in the calculation of the vertical flow exchange is achieved. According to the principle of water conservation, the following relationship between the change in water volume within the ponded area and the flow in the connected conduits of the node exists:

$$V_{new}-V_{old}={∆t}_{1D}({\sum }_{i=1}^k{Q}_i+Q_{lat})$$

(5)

where V_old represents the water volume within the ponded area of the node at the beginning of the time step (m³); V_new represents the water volume within the ponded area of the node at the end of the time step (m³); Q_i is the flow in the ith conduit connected to the node (m³/s), with inflow as positive and outflow as negative; Q_lat is the sum of other inflows, such as surface runoff and dry weather flow (m³/s); and k is the number of conduits connected to the node.

The change in the water volume within the ponded area of the node in the 1D model is the volume of vertical exchange between the 1D and 2D models, so the vertical exchange flow can be calculated as follows:

$$Q_{exf}=\frac{V_{new}-V_{old}}{{∆t}_{1D}}$$

(6)

where Q_exf represents the vertical exchange flow (m³/s), with positive values indicating node overflow and negative values indicating node backflow.

3.1.2. Calculation of Vertical Flow Exchange in Non-Overloaded Nodes

When a node is not overloaded but there is ponding water on the surface cell, the flow calculation for the connected conduits does not need to consider the influence of the surface water level. Therefore, there is no need to update the node’s water level, and the orifice flow formula is used here for the backflow calculation. During the backflow process, the surface water level may change. As a result, the following formula for the backflow calculation considering a variable head is used:

$$t_b=\frac{2S\left(\sqrt{H_1}-\sqrt{H_2}\right)}{c_o{A}_{mh}\sqrt{2g}}$$

(7)

where t_b represents the backflow time (s); S represents the ponded area (m²); H₁ represents the surface water level at the node at the beginning of the backflow time t_b (m); and H₂ represents the surface water level at the node at the end of the backflow time t_b (m).

Assuming that H₂ = 0, Formula (7) is used to calculate the time t_b, which represents the time required for all of the ponding water in the cell to completely backflow. If t_b > Δt_1D, this indicates that within the current time step in the SWMM, not all of the ponding water in the cell backflows. Then, Formulas (8) and (9) are utilized in order to calculate the water level of the surface water and the backflow rate at the end of the current time step in the SWMM.

$${H_2=\left(\sqrt{H_1}-\frac{c_{o}A{∆t}_{1D}\sqrt{2g}}{2S}\right)}^2$$

(8)

$$Q_{exf}=\frac{S\left(H_1-H_2\right)}{{∆t}_{1D}}$$

(9)

If t_b ≤ Δt_1D, this indicates that within the current time step in the SWMM, all of the ponding water in the cell has completely returned into the underground network. Then, the backflow rate at the end of the current time step is calculated with Formula (10):

$$Q_{exf}=\frac{S{H}_1}{{∆t}_{1D}}$$

(10)

If the node is in an overloaded state at the end of Δt_1D, this indicates that it has transitioned from a backflow state to an overflow state within this time step. In this case, the vertical flow exchange of the node is calculated using the principle of conservation of inflow and outflow according to Formula (6).

3.1.3. Correction of Continuity Errors Caused by Node Water Head Updates

Figure 4 illustrates the meanings of various variables during the process of updating the node water head. In the coupling process, the 1D model’s time step $∆t_{1D}^n$ is calculated first. When a node is in an overflow state (shown in Figure 4a), the water head of the node changes from H_old to H_new within the time step. The water volume between these two levels represents the overflow volume of the 1D model. During the process of vertical flow exchange, this volume is transferred from the 1D model to the 2D model. Consequently, the volume of water in the 2D model increases accordingly. In order to ensure water conservation in the coupling process, this volume needs to be removed from the 1D model. Therefore, at the end of $∆t_{1D}^n$, the water head and ponding water volume corresponding to the node in the 1D model should be H_old.

After the overflow volume is added to the 2D model and the 2D model’s corresponding time step is completed, the water level in the 2D cells becomes H_2D. Before the calculation of the next time step in the 1D model ($∆t_{1D}^{n+1}$) begins, it is necessary to update the node water head using H_2D. The water head at the node changes from H_old to H_2D, and the ponding water volume also increases accordingly. The increment in this ponding water volume is equivalent to adding additional virtual water to the 1D model without corresponding to any actual physical process, which may lead to continuity errors.

In order to eliminate this continuity error, it is necessary to eliminate the volume of water between H_old and H_2D during the calculation of the flow exchange in time step $∆t_{1D}^{n+1}$. The elimination method involves obtaining the initial result of the flow exchange Q_exf0 using Formula (6) and then correcting it using Formula (11), in order to eliminate the continuity error.

$$Q_{exf}=Q_{exf0}-\frac{{(H}_{2D}-H_{old})S}{{∆t}_{1D}}$$

(11)

When the node is in an overloaded backflow state (shown in Figure 4b), the water head at the node changes from H_old to H_new. At the end of $∆t_{1D}^n$ in the 1D model, the water head and ponding water volume of this node are calculated according to H_new. The water volume between H_old and H_new represents the vertical backflow volume, which will be removed from the corresponding 2D cells during the respective time step in the 2D model. Updating the node water head using H_2D before the calculation of the next time step in the 1D model $∆t_{1D}^{n+1}$ would cause an increase in the ponding water volume in the 1D model’s nodes. The water head changes from H_new to H_2D, and the ponding water volume changes accordingly. This increase in ponding water volume can be considered as the volume of surface water from the 2D model flowing back to the ponded area of the 1D model’s nodes. This volume of water has actual physical significance and does not lead to continuity errors; therefore, no correction is needed.

3.2. Implementation of Model Coupling

In the coupled model utilizing the improved vertical flow exchange method, the SWMM requires the node ponding functionality to be enabled. The specific calculation steps during the coupling process corresponding to Δt_1D are as follows:

i.

Determining which nodes are in an overloaded state;

ii.

Updating the water heads for overloaded nodes;

iii.

Calculating the backflow for non-overloaded nodes (with the method described in Section 3.1.2) and adding backflow to the corresponding node as inflow;

iv.

One-step evolution of the 1D model (the vertical flow exchange for overloaded nodes is calculated with the method described in Section 3.1.1 by adding the corresponding functions);

v.

Eliminating continuity errors;

vi.

Transferring the vertical flow exchange values to the 2D model;

vii.

Evolution of the 2D model until the end of Δt_1D. This is followed by step i.

Since the SWMM module has the function of calculating the ponding water volume at nodes, the vertical exchange flow under the condition of overloaded nodes can be directly obtained through the change in the ponding water volume at these nodes. Steps (i) to (v) can be implemented within the SWMM module by adding corresponding functions to it.

4. Case Study

4.1. Benchmark Software

In this study, InfoWorks ICM (from here on ‘ICM’) was used as a benchmark. ICM solves the complete 2D Saint Venant equations in a finite volume semi-implicit scheme [35] with a Riemann solver [39]. ICM combines a number of distinctive features such as the analysis and prediction of potential flood extent, depth and velocity, and the modeling of the interaction of surface and underground systems in a fully integrated environment [28].

Generally, in a 1D–2D coupled model created by ICM, the calculation of vertical flow exchange between 1D and 2D models is achieved through nodes of 2D and Gully 2D flood types. For 2D nodes, the exchange of water between the 1D and 2D models is calculated using weir equations, assuming a weir crest level at the node ground level and crest length equal to the node shaft circumference. For Gully 2D nodes, the exchange flow is calculated using a specific (customized) “head/discharge relationship”. The Gully 2D flood-type nodes more accurately describe the water exchange between the 1D drainage system and the 2D surface at different rainwater inlets.

ICM can conduct both Fully Distributed simulations and Semi-Distributed simulations. In this study, only Semi-Distributed simulations were performed and discussed. For theoretical study cases, the nodes with vertical flow exchange do not represent any specific type of rainwater inlet. For practical study cases, the drainage network used was a reduced network, only conduits and computational nodes were included. Runoff volumes were estimated for the sub-catchment and then routed to the corresponding nodes as sub-catchments’ outlets. In this case, only 2D flood type were used for water exchange calculation of both study cases, and the typical discharge value of 0.5 was used for all 2D flood type nodes.

4.2. Theoretical Study Case

The improved method was validated using a theoretical case study, as shown in Figure 5. The drainage system consisted of six nodes and six conduits, and the level surface area was a closed square with a side length of 200 m. Nodes 2 to 5 allowed vertical flow exchange. The elevation of the surface was 0 m, and the roughness coefficient (n value) was 0.025. Water entered the system through node 1, with the inflow gradually increasing from 0 to 1.0 m³/s within the first 10 min of the simulation, remaining constant until the 20th hour and then dropping to 0. Water left the system through node 6, with the outfall set as a free outflow. The system was initially dry. The attributes of the conduits and nodes are shown in

Table 1. Properties of the elements of the drainage system.

ID	Conduit					Node
	Upstream Node	Downstream Node	Dia /m	Len /m	Roughness	X Coordinate	Y Coordinate	Full Depth/m	Ground Elevation/m	Well Dia /m
1	1	2	1.0	100	0.013	−70.711	100.000	1.159	0	0.7
2	2	4	0.6	100	0.013	29.289	100.000	1.259	0	0.7
3	2	3	0.6	100	0.013	100.000	170.711	1.329	0	0.7
4	4	5	0.7	100	0.013	100.000	29.289	1.329	0	0.7
5	3	5	0.7	100	0.013	170.711	100.000	1.400	0	0.7
6	5	6	1.0	100	0.013	270.711	100.000	1.500	0	-

.

In the coupled model proposed in this study, the area of the mesh cell where the node is located was used as the ponded area of the node. Therefore, the impact of the grid mesh on the coupled process is important. Two types of grid meshing structures (2 m × 2 m (S4) and 4 m × 4 m (S16)) were used for the theoretical case.

Figure 6a shows the water balance in the improved model with grid meshing of cell areas of 16 m². During the process of vertical flow exchange, the accumulated overflow volume from the nodes should be equal to the sum of the surface water volume and the accumulated backflow volume from the nodes. Figure 6a shows that the improved model maintained water conservation during the vertical flow exchange process, as the accumulated overflow curve of the nodes perfectly overlapped with the curve of “accumulated backflow + surface water accumulation”. In the process of the coupled 1D–2D simulation, the total water entering the coupled model (accumulated system inflow) needed to be equal to the sum of the water leaving the system (accumulated system outflow) and the water stored within the model (surface water). Figure 6b demonstrates that the overall water conservation of the model was also well satisfied.

Root mean squared error (RMSE) was used as a statistical measure of deviation between “accumulated overflow” with “accumulated backflow + surface water accumulation” for Figure 6a, and “accumulated system inflow” with “accumulated system outflow + surface water” for Figure 6b.

Figure 7 takes node 2 as a typical overflow node, node 3 as a typical backflow node, and the locations of point A and node 3 as representatives for the surface water depth. It compares the calculation results of the improved model and the ICM model with two types of grid meshing. The Nash–Sutcliffe efficiency coefficient (E_NS) was used to quantitatively evaluate the degree of agreement between the simulation results of the two models. It can be seen that the calculation results of the improved vertical flow exchange method with different grid meshes had a high degree of agreement with the results of the ICM model. With the meter-level and ten-meter-level grid meshes, the coupled calculation results obtained with the improved vertical flow exchange method were basically consistent with the results obtained with the ICM model. It can be seen in Figure 7c that the difference in grid meshing had a certain impact on the simulation results for the surface water depth, but the improved method was able to obtain a water depth change curve that was almost consistent with that of the ICM model.

4.3. Practical Study Case

The Dongsha storm drainage system in Wuchang District, Wuhan City, was used as a practical case for the simulation of urban waterlogging; the study area covers a catchment range of 82.73 km². The terrain within the catchment area is relatively flat, with Shahu Lake and Donghu Lake located on the south side. These lakes are the main recipients of runoff. However, as the lake water levels remain relatively high throughout the year, it adversely affects the drainage of the stormwater drainage network, making local low-lying areas prone to waterlogging.

The terrain and drainage system of the study area are shown in the right subfigure of Figure 8. The ICM model and the improved model were both used to establish a coupled 1D–2D model for the study area. An unstructured grid with an average cell size of 450 m² was adopted as the ICM model’s ground grid, and the total number of cells was 210,000. The improved model used a structured grid for the ground, with a cell size of 20 m × 20 m and a total of 240,000 cells. The rainfall data used were from two heavy rainstorms in the area, which occurred on 18 June and 6 July 2016 (hereinafter referred to as the 201606 and 201607 rainstorms). The cumulative rainfall amounts for the two storms were 143.5 mm and 125.6 mm, respectively, and the rainfall intensity curves for both storms is shown in the two left subfigures in Figure 8.

Both coupled models were used to simulate the waterlogging conditions for the two rainstorms. The inundated area of the improved model is shown in Figure 9. The 201606 rainstorm, with its heavy rainfall and short duration, led to severe urban waterlogging. The simulation results indicated many inundated areas and significant water depths on the ground. For the 201607 rainstorm, although it also had heavy rainfall, its longer duration and more uniform distribution of rainfall throughout the event, with a peak intensity approximately half of that of the 201606 rainstorm, resulted in less severe waterlogging and fewer recorded waterlogging points on that day. In both rainstorm simulations, the simulation results of the improved model were basically consistent with the actual recorded waterlogging points on the respective days. The simulation results demonstrated that the coupled model using the improved vertical flow exchange method was able to accurately simulate the inundation range and depth of urban waterlogging.

Figure 10 presents the temporal changes in the surface water volume for both the improved model and the ICM model during the two rainfall events. It can be observed that the vertical flow exchange calculations from the improved model closely aligned with those from the ICM model, as indicated by the high degree of overlap in the curves depicting the temporal variation in the volume of groundwater accumulation. This suggests that the vertical flow exchange calculations at the nodes using the improved method were highly consistent with those from the ICM model.

In the urban waterlogging simulations, the vertical flow exchange at overloaded nodes accounted for most of the overall vertical flow exchange. Figure 11 shows a comparison of the number of different types of vertical flow exchange nodes in the evolution of the 1D model during the coupling process. Based on the statistical analysis of the simulation results, the proportion of vertical flow exchange at overloaded nodes, with respect to the total vertical flow exchange at nodes during the two rainstorm simulations, was 81.43% and 88.66%, respectively. Using the improved method allowed for the direct calculation of the vertical flow exchange results at nodes using the SWMM module, based on the principle of water conservation in most cases. The calculation of the backflow volume using a variable-head orifice flow formula was only needed in a few instances of backflow at non-overloaded nodes. This significantly reduced the amount of vertical flow exchange calculations using the orifice or weir flow formulas, thereby improving the computational efficiency and stability of the coupled simulation.

5. Conclusions and Future Work

Based on the SWMM and an independently developed two-dimensional hydrodynamic model, a coupled hydrological and hydrodynamic model for urban waterlogging was constructed. Within this framework, the vertical flow exchange method was improved using the principle of node water balance. Both a theoretical drainage system and a real urban drainage system were used as case studies in order to validate the improved method in terms of simulation accuracy, reliability and water balance. The main conclusions are as follows.

Calculating the vertical flow exchange using the principle of water conservation can significantly reduce the calculations related to the orifice or weir formulas. Most calculations related to the vertical flow exchange can be implemented with the original functions in the SWMM. Only the backflow at non-overloaded nodes needs to be calculated outside the SWMM by using orifice or weir flow formulas.

The calculation of the vertical flow exchange is integrated into the SWMM calculation module. The SWMM module uses an iterative method to calculate the node water heads and conduit flow rates, with the convergence of node water heads as the stopping condition for each iteration. Thus, numerical harmony is achieved among the node water heads, conduit flows and vertical exchange flows. In addition, numerical stability of the simulation is achieved without any artificially set values in the calculation of the vertical flow exchange.

A comparison of the simulation results from the ICM model with those of the improved model demonstrates that the coupled model using the improved vertical flow exchange method exhibits good accuracy and reliability. Its results align closely with those of the ICM simulation across various grid sizes ranging from meters to hundreds of meters, making this method suitable for the numerical simulation of urban waterlogging.

However, this study still has certain limitations. It only discusses the calculation of vertical flow exchange in the case of Semi-Distributed Urban Stormwater Models. The nodes involved in vertical flow exchange are computational nodes, not specific rainwater inlets. Additionally, due to insufficient monitoring data, there has been no detailed calibration or discussion of the flow coefficients at the exchange nodes in both the ICM and proposed models. These shortcomings need to be addressed in future work.