## ABSTRACT

Drainage modeling that accurately captures urban storm inundation serves as the foundation for flood warning and drainage scheduling. In this paper, we proposed a novel coupling ideology that, by integrating 2D-1D and 1D-2D unidirectional processes, overcomes the drawback of the conventional unidirectional coupling approach that fails to properly represent the rainfall surface catchment dynamics, and provides more coherent hydrological implications compared to the bidirectional coupling concept. This paper first referred to a laboratory experimental case from the literature, applied and analyzed the coupling scheme proposed in this paper and the bidirectional coupling scheme that has been widely studied in recent years, compared the two coupling solutions in terms of the resulting accuracy and applicability, and discussed their respective strengths and weaknesses to validate the reliability of the proposed method. The verified proposed coupling scheme was then applied to the modeling of a real drainage system in a region of Nanjing, China, and the results proved that the coupling mechanism proposed in this study is of practical application value.

## HIGHLIGHTS

A novel coupling approach by combining two unidirectional processes.

Comparison with bidirectional tight coupling scheme in terms of simulation accuracy and applicability.

The proposed coupling idea provides reliable outcomes with enhanced applicability.

## INTRODUCTION

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.

## METHODOLOGY

### Hydraulic models

#### 2D Surface model

*h*is the water depth,

*u*and

*v*are velocity components of the velocity vector in the

*x*and

*y*directions, respectively,

*Z*is the free surface elevation, , , and are the source terms representing the wind, Coriolis force, bottom friction, and a source or a sink of momentum within the domain, and is the momentum diffusion coefficient.

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.

#### 1D Sewer network model

*A*is the conduit cross-sectional area,

*t*is time,

*Q*is the discharge in the pipe,

*g*is the acceleration of gravity,

*H*is the water head in the pipe, and is the friction slope, which can be expressed in terms of the Manning equation used to model steady uniform flow:where

*n*is the Manning roughness coefficient, and

*R*is the hydraulic radius of the flow cross-section.

*H*is the node water depth,

*t*is time, is the net flow rate of the node, is the storage surface area of the node (if it is a storage node), and is the surface area contributed by the links connected to the node.

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.

### Linking methodology

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.

#### Combination of two unidirectional couplings (CTUC)

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.

#### Bidirectional coupling (BC)

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

##### Free weir linkage

##### Submerged weir linkage

##### Orifice linkage

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

## CASE STUDIES AND RESULTS

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.

### Model validation using a laboratory experiment

#### Case study description

^{2}full-scale concrete street section with a sewer system (Figure 4). The street surface consists of a sidewalk with a 15-cm elevation gap and a concrete roadway. The roadway has an approximate transversal slope of 2% up to the sidewalk and a 0.5% longitudinal slope up to an outflow channel. The generated rainfall runoff is discharged into the pipe system through two gully pots (0.3 m × 0.3 m) located along the curb and a lateral outflow channel. The flow is then conveyed to a common outlet of the pipe system. The detailed surface elevation data and conduit parameter data can be obtained in WASHTREET (Naves

*et al.*2019b).

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.

#### Results and discussions

Models . | GP1 . | GP2 . | |||||
---|---|---|---|---|---|---|---|

30 mm/h . | 50 mm/h . | 80 mm/h . | 30 mm/h . | 50 mm/h . | 80 mm/h . | ||

Experiment | 0.0472 | 0.0888 | 0.1414 | 0.1447 | 0.244 | 0.3402 | |

CTUC | Value | 0.05955 | 0.1027 | 0.1535 | 0.1517 | 0.2424 | 0.3465 |

Error (%) | 25.9534 | 15.6532 | 8.5573 | 4.8376 | 0.6557 | 1.8519 | |

BC | Value | 0.06052 | 0.1053 | 0.1584 | 0.1462 | 0.2429 | 0.3426 |

Error (%) | 28.2203 | 18.5811 | 12.0226 | 1.0366 | 0.4508 | 0.7055 |

Models . | GP1 . | GP2 . | |||||
---|---|---|---|---|---|---|---|

