Urban flood modeling with a novel coupling method of surface and sewer hydrodynamic processes

The global climate change and urbanization process have contributed to an increased frequency of flooding, which has become the most frequent weather-related hazard in the last two decades, affecting 2.3 billion people worldwide, more than all other weather-related hazards combined (Wallemacq et al. 2015). At the same time, flooding is one of the world's most destructive natural disasters, not only endangering lives and property but also having a serious negative impact on society and the environment (Han & He 2021; Dong et al. 2022). It is essential to develop an efficient and trustworthy urban drainage model appropriate for various rainfall scenarios and diverse geographical regions to comprehend the development progress of urban flooding, which is the foundation for mitigating flood damage.

The sewer network flow and the overland surface flow constitute two distinct subsystems that jointly make up the urban drainage system. The simulation of overland surface flow is typically carried out with either simplified conceptual model-based hydrological methods or numerical model-based hydrodynamic methods (including one- and two-dimensional hydraulic models). Conceptual model-based hydrological and 1D hydrodynamic approaches for solving the shallow water equations (SWEs), while computationally efficient, are not as accurate as physically based 2D hydrodynamic methods (Djordjević et al. 1999; Teng et al. 2017; Sañudo et al. 2020). Modeling sewer flow is usually accomplished by solving the 1D Saint-Venant equations to examine the ability of a sewer system to cope with intense rainfall (Djordjević et al. 1999; Rossman & Simon 2022). The 1D model simulates the water hydrodynamics in the sewer network with high precision and rapid calculation, but has evident flaws when simulating urban structures (such as buildings) with obvious 2D and 3D characteristics. The 2D model deals primarily with surface water movement, but it is unable to simulate the water flow dynamics in the drainage network as effectively as the 1D model (Dong et al. 2021). Therefore, coupled 1D sewer flow and 2D overland flow models for urban inundation modeling have been extensively investigated over the last few years (Leandro & Martins 2016; Noh et al. 2016; Dong et al. 2022).

The first 1D/2D coupling was developed by Hsu et al. (2000) who used Stormwater Management Model (SWMM) to analyze the flow dynamics in a drainage network and treated its overflow as inputs to a 2D overland flow model. The process of returning overflow to the network is not accounted for with this coupling technique, and therefore, the extent and depth of flooding in the downstream area are often overestimated. In an attempt to overcome this limitation, some academic scholars incorporated the consideration of the reflux process to define the 1D/2D bidirectional interactions (Hsu et al. 2002; Seyoum et al. 2012). These studies can be further categorized into two classes according to whether or not the surface convergence process is represented in 2D. The first category applies rainfall to the 1D hydrologic model, calculating sewer flow first and then interfacing the overload with the 2D surface flow model. The second category directly adds rainfall to the 2D model, where both surface catchment and network overflow are computed from the 2D model. Flood estimation may be erroneous if the rainfall intensity exceeds the intended capacity of the sewer inlets because the hydrologic model's catchment process assumes that runoff from sub-catchments is fully collected and discharged into the sewer network. In contrast, it is a fact that in urban areas, rainfall turns into runoff immediately after it drops to the ground surface and travels across the landscape for a certain distance before it reaches a drainage inlet or manhole. A simplified hydrologic model is unable to adequately clarify this phenomenon. The findings of Chang et al. (2015) further demonstrated that the 1D/2D coupling should apply rainfall directly to the 2D model to describe the precipitation catchment process rather than employing hydrologic routing.

The previously described bidirectional coupling strategy expresses two-way interactions in terms of water exchanges between the sewer system and the surface water system, where the interaction flow is commonly a simplified calculation based on the weir or orifice equations applied to the point-source head gradient. The coefficients within these manhole flow exchange equations bear stochastic properties that inject a degree of uncertainty and potential error into the estimated flow rates (Zhang et al. 2023). And achieving coherence in simulation results necessitates the use of exceedingly small time step, invariably sacrificing the temporal smoothness of the model. Furthermore, the development of the closely coupled model involves a heightened research threshold that not only demands a profound comprehension of the 1D and 2D system dynamics but also requires substantial programming expertise. Although there exists some commercial software with strong generalization and stable performance, such as InfoWorks ICM, MIKE, and TUFLOW (Syme 1992; Innovyze 2019; DHI 2020), they are not open-source and costly to utilize. The current unidirectional coupling methods, which are relatively convenient to implement, tend to only consider the transformation of the overflow from the 1D pipe network to the 2D model, making the results comparatively imprecise. By combining two unidirectional processes – feeding the 2D confluence into the 1D sewer network, and vice versa, channeling the overflow from the 1D network back into the 2D model – we can adopt a straightforward and manageable approach that addresses the conventional unidirectional process's inability to realistically depict the rainfall-runoff process. However, to our knowledge, no studies have yet attempted to compare the effectiveness of this integrated unidirectional method with that of the bidirectional coupling strategies.

