Enhancing urban pluvial flood modeling through graph reconstruction of incomplete sewer networks

This work presents an efficient graph-reconstruction-based approach for generating physical sewer models from incomplete information, addressing the challenge of representing the sewer drainage effect in urban pluvial flood simulations. The approach utilizes graph-based topological analysis and hydraulic design constraints to derive gravitational flow directions and nodal invert elevations in decentralized sewer networks with multiple outfalls. By incorporating linearized programming formulation to solve reconstruction problems, this approach can achieve high computational efficiency, making it suitable for application to city-scale sewer networks with thousands of nodes and links. Tested in Yinchuan, China, the approach integrates with a 1D/2D coupled hydrologic–hydrodynamic model and accurately reproduces maximum inundation depths (R²=0.95) when the complete network layout and regulated facilities are available. Simplifications, such as the adoption of road-based layouts and the omission of regulated facilities, can degrade simulation performance for extreme rainfall events compared to calibrated equifinal methods. However, design rainfall analysis demonstrates that the physical reconstruction approach can reliably outperform equifinal methods, achieving reduced variation and higher accuracy in simulating inundation areas. However, proper configuration of regulated facilities and network connectivity remains crucial, particularly for simulating local inundation during extreme rainfall. Thus, it is recommended that the proposed algorithm be integrated with targeted field investigations to further improve urban pluvial flood simulation performance in data-scarce regions.

1 Introduction

Global climate change and intensified human activities have significantly altered climatic conditions and surface environments, leading to an increase in the frequency and magnitude of urban pluvial floods, particularly where rainfall exceeds the design capacity of drainage systems (Rosenzweig et al., 2018). As key drainage infrastructure, sewer networks can effectively mitigate pluvial flood risks by transporting rainwater to downstream sites such as river channels and detention ponds (Wang et al., 2022). However, their performance can be inherently insufficient and gradually deteriorate due to outdated design and poor maintenance (Tran et al., 2024), potentially resulting in system overflow and worsening surface flood conditions during extreme rainfall events (Schmitt and Scheid, 2020). Thus, accurately modeling the drainage effects of sewer networks is essential for urban pluvial flood simulation and risk attribution (Luan et al., 2024).

However, detailed sewer network information may not always be available, and implementing it in hydraulic models requires laborious work with highly uncertain outcomes (Montalvo et al., 2024). As a result, various simplified approaches have been developed to approximate the effect of sewer drainage. These approaches typically introduce the term sewer-induced mass loss in rainfall–runoff simulations within specific regions, such as roads and inlets (Li et al., 2020; Xing et al., 2021, 2019). However, such approaches fail to explicitly account for the relative differences in hydraulic head between sewer networks and surface water, and thus are incapable of capturing critical processes such as overflow caused by system overload. To address these limitations, Montalvo et al. (2024) proposes a physics-based approach that can be used to design a representative virtual sewer network using publicly available datasets and local design standards. This approach enables the development of a physical sewer hydraulic model with bidirectional interactions between surface and sewer flows, thereby offering a more comprehensive representation of sewer drainage processes. The virtual network design (VND) consists mainly of two steps:

• Layout definition. This step determines the geometric connections and hydraulic directions between nodes in the sewer network. Typical algorithms include surface-elevation-driven delineation (Blumensaat et al., 2012; Duque et al., 2022), graph-theory-based modeling (Bakhshipour et al., 2019; Kim et al., 2021), and street-network-guided configuration (Chegini and Li, 2022; Reyes-Silva et al., 2023). Most existing work assumes the target system to be a centralized layout with a single outfall (Hesarkazzazi et al., 2022), facilitating the application of tree-based solutions such as the minimum spanning tree (Reyes-Silva et al., 2023) and the Steiner minimal tree (Machine Hsie and Huang, 2019). However, this simplification ignores the rich topological structures of real-world decentralized networks with potential cycles and forests (Bakhshipour et al., 2019). Furthermore, most approaches assume that gravitational flow direction aligns with surface gradients – an assumption that may not hold in flat urban areas or in the presence of inaccurate elevation data (Hesarkazzazi et al., 2022; Dunton and Gardoni, 2024).

• Hydraulic design. In this step, the sizes and slopes of pipes are determined, as well as the configuration of rainwater treatment facilities, with the objective of satisfying the design discharge requirements estimated by the Rational Method (Wang and Wang, 2018; Reyes-Silva et al., 2023). To achieve this, most approaches formulate the design task as an optimization problem and employ mathematical programming techniques such as linear programming (Swamee and Sharma, 2013), mixed-integer linear programming (Safavi and Geranmehr, 2017), nonlinear programming (Li and Matthew, 1990), and multi-objective programming (Bakhshipour et al., 2021). However, it requires significant computational efforts to obtain global optima in nonlinear or integer cases. Thus, approximate approaches such as piecewise linearization (Elimam et al., 1989), evolutionary algorithms (Wang et al., 2017; Bakhshipour et al., 2021), and topological-driven prediction-correction (Sitzenfrei et al., 2020; Chegini and Li, 2022) are proposed to alleviate the computational burden associated with optimization procedures. Nevertheless, a trade-off between approximation-based solution quality and computational efficiency still exists for large-scale cases (Yu et al., 2024).

Considering previous achievements in VND, this study further proposes a virtual network reconstruction (VNR) approach to improve the physical realism of sewer representations in urban pluvial flood simulations. While VND optimizes network design for economic costs and drainage performance, VNR estimates the most possible network properties constrained by partially known information. Although existing tools such as SWMManywhere (Dobson et al., 2025) offer practical solutions for high-quality VND using publicly available datasets, certain discrepancies between virtual and real networks exist due to inherent assumptions and simplifications, such as predefined road alignments for network topology, river proximity for outfall placement, and ideal hydraulic design using the Rational Method. These discrepancies can potentially alter hydrological connectivity and lead to the poor performance of flood simulations (Tran et al., 2024). Thus, to enhance the physical realism of virtual networks, additional efforts should be devoted to incorporating partially known information as constraints for network generation. This leads to a critical challenge: while basic network attributes, such as terminal nodes, spatial layout, and pipe sizes, can be identified and easily adjusted through existing VND toolsets, field surveys, and engineering drawings, key hydraulic parameters, such as flow direction and invert elevation, often remain unknown or are hard to estimate efficiently at city scale. To bridge this gap, we present our VNR approach, which integrates partial network information with surface elevation data to derive the most possible hydraulic parameters, thereby reducing the simulation uncertainty associated with sewer drainage effects in urban flood modeling.

