Modeling Urban Flood Inundation and Recession Impacted by Manholes

Abstract

Urban flooding, caused by unusually intense rainfall and failure of storm water drainage, has become more frequent and severe in many cities around the world. Most of the earlier studies focused on overland flooding caused by intense rainfall, with little attention given to floods caused by failures of the drainage system. However, the drainage system contributions to flood vulnerability have increased over time as they aged and became inadequate to handle the design floods. Adaption of the drainages for such vulnerability requires a quantitative assessment of their contribution to flood levels and spatial extent during and after flooding events. Here, we couple the one-dimensional Storm Water Management Model (SWMM) to a new flood inundation and recession model (namely FIRM) to characterize the spatial extent and depth of manhole flooding and recession. The manhole overflow from the SWMM model and a fine-resolution elevation map are applied as inputs in FIRM to delineate the spatial extent and depth of flooding during and aftermath of a storm event. The model is tested for two manhole flooding events in the City of Edmonds in Washington, USA. Our two case studies show reasonable match between the observed and modeled flood spatial extents and highlight the importance of considering manholes in urban flood simulations.

1. Introduction

Flooding is one of the most frequent weather-related natural disasters and affects many people around the world every year [1]. Major floods often cause significant impacts to communities and economies [1,2]. For example, in the United States alone, flood damages cost $260 billion (USD) per year from 1980 to 2013 [3]. The National Flood Insurance Program (NFIP) paid on average $2.9 billion a year between 2000 and 2018 [4]. Similarly, flooding caused more than 700 fatalities and at least €25 billion economic losses in Europe between 1998 and 2004 [5]. Global warming is expected to lead to more frequent extreme precipitation events, increasing flood hazards in many cities around the world [2,6,7].

Despite the severe damages that flooding can induce, floods are an important part of life in various regions of the world where people have been adapting for centuries. They support riparian ecosystems dependent on flood inundated zones [8], and they are the principal source of groundwater recharge in many arid and semi-arid settings [9,10,11,12,13,14]. The negative impacts of floods, however, are considerable, and are mostly associated with the unexpected magnitudes and frequencies of floods, which are connected to climate change, rapid expansion of urbanized areas, and inadequate and aging urban drainage systems [2,15,16,17,18,19,20].

Urban flooding and associated damages to properties account for 73% of the $107.8 billion total damages caused by floods from 1960 to 2016 in the United States [21]. Thus, accurate flood monitoring and estimation are essential to reduce flood impacts and vulnerabilities while supporting urban planning and ecosystems.

The spatial and temporal characteristics of floods in urban areas are complex due to the widespread change to the land uses [22], which introduces micro-urban features such as buildings, roads and drainage networks [23,24]. The specific urban infrastructures that affect flood form a storm event includes the type and geometry of building, garage ramps, light wells, pillars, and yards at or just beneath the ground surface [25,26]. Overall, it is an established concept that an increase in impervious areas and connected conveyance systems in urban areas increase peak discharges and volumes [27,28]. However, flood mitigation and management requires detailed information about the spatial extents of floods, their water levels (i.e., flood depth), and flow velocities [29]. Such information is often impractical to measure directly. Consequently, empirical and hydraulic models are widely used to estimate these parameters [29,30,31] and assess associated flooding risk [18,32,33,34,35,36].

Although hydrological modeling can simulate both surface and subsurface processes adequately at a watershed level, most of them are unable to simulate urban flooding accurately [37]. This is partly due to the difficulty of defining model boundary conditions and the complex nature of flood propagation in urban areas. In addition, most flood inundation models lack the potential contribution of manholes overflow to the flooding. The majority of these models are commercial (e.g., MIKE FLOOD, XPSWMM, and FLO2D), and are not accessible for most users. They are complex, require a large number of datasets, and computational resources, which make them ineffective for most applications [29,38]. Thus, simplified flood inundation modeling techniques are commonly used to determine the spatial extent of flooding in urban areas [37,39]. These models only require a digital elevation map and mathematical representation of floodwater propagation in a given area [40]. These models were used to identify the potential flood-prone regions during a given storm event, but not the aftermath of the flood event. Furthermore, the simplified models often do not simulate the watershed and the drainage system. Consequently, coupling hydrologic models, one-dimensional (1D) hydrodynamic models, and simplified flood inundation models are gaining attention as an alternative approach to simulate flood inundations in urban areas [29].

Despite well-documented effects of manholes on urban flooding [19,20,37,41,42], there exists relatively limited research [43,44,45,46,47] that directly incorporate manhole overflow in simulations of urban flooding. Chen et al. [48] used coupled surface and sewer flow modeling and found approximately 2 m surge flow from manholes. Leandro and Martins [43] used a two-dimensional (2D) flood inundation model to simulate bi-directional flow interaction between sewer and overland flow, to estimate the volume of sewer serge from manholes ranging from 15,992–18,404 m³ and a maximum possible depth of 0.8 m. Son et al. [44] coupled a one-dimensional storm water management model with a 2D overland flow model from a manhole overland flow with a maximum overland flow, ranging from 2–5 m³/s at different manhole locations, and estimated a 0.9 m flood depth and 2.5 m/s flood velocity. Jang et al. [46] used a coupled one-dimensional sewer flow with the two-dimensional overland flow and estimated manhole overflow depth up to 2 m. Seyoum et al. [47] coupled 1D sewer and 2D dimensional flood inundation models to simulate urban flooding and showed the combined flood depth variations from 0.3 to 0.8 m in their study area. Manhole overflow is a critical issue in urban areas, yet their contribution to flooding is not well understood. Most hydrodynamic models assume that excess water will pond around the manhole and return back or will be lost from the system after the flood recedes [49]. Major cities have aging infrastructures and drainage systems designed and built more than a decade years ago [50]. These systems are increasingly underperforming due to their design assumption of stationary storms and flood events [51,52]. This assumption often causes inadequacy to handle the rising flood risk caused by increased storms and impervious layers.