To address the above issues, this paper investigated the coupling technique to merge the outcomes of two unidirectional processes and contrasted it with the bidirectional coupling method in terms of the result correctness, calculated stability, and implementation requirements. The 2D model was applied to mimic the rainfall catchment and surface diffusion process, whereas the 1D model solely analyzes the internal flow of the sewer network. The 2D model adopted Telemac-2D and the 1D model employed SWMM (Moulinec et al. 2011; McDonnell et al. 2021), both of which are open source and hence free and accessible. The details of the two coupling approaches are covered in the section that follows.

The paper is organized as follows: first, a detailed procedural description is presented for both the new coupling scheme proposed in this study and the existing bidirectional coupling scheme; then, two case discussions are introduced to evaluate the effectiveness of the proposed methodology, where the first one cites a large-scale laboratory experiment from the literature to validate the correctness and reasonableness of the proposed model (Naves et al. 2020), and the second one applies the method to simulate a real drainage system in a specific area of Nanjing, China, to validate the applicability of the methodology to large-scale drainage systems; and finally, the major conclusions are highlighted.

Among several numerical solvers, Telemac-2D is preferred in this paper since it simulates the moving boundaries of flows on a horizontal plane by the finite element method, which offers the strength of faster speed, lower memory requirements, better robustness to absorb the drastic variation in the transition between wet and dry regions, and especially the availability of the open source (Heniche et al. 2000). Hierarchical mesh resolution is possible with Telemac-2D finite element grids (Hervouet 2007). The average vertical water depth and velocity at each point in the resolution grid are obtained by solving the SWEs on triangular or quadrilateral elements.

The equations are discretized in finite difference form with spatial dispersion equal to the conduit length and are solved implicitly and iteratively using Picard's method over a given time step.

The majority of surface water that contributes to urban flooding is collected in the subterranean sewer system through manholes and artificial outlets (Dong et al. 2022). When the precipitation intensity surpasses the drainage capacity of the network, it causes rainwater overflow and the emergence of surface flooding.

Eventually, the 2D model results of the two processes can be aggregated with the water depths at identical spatial locations within the same time step to yield a comprehensive simulation of the entire rainfall-runoff overflow flooding procedure. The fact is that if variable time steps are adopted, it is impossible for two different 2D models to be synchronized in time, and the grid vertex locality cannot correspond accurately. Therefore, in this study, the 2D model results of the 2D-1D process are taken as the base results, with the 2D model results of the 1D-2D process as the additional results. The bathymetry values of the additional results are linearly interpolated on the time scale to correspond to the time steps of the base results, and on the spatial scale the KD-Tree algorithm (Bentley 1975) is employed to search for additional result nodes adjacent to the base result nodes at a certain distance (the mesh default size in this paper), and their average water depths are computed to be appended to the base result nodes' water depths. The treatment can be parallelized by partitioning the time series to improve computational efficiency.

The two straightforward and simple unidirectional coupling processes that make up the CTUC method's implementation eliminate the need to develop a customized 2D computational engine and coupling module and enable the direct adoption of pre-existing mature software (whether open source or not, e.g., LISFLOOD and HEC-RAS (Shustikova et al. 2019)) that focuses primarily on 2D simulation capabilities, offering a broader spectrum of application.

A precise characterization of the flow interaction between surface and sewer systems is important in the bidirectional coupling method. It synchronizes the simulation process of 1D and 2D models, and the water level difference between the node of the sewer network and the corresponding land surface is assessed at each time step to determine whether it is a nodal overflow or an inflow. Nodal inflow can be subdivided into free and submerged inflow according to the relation between the nodal water level and the adjacent surface elevation. Therefore, bidirectional coupling can be divided into the following three scenarios (Chen et al. 2007).

In the above equations, is the interacting discharge, is the weir discharge coefficient, is the orifice discharge coefficient, is the weir crest width, which is the manhole perimeter in this model, is the manhole area, g is the gravitational acceleration, is the hydraulic head at the manhole, is the surface elevation, is surface water level, and is the surface water depth, thus, ⁠.

A full-scale dual-drainage laboratory experimental case from the literature (Naves et al. 2020), which emulated the hydraulic and sediment transport process using a realistic rainfall simulator, was adopted in this work and the differences between the two coupling schemes concerning accuracy and implementation were analyzed and compared as the first step in verifying the dependability of the CTUC method. Then, the CTUC method was utilized to simulate the performance of a real large-scale drainage system in NJ, China during a typical local rainfall circumstance.