In the remainder of this paper, we first describe the proposed workflow for VNR (Sect. 2) along with the corresponding 1D/2D coupled hydrologic-hydrodynamic model for urban pluvial flood simulation (Sect. 3). We then conduct a comparative analysis of flood simulation performance given varying levels of information completeness during reconstruction, taking into account factors such as regulated facility removal and road-based layout simplification (Sect. 4).

2 Graph-based sewer network reconstruction

2.1 Graph representation of sewer networks

The necessary inputs for sewer network reconstruction include the sewer network layout with specified outfalls and the corresponding digital elevation model (DEM). For convenience of topological analysis, we represent the sewer network as an undirected graph $G=(V,E)$, where V={v_i} and $E=\mathit{\left\{}{e}_{i,j}\mathit{\right\}}$ represent the sets of nodes and links, respectively (Shi et al., 2023). Then we can specify the link directions as related gravitational flow directions and convert the original undirected graph to the corresponding directed form $G_{\mathrm{D}}=(V,E_{\mathrm{D}})$, where $E_{\mathrm{D}}=\mathit{\left\{}{\overrightarrow{e}}_{i,j}\mathit{\right\}}$ represents the set of directed links evaluated as follows:

$$\begin{array}{}\text{(1)}& {\overrightarrow{e}}_{i,j}=\left\{\begin{array}{ll}\mathrm{1},& v_j\in {\mathcal{N}}_{G_{\mathrm{D}}}\left(v_i\right)\\ \mathrm{0},& v_j\notin {\mathcal{N}}_{G_{\mathrm{D}}}\left(v_i\right)\end{array}\right.,\forall v_j\in {\mathcal{N}}_G\left(v_i\right)\end{array}$$

where ${\mathcal{N}}_{G_{\mathrm{D}}}\left(v_i\right)$ and 𝒩_G(v_i) represent the neighboring nodes of node v_j in G_D and G.

For one node v_i, its out-degree d_out,i and in-degree d_in,i can be defined as the number of directed links pointed from/to it:

$$\begin{array}{}\text{(2)}& d_{\mathrm{out},i}=\sum _{v_j\in {\mathcal{N}}_G\left(v_i\right)}{\overrightarrow{e}}_{i,j},d_{\mathrm{in},i}=\sum _{v_j\in {\mathcal{N}}_G\left(v_i\right)}{\overrightarrow{e}}_{j,i},\end{array}$$

In addition to these basic attributes, some useful structures can be derived from the graph representation. Considering the potential cycles in the geometric layout of G, we identify the set of simple cycles 𝒞={C_k} defined as the set of closed paths where no node appears twice (Johnson, 1975), such as C₁ and C₂ in Fig. 1. For a cycle C_k∈𝒞, we group its internal nodes into the corresponding set, denoted by CV_k, where the ith node is indexed by ${\mathrm{CV}}_k^{\left(i\right)}$. Then its out-degree dC_out,k and in-degree dC_in,k can be defined as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}d{C}_{\mathrm{out},k}& =\sum _{v_i\in {\mathrm{CV}}_k}\sum _{v_j\in {\mathcal{N}}_G\left(v_i\right)\backslash {\mathrm{CV}}_k}{\overrightarrow{e}}_{i,j},d{C}_{\mathrm{in},k}\\ & =\sum _{v_i\in {\mathrm{CV}}_k}\sum _{v_j\in {\mathcal{N}}_G\left(v_i\right)\backslash {\mathrm{CV}}_k}{\overrightarrow{e}}_{j,i}\end{array}\end{array}$$

where 𝒩_G(v_i)\CV_k represents the difference set between 𝒩_G(v_i) and CV_k, i.e., the set of nodes that are neighboring nodes of v_i but not in C_k.

Figure 1 Cycles, islands, outfalls, and sources in the graph-based sewer network representation.

According to Eq. (3), we define a cycle C_k as an “island” if the sum of its in-degree and out-degree is equal to 1, such as C₂ in Fig. 1. We further denote the set of nodes located in the islands of G as ISD(G):

$$\begin{array}{}\text{(4)}& \mathrm{ISD}\left(G\right)=\bigcup _{C_k}{\mathrm{CV}}_k\mathrm{s}.\mathrm{t}.d{C}_{\mathrm{out},k}+d{C}_{\mathrm{in},k}=\mathrm{1}\end{array}$$

In addition to general properties, we also consider some sewer-related hydraulic characteristics for the convenience of the following analysis. For a node v_i∈V, we include it in the set of outfalls (denoted as 𝒪) if its out-degree is 0 and in-degree is 1, such as v₁₃ in Fig. 1, or in the set of sources (denoted as 𝒮) if its in-degree is 0 and out-degree is 1, such as v₁ in Fig. 1:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathcal{O}=\left\{v_i\mid v_i\in V,d_{\mathrm{out},i}=\mathrm{0},d_{\mathrm{in},i}=\mathrm{1}\right\},\\ & \mathcal{S}=\left\{v_i\mid v_i\in V,d_{\mathrm{out},i}=\mathrm{1},d_{\mathrm{in},i}=\mathrm{0}\right\}\end{array}\end{array}$$

During practical usage, the set of sources can be derived as the difference set between the set of predefined outfalls and the set of nodes whose in-degree and out-degree sum to 1.

To allow for hydraulic simulation, we attach the invert elevation attribute to v_i (denoted as z_i), which can be initialized as follows:

$$\begin{array}{}\text{(6)}& z_i^{\left(\mathrm{0}\right)}=z_{g,i}-D_{\min }\end{array}$$