The main objective of the study is to develop and test a new flood inundation and recession method (FIRM) that can readily be used to simulate excess floodwaters generated from manholes. The United States Environmental Protection Agency (EPA)-Storm Water Management Model (SWMM) is commonly used to simulate the complex hydrological and hydraulic processes in urban areas, but it does not have the configuration to simulate flood inundation and recession from manholes overflows. Our method uses manhole overland flow volumes and depths from SWMM output and digital elevation data of the study area to simulate flood inundation, recession, and depth. The flood inundation model uses the flat-water assumption to distribute the overflow. The recession model uses the location of manhole and surrounding topographic variation to determine whether the inundated region is going to drain to the manhole or pond in localized regions. The model was tested using a synthetic case study and a flooding event in Edmonds, an urban region in Washington State. Since the inputs (elevation and overflow) and outputs (flood area and depth) are known, the synthetic case study was used as a proof-of-concept to validate the model accuracy before applying it to a real-world problem. It allowed us to test the model performance using different scenarios under variable hypothetical case studies. For our Edmonds case study, previous study reports, social media, and news reports were used to delineate the flood boundary (i.e., used as observational data). The flooded area estimated by FIRM was compared against the reconstructed flood area visually and using statistical measures to evaluate our model’s ability to identify areas that were flooded versus dry during the actual flood event.

2. Methods

This section details (a) the hydrodynamic model “SWMM” applied to Edmond case study (Section 2.1), (b) the model domain, calibration, and validation (Section 2.2), (c) the simulation of manholes’ overflow and recession using FIRM (Section 2.3), and (d) the evaluation of the FIRM performance using statistical measures (Section 2.4).

In addition to the synthetic case study used for proof-of-concept, the FIRM was applied to simulate manhole flooding in the city of Edmonds. The SWMM model domain was discretized into detailed sub-catchments to incorporate key hydrologic and hydraulic components. Three years of sub-daily rainfall and lake level data were used to develop and calibrate the model using the automated differential evolution optimization methods [53]. The SWMM model was validated against a one-year simulation. Following the SWMM simulation, we simulated the spatial extent and depth of flood inundation and recession from flooded manholes. A high-resolution digital surface model (DSM) was used due to its detail information of surface features, including infrastructures. Unlike to the digital elevation map (DEM), the DSM contains surface elevation and surface features [54].

2.1. Hydrodynamic Modeling Using SWMM

SWMM is one of the most widely used hydrologic and hydraulic models in urban settings [55]. It is capable of simulating event-based or continuous rainfall-runoff processes that are useful for both water quality and quantity analyses in urban areas [56].

SWMM uses spatially distributed and temporally discrete processes to simulate the hydrological and hydraulic state variables [49]. As the simulation progresses, the state variables will be updated and stored as follow.

$${\mathrm{X}}_{\mathrm{t}}=\mathrm{f}\left({\mathrm{X}}_{\mathrm{t}-1},{\mathrm{I}}_{\mathrm{t}},\mathrm{P}\right),$$

(1)

$${\mathrm{Y}}_{\mathrm{t}}=\mathrm{g}\left({\mathrm{X}}_{\mathrm{t}},\mathrm{P}\right),$$

(2)

where f and g represent the functions that calculate and update the state and output variables respectively. X_t represents state variables (such as flow rate and depth in a drainage network link), Y_t represents output variables (such as runoff flow rate at each sub-catchments and outlet), P represents the constant parameter, I_t represents input variables (such as rainfall and temperature) at a given time.

The hydraulic simulation of SWMM involves water and contaminants transport through the conveyance portion of the drainage network. External flow sources entered the drainage network using inflow nodes and transported through pipes and storage components and finally exit at outflow nodes [49]. The flow equation through links can be solved using the dynamic and kinematic wave routing equation. The dynamic wave analysis uses the Saint-Venant flow equation, whereas the kinematic wave method uses the simplified form of the momentum equation to estimate the flow condition though the conveyance system [49]. The dynamic wave analysis has advantages in simulating gradually varied flow conditions, such as surge, in the drainage systems. This results in having the dynamic wave method, to have different pressure and friction force part of the momentum. We used a dynamic wave routing method to simulate the gradually varied flow in urban drainage system.

2.2. Manhole Overland Flow Inundation and Recession Modeling

