In the United States and around the world, flood maps are used to communicate flood risk, set insurance rates, and determine infrastructure policy. The extent and depth of flooding are commonly simulated by using either physics-based hydrodynamic models or a more simplistic static approach. In addition to models, empirical methods delineate historical floods using observations such as on-ground measurements, aerial photographs, and satellite imagery. However, the limited availability of observed flooding data and their applicability solely to past events are significant limitations of this approach. Empirical methods are commonly used to calibrate and validate hydrodynamic models. A comprehensive description and comparison of the three methods can be found in Teng et al. (2017).
Physics-based hydrodynamic models allow for an accurate simulation of flooding since they include physical parameters such as wind speed, atmospheric pressure, bottom friction, etc. Also, hydrodynamic models can predict the evolution of a flood event, providing a temporal component. The use of these models requires expert skills, specialized software, and high-spec computers that often limit their use to specialized engineering sectors. Furthermore, due to their computational demand, they are not suitable for simulating an ensemble of multiple flooding conditions such as sea-level rise projections and the corresponding uncertainties, storm surge scenarios, and digital elevation model realizations (Amante, 2019;Barnard et al., 2019). For the same reason, using high-resolution elevation data in hydrodynamic models is often limited to small-scale study areas. Coarse spatial resolution is used in regional to larger spatial scales, which can lead to inaccuracy in flood maps. A coarse resolution tends to average out topo-bathymetry information, such that some regions are incorrectly modeled as wet or dry; this is particularly a problem when coarse resolution prevents hydraulic connectivity to a flood-prone region. Previous research has emphasized the significance of resolution in achieving accurate simulations of inundation events (Williams and Lück-Vogel, 2022). Nevertheless, the demand for high spatial scale in inundation modeling poses a substantial challenge for hydrodynamic models when applied to relatively large areas.
An alternative approach for simulating flooding is the use of machine learning (ML) algorithms. In recent years, the use of ML has been shown promising in the field of flood mapping. ML algorithms have been mainly used to delineate the flood during or after its occurrence rather than to predict/simulate the flood extent (Bentivoglio et al., 2022). They face two significant drawbacks: the requirement for extensive training data (recorded flooded events or physics-based hydrodynamic model simulations) and, while they can exhibit computational efficiency with small datasets, they can also prove computationally demanding. The computational intensity of ML hinges on several factors, including the size of the training dataset, hyperparameter tuning (which often involves numerous iterations), and the model's architectural complexity, such as the number of hidden layers and neurons employed. In addition, the multitude of available ML algorithms often leads to the need for ensemble techniques, where their outputs are combined to enhance predictive accuracy, escalating the computational demand.
An alternative to computationally costly hydrodynamic models and the substantial data demands ML algorithms is the use of static approaches, also known as the "bathtub", in which all regions below a flood water level are considered to be flooded. Overland inundation estimates through static methods have traditionally been based on a GIS (geographic information systems) framework (e.g., Amante, 2019;Li et al., 2009;Torresan et al., 2012;Yunus et al., 2016), often processed in ESRI's ArcGIS (e.g., Breilh et al., 2013;Maloney and Preston, 2014;Patrick et al., 2015;Perini et al., 2016;Seenath et al., 2016;Williams and Lück-Vogel, 2020;Zachry et al., 2015) and to a lesser extent in QGIS software (e.g., de Lima et al., 2021;Paulik et al., 2021;de Leo et al., 2022). Other static models include the Height Above Nearest Drainage (HAND, Nobre et al., 2011), Teng Vaze Dutta (TVD, Teng et al., 2013), and Floodwater Depth Estimation Tool (FwDET, Cohen et al., 2018) models. HAND is a drainage normalized version of a digital elevation model (DEM), where flooding extent and depth are identified as those cells that belong to a mutually connecting flow path. The model needs a water level and a DEM as input to calculate the flood depth and extent. In contrast, TVD and FwDET can simulate flood depth but require an input for the flood extent. Teng et al. (2022) offers a comprehensive comparison among these three static approaches and compare the results against hydrodynamic models.
Static approaches are easy to use and computationally very efficient in comparison to hydrodynamic models and machine learning algorithms. Static models allow the rapid simulation of a large set of realizations and they necessitate minimal input data for implementation, mainly terrain elevation data, such as DEMs, and a flood water level. However, results from static models are expected to be less accurate than hydrodynamic models because static methods do not account for the physical processes involved in an inundation event, such as bottom friction, which effectively reduces the flooding depth and extent. Therefore, static methods typically overestimate the inundation (Ramirez et al., 2016;Teng et al., 2022). Similarly, static methods that disregard hydrological connectivity are anticipated to exacerbate flood overestimations. This is due to the fact that areas below the flood water level, which are not connected to the sea or river, may erroneously be identified as flooded regions.
Another drawback of static approaches is that the flood water level is typically assumed to be uniform over the study area, however, water level is spatially variable. Some exceptions exist in the literature, where a spatially varying water level can be applied in static methods, such as in the HAND model. There is an ongoing discussion in the literature on the use of static flooding methods. Some authors (Ramirez et al., 2016;Seenath et al., 2016) suggested that static models should be avoided when computational speed is not of primary concern. Others have compared flood extents and depths resulting from physics-based models versus static approaches, finding relatively good agreement (Orton et al., 2015a). Teng et al. (2017) argues that the disparity between static approaches and hydrodynamic models varies depending on the topographical complexity. Increased topographical intricacy results in a greater level mismatch. On the other hand, hydrodynamic and static models exhibit good agreement in clear flow paths.
Here, we implement a static model to determine the flood extent and depth in MATLAB (MatFlood). MatFlood takes into account hydrological connectivity and, acknowledging the shortcomings of static methods, includes the attenuation of the inundation through a reduction factor that accounts for the distance to the source of flooding. The reduction factor can be spatially uniform (i.e., constant across the study area) but also spatially varying, using an optimal interpolation method. Few studies have previously applied static methods using reduction factor, but it has been assumed to be constant across the study area (e.g., Williams and Lück-Vogel 2020;Ward et al., 2020). In addition, Mat-Flood allows for the calculation of a spatially varying flood water level, based on gauge's observations or model points distributed across the study area. The simplicity, minimized data requirements, and computational efficiency of MatFlood set it apart from hydrodynamic models and ML algorithms. These characteristics render MatFlood well-suited for generating an ensemble of flood simulations and estimating the extent of inundation under diverse conditions, including different mean sea level scenarios, uncertainties, extreme events, and DEM realizations. In addition, MatFlood could be used in flooding emergency systems by facilitating the generation of high-resolution flooding maps over extensive regions in near real-time.
Conceptually, MatFlood shares similarities with other static approaches, such as the HAND model, that simulate flooding by accounting for hydrologic connectivity. In addition, the HAND model accommodates spatially varying flood water level, although does not include a reduction factor. MatFlood extends the literature by offering a MATLAB-based algorithm that incorporates optional features (spatially varying reduction factors and flood water levels). MatFlood's high computational efficiency allows for the generation of multiple flood maps using different inputs or scenarios, as well as the use of highresolution DEMs.
In the following sections, we describe the algorithm. We then apply MatFlood to assess the flood extent and depth of Hurricane Sandy (2012) in the New York/New Jersey area to illustrate its use. Results are compared against observed high water marks (HWMs) and a hydrodynamic model.
MatFlood comprises two main components. The first component maps the flooding depth based on the water level during an inundation event. The second component, which is optional, employs a reduction factor to attenuate the extent and depth of flooding. A visual workflow of MatFlood is provided in Fig. 1. MatFlood requires three primary inputs:
(1) A digital land elevation model (DEM) in a raster format, including the elevation (z) and the coordinates. DEMs include Digital Surface Models (DSMs) or Digital Terrain Models (DTMs). Previous analyses have demonstrated that the choice between using DSM or DTM has a considerable impact on the accuracy of simulated flood scenarios (Williams and Lück-Vogel, 2020). As previously stated, utilizing a high-resolution DEM is advisable to obtain more accurate results. However, it is important to highlight that MatFlood is compatible with both DSM and DTM, and it can accommodate various resolutions of DEM datasets. (2) A flood water level of interest. The flood water level refers to the vertical water level producing the inundation and should be based on the same geodetic reference frame vertical datum as the DEM; and (3) A single point which identifies the main water body from which flooding occurs, hereinafter referred as lp (location point). We use the term "main water body" because MatFlood can be applied to any interface of water-land (lakes, rivers, and the ocean). In the following, we will refer to the main water body as "sea", for simplicity. The location point, lp, is a single longitude and latitude point, indicated by the user, and located within the sea (the sea polygon is referred as the principal water polygon, PPW1 in Fig. 2).
In the following sections we introduce the algorithm for simulating inundation (Section 2.1), the description of the additional features to obtain spatially varying flood water level (Section 2.2), and to attenuate the flood depth and extent by applying a reduction factor (Section 2.3). The methodology is illustrated by considering a basic topography z of land/sea (Fig. 2), and a flood water level of 1 m. The horizontal resolution is 1 square meter (each cell in Fig. 2). Within the sea (the principal water polygon, PPW1), there are two polygons that represent land areas of different elevations (children polygon in water 1 and 2, CPW1 and CPW2). In the land area (PPL1), there are four embedded areas of lower elevations (CPL1 to CPL4).
A first estimation of the flooded area is obtained by subtracting the flood water level (fwl) from the land elevation (z), obtaining a flooded elevation grid (z F ) of size Y x.
Since both the flood water level (fwl) and the DEM (z) are referred to the same vertical datum, positive values in z F represent land/dry areas while negative and zero values indicate water. In z F , elevations equal or lower than the flood water level are inundated. In the example, the sea polygon now includes PPW1, CPL2, CPL3 and CPL4 (Fig. 3). In this step, hydrological connectivity has not been considered yet and therefore, z F also includes flooded areas not connected to the sea (CPL1 and CPW2 in Fig. 3). In order to account for hydro-connectivity, the algorithm must pinpoint and separate the sea polygon. To do so, z F is converted into a binary black and white matrix of ones (water) and zeros (land), z BW :
Producing a binary matrix allows us to identify each polygon of water and land, regardless of the elevation. To isolate the polygons in z BW , we modify an existing MATLAB image processing function, bwboundaries (Gonzalez et al., 2004), which identifies the outermost objects ("parents") as well as the successive polygons enclosed by them ("children"). The modified function, bwboundaries_SFM, enables for a faster way to delineate all polygons in z BW . We use bwboundaries_SFM accounting for 8-pixel connectivity, i.e., pixels are connected if their edges or corners touch either along the horizontal, vertical or diagonal direction.
The output, P, is a Y x matrix containing an identification number (id) for each polygon detected (Fig. 4). The polygon including the sea and the flooded areas, P mwb,Fl , is identified by the algorithm since it embeds lp. The depth of the inundated areas is calculated as the difference between the topography (z) and the flood water level (fwl) for those cells within P mwb,Fl:
where z FHC (Flooded topography accounting for Hydrological Connectivity) is a Y x matrix where only areas below the flood water level and hydrologically connected to the sea are flooded (Fig. 5).
Typically, in static flood modelling, the flood water level is assumed to be spatially uniform over the entire domain (as done in Section 2.1). By doing so, the variability of the water level along the shoreline is neglected, potentially leading to overestimate or underestimate flooding in some areas. Consideration of varying spatial water levels becomes particularly pertinent in the simulation of flooding across extensive areas, where larger water level variability is anticipated. Conversely, smaller areas along the coast may exhibit less pronounced variations in water levels, making the distinction between employing a spatially uniform versus varying flood water level less significant.
The spatial variability of the water levels can be captured by a set of observed/modeled water levels scattered over the study area. For simplicity, we refer to these as monitoring stations; W N i=1 = {w i ,…,w N }, being N the number of monitoring stations where water level information is available. Following Agulles et al. (2020), we apply an optimal interpolation method to interpolate the water level (W) information at the monitoring stations over the entire study area. The resulting spatially varying flood water level ( Ûwl) consists of a linear combination of the water levels at the monitoring stations (W) and a background field (i.e., first guess), with weights determined from the covariances of W and the background. In our case, the background is determined as the average between the water levels contained in W. In other words, Ûwl is a 2-dimensional grid containing water level anomalies above/below the average, where cells tend to values similar to the closest monitoring station. Cells far from the monitoring station tend to the average of W. Fig. 6 shows the spatially varying flood water level for the illustrating example. Note that the topography, and therefore the flood path, is not considered when calculating the spatially varying water level. The original optimal interpolation algorithm formulation can be found elsewhere (e.g., see Bretherton et al., 1976;Daley, 1993); here, only the major features of the algorithm are summarized.
Ûwl, of size Y x, is obtained by:
where BK is a Y x matrix containing the background information. S is a M x N matrix containing the covariance of the Euclidean distances between the digital elevation model and the monitoring stations (W), which are then modeled using a Gaussian function: Fig. 4. Polygons (parents and children) found in z BW . Five different polygons were identified: parent polygons are id = 1, 2. Inside the parent polygons, two children polygons (id = 3, 5) are located within the sea polygon, P mwb,Fl , and one children polygon is in the main land area (id = 4). Red dot shows the location point, lp. PS = e -d1 2 ij/ 2L 2 (6) d1 ij is the distance between each j point of the digital elevation model and each monitoring station contained in W.
L is the spatial correlation length scale. Smaller values of L result in a lower spatial variability of flood water levels, constraining the influence of each individual water level value to the vicinity of the corresponding monitoring station. Simultaneously, smaller values of L better preserve the water level magnitude at the monitoring station. In contrast, higher values of L generate a broader spatial footprint for each individual flood water level. However, the water levels at the monitoring stations become slightly smoothed out when using higher values of L. The objective of using the optimal interpolation method is to simulate a spatially varying flood water level that accurately represents the maximum variability across the study area while still preserving the flood water level data at the monitoring stations. Therefore, the algorithm estimates the optimal value of L as the maximum value of the ratio between the standard deviation of the interpolated flood water level (STD( Ûwl)) and the mean difference between the interpolated grid water level at the monitoring stations and the flood water level recorded by the monitoring stations ( Ûwl w -W):
Values of L from 0.5 to 100 with 0.5 increments are tested. D* is a N x N matrix containing the covariance of the Euclidean distances across W. Again, D* is modeled using a Gaussian function:
In this case, d2 ij is the distance between the monitoring stations i and j. We assume that the water levels (from simulations or observations) are perfect; however, water level gauge data and model outputs each have uncertainties. We include an error to the covariance matrix D*. ε is a N x N diagonal matrix containing an observational error of 0.01. The covariance matrix D* is then modified as: 5), is a N x 1 vector containing the anomalies between the water levels at the monitoring stations (W) and the background:
The area flooded by a spatially varying water level (z FHC,V ) is obtained by applying the method outlined in Section 2.1 using the spatially varying water level instead of a uniform flood water level.
Following the example from Section 2.1, Fig. 6 shows the resulting spatially varying flood water level obtained from five different monitoring stations, in meters, W = {0.2,0.7,1,0.7,0.7,0.5}. Fig. 7 shows the resulting flood map when using the spatially varying flood water level (z FHC,V ). Note, both flood depth and extent are reduced in comparison to using a uniform flood water level (z FHC ). Compared to Fig. 5, the use of a spatially variable water surface has prevented flooding of the CPL3 region and therefore, CPL4. Also, closer cells to the W = 0.2 m (Fig. 6), are not inundated, since those grid cells have an elevation of 0.5 m.
The static approach neglects the hydrodynamics of flooding, and effectively assumes that all potentially flooded regions will be instantaneously flooded when a flood reaches its peak water level. However, a combination of hydraulic controls, form drag, and hydraulic roughness impedes flood currents, and imposes a timescale for flooding to occur at each point. Thus, a static approach can greatly overestimate the extent and the depth of flooding, as mentioned in the Introduction.
To address the expected overestimation of the inundation, and following previous studies (Ward et al., 2020;Williams and Lück-Vogel, 2020), MatFlood incorporates a reduction factor that diminishes both the flood depth (representing the amount of water in the vertical axis) and the flood extent (referring to the flooded area). The attenuation of the inundation depends on the distance to the coast and a user-defined rate, so it can be adjusted to the particularities of the study area. The reduction of the flood is applied as a post-processing step after generating the flood map (Fig. 1). It should be noted that the use of the reduction factor is optional. Additionally, the reduction factor can be either spatially uniform or spatially varying, depending on the specific requirements of the analysis. Detailed explanations are provided in the subsequent sections.
The methodology consists of (a) delineating the shoreline, (b) calculate the distance of each point in the elevation data (DEM) to the shoreline and, (c) applying the user-defined rate. In this study, the shoreline refers to the intersection between the sea and land under calm conditions, representing the pre-flooding state.
To delineate the shoreline, we use the MATLAB function edge, which takes a binary image as input and returns a binary image of the same size, with 1's where the edge is found and 0's elsewhere. Note, the shoreline is identified using the initial elevation data (z) and not the flood map. The input binary matrix (z BW0 ) is defined to take values of 1 for the cells located in the sea and 0 otherwise. In order to obtain z BW0 , Equations ( 1)-( 3) are applied to z using a uniform flood water level of 0 m. Then, Equation ( 4) is modified to:
Then edge is applied to z BW0 :
SH contains the coordinates of the shoreline. The user-defined flood reduction rate (RF) is defined as the vertical water depth reduced (v rwd ) over a horizontal distance (h d ):
This rate is then applied to the flooded map (z FHC ) as:
D* contains the covariance of the Euclidean distances between each point in the computational grid and the shoreline. This calculation can be computationally expensive in large areas or when using high spatial resolution. To improve the computational demand, the distances are calculated over a reduced-resolution computational grid and then interpolated back again to the original resolution. Note that only the distances are calculated over a reduced-resolution grid, the flooded map retains its original resolution so the accuracy of the flooding maps is not impacted. Areas that are no longer flooded due to the application of the reduction factor take values of the initial elevation data (z):
Therefore, based on Equation ( 14), the amount of flooding attenuated by the reduction factor depends on the distance to the coast and the flood depth, which is ultimately determined by topography elevation. In some instances, the hydraulic connectivity might not hold after applying the reduction factor. Hydrological connectivity is enforced by running the static flooding algorithm again (Section 2.1) using a flood water level of 0 m. By doing so, we obtain the flood depth and extent attenuated by the reduction factor and accounting for hydro-connectivity (noted as z FHC,R ). Note, in this case we have applied the reduction factor to the map flooded by a uniform flood water level (z FHC ), since we have showed that the spatially varying flood water level already reduces the extent and depth of the inundation. However, the map flooded by a spatially varying flood water level can be used instead.
Following the optimal interpolation method used to obtain spatially varying flood water levels (Section 2.2), the algorithm allows for the calculation of a spatially varying reduction factor given a set of reduction factors distributed over the study area. The approach follows the same steps as in Section 2.2 using a set of user-defined flood reduction rates instead of flood water levels.
For demonstration proposes, we have employed a rate of 0.0417 to the previously simulated flooding area. Note, this value is chosen to illustrate the methodology but it has no physical meaning. In this example, 0.0417 reduces the flood depth by 0.5 m every 12 m in the horizontal. Fig. 8 shows the results after accounting for water connectivity. The flood depth and extent have been attenuated as can be seen when comparing the flood map before (Fig. 5) and after (Fig. 8) applying the reduction factor: CPL3 and CPL4 are not flooded anymore. In addition, the flood depth was reduced by 0.5 m in the landward area in CPL2.
We use MatFlood to assess the flood depth and extent of Superstorm Sandy, which caused an extensive storm surge in the states of New York and New Jersey on October 30th, 2012, resulting $50 billion in damages and 147 direct deaths (NWS, 2013). The results of MatFlood are compared to observed high water marks (HWMs) and to a physics-based hydrodynamic model (referred as PBM for simplicity).
MatFlood's flooding simulation is performed using the Continuously Updated Digital Elevation Model elevation data, referenced to NAVD88 (NCEI, 2022). The spatial resolution of the elevation data is 1/9 of an arc-second (~3.43 m). The study area covers the southwestern part of the state of New York and the eastern part of New Jersey (Fig. 8). The elevation data (and therefore computational grid, see Section 2.1), contains 24,312 (Y) × 16,212 (X) cells. Maximum water levels (above NAVD88) from 29/October 30, 2012 are retrieved from 6 tide gauges (monitoring stations) from NOAA and the USGS: Bergen Point, Jamaica Bay at Inwood, Kings Point, Newark Bay, Rockaway Inlet, and The Battery.
A total of 312 HWMs are retrieved from the U.S. Geological Survey Short-Term Network Data Portal (USGS, 2022). HWMs consist of flood depth measurements (above NAVD88) based on, for instance, debris and mud washed ashore, stain lines, and seed lines. Generally, HWMs are situated at the inland boundary of the flood. However, this may not be the case for all observations.
The PBM consisted in the coupled ADvanced CIRCulation (ADCIRC)/ Simulating Waves Nearshore (SWAN) models, which include the effect of storm tides and significant wave heights. The model used the US FEMA Region II operational unstructured numerical grid, where the smallest distance between nodes was 70 m. It should be noted that the ADCIRC model was edited to "close the tide gates" in order to improve the results. The simulations used wind and atmospheric pressure reanalysis from OceanWeather, Inc. Details on the numerical modeling can be found in Brandon et al. (2016) and the Superstorm Sandy simulation is described in Orton et al. (2015b). In order to compare the PBM with MatFlood, the flood depths from the PBM where gridded and interpolated to the MatFlood resolution (we recall that the PBM was run on an unstructured computational grid).
Fig. 9 shows the study area, including the topography, the locations of the tide gauges and HWMs. Water level observations reveal limited variability (Fig. 10): three tide gauges, located on the western study area, show higher than average water level (from 3.4 to 3.6 m) while the other three tide gauges, on the eastern, show lower water level (from 3.08 to 3.2 m).
We compare the area flooded under different model setups: (1) Neglecting hydrological connectivity. In this case, the inundation is calculated six times, one for each tide gauge water level. For each of this simulations, the flood water level is spatially uniform. (2) Accounting for water connectivity. As in 1), the inundation is calculated six times, one for each tide gauge water level and the flood water level is spatially uniform. (3) MatFlood using a spatially varying flood water level. Observations from the six tide gauges are utilized to obtain the varying flood water level, as in Section 2.2. The resulting spatially flood water level is shown in Fig. 10. (4) Applying the reduction factor: the inundation resulting from 3) is reduced by a 3.75 × 10 -4 rate.
(5) Applying a spatially varying reduction factor: the inundation resulting from 3) is reduced by a rate spatially varying from 0 to 3.75 × 10 -4 .
MatFlood results are compared against the PBM and HWMs using three metrics: the percentage of coverage, the root mean square error (RMSE) and the critical success index (CSI). The percentage of coverage indicates the number of HWMs that fall within the simulated flood area, out of the total of the 312 HWMs. The root mean square error (RMSE) shows the mean difference in flood depth between the HWMs and the simulated area, in meters. Finally, the critical success index shows the overall overlap between MatFlood and the PBM, accounting for both under and over prediction being 0 = no skill and 1 = perfect skill. The critical success index is calculated as follow:
Being cw the simulations of wet points, fd the false dry simulations and fw false wet grid simulations.
Table 1 presents the results obtained from comparing MatFlood, HWMs, and the PBM, including the respective flooded areas for each model setup. Setups 1 and 2 provide the mean values derived from six simulations, along with the minimum and maximum values indicated within brackets. The PBM model exhibits a percentage of coverage of 47.32%, accompanied by a root mean square error (RMSE) of 1.09 m. This percentage of coverage is consistently lower than MatFlood (from 77.29 to 78.55%). Also, the RMSE is larger than using MatFlood (from 0.73 to 0.75 m).
As anticipated, neglecting hydrological connectivity leads to larger flooded areas (from 298.24 to 364.80 km 2 ) than when hydrological connectivity is considered (ranging from 279.75 to 354.94 km 2 ). Due to the reduced flooded area, the percentage of coverage also decreases when accounting for hydrological connectivity, from 77.55% to 76.39%. On average, the critical success index improves by 0.01 when hydrological connectivity is considered, and substantial differences are observed when examining the minimum and maximum values (Table 1). By employing spatially varying flood water levels (model setup 3), the flooded area increases in comparison to using a uniform flood water level as well as the percentage of coverage. The critical success index reduces by 0.01 m. The utilization of a spatially uniform reduction factor significantly attenuates the inundation to 287.52 km 2 and improves the RMSE by 0.02 m. The spatially varying reduction factor produces a larger attenuation of the inundation, to 209.35 km 2 . However, this reduction also results in a decrease in the critical success index to 0.45. The spatially varying reduction factor improves the RMSE by 0.01. The flood depth obtained from MatFlood using a spatially varying flood water level and a uniform reduction factor of 3.75 × 10 -4 (model setup 3) is shown in Fig. 11. Fig. 12 compares the flood extent of Mat-Flood (model setup 3) and the PBM. The differences in flood depth between MatFlood, the PBM and the HWMs are shown in Fig. 13. The PBM generally result in higher flood depths than those observed by the HWMs (Fig. 13a), particularly along enclosed areas (Fig. 13d). However flood depths are larger in the PBM at the open coast. The distribution of the differences between MatFlood and observations are more centered toward 0 (Fig. 13b). Similarly to the PBM, MatFlood tends to underestimate the flood depth in comparison with HWMs, especially in the enclosed areas (blue dots in Fig. 13e). For completeness, we include a comparison between MatFlood and the PBM at the HWM's locations (Fig. 13c and f). MatFlood tends to overestimate flood depth in comparison to the PBM (Fig. 13c) although not spatial coherence is found (Fig. 13f).
The main advantage of static methods against numerical models is the computation efficiency. The computational cost of MatFlood is evaluated against a static flooding method performed on a GIS framework, since this is one of the most widely used approach. The GIS simulation is processed in the ArcGIS Pro software. The approach followed in ArcGIS Pro mimics the methodology in MatFlood. Both, GISbased and MatFlood static simulations were performed using the same DEM (Section 3.2) and flood water level. Both models account for hydrological connection, using a uniform varying water level and not applying the reduction factor (further details in how the GIS-model was implemented in ArcGIS Pro are provided in the Supplementary Information). We run the tool using an AMD Ryzen 9 5950X machine with a 16-Core Processor of 3.40 GHz and 128 GB of RAM. When using the most generic set up (i.e., without calculating the spatially varying flood water level and without applying the reduction factor, Section 2.1) the calculation takes 28.1 s for this case study. Obtaining a spatially varying flood water level and applying the reduction factor takes longer: 5.5 and 5.9 min, respectively. Note, the computational cost when applying the reduction factor mainly results from the calculation of the distances to the shoreline (Section 2.3). Thus, once the distances are calculated, several values of reduction factor can be tested with little computational cost (in the order of seconds).
We simulate the flood area by applying a GIS-based inundation routine through ArcGIS Pro 2.0.0. We use a uniform flood water level of 2 m both in ArcGIS and MatFlood. Williams and Lück-Vogel (2020) developed an external ArcGIS tool aimed at reducing the flood area using a roughness coefficient to the flooded area. However, the implementation of this roughness coefficient greatly differs from the reduction factor applied in the present work. Therefore, we limit this analysis to the most common GIS-based setup, which does not include any reduction factor of the inundation. Likewise, obtaining a spatially varying flood water level is not possible in ArcGIS, and is therefore not included in MatFlood for this example. The resulting computational time in ArcGIS Pro 2.0.0 is 3.2 min on average (standard deviation of 0.73 min). Using the same initial conditions, MatFlood takes 26.27 s on average (standard deviation of 0.63 s). The steps taken in ArcGIS Pro 2.0.0 to obtain the flood depth are outlined in the Supplementary Information.
This study presents a rapid and user-friendly MATLAB-based algorithm, MatFlood, designed for mapping flood depth and extent, which improves on the widely used GIS-based static models. MatFlood considers water connectivity and includes a reduction factor, in order to more realistically represent flooding conditions. In addition, MatFlood allows for spatially varying flood water level. This feature might improve the representation of the flood water levels along the coast, particularly in larger areas, where larger variability across water levels is expected. Results show that using a spatially varying flood water level impacts the flood depth and extent in comparison to using a spatially uniform flood water level, as it is typically done when using static approaches to assess flooding. MatFlood was compared with a GIS-based system using ArcGIS Pro 2.0.0. Results showed that MatFlood is computationally nearly an order of magnitude more efficient than traditional GIS algorithms.
After describing the algorithm, we applied MatFlood to assess the inundation of Superstorm Sandy in the New Jersey/New York area. The simulation of flood depths was conducted using various model configurations, all of which demonstrated comparable performance metrics when compared to High Water Marks (HWMs) and a physics-based numerical model (PBM). However, additional analyses are required to assess the model's capabilities since the performance of models can be influenced by the quality of the observations they are compared against, the magnitude of the flood event, the complexity of the topography, and the resolution of the DEM (Teng et al., 2017(Teng et al., , 2022)). Also, the model' performance should be further analyzed in larger study areas, where the