Coupled sink and flow accumulation analyses with single flow direction and multiple flow direction algorithms

Flooding is a widespread natural disaster, the occurrence of which is trending upwards due to climate change, population growth and urbanization. Flooding incurs significant monetary costs due to its damage to infrastructure and endangers human lives (Formetta & Feyen 2019). From the different types of flood events, flash floods are particularly deadly due to their rapid and unpredictable onset (Jonkman 2003).

As a response to these threats, flood risks must be properly determined, evaluated, and communicated, followed by the implementation of well-informed mitigation measures. These are crucial steps in pluvial flood risk management (PFRM) (DWA 2016). They can be applied at an urban (local) scale as well as implemented at the regional level. Regional-scale mitigation plans provide a more comprehensive analysis and can be more effective when determining the optimal placement of mitigation measures (Webber et al. 2021; Hamers et al. 2023). They are, however, more demanding in terms of preparation and computation. Therefore, it is necessary to provide decision-makers with the ability to rapidly determine risks and test different mitigation plans.

There is currently a wide range of tools available for urban inundation modeling that can be categorized according to their computational requirements. On the higher end of the spectrum there are methods solving the 2D Shallow-Water Equations (SWES), which can be additionally coupled with a 1D drainage network model. The resulting model output provides high quality results, albeit at a high computational cost. A simulation typically requires a solution of the SWE for each timestep to obtain water depths and velocities (Teng et al. 2017).

Cellular Automata (CA) flood models provide another alternative, whereby the model domain is conceptualized as a grid in which each cell has a state and interacts with grid neighbors according to predefined transition rules. For example, Dottori & Todini (2010) and Jahanbazi et al. (2017) both developed CA methods solving the diffusive wave equations over multiple timesteps. Guidolin et al. (2016) and Jamali et al. (2019) took alternative approaches and defined simpler transition rules rather than using the SWE. The advantages of this method class are the reduced computational effort and improved ease-of-use.

Data-driven flood modeling is an emerging field that relies on machine learning algorithms that are trained on data from flood inventories to quickly produce flood depth and risk maps. The performance of a model, however, is limited by the amount of training data available. In addition, the model transferability to other areas outside the training domain is not reliable and continues to be a major challenge (Bentivoglio et al. 2022; Seleem et al. 2022, 2023).

A category of methods, referred to from here on as Static Flood Analysis (SFA), can be used to model flooding with more simplified approaches. They can be understood as a 1D flood routing technique by which surface excess water is routed along a network of topographic depressions connected by flowpaths. Following the method deviced by Djordjević et al. (1999), sink areas are first extracted from a Digital Elevation Model (DEM) along with their attributes, such as their elevation extents, as well as their storage volume curve and overflow point. The sinks are then linked by a flowpath network that is derived from flow accumulation (FA) results. Notable examples from this modeling domain include Maksimović et al. (2009), Jahanbazi & Egger (2014), Balstrøm & Crawford (2018) and Samela et al. (2020).

FA algorithms have been used for decades since their inception by O'Callaghan & Mark (1984). Broadly, they utilize the DEM topography to determine the number of cells that drain into each cell. Since then, additional FA algorithms have been developed with the aim of improving accuracy and performance, as well as extending functionality. FA methods can be roughly subdivided into two categories: Single Flow Direction (SFD) and Multiple Flow Direction (MFD).

SFD methods will only drain the flow from each cell into one of its neighbors, namely the one possessing the steepest downslope gradient. The D8 method by O'Callaghan & Mark (1984) is the most prevalent among the SFD algorithms. Since then, improvements aimed toward lowering its time complexity have been carried out by: Yao & Shi (2015), Wang et al. (2011) and Jiang et al. (2013). More recently, Zhou et al. (2019) developed an SFD algorithm that yielded the fastest results when compared to a variety of its counterparts.

MFD methods are formulated to distribute the flow of each cell to all of its downstream neighbors. Freeman (1991) first calculated divergent flow by applying a flow partition coefficient for each downstream flow direction. Successive improvements have been made by Quinn et al. (1991), Quinn et al. (1995), and Qin et al. (2007).

A selection of methods to model flash floods have been presented ranging from advanced to simple in terms of their implementation, practical use and performance. Following these assessments, there is a clear need for tools that allow for rapid computation of large, regional-scale models. In addition, modeling outputs must provide detailed quantitative information in a variety of formats that is relevant for PFRM. The aim of this study is to develop a deterministic SFA toolset that will build upon the current state-of-the art methods described previously. The toolset will take the following inputs:

• A raster layer DEM.

• The total rain amount [mm] for the event.

• A vector layer containing building polygon features with height attributes (optional).

• A vector layer containing line features representing known streams (optional).

• A vector layer containing polygon features representing the landuse type with runoff coefficient attributes (optional).

Also, be capable of producing the following outputs:

• A raster layer representing the overland FA volumes.

• A topographic sinks vector layer enriched with hydrological attributes.

• A vector layer representing the flood extent within each sink.

• A raster layer representing the flood depths within each sink.