SWMM considers overflow from manholes as either surface ponding for a certain period and return back to the drainage or lose from the system [49]. However, these approaches do not allow for direct estimation of the associated flood spatial and temporal extents. Besides, all the overland flow generated from the node does not necessarily return to the drainage, as some will recess, and others might be isolated from the manholes and ponded in depressions. In this study, we developed a simplified grid-based flood module to propagate and recess overflow from and to a manhole spatially. The module requires a gridded surface elevation map, locations of manholes, and total overflow volume and depth at the manholes. The flood inundation computation assumes that (i) water flows from a higher elevation to lower elevation because of gravity, (ii) water spreads spatially by maintaining its level surface (‘flat’ water assumption), and (iii) flood fills first the nearest and connected cell with the lowest elevation. The module estimates both the flood depth and spatial extent by propagating the excess pressurized water from the manholes to the surrounding areas according to the topographic variations. A high-resolution (1 m × 1 m) digital surface terrain from Light Detection and Ranging (LiDAR) was used to capture the flooding extents accurately.

The flood recession module uses the location of manholes, the areal extent of the flood, as well as topography of the region to determine the aftermath of the flooded area. The flood recession method assumes that all flooded cells drain back to the manholes if there exists a flow path connecting them to the manholes; otherwise, it will remain ponded. Some of the flooded water in the local depressions can be disconnected from the manhole and will not fully be drained. Otherwise, the ponded locations can be distant from the manholes due to the topographic barrier.

2.2.1. Manhole Overland Flow Inundation Modeling

Our simulation of flood inundation is based on the elevation variations in neighboring cells. Starting from the cells containing the flooded manhole, floodwaters propagate to neighboring cells if the cell elevation is lower-than and connected-to a flooded cell. We used the D8 neighboring algorithm, which evaluates the elevations of the eight adjacent neighboring cells to each flooded cell, to determine the preferred flow direction (Figure 1). The flooded cells will maintain level flood surface and expand spatially as the lowest elevations are filled with available excess water. For multiple manholes overland flow, iterations are used to distribute all the overflows by maintaining the level flooded surface.

Before the flooding starts, the area is assumed dry, or the flood depths at each grid cell are known. Once the manhole overflow starts, the floodwater inundates neighboring and connected grid cells that have lower elevations (Figure 1). This inundation process will stop when the excess overland flow is completely allocated, and no water left to inundate the next dry cell. The flood inundation calculation is performed based on three main cases. In the first case, if the manhole cell elevation is smaller than that of the neighboring cells, the overflow will accumulate in the manhole cell until it reaches the minimum elevation of the neighboring grid cell. Once the overflow surface elevation exceeds the minimum elevation of a neighboring cell, the water will propagate spatially as long as there is enough overflow. In the second case, if the elevation of the manhole is higher than the elevations of neighbor cells, the overflow will be allocated automatically to the neighboring cells regardless of their slopes but maintaining level flooded surface. In the third case, if the elevation of the manhole is the same as the elevation of any neighboring cell, the excess water will be distributed equally among cells with the same elevations.

Figure 1 and Figure 2 show a schematic representation of flood inundations caused by overflow from one manhole and two manholes respectively. When the flood surface (surface elevation plus inundation depth) exceeds the elevation of the surrounding cells, the water will start to flow to the cells with the lowest elevations. As the overland flow increases, the neighboring cells will be gradually flooded while maintaining level flood surface. For example, when the overland volume is 49 units (Figure 1d), the extent of the flood coverage is first determined, based on the amount of the overland flow and elevation, and then the floodwater is distributed iteratively among the grids within the flooded area to maintain level flood surface. This iterative process considers the available flood volume, and recursively fills the next small neighboring cells.

Not explicitly considering the slope and land covers are the main limitation of the presented flood inundation approach. Similar to other grid-based hydrological and hydraulic modeling, the DSM resolution poses a high level of uncertainty in the model representation. Thus, the use of fine-resolution DSM and surface maps are critical for estimating the floods from manholes accurately [37]. For this study, we have used a one-meter resolution Light Detection and Ranging (LiDAR) DSM obtained from Washington State Department of Natural Resources. The model and the test cases do not include manhole inundation caused my previous storm and overland flow the drainage. Thus, we assumed dry terrain prior to the overflow from the manholes. However, in most cases, as the reviewer correctly stated, the manholes and their surrounding areas might already be flooded before the overflow. One simple modification is to adjust the terrain elevations for the inundation depth and then simulate the manhole overflow on top. Such an approach may work if the inundations surrounding the manholes that resulted from the previous storm were stationary, and have known flood depths. However, in most cases, the overland and overflow inundations from the drainage area and manholes, respectively, happen simultaneously, requiring a coupled simulation of both inundations.

2.2.2. Recession Modeling Associated with Manholes

Similar to the overland flow inundation, the flood recession to a given manhole considers the elevations of the surrounding grid cells and their connectivity to the manhole. For any given cell, the flooded water above the elevations of the adjacent grid cells will drain to the manhole if a flow path exists. Otherwise, it will remain ponded in the depressions. For example, when a manhole elevation is higher than the neighboring cells, only partial recession will occur. Starting from the manhole cell and expanding outward using the D8 algorithm, the recession model identifies flooded cells and drain them completely if their elevations are above one of the eight neighboring cells. Else drains part of the flooded water above the minimum elevation of the eight neighboring cells. As the grid cell expand spatially by searching adjacent eight neighboring, the connectivity of the target cell and the neighbor is traced using a flow path algorithm. For each flooded cell, the flow path algorithm detects the presence of any possible flow path to the manhole.