Rainfall simulators are capable of producing rainfall with intensities of 30, 50, and 80 mm/h. The simulation employed measured rainfall intensity maps to approximate the spatial distribution of actual rainfall. To compare the results of the CTUC and BC methods with the experimental data, the discharges from the two gully pots were selected as flow indicators, the measuring point S1 as surface water depth indicator, and the measuring points S2 and S3 as inside-pipe water depth indicators (Figure 4) for reference.

Discretizing the surface domain with an unstructured triangular mesh of default size 0.05 m, the CTUC method generated a total of 29,508 elements, while the BC method created more elements irrespective of the boundary geometry, with a total of 30,186. The two gully pots and the outlet channel are set up as open boundaries for free outflow. Following previous studies conducted in this laboratory facility (Naves et al. 2019a), the Manning roughness coefficients for the surface and pipes were set at 0.016 and 0.008 s.m^−1/3, respectively. For the flow exchange formula of the BC method, the weir and orifice coefficients for the inlets were set at 0.52 and 0.62 after model testing, respectively.

Since the BC idea intervenes in the variable states of the two subsystems at each time step, it is critical to restrict the magnitude of the interactive fluxes to ensure the stability of their results. In this paper, the tuning of this limitation is manifested mostly in the coefficient α of Equation (13), where α is too large for the results to oscillate significantly, and α is too small to describe the exact amount of the interaction flow. The determination of α value is associated with the overall runoff quantity size and time step size, which is an arithmetic and time-consuming task. In addition, the BC concept cannot be implemented in a mature calculating engine without open source, which necessitates the developers to possess a profound mastery of both hydraulic mechanisms and computer programming. The CTUC approach, on the other hand, integrates two continuous hydrological processes and eliminates the necessity for additional stabilization safeguards. Moreover, the CTUC approach can be realized by utilizing mature commercial software for 2D and 1D computing that is widely available in the market, which allows for wider applicability.

The unstructured triangular mesh with a default resolution of 6 m was generated, and the CTUC method sets each manhole as a square open boundary. A finer grid resolution of 1 m was implemented at the boundaries, and the square sides were 2 m long. The total number of grids for the area was 1,192,979. A precipitation hydrograph with a duration of 2 h is shown in Figure 6. The rainfall distribution is spatially homogeneous. The overall simulation duration was 4 h, and the time step was set to 5 s, satisfying the Courant–Friedrichs–Lewy (CFL) stability constraint (Courant et al. 1967).

We first explored the influence of different values in the CTUC method on the results (Table 2). Because of the relatively large size of the case area and the simplified drainage network model, which may not accurately reflect the positioning of the actual drainage inlets, the surface water volume rather than water depth was considered as the indicator to capture the overall inundation degree. It is clear from Table 2 that the lower the ⁠, the more pronounced restriction on the rainwater entering the network occurs, and more stagnant water remains in the surface domain, but due to the finite capacity of the network itself, the water quantity available to enter the network does not boost endlessly by increasing the ⁠. Further, the volumetric differences between the various settings were not evidently reflected over the full bathymetry map due to the large study area. Consequently, we disregarded the threshold limitation (i.e., ⁠) for the following discussions.

Among the 1D/2D coupled models that are extensively employed to simulate urban stormwater inundation developments, one generally acknowledged viewpoint is that the unidirectional coupled model from 1D network overflow to 2D surface diffusion cannot accurately characterize the rainfall surface confluence process, and hence is less dependable than the bidirectional coupled model. However, this paper proposed a new coupling scheme CTUC that combines the two unidirectional coupling processes of 2D surface confluence to 1D network and network overflow to surface flooding, and compared it with the contemporary bidirectional coupling idea BC in terms of the result accuracy and implementation requirements. This study adopted the 1D engine SWMM and the 2D engine TELEMAC, both of which are proven software that are open source and stable in operation and maintenance.

Through a laboratory experiment case designed to mimic urban street flooding, we first verified the reliability of the proposed coupling approach. The simulation results obtained by the CTUC and BC methods were compared with the experimental findings, separately, indicating that the CTUC and BC methods produced comparable results in terms of accuracy, both of which can satisfy the engineering requirements. The CTUC solution outperformed the BC marginally in the surface water depth, while it was slightly inferior to the BC strategy in the pipe water depth. The CTUC concept has a better calculation stability and applicability range than the BC method. Afterward, we applied the CTUC method to establish an actual urban drainage system model in a region of Nanjing, China, to synchronize the spatial and temporal variation of waterlogging under a 3-year rainfall condition, and the results were consistent with the physical topography and engineering characteristics of the urban area, which provided an efficient theoretical basis for the urban flood control design and renovation.

To summarize, the CTUC method proposed in this paper can generate relatively satisfactory and accurate simulation results for large-scale drainage systems, with improved computational stability and lower implementation threshold.

This work was financially supported by the National Natural Science Foundation of China (NSFC) [Grant 51978493] and [Grant 51778452].

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.