The proposed method will implement a dynamic coupling procedure wherein sink and flowpath interactions are solved via an exchange of volumes. The effects of sink storage will also be reflected in the FA results. Besides using an SFD algorithm, this study will also develop a new variant of the MFD algorithm and compare any differences between these two approaches in terms of accuracy and performance.

Programming was done in C + +. The following Free and Open Source Software (FOSS) was used to implement the following methods described in this paper:

• (GDAL/OGR) Geospatial Data Abstraction Library, version 3.6.1. (GDAL/OGR contributors 2023).

• QGIS Geographic Information System, version 3.18.3. (QGIS Development Team 2023).

• WhiteboxTools (WBT), version v1.5.0. (Lindsay 2014).

The proposed SFA method can be roughly divided into three consecutive steps as shown in the graphical abstract. The terrain and sink preparation form part of the model pre-processing step. The main simulation procedure is the coupled analysis where FA and bidirectional volume exchange with sinks takes place. The sink enrichment and flood depths/extents generation are the final procedures for output generation.

Due to its computational efficiency, the SFD algorithm by Zhou et al. (2019) was chosen to compute FA and was modified to consider volume exchanges between sink objects and FA volumes. Supplementary material, Algorithm 1 provides pseudocode for the developed method.

The flow direction grid codifying the cardinal direction of flow according to the steepest descent gradient around each cell is calculated from the filled DEM by using the WBT D8Pointer tool. It is loaded into memory as an array of byte-sized elements containing one of eight possible flow direction values found within the set ⁠. A value of zero is considered invalid and represents no flow direction.

The steps for the calculation of NIDP values are provided in the supplementary material (Supplementary material, Algorithm 2). nextCell is a function returning a Boolean value. For an input cell that drains into another cell, the function returns a true value and the cell number is updated to the next cell number. If the cell has no neighbors to drain into or if the neighbor is out of bounds, the function returns a false value.

The SFD algorithm (1) was extended to enable MFD computations. The pseudocode is provided in the Supplementary material, Algorithm 3. The flow direction variable must be capable of encoding all possible directions for each cell. To achieve this, an 8-bit unsigned integer is utilized as described by Ortega & Rueda (2010), wherein each bit represents a flow direction. The range of valid flow directions is thus contained within the set ⁠.

NIDP values are assigned to each raster cell following the steps in Supplementary material, Algorithm 4. The nextCell function was modified to accept an additional flow direction parameter.

Once the FA procedure is finished, each sink is enriched with the following hydrological attributes listed in Table 1. 

The attributes listed in the first row are obtained during the sink preparation step. The other attributes are derived from the coupled FA results. The flow ratio of each sink is the inflow divided by the sink volume. The flood elevation is linearly interpolated according to the inflow and the stage curve of the sink.

The extent of flooding is extracted from the DEM by filtering all raster elevation cells that are higher than the flood elevation attribute of each sink and constructing polygon features to represent the flooding extent. In a similar fashion, rasterized flood depths within each sink are calculated by subtracting the DEM elevation from each cell within a flooded sink area from the sink flood elevation attribute.