Figure 3 shows examples of flood recession based on single and multiple manholes. The result from the flood inundation estimate is given in Figure 3a. It is assumed that when rainfall ceased, the holding capacity of the drainage system decreases, and overland flow starts to recess toward the manholes. Figure 3b,c illustrate how flooded surfaces are recessing for a single and two manholes scenarios. The recession associated with a single manhole (Figure 3b) indicates two regions of ponding—one next to the manhole and another away from manhole caused by local topographic barrier. The ponded regions away from the manhole are not connected to the manhole hydraulically. The second case introduced an additional manhole in the depression region, and thereby the flood can be drained by the two independent manholes. The left manhole recess similarly as Figure 3b. The second manhole drains the ponded floodwater caused by a topographic barrier (Figure 3c).

2.3. Study Site: The Hall Creek Watershed

The Hall Creek watershed is an urban watershed (Figure 4), containing four major cities near Seattle, WA. These cities are Edmonds, Esperance, Lynnwood, and Mountlake Terrace, which are part of the northern Seattle–Tacoma–Bellevue metropolitan region. The predominant land cover type in the watershed is a developed urban region, which accounts for 96% of the land cover. The rest is covered by forest and water bodies (Figure 4b). The watershed is frequently affected by prolonged storm and flooding events.

The Hall Creek is intermittent and drains toward Lake Ballinger. It is the main tributary to the Lake Ballinger [57,58], which discharges to the downstream McAleer Creek. The creek does not have a monitoring discharge and stage data. Since the lake level is highly influenced by storm events (Figure 5) and the flow from the Hall Creek, the available level data were used to calibrate the SWMM model.

The flooding in the study area is mostly caused by storm events, while the urban area and street intersections are often affected by manhole overflow. For example, during the 19 September 2016 storm events, two manholes overflowed and caused flooding in the surrounding area. This flood event was used to validate our flood inundation and recession methodologies. One of the manholes flooded areas across a highway (Case 1), while the other caused flooding alongside the highway (Case 2) (Figure 6).

2.3.1. Data

High resolution observed meteorological data (obtained from the King County’s watersheds and rivers database) and sub-hourly lake level data (obtained from the city of Edmonds) were used to develop the SWMM model. Surface elevation from a one-meter resolution LiDAR data (obtained from Washington State Department of Natural Resources) were used for the flood inundation and recession modeling. We have extracted the conveyance network system of the city of Edmonds and Mountlake Terrace from their respective Geographical Information System (GIS) Department. We have identified the type, location, and possible flow direction of the sewer system. The land cover data were obtained from the US Department of Agriculture and were used to estimate the percentage of impervious layer. For the purpose of calibration and validation of the FIRM model, the observed flood boundaries are delineated using related images (Figure 6) and texts from social media users (such as locations). We used Google Earth and high-resolution LiDAR data to determine the relative topographical variation and delineated the inferred flood boundary by considering 360-street and panoramic view. We also used the city of Edmonds GIS-dataset to correlate the inferred flood boundary and building footprints to check if there exists a mismatch between the building blocks and the street boundaries.

2.3.2. SWMM Model

The SWMM model was discretized into 32 sub-catchments based on the hydrological and drainage network criteria. These criteria include percent of land cover, slope, availability of conveyance network system, and percent of impervious layers. A 5% of land use, slope, and soil type is used to subdivide the catchment. The sub-catchment layer, conduit layer and node layers are used for urban watershed discretization. After the watershed was discretized, the external inputs such as precipitation, temperature, and evapotranspiration were extracted and applied for each sub-catchment. The hydraulic properties of manholes, storage, ditches, culverts, and other structures were incorporated. There are total of 106 manholes, one storage, and 108 nodes connected by conduits in the study area. Only 40% of the manholes were considered in our study based on consideration of data availabilities and the computational requirement of the SWMM model. The Hall Creek flux are represented using SWMM’s inflow package. Despite the reasonable performance of our model in representing the hydrology and hydraulics of the study area (Section 3.1), not considering all the hydraulic structures in our model might have introduced some level of uncertainty that need for a further study. SWMM simulations can take hours if the model domain is large and accompanied by detailed complex hydraulic structures and sub-daily meteorological and hydrological input variables. The model simulation was conducted using the “swmmr” R-package [55].

Model calibration can be performed using either manual or automated method [59]. In this study, we used both methods of calibration to ensure usage of their advantages. First, we used manual calibration to identify sensitive parameters. The sensitive parameters was then further calibrated using differential evolution (DE) method, which finds the global optimum parameters values for continuous and differentiable functions [60], based on successive generation and transformation of the parameters values under a given fitness-measured criteria [53]. The DE requires defining parameters upper and lower bounds, objective function, and the lower or the upper optimal solution goal. The algorithm starts by randomly dividing the parameters values in to three distinct populations. The parameters values from each population are then combined to generate the next sets of populations that minimize the objective function. To ensure global optimal solution, the algorithm uses mutation to include non-optimal parameters values in the new populations [53]. The evolution continues until it meets the objective function criteria.

Due to the lack of discharge observational data for the Hall Creek, the model calibration and validation were performed based on the lake level fluctuation of Lake Ballinger, which is located at the outlet of the creek. Since the creek is a main feed to the lake, the lake level fluctuations reflect the changes in the creek discharge. The model calibration includes the initial condition of the model, which enables to determine the calibration parameter ranges, and optimization of the parameters using differential evolution optimized method [61]. Nearly three (2.7) years of data were used to calibrate the model, and one-year of data were used for validation.