where $z_i^{\left(\mathrm{0}\right)}$ is the initial value of z_i; z_g,i is the surface elevation at the location of v_i; D_min is the minimum allowable cover depth and is set as D_min=0.6 m by default (MOHURD, 2021).

By comparing the invert elevations of neighboring nodes, link directions can be initialized as follows:

$$\begin{array}{}\text{(7)}& {\overrightarrow{e}}_{i,j}^{\left(\mathrm{0}\right)}=\left\{\begin{array}{ll}\mathrm{1},& z_i^{\left(\mathrm{0}\right)}\ge z_j^{\left(\mathrm{0}\right)}\\ \mathrm{0},& z_i^{\left(\mathrm{0}\right)}<z_j^{\left(\mathrm{0}\right)}\end{array}\right.\end{array}$$

2.2 Gravitational flow direction derivation

Since flow direction initialization using Eq. (7) does not consider the necessary topological constraints for a feasible layout of the sewer network (Fig. 2), it may produce an unfeasible or unreasonable layout that requires further correction (Reyes-Silva et al., 2023).

Figure 2 Topological constraints for a feasible sewer network.

To ensure flow directions are preserved as much as possible during correction, we formulate the following 0–1 programming problems with linear constraints summarized in Table 1:

$$\begin{array}{}\text{(8)}& \underset{e_{ij}}{min}\sum _{{\overrightarrow{e}}_{i,j}}\left|{\overrightarrow{e}}_{i,j}-{\overrightarrow{e}}_{i,j}^{\left(\mathrm{0}\right)}\right|\end{array}$$

Table 1 Optimization constraints for gravitational flow direction derivation. $\left|C_k\right|$ is the number of nodes in C_k; ${\mathcal{E}}_{\mathrm{D}}^{*}$ is the set of links with known directions (denoted by ${\overrightarrow{e}}_{i,j}^{*}$); ${\overrightarrow{e}}_{{\mathrm{CV}}_k^{\left(\right|{\mathcal{C}}_k\left|\right)},{\mathrm{CV}}_k^{\left(\mathrm{1}\right)}}+\sum _{i=\mathrm{1}}^{\left|{\mathcal{C}}_k\right|-\mathrm{1}}{\overrightarrow{e}}_{{\mathrm{CV}}_k^{\left(i\right)},{\mathrm{CV}}_k^{(i+\mathrm{1})}}$ is the sum of consecutively connected direction values along C_k.

Prior to optimization, pipes are discretized into 50 m segments following common inlet spacing (MOHURD, 2021). This ensures uniform link lengths for objective evaluation and provides computational elements with inlets at segment endpoints, enabling surface–sewer exchange in the hydraulic simulation.

2.3 Nodal invert elevation derivation

Correcting gravitational flow directions enable the derivation of nodal invert elevations required for downstream hydraulic simulation. Based on invert elevation initialized by Eq. (6), further adjustment is needed to keep the slopes of links within the range of minimum and maximum allowable slopes (Safavi and Geranmehr, 2017).

To achieve the least deviation from the initial results considering sewer construction costs proportional to excavation depth (Swamee and Sharma, 2013; Machine Hsie and Huang, 2019), we introduce the following nonlinear optimization problem to estimate the invert elevation of nodes:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& \underset{z_i}{min}\sum _{z_i}\left|z_i-z_i^{\left(\mathrm{0}\right)}\right|\\ & \text{s.t.}\left\{\begin{array}{ll}{i}_{i,j}\le -i_{\min ,i,j},& \forall {\overrightarrow{e}}_{i,j}=\mathrm{1}\\ i_{i,j}\ge -i_{\max ,i,j},& \forall {\overrightarrow{e}}_{i,j}=\mathrm{1}\end{array}\right.\end{array}\end{array}$$

where i_i,j is the slope of ${\overrightarrow{e}}_{i,j}$; $i_{\min ,i,j},i_{\max ,i,j}$ are the minimum and maximum allowable slopes of ${\overrightarrow{e}}_{i,j}$, summarized from the Standard for Design of Outdoor Wastewater Engineering (MOHURD, 2021), considering both explicitly stated slope limits and those derived from velocity constraints (Swamee and Sharma, 2013).

For simplicity, but without loss of generality, i_i,j is estimated using the invert elevation of two consecutive nodes, ensuring the linearity of constraints for the above optimization problem:

$$\begin{array}{}\text{(10)}& i_{i,j}=-{\displaystyle \frac{\mathrm{1}}{d_{i,j}}}{z}_i+{\displaystyle \frac{\mathrm{1}}{d_{i,j}}}{z}_j\end{array}$$

It should also be noted that Eq. (10) can be extended to multiple consecutive nodes. The corresponding slope estimator remains a linear combination of nodal invert elevation (see Sect. S1 in the Supplement), which maintains the linearity of constraints in the original optimization problem.

Table 2 Slope constraints with respect to different sizes according to the Standard for Design of Outdoor Wastewater Engineering (MOHURD, 2021).

2.4 Computational efficiency considerations

The objective functions of Eqs. (8) and (9) are expressed as the sum of absolute deviations and can be further reduced to an equivalent linear form following (Wagner, 1959; Giangrande et al., 2013), thus ensuring their solution efficiency for large-scale applications.

Without loss of generality, for minimizing the summation of absolute values:

$$\begin{array}{}\text{(11)}& L=\sum _{i=\mathrm{1}}^N\left|k_i{x}_i-b_i\right|\end{array}$$

where k_i,b_i are fixed numbers with respect to N decision variables x_i, we can introduce auxiliary decision variables $z_i=\left|k_i{x}_i-b_i\right|$ constrained by:

$$\begin{array}{}\text{(12)}& \left\{\begin{array}{l}{z}_i\ge k_i{x}_i-b_i\\ z_i\ge -k_i{x}_i+b_i\\ z_i\ge \mathrm{0}\end{array}\right.\end{array}$$

Thus, the original unconstrained nonlinear problem can be reformulated into an equivalent linear one with additional constraints:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \underset{x_i,z_i}{min}\sum _{i=\mathrm{1}}^N{z}_i\\ & \text{s.t.}\left\{\begin{array}{l}{z}_{\mathrm{1}}\ge k_i{x}_{\mathrm{1}}-b_{\mathrm{0}},z_{\mathrm{1}}\ge -k_i{x}_{\mathrm{1}}+b_{\mathrm{0}},z_{\mathrm{1}}\ge \mathrm{0}\\ \mathrm{\dots },\mathrm{\dots }\\ z_N\ge k_i{x}_N-b_N,z_N\ge -k_i{x}_N+b_N,z_N\ge \mathrm{0}\end{array}\right.\end{array}\end{array}$$

Furthermore, if additional constraints are present with Eq. (11), they can be added directly to the reformulated problem (Eq. 13) without loss of equivalence. This reformulation strategy has been used in previous research (Giangrande et al., 2013) and thus can be easily extended to our work on efficient sewer network reconstruction (Eqs. 8 and 9). In terms of solving the reformulated linear programming problems, we select the open-source Computational Infrastructure for Operations Research (COIN-OR) branch-and-cut solver (CBC) (Forrest and Lougee-Heimer, 2005) through its Python interface in PuLP (Mitchell et al., 2011).

3 AUTOSHED: A 1D/2D coupled hydrologic-hydrodynamic model

Using the sewer network reconstructed from incomplete information, we further adopt a 1D/2D coupled hydrologic-hydrodynamic model, AUTOSHED, to simulate pluvial flood process in an urbanized area. AUTOSHED delineates the whole simulation domain into multiple triangular-shaped units (TSUs). For each unit, AUTOSHED uses the rain-on-grid approach, also known as the fully hydrodynamic approach (Hall, 2015; Chen and Huang, 2024; Perrini et al., 2024), and performs sequential hydrologic-hydrodynamic computations for surface runoff generation and routing, respectively.

3.1 Surface runoff generation and routing

The surface runoff generation scheme is mainly built upon the core module of the Surface Urban Energy and Water Balance Scheme (SUEWS) (Järvi et al., 2011) and Tsinghua Integrated Hydrological Modeling System (THIHMS) (Ni et al., 2008), except for overland and channel flow routing. Each TSU is divided into impervious and pervious parts according to the empirical impervious ratio α_imp derived from the unit land use statistic by area-weighted averaging (Cao et al., 2020; Krebs et al., 2014). And the total surface runoff ΔS_sr is calculated as the area-weighted sum of runoff from the impervious and pervious parts, considering typical hydrological processes such as vegetation interception, building interception, and soil infiltration (Fig. 3). Details of methodologies for rainfall–runoff simulation can be found in the work of Ni et al. (2008) and Järvi et al. (2011).

Figure 3 Flux calculation within a typical triangular-shaped unit and across the surface–sewer system.

After obtaining the total surface runoff, a well-balanced formulation of 2D shallow water equations (SWEs) (Song et al., 2011) is adopted to simulate the subsequent surface routing process:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\displaystyle \frac{\partial }{\partial t}}\left[\begin{array}{c}h\\ hu\\ hv\end{array}\right]+{\displaystyle \frac{\partial }{\partial x}}\left[\begin{array}{c}hu\\ h{u}^{\mathrm{2}}+\mathrm{g}(h^{\mathrm{2}}-z_b^{\mathrm{2}})/\mathrm{2}\\ huv\end{array}\right]\\ & +{\displaystyle \frac{\partial }{\partial y}}\left[\begin{array}{c}hu\\ huv\\ h{v}^{\mathrm{2}}+\mathrm{g}(h^{\mathrm{2}}-z_b^{\mathrm{2}})/\mathrm{2}\end{array}\right]\\ & =\left[\begin{array}{c}\mathrm{\Delta }{S}_{\mathrm{sr}}-q_{\mathrm{sd}}\\ -\mathrm{g}(h+z_{\mathrm{b}})\frac{\partial z_b}{\partial x}-\mathrm{g}h\frac{n^{\mathrm{2}}u\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}}{h^{\mathrm{4}/\mathrm{3}}}\\ -\mathrm{g}(h+z_{\mathrm{b}})\frac{\partial z_b}{\partial y}-\mathrm{g}h\frac{n^{\mathrm{2}}v\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}}{h^{\mathrm{4}/\mathrm{3}}}\end{array}\right]\end{array}\end{array}$$

where h is the water depth; u,v are the depth-averaged velocity components in the x- and y-directions; z_b is the bed elevation; g is the gravitational acceleration; q_sd is the mass loss term caused by sewer drainage; n is the empirical Manning coefficient initialized as n₀ by land-use type and further increased according to the intersected area ratio of TSUs with buildings (denoted as b_f) by the building roughness method (Schubert and Sanders, 2012), where the maximum roughness is set as 0.5 s m${}^{-\mathrm{1}/\mathrm{3}}$. Details on the methodologies used for solving 2D SWEs can be found in the work of Song et al. (2011).

3.2 Sewer flow simulation

In TSUs with rainwater inlets, vertical linkages enable bidirectional coupling between surface and sewer flows. A modified form of the 1D Saint-Venant equations is introduced to simulate the sewer flow process using a node-link approach (Rossman and Huber, 2017):

$$\begin{array}{}\text{(15)}& \left\{\begin{array}{l}\frac{\partial H}{\partial t}=\frac{Q_{\text{1D-2D}}+\sum Q_{\mathrm{nl}}}{A_{\mathrm{as}}}\\ \frac{\partial Q}{\partial t}=\mathrm{2}U\frac{\partial A}{\partial t}+U^{\mathrm{2}}\frac{\partial A}{\partial x}-gA\frac{\partial H}{\partial x}-gA{S}_{\mathrm{f}}\end{array}\right.\end{array}$$

where H is the total nodal hydraulic head of water; Q is the flow rate; U is the flow velocity; A is the cross-sectional area; A_as is the node assembly surface area consisting of the node's storage surface area and half the surface area of connected links (Rossman, 2015); S_f is the friction slope parameterized by the Manning formula; Q_nl is the flow contributed by neighboring links in the sewer system; Q_1D-2D is the mass source term exchanged between the surface flow and the sewer flow. Further details regarding the methodologies used for solving 1D Saint-Venant equations under regulation can be found in Rossman (2015). The corresponding surface–sewer flow coupling scheme is given in the Sect. S2.

4 Performance analysis: from graph reconstruction to urban flood simulation

4.1 Study area and data

The study area is the main city zone of Yinchuan (CYC), Ningxia Hui autonomous region, situated in the northwestern part of China, covering an area of 137.15 km². Despite being located in an arid region with an average annual rainfall of 189 mm and annual evaporation of 1825 mm, CYC has experienced several flood disasters due to increasing extreme rainfall events in recent years (Lu et al., 2024). According to local documents, most of the drainage system in CYC is a combined sewer system. Its drainage capacity is largely constrained by downstream wastewater treatment plants (WTPs) due to restrictive environmental regulations, which can cause overflow at rainwater inlets and further increase surface inundation risks.

To establish the urban pluvial flood model for scenario analysis, we collected the following multi-source geographical datasets :

• 5 m digital elevation model (DEM, Fig. 4a) generated from the stereo images of ZiYuan-3 satellite. Following common DEM preprocessing practices for urban flood simulation, we delete the inaccurate building roof elevation from the original product (Chen and Huang, 2024) and lower the surface elevation of grids containing roads by 0.2 m (Liu et al., 2022) to improve the representation of hydrological connectivity in DEM.

• 10 m land-use map (Fig. 4b) from the European Space Agency (ESA) WorldCover 2021 product (Zanaga et al., 2022), which covers tree, grassland, cropland, built-up and bare-land area in the study area.

• Building footprint (Fig. 4a) and road network layers (Fig. 4b) downloaded via the application programming interface of the Baidu map. Considering the potential underestimation of the impervious area identified from remote sensing datasets (Weng, 2012), we further reclassified pixels intersecting with buildings and roads as the built-up land-use type.

• Sewer network layout and the corresponding attribute of pipeline size (Fig. 5a) were digitized from investigation and design documents authorized by the local government. It should be noted that when detailed information was unavailable, previous research would generate a virtual sewer network from the street network (Montalvo et al., 2024). To investigate the potential effects of such simplification, we further created a road-based sewer network by identifying pipes that intersected with the major roadway zone and then performed the necessary reconnections at pipes broken intermediately (Fig. 5b). Here, we define the term “major road” (Fig. 4b) as including primary, secondary, tertiary and trunk roads. Roadway zones are generated by buffering road centerlines with widths determined according to The Design of Urban Road Engineering (CJJ 37-2012).

Figure 4 Geographical datasets in the study area. (a) Digital elevation model. (b) Land-use map and road networks. Basemaps are derived from ESRI World Imagery (Credit: Esri, TomTom, Garmin, FAO, NOAA, USGS © OpenStreetMap contributors, and the GIS User Community).

Figure 5 Sewer network in the study area. (a) Complete sewer network digitized from the local document. (b) Simplified sewer network filtered by the major roadway zone. Basemaps are derived from ESRI World Imagery (Credit: Esri, TomTom, Garmin, FAO, NOAA, USGS, © OpenStreetMap contributors, and the GIS User Community).

The design rainfall formula in CYC is calculated based on local rainfall statistics:

$$\begin{array}{}\text{(16)}& q={\displaystyle \frac{\mathrm{551.4}(\mathrm{1}+\mathrm{0.584}\lg P)}{(t+\mathrm{11}{)}^{\mathrm{0.669}}}}\end{array}$$

where q is the rainfall intensity, L (s ha)⁻¹; P is the return period; t is the rainfall duration in min.

Seven designed rainfalls with return years of 1a, 5a, 10a, 20a, 50a, 100a, and 200a are selected for the scenario analysis, where each rainfall event lasts for 2 h (see Fig. S1 in the Supplement), followed by an additional 2 h period without any precipitation.

4.2 Baseline configurations

Five baseline configurations are selected for comparing the performance of the proposed approach (denoted as the FSR approach) and other approximation methods:

• Full-Sewer-reconstruction-No-controls (FSN), which uses the same reconstruction approach but ignores the existence of regulated facilities such as pump stations and wastewater treatment plants. To be specific, regulated facilities located at outflow and intermediate points are configured as free outfalls and uncontrolled junctions, respectively.

• Road-based-Sewer-reconstruction-No-controls (RSN) which uses the same reconstruction approach but adopts the layout of sewers filtered by the major roadway zone (Fig. 5b) and ignores the existence of regulated facilities.

• Road-based-Sewer-reconstruction-with-Controls (RSC) which uses the reconstructed sewer network from RSN and further incorporates regulated facilities at corresponding locations.

• Road-Drainage Approximation (RDA) (Xing et al., 2021; Li et al., 2023), which removes water mass at a uniform rate i_d within roadway zones buffered from road centerlines, with widths evaluated using The Design of Urban Road Engineering (CJJ 37-2012).

• No-Drainage Approximation (NDA) which ignores the drainage effects of the sewer network.

According to whether the sewer hydraulic module is reconstructed for simulation with AUTOSHED, we can further classify the above approximation methods into two categories: physical approaches (including FSR, FSN, RSN, and RSC), and equifinal ones (including RDA and NDA). This means that similar flood simulation results can be achieved by different methods with appropriate parameters of drainage effects (Dobson et al., 2025).

To quantify the simulation performance of different approximation methods under various design scenarios, we calculate the probability of detection (POD), false alarm rate (FAR), critical success index (CSI), and bias (BIAS) of inundation area simulations (McGrath et al., 2018), using results from the FSR approach as the reference (see Sect. S3).

4.3 Sewer network reconstruction performance

According to the layout of the digitized sewer network (Fig. 5), gravitational flow directions are initialized by surface elevation (Eq. 7) and then estimated by solving 0–1 programming problems (Eq. 8). For simplicity, Fig. 6 summarizes the results for the initialized and reconstructed directions of the complete sewer network with 11 outfall points consisting of wastewater treatment plants and pumping stations. The directions of 359 among 908 pipes are reversed to satisfy topological constraints. The majority of the modified pipes are located in the southern part of the study area and exhibit a pattern of cyclic structures over relatively flat areas, posing challenges for the direct identification of dominant flow directions and thus necessitating the incorporation of additional topological constraints (Table 1) for rational reconstruction.

Figure 6 Comparison of the initial and reconstructed gravitational flow directions of the complete sewer network in the study area. (a) Flow directions initialized by surface elevation. (b) Flow directions reconstructed from topological constraints. Basemaps are derived from ESRI World Imagery (Credit: Esri, TomTom, Garmin, FAO, NOAA, USGS, © OpenStreetMap contributors, and the GIS User Community).

According to the reconstructed flow directions, nodal invert elevation is further derived from values initialized by constant cover depth (Eq. 6). The slope distributions before/after adjustment are summarized in Fig. 7. The initial slopes of both the complete and simplified sewer network exhibit a relatively symmetric distribution centered on 0, indicating adverse slopes in nearly half the links due to the reversal of nearly 40 % of the pipes in the previous step. After reconstruction using Eq. (9), slope distributions shift to strictly positive ones, where the majority of slopes are around 0.001. Additionally, the simplified sewer network shares a similar pattern of slope distribution as the complete one, except for minor increases around 0 because its layout is slightly modified due to the aforementioned reconnection of broken pipes in Sect. 4.1. Thus, it can be concluded that the proposed reconstruction approach is robust to empirical simplification of the sewer network layout following road networks and can be stably extended to applications using street-based network generation.

Figure 7 Comparison of the initial and reconstructed sewer slope distributions. (a) Complete sewer network used in the FSR and FSN approaches. (b) Road-based sewer network used in the RSN and RSC approach.

Owing to the linearized programming formulation used in the process of reconstructing gravitational flow directions and nodal invert elevation (Sect. 2.4), 1D sewer hydraulic models can be built from prepared geographical datasets in 33 s on a local desktop using a single core of Intel(R) Core(TM) i7-14700KF CPU and 400 MB RAM occupancy. For reconstructed sewer networks initialized with a uniform spacing of 50 m (cf. Sect. 2.2), the FSR/FSN and RSN/RSC models have 8242 links and 8131 inlets, and 5268 links and 5188 inlets, respectively. Furthermore, the entire study area is discretized into 7222150 TSUs with a median area of 18.2 m² for surface runoff generation and routing computation.

4.4 Model calibration and evaluation

To compare pluvial flood simulation performance across different approaches, we first calibrate the model parameters using a severe storm that hit CYC on 11 July 2022 (denoted as the July 2022 storm). This storm caused an average total rainfall of 57.8 mm and a peak intensity of 44.5 mm h⁻¹, resulting in 31 inundated points, as reported by the local government, of which 18 were recorded to have the corresponding maximum waterlogging depth. Considering its representativeness as a disastrous torrential event and data availability, we use the recorded maximum waterlogging depth to calibrate the parameters in the aforementioned model and the related approximation approaches for the sewer drainage effect.

For AUTOSHED, the main physical parameters include the empirical impervious ratio α_imp and the Manning coefficient n estimated from land-use types. Considering the dominance of built-up area, and to avoid possible parameter equifinality between physical processes such as soil infiltration and sewer drainage, we tested the model performance using the recommended empirical values of the corresponding parameters from the Technical Specification for Construction and Application of Mathematical Model of Urban Flooding Prevention and Control System (DB11/T 2074-2022) and finalized the values of parameters without exhaustive search (see Table S1 in the Supplement for details). Figure 8 shows the inundation map simulated by AUTOSHED using the sewer hydraulic module reconstructed from the FSR approach for the July 2022 storm.

Figure 8 Maximum inundation depth map simulated by the FSR approach for the July 2022 storm. Basemaps are derived from ESRI World Imagery (Credit: Esri, TomTom, Garmin, FAO, NOAA, USGS, © OpenStreetMap contributors, and the GIS User Community).

The coefficient of determination (R²) is used to quantify the overall performance of the models using different parameters. According to Fig. 9a, the FSR approach achieves an R² of 0.95, indicating strong agreement between the simulated and observed values and establishing the physical validity of the model for further comparison. In addition, we further examine the performance of the FSR approach with varying inlet spacings (denoted as ILS) from 30 to 100 m. According to Fig. 9b, the FSR approach demonstrates relatively robust performance with R² ranging from 0.85 to 0.95 where ILS of 50 m achieves the highest R². It is worth noting that increasing ILS leads to both underestimation and overestimation, such as P5 and P2 in Fig. 9b, highlighting the dual role of rainwater inlets in either alleviating surface inundation through drainage or contributing to it when overloaded.

Figure 9 Simulation performance of the maximum inundation depth at inundated points using different approximation methods and parameters. The index of inundated points follow Fig. 8.

Based on calibrated parameters of the physical model, we further use the grid search strategy to optimize the most possible equivalent drainage rate i_d for use in the RDA method when detailed sewer network information is unavailable. Specifically, we set the possible range of i_d according to the local design drainage standard and perform simulations with i_d = 2.5, 5.0, 7.5, 10.0, 12.5, 15.0 mm h⁻¹. According to Fig. 9c, i_d = 2.5 mm h⁻¹ achieves the highest R² and is thus selected as the final value for the drainage rate for further scenario analysis.

4.5 Extreme rainfall analysis: both local inlets and global regulation matter

After calibrating the AUTOSHED model integrated with the reconstructed sewer hydraulic modules, we now delve deeper into how the proposed graph-reconstruction-based approach performs compared to simplified alternatives across different rainfall scenarios and varying levels of available information. The aforementioned July 2022 storm event, with an average total rainfall exceeding the 500-year return period threshold (40.8 mm), highlights its extremity in an arid region. According to Fig. 9, when information on regulated facilities and detailed layouts is unavailable, simplified physical approaches, such as FSN, RSC and RSN, exhibit relatively lower accuracy than calibrated equifinal ones, which approximate limited drainage capacities during extreme rainfall events using a small equivalent drainage rate.

As physical approaches incorporate drainage effects at rainwater inlets, we further investigate the relationship between simulation errors in maximum inundation depth and the distance to the nearest inlet, denoted as D_f for the complete sewer network and D_s for the simplified one (see Sect. 4.1). According to evaluation results (see Fig. S2), simulation error variability decreases as D_f and D_s increase, with the maximum absolute value of deviations remaining less than 5 cm for distances exceeding 30 m, where local microtopography predominantly influences inundation. Furthermore, NDA can underestimate the maximum inundation depth at points such as P8 and P18, suggesting the occurrence of potential overflows around nearby rainwater inlets due to insufficient drainage capacity during extreme rainfall.

To assess differences in simulated sewer drainage performance, we further summarize accumulated discharges at rainwater inlets and conduct two steps of comparisons based on the available information.

• Effects of neglecting regulation: when regulated facilities are excluded, FSN can significantly underestimate the maximum inundation depth at points such as P2 compared to FSR. This mainly results from overestimated inflows ($\mathrm{\Delta }\sum Q_{\text{1D-2D}}>\mathrm{0}$) at surrounding inlets, as seen at P2 in Fig. 10a and e. Similar overestimation can also be observed when comparing RSN with RSC, but with greatly decreased magnitude, since road-based simplification reduces the rainwater collected into the sewer system and thus alleviates the drainage burden. However, when regulated facilities are deleted, FSN also overestimates the maximum inundation depth at certain points such as P1 and P3, which may seem counterintuitive. Analysis of P1's nearby inlets reveals lower inflows using FSN compared to FSR, particularly at IL1 (Fig. 10b and f), causing model overestimation of surrounding maximum inundation depth. After tracking accumulated discharges at downstream inlets with respect to IL1, we further locate the drainage bottleneck at one pump station, PM1. As illustrated in Fig. 10d, its removal can significantly increase inflows at upstream inlets, further increasing water heads at adjacent sewer nodes and consequently reducing upstream inlet drainage capacities due to backwater effects.

• Impact of network simplification: when approximating the sewer layout using major roads, RSN generally overestimates the maximum water depth at inundated points within 30 m of rainwater inlets compared to FSN, driven by intensified overflows from excessive surface water flowing through a single drainage pathway, as seen at P1 in Fig. 11a and c. However, some inundated points, like P13 and P18, where surrounding inlets are removed, may also underestimate maximum water depths. This occurs because these inlets initially contribute to surface overflows in FSN-based simulations, as evidenced by the underestimation of maximum water depth when NDA is adopted. In particular, we can also observe the overestimation of maximum water depths at certain flood points without removing local inlets during simplification, such as in P2 and P4 with ΔD=0. Examining the simulated total discharges at rainwater inlets near P2 (Fig. 11b and d) reveals that increased local inflows are outweighed by amplified upstream overflows, thus exacerbating downstream inundation conditions.

Figure 10 Comparison of accumulated discharges (∑Q_1D-2D) at rainwater inlets during the July 2022 storm, where positive values indicate inflows. (a–b) FSN simulated discharges. (c–f) Simulated discharge differences (∑ΔQ_1D-2D) between FSN and FSR, defined as FSN minus FSR. Basemaps are provided by Tianditu (https://map.tianditu.gov.cn/, last access: 14 October 2025).

Figure 11 Comparison of accumulated discharges at rainwater inlets during the July 2022 storm. (a–b) RSN simulated discharges. (c–d) Simulated discharge differences between RSN and FSN defined by RSN minus FSN. Basemaps are provided by Tianditu (https://map.tianditu.gov.cn/).

In summary, while simplified physical approaches based on major roads can yield a similar pipe slope distribution (Sect. 4.3), achieving an accurate simulation of maximum inundation depths requires the proper configuration of regulated facilities and drainage pathways. Otherwise, equifinal methods, which offer flexible parameter calibration, can offer better performance, especially when observational data are available.

4.6 Design rainfall analysis: physical approaches outperform but require more complete information

To assess whether equifinal approaches exhibit comparable or even better performance than physical counterparts, we further compare inundated areas simulated by different approaches under design rainfalls where the FSR-based results serve as the reference.

Figure 12 Comparison of simulated inundation area under design rainfall scenarios using different approximation methods.

In terms of detecting potential inundated areas, all approaches show only minor variations with high POD values ranging between 0.85 and 0.95. The NDA and RDA approaches even slightly outperform FSN, RSC, and RSN under extreme rainfall conditions, aligning with previous results in Sect. 4.5. As rainfall return periods decrease, the FAR values of the RSC, RSN, NDA and RDA approaches gradually increase, with the latter two showing much more significant variations. This indicates that the performance of the equivalent drainage rate i_d, calibrated using a single event, still has high uncertainty when applied to rainfall scenarios of varying intensities. Therefore, while differences in PODs are minimal, the dramatically higher FARs of equifinal approaches widen the overall performance gap, as measured by CSI, especially for rainfall events with return periods below 50 years. Regarding bias, physical approaches that adopt reconstructed networks but without controls consistently underestimate the total inundated area with BIAS values less than 1. In contrast, other equifinal approaches tend to overestimate the inundated area. Specifically, RDA and NDA can overestimate the total inundated area by more than 8 % and 10 %, respectively, for rainfall events with return periods less than 100 years. Thus, a constant equivalent drainage rate is insufficient to achieve an accurate simulation of sewer drainage effects under diverse rainfall conditions. However, physical approaches can also yield biased estimations of the inundated area if regulated facilities are not properly configured.

From the perspective of water balance, the total inundated area (or volume) is closely related to the volume of rainwater drained outside the system. Thus, we further compare the discharge processes at outfalls, simulated using different physical approaches (see Fig. S3). Being the solution reconstructed using the most detailed information, FSR clearly reveals insufficient drainage capacities at most outfalls, such as OF1. When regulated facilities are excluded but the sewer layout remains unchanged, FSN demonstrates significantly higher discharges, especially at major outflow points such as OF2 and OF3, recording peak flows that are over 300 % higher, thereby leading to an underestimation of total inundated areas (Fig. 12). When the sewer layout is further simplified using the roads layout, the simulated outflows decrease by approximately 27 % due to reduced inlets that collect rainwater into the sewer system. This attribution can be further supported by the improved alignment of outflow hydrographs between RSN and RSC, which exhibit fewer truncations during peak flows under regulated conditions when compared to FSR and FSN. This indicates an alleviation of drainage deficiency due to reduced rainwater inflow. However, certain outflow points, such as OF3, still yield substantially higher discharges with simplified networks. A comparison of the layouts used in FSN and RSN reveals that a main pipe, measuring 6.8 m × 2.6 m (near OF1 and OF8, cf. Fig. 13) is removed during simplification, blocking the drainage pathway that would otherwise divert upstream rainwater to the northern outflow points such as OF1 and OF8, eventually magnifying the inflow to OF3. Thus, for real urban environments with complicated sewer layouts due to practical factors such as multistage construction, direct geometric simplification through major roads may remove some critical connectivity inside a graph-like drainage system and thus fail to accurately represent the drainage effects of the sewer network due to topological incompleteness.

Figure 13 Comparison of simulated flows at 90 percentiles (Q₉₀) of pipes under design rainfall scenarios using physical approaches with/without network simplification. (a) FSN approach. (b) RSN approach. Basemaps are provided by Tianditu (https://map.tianditu.gov.cn/).

4.7 Outlook and limitations

The proposed graph-reconstruction-based approach offers an efficient solution for generating physical sewer models using available data on sewer layouts and surface elevation, which can be obtained from local surveys (Xing et al., 2022), engineering drawings (Lyu et al., 2018), or open source software such as SWMManywhere (Dobson et al., 2025). Given the significant influence of network topology and regulated facilities on inundation simulation results (as discussed in Sect. 4.5 and 4.6), in addition to the aforementioned data sources, it is imperative to further implement intelligent surveying systems targeting relevant critical components such as inlets and outfalls. These data can then be integrated with the proposed reconstruction algorithm to achieve more reliable pluvial flood modeling. In particular, outfall identification still remains a difficult task, especially for decentralized sewer networks in relatively flat urban areas, challenging the applicability of existing topography-driven approaches (Dobson et al., 2025) and highlighting the need for user-defined configurations based on prior knowledge of the area (Reyes-Silva et al., 2023). In addition, it should be acknowledged that given the inherent quality issues of multi-source geographical datasets, the mathematical programming procedures in the proposed approach may occasionally fail because of infeasible configurations such as over-elevated outfalls and dangling pipelines. While introducing data cleaning procedures may help alleviate such issues and streamline the downstreaming reconstruction process, it also necessitates a trade-off between algorithm robustness and accuracy owing to embedded assumptions during quality control, which deserves further investigations. In addition to related cleaning functions implemented in existing packages like SWMManywhere (Dobson et al., 2025), we propose leveraging the built-in functionality of irreducible inconsistent subset (IIS) computation within our optimization-based framework, which can be steadily accessed in solvers such as Gurobi (https://docs.gurobi.com/projects/optimizer, last access: 14 October 2025). This tool can help identify the specific constraints responsible for infeasibility and thus support the spatial pinpointing of ill-conditioned sewer elements that over-constrain the optimization and may require further correction or elimination. In summary, additional integration with other existing tools will further enhance the practicality of the proposed reconstruction method and facilitate a more realistic representation of drainage effects in urban pluvial flood simulation.

5 Conclusions

In this work, we present a graph-reconstruction-based approach for building a physical sewer hydraulic model, leveraging incomplete information and a surface elevation model, while adhering to design guidelines. In contrast to conventional optimal-design-based approaches, our approach seeks to reconstruct the most possible sewer network with related hydraulic properties, such as gravitational flow direction and nodal invert elevation, constrained by minimal deviation from the initial status based on available information. The resulting mathematical programming problems can be equivalently transformed into corresponding linearized formulations, ensuring their computational efficiency for large-scale applications. We apply the proposed approach to Yinchuan's main city zone, characterized by multiple outfalls consisting of wastewater treatment plants and pump stations, yielding the following key insights:

• The reconstructed sewer hydraulic module, integrated into a 1D/2D coupled hydrologic-hydrodynamic model, achieves high accuracy in simulating maximum inundation depths, with an R² of 0.95. A comparison with road-layout-based simplifications further demonstrates the robustness of our approach in reconstructing sewer hydraulic properties, as the simplified model maintains a similar distribution of pipe slopes.

• While equifinal approaches ease parameter calibration by encapsulating complex sewer drainage processes into equivalent drainage rates, the proposed reconstruction approach can reliably achieve superior performance in simulating inundated areas across various rainfall scenarios. This improvement stems from its physically grounded structure in simulating the bidirectional interaction between the ground surface and the sewer system, indicating that physical reconstruction is preferable to approximate sewer drainage effects in operational applications.

• The completeness of the available information has non-negligible effects on the simulated hydraulic performance of reconstructed sewer networks. Simplifications involving road-based layouts and regulated facility removal lead to complex deviation patterns in local maximum inundation depths and outflow hydrographs, with both overestimations and underestimations arising from the combined effects of surface microtopography, sewer connectivity modification, and regulated facility operation.

In conclusion, we provide an efficient sewer reconstruction approach that improves simulation performance and understanding of sewer drainage effects during urban pluvial floods, accounting for varying levels of information completeness. However, uncertainties persist in reconstructed physical models when detailed information on regulated facilities and drainage pathways is lacking. Thus, future research will focus on integrating advanced techniques, such as data assimilation and intelligent surveying systems, with reconstructed physical models to augment the realism of sewer drainage process representation.