The DEM data used in this study are derived from ©GeoBasis-DE/LGB, dl-de/by-2-0 (https://data.geobasis-bb.de/geobasis/daten/dgm/xyz). Landuse, streams, buildings, and map data were derived from OpenStreetMap (OSM) data distributed under the Open Database License (ODbL).

All computations were carried out on a Windows 11 machine running on a 12th Gen Intel(R) Core(TM) i7-1260P processor with a clock rate of 2.10 GHz and 16 GB of RAM.

The FA volumes obtained from the SFD show a network of clearly defined flowpaths, while those obtained from the MFD bear a closer resemblance to a network of interconnected ponds. Upon visual inspection, it was noted that the T2D results also show flooding outside of sink areas, while the flood depth results from the SFA methods do not.

The statistical comparison results from Equations (5) and (6) are shown in Table 3. 

Both methods achieved similar results, with the MFD performing marginally better than the SFD method in the NSE and MCC metrics, albeit at a slightly higher computational cost. Both models were each computed in under 1 s while showing a high level of agreement with the hydrodynamic results. In contrast, the mesh generation for the T2D model required 65 s and the parallel simulation on 12 CPU cores finished in 1 h 18 min. The lower MCC values suggest that a significant portion of incorrectly categorized cells do not pose much statistical significance because the flood depth differences at these cells are relatively low.

A similar SFA method developed and tested by Samela et al. (2020) also showed good accuracy for events with higher recurrence periods. Their study was conducted for a small (⁠4 km²) and large (⁠72 km²) model and computed in 99 and 322 s, respectively. They showed increasing precipitation events yielded higher accuracy with a maximum MCC value of 0.74 for a 130 mm precipitation scenario. This assessment also coheres with the results of a model benchmarking series for a similar SFA method carried out by Seleem et al. (2021), achieving MCC values similar to those obtained in this study (⁠ 0.7). However, results were shown to be model-dependent and the flood prediction accuracy can decrease for events where flooding begins to spread beyond sink features. This could lead to an underreporting of flooding extents for unusually extreme precipitation events.

A more sophisticated CA method with transition rules developed by Jamali et al. (2019) was tested on a pluvial flood model of similar size (0.3 km²) and yielded a maximum MCC of 0.86 when validated against a hydrodynamic model. The increased accuracy came at a higher computational expense, as the model ran in 128 s without parallelization, while the SFA calculation in this study was nearly instantaneous.

The regional-scale case computation yielded hydrological data for a total of 44.611 topographic sinks comprising a total area of 47.7 km². The difference of flood extent area was 6,769 m² between the SFD method and the MFD method, which showed flooded areas totaling 15.639 and 15.646 km², respectively. The model finished execution in 428 s for the SFD method, and 688 s for the MFD method. The FA volumes as well as the topographic sink and flood extent polygons are shown in Figure 12 for the SFD coupled method and in Figure 13 for the MFD coupled method.

As with the validation case, the differences between the SFD and MFD variants were barely perceptible. It is thus unlikely one method provides a clear advantage in terms of flood extent calculation over the other. After evaluating the FA raster results, however, it can be concluded that flow dispersion effects and flow exchange between flowpaths is better represented in the MFD models compared to the SFD method. Areas with higher flow dispersion are generally considered to have a higher propensity toward flooding. MFD results can thus better enable users to visually extract more information regarding flow patterns in the model. These assessments are highly subjective, however. Numerical computation of FA volumes in a hydrodynamic simulation (e.g. TELEMAC-2D) should be carried out in future studies to better compare FA results outside topographic sinks. This could help to further highlight differences between the SFD and MFD methods using a quantitative approach.

In this study, an SFA method capable of rapidly generating flood inundation maps for local and regional-scale models was developed, tested and validated for a small urban flash flood model. Finally, a regional-scale model was computed using the proposed methods. The SFA method computes FA volumes while accounting for depression storage in the terrain by coupling topographic sink objects with terrain flowpaths. Two variants of the method were proposed: a method utilizing an SFD FA algorithm by Zhou et al. (2019) and another utilizing an MFD algorithm. To achieve this, a novel extension of the SFD algorithm by Zhou et al. (2019) was developed that enables MFD computations.

As a reference for the evaluation, the water depth results from a hydrodynamic model (TELEMAC-2D) were used. The obtained NSE index and MCC values were 0.8 and 0.7, respectively. The proposed methods were thus shown to effectively predict the extents and depths of flooding. The model computation was nearly instantaneous and clearly outperforms the hydrodynamic method in this regard. According to both metrics and a visual analysis of the inundation maps, the differences between the flood depths and extents between SFD and MFD methods were barely perceptible. Overall, the MFD method yielded slightly larger flood areas in total for both the validation and regional-scale cases. Though the MFD method slightly outperformed the SFD method, this difference cannot be considered to be statistically significant.

While the comparison of inundation depths between the TELEMAC-2D and both proposed SFA methods showed good agreement, a wider variety of models should be validated for rainfall events of different intensities. The authors of this study propose conducting a benchmarking series for the proposed SFA methods. The series could determine result differences based on model size, rain event, landuse as well as explore the possibility of calibrating SFA models of the MFD variant with the upper and lower exponent bounds (⁠ and ⁠) from Equation (4).

The main limitations of the proposed methods are as follows:

• The methods are inherently incapable of detecting flood areas outside topographic sinks.

• The temporal evolution of field parameters is not considered. Events with the same precipitation amount but different durations will always produce identical results, which may be incorrect.

• According to Samela et al. (2020) and Seleem et al. (2021), the accuracy of SFA methods is significantly impacted by the amount of precipitation as well as the retention capacity of sinks in a model. Flood area classification generally improves with increasing precipitation amounts. However, accuracy can also decline for events where flooding spreads beyond the sink extents. These are issues yet to be addressed.

The proposed toolset could benefit from further additions:

• The runoff depths outside of sink areas could be estimated using a kinematic wave approximation provided by Overton & Meadows (1976). This could help increase the accuracy of the SFA method while counteracting one of its main drawbacks.

• The consideration of infiltration could be extended to include other methods such as Green-Ampt or Curve Numbers (CN), such as implemented by Samela et al. (2020) and Wang & Jin (2001), respectively.

• The use of spatially distributed rainfall could be enabled following the work done by Samela et al. (2020). This would also allow the methods to be applied for real-time flood predictions based on radar precipitation forecasts.

• Sediment and/or tracer transport could also be implemented with minimal effort.

In conclusion, the SFA method developed provides users with the ability to rapidly generate flood maps for risk assessment and aiding decision-making during urban planning. The enriched sink attributes provide useful metrics for the categorization of flood risks. They also provide concrete values that can be used by planners to incorporate sinks into mitigation measures by modifying their depths and/or converting sinks into measures such as swales, trenches or retention basins.

This work was funded by Ingenieurgesellschaft Prof. Dr Sieker mbH. The authors would like to thank all experts for their contributions in the field of rainwater resources management.

All relevant data are included in the paper or its Supplementary Information.

The research presented in this article was financially supported by Ingenieurgesellschaft Prof. Dr Sieker mbH.