The initialization was used to identify sensitive parameters and their respective parameter ranges. The model was then optimized using the differential evolution method, namely the “DEoptim” packages [61]. Differential Evolution (DE) is a genetic algorithm that finds global optimum values for continuous and differentiable functions [60] based on successive generation and transformation of the parameter sets [53]. The DE requires parameters for upper and lower bounds and an objective function for the optimization. During each evolution, in addition to identifying the better parameter sets (or population), the algorithm also introduces a random change to those parameter sets to ultimately get the global optimal parameter values.

The model simulation includes spin up, main (calibration), and post-audit (validation) simulations. Fitness measure statistics, including Nash–Sutcliffe efficiency (NSE), percent bias (PBIAS), root-mean-square error (RMSE), ratio of the RMSE to the standard deviation of measured data (RSR), and Kling–Gupta efficiency (KGE) were used to evaluate the model calibration results. The NSE compares the variance of the residuals (or fitting difference) with the variance of observed lake levels [62]. The PBIAS measures the average residuals or deviations of model results from observed lake levels [63]. The RMSE measures how spread out these residuals from the model results. The RSR is a normalized RMSE by the standard deviation of the observed lake level [64]. The KGE uses the idea of diagnostic decomposition, where the NSE is breakdown into three components (the relative importance of correlation, bias, and variance difference) [63,65]. The KGE ranges from negative infinity to one, with the optimal model prediction having the KGE value close to one.

Where Y^obs is observed, Y^sim is simulated, Y^mean is the mean of the observed lake level change, r is correlation coefficient between the modeled and observed lake levels; γ is a ratio between the standard deviation of modeled and observed lake levels and β is a ratio between the standard deviation and mean of the modeled and observed lake levels (

Table 1. Representing statistical model performance criteria.

Statistics	Ranges	Optimal Value
$\mathrm{NSE}=1-\left[\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Y}}_{\mathrm{i}}^{\mathrm{obs}}-{\mathrm{Y}}_{\mathrm{i}}^{\mathrm{sim}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Y}}_{\mathrm{i}}^{\mathrm{obs}}-{\mathrm{Y}}^{\mathrm{mean}}\right)}^2}\right]$	−∞ to 1	1
$\mathrm{PBIAS}=\left[\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{Y}}_{\mathrm{i}}^{\mathrm{obs}}-{\mathrm{Y}}_{\mathrm{i}}^{\mathrm{sim}}\right)\times 100}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{Y}}_{\mathrm{i}}^{\mathrm{obs}}\right)}\right]$	0 to 100	0
$\mathrm{RSR}=\frac{\mathrm{RMSE}}{\mathrm{STDEVobs}}=\frac{\left[\sqrt{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Y}}_{\mathrm{i}}^{\mathrm{Obs}}-{\mathrm{Y}}_{\mathrm{i}}^{\mathrm{sim}}\right)}^2}\right]}{\left[\sqrt{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Y}}_{\mathrm{i}}^{\mathrm{Obs}}-{\mathrm{Y}}^{\mathrm{mean}}\right)}^2}\right]}$	0 to 1	0
$\mathrm{KGE}=1-\sqrt{{\left(\mathrm{r}-1\right)}^2+{\left(\mathsf{\gamma }-1\right)}^2+{\left(\mathsf{\beta }-1\right)}^2}$ $\mathsf{\beta }=\frac{{\mathsf{\mu }}_{\mathrm{s}}}{{\mathsf{\mu }}_{\mathrm{o}}}$ $\mathsf{\gamma }=\frac{{\mathrm{CV}}_{\mathrm{s}}}{{\mathrm{CV}}_{\mathrm{o}}}=\frac{{\mathsf{\sigma }}_{\mathrm{s}}/{\mathsf{\mu }}_{\mathrm{s}}}{{\mathsf{\sigma }}_{\mathrm{o}}/{\mathsf{\mu }}_{\mathrm{o}}}$	0 to 1	1

).

Manhole overflows were extracted from the calibrated SWMM model and were used as input for the FIRM simulation. The spatial extent of the model is compared to the observed flood regions. Since there was no direct measurement of flood spatial extent and depth, the observed flood area was reconstructed based on pictures of the area taken during the flood event and obtained from social media twitter. The social media information includes both pictures, texts, and street names, allowing us to identify the exact locations and extent of the flood. The previous reports in the cities also indicated that there had been multiple incidences of pluvial flooding near side roads and along intersections. This information was used to delineate the observed flood region, which was then used to validate the flood inundation model. The process of flood boundary delineation is depicted, in Figure 6, which show how the flooded region in the study area was extracted from social media outlets. The images and the text by users are used to identify the exact location of the flooded regions. We used google-earth and high-resolution LiDAR data to determine the relative topographical variation and delineated the inferred flood boundary by considering 360-street and panoramic view. We also used the city of Edmonds GIS-dataset to correlate the inferred flood boundary and building footprints to check if there exist mismatch between the ground and the street boundaries.

2.4. Model Inundation Accuracy