30 mm/h . | 50 mm/h . | 80 mm/h . | 30 mm/h . | 50 mm/h . | 80 mm/h . | ||

Experiment | 0.0472 | 0.0888 | 0.1414 | 0.1447 | 0.244 | 0.3402 | |

CTUC | Value | 0.05955 | 0.1027 | 0.1535 | 0.1517 | 0.2424 | 0.3465 |

Error (%) | 25.9534 | 15.6532 | 8.5573 | 4.8376 | 0.6557 | 1.8519 | |

BC | Value | 0.06052 | 0.1053 | 0.1584 | 0.1462 | 0.2429 | 0.3426 |

Error (%) | 28.2203 | 18.5811 | 12.0226 | 1.0366 | 0.4508 | 0.7055 |

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.

### Application of a real urban drainage system

#### Case study description

^{2}, and its topography and network layout are shown in Figure 6. As there was insufficient field data to calibrate the model parameters, the Manning roughness coefficients for the ground surface and the network pipes were consistently specified as 0.0125 and 0.0128 s.m

^{−1/3}, respectively. The system is a stormwater network that does not receive domestic sewage and consists of 358 manholes, 357 pipes, and 76 outfalls. The pipe diameters are generally large, ranging from 500 to 2,000 mm, as the network has been simplified.

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

#### Results and discussions

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.

Q_{max} (m^{3}/s)
. | Final volume (m^{3})
. | Rainfall (m^{3})
. | Entering the network (m^{3})
. | Outflow (m^{3})
. | Relative error . | Average water depth at peak ponding time (m) . |
---|---|---|---|---|---|---|

0.5 | 658475.40 | 1225310.89 | 332673.74 | 234225.07 | 9.62 × 10^{−5} | 0.05033 |

1 | 605790.80 | 1225310.89 | 385598.65 | 233984.84 | 1.05 × 10^{−4} | 0.04915 |

2 | 584575.60 | 1225310.89 | 406911.65 | 233887.05 | 1.08 × 10^{−4} | 0.04862 |

5 | 571293.10 | 1225310.89 | 420194.15 | 233887.05 | 1.10 × 10^{−4} | 0.04851 |

∞ | 571236.00 | 1225310.89 | 420194.15 | 233944.21 | 1.11 × 10^{−4} | 0.04850 |

Q_{max} (m^{3}/s)
. | Final volume (m^{3})
. | Rainfall (m^{3})
. | Entering the network (m^{3})
. | Outflow (m^{3})
. | Relative error . | Average water depth at peak ponding time (m) . |
---|---|---|---|---|---|---|

0.5 | 658475.40 | 1225310.89 | 332673.74 | 234225.07 | 9.62 × 10^{−5} | 0.05033 |

1 | 605790.80 | 1225310.89 | 385598.65 | 233984.84 | 1.05 × 10^{−4} | 0.04915 |

2 | 584575.60 | 1225310.89 | 406911.65 | 233887.05 | 1.08 × 10^{−4} | 0.04862 |

5 | 571293.10 | 1225310.89 | 420194.15 | 233887.05 | 1.10 × 10^{−4} | 0.04851 |

∞ | 571236.00 | 1225310.89 | 420194.15 | 233944.21 | 1.11 × 10^{−4} | 0.04850 |

*t*= 2,400 s. The model reached the maximum water accumulation at

*t*= 5,100 s after the simulation started, producing the largest water volume of 874,157.61 m

^{3}, thereafter the inundation range progressively decreased in the presence of the sewer network. The surface ponding volume reduced and the reduction rate got slower with time, so did the increasing rate of volume entering the network, indicating that the sewer system gradually reached the maximum drainage capacity. After

*t*= 12,000 s, the surface runoff around the street entrances was almost drained, so the sewer system contributed little to the water depth distribution over this period, with only partial waterlogging in the low-lying areas.

## CONCLUSIONS

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.

## FUNDING

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

## DATA AVAILABILITY STATEMENT

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

## CONFLICT OF INTEREST

The authors declare there is no conflict.