The model’s ability to detect the spatial accuracy of the flood extent was evaluated based on the true positive rate (TPR), the positive predictive value (PPV), the modified fit (MF), and the modified bias (MB) methods. The TPR and PPV are derived from the confusion matrix [17], which is a 2-by−2 matrix containing the TPR and PPV for gridded simulated and observed flood conditions (

Table 2. A confusion matrix for flood performance.

	Flooded in Observed Boundary	Dry in Observed Boundary
Flooded in FIRM	True flood (TP)	False flood (FP)
Dry in FIRM	False dry (FN)	True dry (TN)

). These statistics were used recently to assess the flood inundation model performance in [39,66,67]. TPR measures how well the modeled flood region replicates the observed flood boundary. The maximum TPR (100%) represents that the model fully captures the observed flooded regions. The TPR indicates the model tendency to under-predict the flood hazard [39,66,67]. PPV measures how well the model captures flooded in the model but dry in observation. The value ranges from 0%, indicating over prediction of the flood extent, to 100% for accurately captured the observed boundary.

$$\mathrm{TPR}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}\times 100,$$

(3)

$$\mathrm{PPV}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}\times 100,$$

(4)

Flood inundation and recession model (FIRM); true positive rate (TPR); the positive predictive value (PPV); true positive (TP); false negative (FN); false positive (FP); true negative (TN).

Where TP represents the flooded regions in both the observation and model simulation, FP represents the flooded region in the model but dry in the observation, and FN represents regions flooded in reality (i.e., in observation) but simulated as dry. The TPR and PPV percentages represent the overlapping rate between simulated and observed flood areas. Higher percentage values of TPR and PPV indicate higher accuracy of the flood inundation model. The two statistics must be used in combination to measure the accuracy of the model since they each evaluate the different performance of the model. Specifically, the TPR and PPV measure how well the model captures the observed flood and flooded pixel that are dry in observation, respectively. For example, the TRP value can be 100% if the model captures all the observed flooded cells even though it may also consider some dry cells as flooded (refer to Equation (3)). Similarly, the PPV value can be 100% if the model captures all the dry cell even though it may also consider some flooded cell as dry.

Other methods of model performance evaluation are fit and bias indicators, which are also commonly used for flood inundation modeling [40,44,68,69,70]. Previous studies used both fit and bias indicators to compare flood inundation extents between different models. For this research, we infer observed flood extents based on street photos taken during the actual flood event and compare them with the flood inundation results from our model. The modified fit indicator is calculated based on overlapping areas between observed and simulated inundated areas. The indicator ranges from 0% to 100% for poor and ideal model performances, respectively. While the modified bias indicates the overall difference between the simulated and observed flood extents. Positive and negative modified biases indicate an overestimation and underestimation of flood extents by the model, respectively.

$$\mathrm{Modified}\mathrm{fit}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}+\mathrm{FN}}\times 100,$$

(5)

$$\mathrm{Modified}\mathrm{bias}=\left(\frac{\mathrm{TP}+\mathrm{FP}}{\mathrm{TP}+\mathrm{FN}}-1\right)\times 100,$$

(6)

3. Results and Discussion

3.1. SWMM Model Calibration and Validation

We compared simulated versus observed lake level changes, and show our model captures the lake level fluctuations reasonably well (Figure 7, Figure 8 and Figure 9). The main objective of the SWMM simulation was to estimate the flood condition in Edmonds from the 19 September, 2016 storm event. For model stability and to represent preexisting hydrological conditions such as soil moisture content, we considered ten months (from 1 August, 2015 to 31 May 2016) of simulation as a spin up. The model calibrated using data from 1 June 2015 to 31 January 2018; and validated using data from 1 February, 2018 to 10 January 2019. After the sensitivity analysis using the manual calibration, we identified five parameters for model calibration. These includes the imperviousness percentage, width, roughness coefficient, depression storage, and the hydraulic conductivity of the soil. The calibration is performed in the DEoptim R-package.

Table 3. Model parameters upper and lower limits, as well as their optimal values.

Parameters	Lower–Upper Bound	Optimal Values
Impervious (%)	25–90	70
Width (m)	150–300	152
Roughness (−)	0.01–0.03	0.012
Depression Storage (mm)	1.2–5.2	1.78
Hydraulic Conductivity (mm/h)	0.1–3	0.11

indicates the upper and lower limits of the model parameters used in the calibration process.

Figure 7 and Figure 8 shows the comparison of the model simulation and observed time series plots for daily and monthly average lake level change for the calibration and validation simulation periods, respectively. The figures demonstrate that the lake level changes as a result of storm events were reasonably captured for the calibration period. The validation results for both daily and monthly observed and simulated lake levels show the ability of the model to predict beyond the calibration period.

Figure 9 indicates the correlation between the observed and simulated lake level change for spin-up, the main model simulation, and the validation period. The regression coefficient (R²) is 0.42 for the spin up period, 0.83 for the calibration period, and 0.77 for the validation model simulation period. The correlation coefficients for the calibration and validation period also confirm the model captured the observed lake level reasonably well. The simulation was also evaluated using the NSE, KGE, RSR, and PBIAS.

The statistical summary of the spin up, calibration and validation simulations results are presented in

Table 4. Model performance statistics to evaluate the Storm Water Management Model (SWMM) daily and monthly lake water level simulations.

Simulation	KGE		NSE		RSR		PBIAS		Performance Rating [71]
	Daily	Mon	Daily	Mon	Daily	Mon	Daily	Mon	Daily	Mon
Spin Up	0.64	0.61	−0.31	−1.15	1.14	RSR	−0.10	−0.10	Unsat *	Unsat *
Calibration	0.91	0.96	0.82	0.94	0.43	1.39	0.00	0.00	V. good ^	V. good ^
Validation	0.88	0.95	0.67	0.81	0.57	0.24	0.00	0.00	Good	V. good ^

Kling–Gupta efficiency (KGE); Nash–Sutcliffe efficiency (NSE); ratio of the RMSE to the standard deviation of measured data (RSR); percent bias (PBIAS); * Unsatisfactory; ^ Very good.

. For the calibration and validation periods, the KGE are 0.91 and 0.88, respectively, while the NSE is 0.82 and 0.67. These indicate satisfactory model performance. The RSR, which indicates the variation in the residuals, is between 0 and 0.5. This is considered a very good performance [71]. The PBIAS values are also close to 0, which confirm that the model simulates the observed water level with minimal bias. Compared to the daily model performance, the monthly model performance is improved.

3.2. Flood Inundation and Recession

The spatial extent and depth of the flood inundation were simulated for a pluvial flood event that happened on 19 September 2016 using FIRM. A detailed conveyance network system with sub-daily meteorological data, such as rainfall data, were used to capture the flood event using SWMM continuous simulation. The FIRM simulate the floods caused by overflow from two manholes, and its results were compared with the inferred flood boundary (Figure 10). Figure 10 represents the areal extent and depth of the flood inundation at two locations. The red dot-line represents the inferred flood boundary, and the black point represents manholes. The color represents flooding depth. The fine-resolution LiDAR data were able to identify detailed urban infrastructures, such as building food prints, and streets. The FIRM was able to identify elevated urban infrastructures and with low-lying streets. Overall, coupling both the 1D SWMM model with the 2D FIRM model was able to delineate the spatial extent and depth of the flood generated from manholes overflow in the study area.

The flood started from a given manhole and propagated spatially by filling any neighboring grid cells with lower elevations. As shown in Figure 10a, representing Case 1, the flood is concentrated across a highway going west to east. The topographic variations around the manhole are relatively small, with the mean slope of the flooded region being 1.2 degrees (Figure 11a). The flat slope, particularly along the street, enables the overland flow to inundate along the street perpendicular to the main highway. Based on the FIRM result, the areal extent, maximum flood depth and volume of the inundation region is 2129 m², 0.6 m, and 710 m³, respectively. The depth of the flood is controlled by the local elevation and the amount of the excess overland flow. The low-lying part of the street generally has deeper flood depth compared to the peripheral part of the flood extent. Hence, the pavement is often elevated compared with the street elevation. Another manhole flooding in our study (Figure 10b, Case 2) was used to evaluate the flood inundation simulation. In this case, the manhole is located in the steeper part of the street, where the average slope of the street is 5.2 degrees from the west to east (Figure 11b). Consequently, the flood is mostly concentrated along with the highway in the north and south direction, and did not spread much laterally. The volume, areal extent, and maximum depth of the flooded region are 38 m³, 426 m², and 0.15 m, respectively.

Based on the two cases, we observe that the flood inundation model considers the spatial heterogeneity of the surface feature to inundate, for example, flood inundation algorithm was able to differentiate building with streets (Figure 10a) and intersections of streets with varying slope and elevation (Figure 10b). FIRM inundates the lowest level and cover wider area for a relatively gentle region. Conversely, for the manhole located in a steeper region, the algorithm follows the preferred flow direction and inundate relatively smaller area.

To determine the flood recession, we assumed that once the storm ceased and the pipe full capacity decreased, the water eventually drains back to the conveyance system unless it is isolated from the manhole because of depression storages. Accordingly, some of the inundated floodwaters may recess, while the remaining waters are left as ponding associated with local topographic barriers. The result for Case 1 shows that most of the floodwaters drain since the manhole is located at a lower elevation. There is some ponded water away from the manhole due to possible topographic barriers between the manhole and the flooded areas (Figure 12a). Since the manhole is located at a lower elevation, most of the flooded water drained to the manhole. The volume, areal extent, and maximum depth of the ponded region are 15 m³, 89 m², and 0.17 m, respectively. For Case 2, the ponded water is generally concentrated near the manhole due to the local topographical depression around the manhole (Figure 12b). The flood volume, the areal and depth of the ponded water decreases to 3 m³, 85 m², and 0.06 m, respectively. The results showed the FIRM abilities to determine the maximum flood extent and the extent aftermath of a given storm. Each information is important to assess the flooding risk and associated potential short-term (flood inundation) and long-term (flood recession) impacts.

3.3. Model Inundation and Recession Accuracy

In addition to the visual comparison of the observed and simulated flood areas, the model performance evaluation for the inundation was carried out using statistical measures that compare the simulated and observed gridded flood areas. This enables us to identify and assess the model performance based on how accurately the model predicted the observed flooded and dry regions. We adopted the TPR and PPV from [66] and used the modified version of the fitting (MF) and bias (MB) indicators from [44,70]. The model inundation extents are compared with the inferred observed flood boundaries that were extracted from photos of the flood boundaries.

Table 5. Statistical evaluations of the flood inundation model based on inferred flood area at the two manhole locations (Case 1 and Case 2).

Inundation Model Performance	Case 1	Case 2
True positive rate, TPR (%)	89.04	71.31
Positive predictive value, PPV (%)	95.44	97.25
Modified fit, MF (%)	85.04	69.90
Modified bias, MB (%)	−6.71	−26.68

shows the model performance indicators used to evaluate the flood inundation model based on the inferred flood regions. The TPR for Case 1 (89%) indicates the model’s ability to capture the flooded grid cells, while for Case 2, the TRP is 71%, indicating the model predicted a relatively more flooded region as dry land. Thus, the model under predict the flood hazard in Case 2. The PPV for Case 1 and Case 2 are found to be 95.4% and 97.25%, respectively. These indicate the model’s ability to capture the observed none flooded cells. The MF of 85% for Case 1 indicates that the model has better agreements with the observed flood boundary compared to that of Case 2, which has MF of 69.90%, indicating the existence of relatively large variation between the predicted and observed flooded regions. The negative values of the MB for both cases indicates that the flooded regions in both cases were underestimated. Overall, the relative errors are relatively higher for Case 2 compared to Case 1. This is due to underestimating the flood hazard in Case 2 compare with the observed boundary, and possibly due lack of the FIRM to simulate the impact of direct rainfall during the flood events or lack to represent flood water loses into the buildings.

The low TPR values compared to the PPV values suggest that the flood inundation model underestimated the total flooded regions for both cases; but predicted well the flooded region within the observed flood boundaries (Figure 10). The underestimation of the flood areas might be due to the model inundation algorithm not incorporating the additional flooding resulted from direct precipitation or generated surface runoff. The relatively poor performance of the model for Case 2 might also be due to the inundation algorithm limitation to incorporate direct additional flood flux from upstream overland flow into the flooded regions. In addition, the relatively higher slopes in the area, which facilities rapid overland flow from the manhole toward the low elevated regions, may have impacted the model performance. Moreover, the relatively better model performance for Case 1 might be because of the relatively homogenous topography in the area, which is well represented by the 1-m LiDAR data.

Coupling the FIRM model and hydrodynamic model with projected future storm scenarios may help to identify areas that may experience future manhole flooding. This modeling capability can help to better assess flooding risk, and improve designs of storm water drainage systems in flood-prone urban areas. To further improve the work, it is important to consider direct rainfall during storm events and other sources or losses of floodwaters (e.g., possible loses of floodwater by draining into the building or addition of excess runoff from rooftops), as well as the land cover and slopes.

4. Conclusions

We have presented effective flood inundation and recession methodologies that use overflow from given manholes and topography of an urban region. We used the SWWM model to estimate the volume of overflow from manholes. In order to determine the associated flood depth and extent during and after storm events, we developed a flood inundation and recession model (FIRM) that uses high-resolution LIDAR elevation data. SWMM was developed on the basis of watershed characteristics and the drainage conveyance network in the area. SWMM was calibrated using a differential evolution optimization method and validated based on observed lake level data at the outlet of the watershed. The manhole overflows were extracted and used in FIRM to delineate the spatial extent and flood depths.

The spatial extent of the simulated flood area was compared with the observed flood boundary, which was derived from social media pictures and reports from the cities. Two case studies, based on flood events in Edmonds, WA, were considered to evaluate the flood inundation and recession model. In these case studies, the flood occurred across and along a main highway under different topographical characteristics. The results showed that the spatial extent of the flood regions is highly influenced by local topography and the position of the manholes. Particularly, the spatial arrangement of the manhole and the slope of nearby areas are crucial for determining the spatial extent, spatial heterogeneity of the flood depth, and selecting preferential flow paths to inundate low-lying areas. The model is able to capture the flood extent for manhole overland flow in fluvial flood events. Incorporating the direct impact of rainfall on the fluvial flood event can improve the representation of the physical process and the accuracy of the model. As the flood recession observation data are scarce, the performance of the flood recession model result was difficult to quantify. Finally, proper understanding and representation of the study area, the boundary condition, and engineering structures are important for the flood inundation and recession modeling associated with manhole overland flow.

Regional authorities can utilize the presented model (FIRM) by coupling with existing hydrodynamic modeling (e.g., SWMM) to quantify flood hazard based on pluvial generated overland flooding and manholes induced flooding in urban areas, where flood mechanism is complex and modified by local infrastructures. FIRM can be used to estimate the areal extent and depth of flood caused by manholes overland flow during a flood event (flood inundation) and after the event is over (flood recession). Because of the relative simplicity of the model and its uses of readily available data, the model can be used for a real-time assessment of flood progression and to identify potential impact areas. The model’s ability to simulate flood recession will also allow identifying areas where the flooded water will remain ponded for days after the floodwater subsides. The ponded water, or the floodwater not drained, can impact human health and properties. In addition to the real-time forecast of flood inundation and estimation of the aftermath ponding condition for the existing drainage system, the model can be used to design better a new or retrofit the current drainage system to minimize the overflow and ponding after the flood events. The model can be used to assess the flood condition under multiple storms, watershed conditions, and drainage scenarios. It can contribute to our understanding of climate change and appropriate engineering designs for mitigation.