A Spatial Analysis Approach for Urban Flood Occurrence and Flood Impact Based on Geomorphological, Meteorological, and Hydrological Factors

Abstract

Urban flooding can cause significant infrastructure and property damage to cities, loss of human life, disruption of human activities, and other problems and negative consequences on people and the local government administration. The objective of this research work is to investigate the relation between urban flood occurrence and potentially flood-triggering factors. The analysis is performed in the western part of Athens Basin (Attica, Greece), where over the past decades several flood events caused human losses and damages to properties and infrastructure. Flood impact is measured by the number of citizen calls for help to the emergency line of the fire service, while potentially influencing factors are several geomorphological characteristics of the area and hydrometeorological indices regarding storms, which were determined with the aid of GIS techniques. The analysis is based on the investigation on binary logistic regression and generalized linear regression models that are used to build relationships between the potentially flood-influencing factors and the flood occurrence/impact for three events that were selected for reasons of comparison. The entire analysis highlights the variations attributed to the consideration of different factors, events, as well as to the different cell size of the grid used in the analysis. Results indicate that, the binary logistic regression model performed for flood occurrence achieves higher predictability, compared to the ability of the model used to describe flood impact.

1. Introduction

Floods are the most frequent type of natural disasters and they affect the highest number of people globally [1]. Especially in urban areas, flood risk is particularly high [2] and a number of trends suggest that the problem of urban flooding is likely to increase, notably after considering the possibility for climate change to lead to more extreme rainfall. As more people move to the cities, they inevitably turn green areas into impervious areas, increasing urban runoff, and, as more people live in densely populated urban areas, their exposure to flood hazards is increased [3]. Urbanization has also an adverse impact on the urban hydrological processes, such as accelerating runoff flow velocity and increasing peak flow [4,5,6], thus, increasing the urban flood risks.

The improvement of urban flood risk management practices has become a global priority because the world’s population is becoming increasingly urban. As more and more of the population is often situated either on areas with high density buildings, or on flood plains and low-lying coastal areas that are already highly urbanized areas, engineering measurements for flood mitigation are difficult to obtain. Thus, contemporary practices focus on non-engineering measurements, such as flood risk warning, as the main way for managing urban flood risk [2]. However, when trying to predict the susceptibility of different parts of a city to flooding and to determine the vulnerability of individual structures, the complexity of urban hydrology together with the influence of the multiple and continuously changing factors lead to significant difficulties [7,8,9].

New approaches regarding flood impact modeling, mitigation, and management require accurate modeling to estimate inundation extents and other hydraulic parameters in urbanized flood-prone areas, such as water depths, discharge and flow velocity, and the mapping and management incorporate more details regarding various hydro-meteorological and geomorphological conditions that characterize the vulnerable areas. There are also many approaches that arise from the consideration of various flood-related factors; for instance, the damage relationships can be functions of a number of damage-influencing factors [10]. A distinction has been made by Merz et al. [11] between impact parameters, i.e., the characteristics of the flood that causes the damage, and resistance parameters, i.e., the parameters that relate to the resisting object.

The knowledge and study of the factors that control flooding and their intensity may decisively contribute to a more realistic assessment of flood risk, while they provide important information for understanding the flood-producing mechanisms. There is a variety of factors that influence the occurrence and the severity of flood events; these factors can be described by meteorological, geological, geomorphological, hydrological, and land use characteristics of flood-prone areas that influence the spatial distribution of flood occurrence [12,13]. Previous works have highlighted the effect of the meteorological characteristics of storms (e.g., [14,15,16]), the catchment’s geomorphological characteristics (e.g., [17,18,19,20]), and the building components (e.g., [21,22]). Finally, research regarding parameters that affect the spatial distribution of flooding contributes to the determination of areas with increased propensity to experience harm [12] and denotes the most significant factors that are incorporated in a GIS-based multi-criteria analysis regarding flood susceptibility mapping (e.g., [23,24,25]). The relationship between the aforementioned factors and the spatial distribution of floods has been examined using regression models in various studies [12,26,27], all of which contribute to predict flood susceptibility based on the characteristics of known flooded and not-flooded locations. All these approaches utilize GIS techniques to represent the values of flood-influencing factors and for the creation of thematic maps as raster datasets. Finally, several approaches (e.g., [28,29]) take into account factors related to topography, land use distribution, and rainfall–runoff event characteristics (single scenario event) and implement GIS techniques for flood damage mapping; a framework that is particularly important also for large-scale mapping of direct economic flood damage in data-scarce environments.

The current research aims to investigate regression models and, particularly, binary logistic regression and generalized linear regression for flood occurrence and flood impact identification, respectively, based on the consideration of 16 influencing factors that were estimated using GIS techniques. The analysis is performed for the western part of Athens basin, Attica region, which is frequently hit by heavy flooding. The socioeconomic development in this area over the past few years is linked to the expansion of human activities and infrastructure into the floodplain in high-risk areas, leading to consequential effects on properties and infrastructure, even to loss of human life. Three rainfall-induced flood events are analyzed to compare the corresponding results. These events occurred on 22 February 2013 (maximum and average total precipitation depth: 98 and 78 mm), 24 October 2014 (maximum and average total precipitation depth: 49 and 42 mm) and 22 October 2015 (maximum and average total precipitation depth: 71 and 50 mm). The latter was a more localized event, the event of 2014 is linked to a storm with a total maximum rainfall depth of about 130 mm in five hours in the northwestern part of Athens that caused 1223 damage locations, while the first one affected the entire Athens basin (Figure 1). The proposed methods have not applied before in a comparative way for a specific region and this analysis aspires to contribute to flood susceptibility mapping over urban and peri-urban areas. Therefore, the influence of the cell size used in the analysis is further examined, and for this reason, results are given for two indicative scenarios. To the best of our knowledge, this is the first application that considers this number of factors and provides a comparison that highlights the uncertainties in following a standard approach for assessing and mapping of urban susceptibility to flood hazard. The next section provides a description of the data used and the methodological scene. Section 3 presents the corresponding values from the model applications, and the discussion illustrates the results’ variability and provides the necessary comparisons based on previous research. Concluding remarks are presented in Section 4.

2. Materials and Methods

2.1. Methodological Scene

The statistical method that is implemented to examine the influence of various factors on the spatial distribution of urban flooding requires the introduction of the recorded flooded locations as point features in a GIS environment and the determination of several flood-related factors’ values. The first dataset, derived from a geocoding process applied to data, was provided by the Integrated Emergency Coordination Centre of the Hellenic Fire Service and consists of an exhaustive list of all the emergency calls the Centre received the period between 2005 and 2016 including the type, exact time, and the exact location (address) of the flooded property for which the call was made.

To quantify the impact on space, the study area is divided into cells by creating a staggered grid of a given cell size. In this analysis, two approaches were followed; a 400 × 400 m² and a 1000 × 1000 m² grid, to estimate the corresponding variations in the model output. By representing the actual flooded locations as points on GIS and then overlaying this grid, each cell is then given a value, describing the number of flooded locations (i.e., number of points within each cell), as a means of expressing the intensity of the event, that is then called as flood impact ‘FI’, attributing a value equal to ‘0′, ‘1′, …, n. Then, a Boolean map is created that denotes whether at least one incident is observed (‘1′) or no flooded locations observed (‘0′), as a way to introduce two values that are used as indicators for flood occurrence (hereinafter ‘FO’) in space, representing, in essence, the spatial distribution of flooding across the study area on a grid basis. FI and FO are the dependent variables in the regression models (generalized linear regression and binary logistic regression, respectively).

To perform a multi-variate linear regression analysis, one or more predictors or independent variables should be determined. These are the factors chosen on the basis of relevant literature regarding parameters controlling flooding in urban environments [12,23,25,30]. The factors discussed in detail in the following paragraph are estimated as spatially distributed values across the study area and are represented through raster datasets in GIS, in order to finally assign a value in each cell of the two grids used in this analysis. As all factors are initially estimated in a cell size equal to 5 m, which is also the resolution of the digital elevation model, then each factor’s values are calculated for the analysis grid by averaging the respective values of the raster that fall within each cell of the grid, i.e., mean value of 6400 and 40,000 × 25 m²-cells for the case of a 400 × 400 m² and a 1000 × 1000 m² grid, respectively.

The statistical correlation between the factors and floods is explored using a Generalized Linear Regression Model (GLM) and Binary Logistic Regression (BLR), for FI and FO, respectively. In both models, the independent variables (predictors) are the 16 potentially influencing factors and the dependent variable is either the number of flooded locations in a cell (GLM case) or the dichotomous value of whether or not any flooded locations are recorded in a cell (BLR case). In order to denote the potential influence of each parameter, it was assumed that the higher statistical correlation between flood distribution and a factor’s value is observed, the higher the influence (positive or negative) it has on flood occurrence, as this assumption has been used in previous relevant studies [12,27,31,32,33].

Finally, the determination of the most important factors regarding flood susceptibility in a built environment requires the investigation of various combinations among independent variables, which were also tested with the purpose of suggesting the best fit regression model and, if possible, reducing the number of independent variables. The main steps followed for the preparation of the variables are summarized in Figure 2. The factors’ values and variable preparation was completed in ArcGIS Desktop 10.8.2 [34] and the statistical analysis was carried out using IBM SPSS Statistics for Windows, Version 28 [35].

2.2. Factor Selection and Visualization

The selection of influencing factors was undertaken after conducting a deep literature review, with the aim of including in the present analysis factors related to geomorphological and hydro-meteorological characteristics, as well. Flood-related factors can be divided into two main categories (Figure 2): (i) the static factors, which are the factors that remain stable for all events examined and (ii) dynamic, which are the rainfall-related factors and differ among events.

To begin with the factors from the first category, the Digital Elevation Model (DEM) is used to represent the first factor named Elevation (hereinafter ‘ELV’). This dataset is provided by the National Cadastre and Mapping Agency S.A. (Copyright © 2022, National Cadastre and Mapping Agency S.A., Athens, Greece); the format of the files is ERDAS Imaging, the pixel size on the ground is 5 m, the geodetic system is the Hellenic Geodetic Reference System 1987 (known as Greek Grid-GGRSA87), the geometric accuracy of the product is RMSEz ≤ 2.00 m, and the absolute accuracy ≤ 3.92 m for a 95% confidence level. The same layer is used to create other geomorphology-related factors: the layer of Surface Slope (SSL) in degrees, the flow accumulation (FLA) that expresses the contributing area per cell and this layer was used for the stream delineation to finally create another layer for the Euclidean distance from streams (STR). Regarding the distance from the hydrographic network, two additional factors are introduced based on the approach of Feloni et al. [25]. The first factor stands for the Horizontal Overland Flow Distance (HOFD) that expresses the actual movement of water from cell to cell and not the Euclidean distance. This criterion is also a product of the stream network of the area or, more generally, of the DEM [36]. The smaller the flow path from the cell to the riverbed, the more flood-vulnerable the cell is. The second factor stands for the Vertical Distance of Channel Network (VDCN) and expresses the vertical distance between cell elevations and the elevations calculated for the channel network in that cell. Generally, the cells which are out of the stream network will be assigned a value that represents the elevation difference between those cells and the channel that flows through them, in case it existed; the smaller the distance, the more vulnerable the flood area.

There are two additional static factors that are not geomorphology-related, however, they express the hydrological characteristics of the area. The first expresses the degree of surface imperviousness and, to express this attribute, a freely available dataset regarding the degree of soil sealing (SOS) is used that is provided by the European Environmental Agency Fast Track Service Precursor on Land Monitoring [37]. This was assumed as a particularly important factor as anthropogenic soil sealing is one of the most important types of soil degradation. It results in the loss of fertile soils, reduced water infiltration capacity, diminished recharging of aquifers and, when high-intensity rainfall occurs, increased runoff and flooding [10]. The second factor named Composite Curve Number (CNcomp) depicts the Curve number [38], an empirical parameter used in hydrology to estimate infiltration, which is calculated according to the hydrologic soil group, the land cover type, and the soil moisture conditions of a catchment. In this study, the same approach [37] was followed, as the distributed CN values are estimated for average soil moisture conditions and are based on the CORINE Land Cover dataset of the year 2018. However, as a significant part of the region is urban, the information on the percentage of Imperviousness (also see SOS factor) is considered for this criterion, as it captures the percentage of soil sealing and it allows definition of more detailed values of CN inside the urban polygons of the CORINE. Thus, CNcomp is a modified version of the CN layer and it is calculated according to the formula:

$${\mathrm{CN}}_{\mathrm{comp}}=\frac{\mathrm{Imperviousness}}{100}\left(\mathrm{CN}+\left(99-\mathrm{CN}\right)\right)$$

(1)

Finally, dynamic factors are the factors stem from the consideration of the corresponding rainfall measurements per event. In this analysis, eight factors that are directly rainfall-related and one that describes the hydrological conditions (soil moisture level-related) are introduced. Particularly, the meteorological characteristics that can be visualized as raster layers in a GIS environment are the factors of total precipitation depth per event (PCP) and the maximum rainfall intensities of various durations, i, (RIi) that were also considered in previous research for the same part of Athens basin [12,39]. The additional rainfall-related factor controlling flood occurrence that is relevant to the previous conditions of the basin, i.e., the soil moisture level, and is usually expressed by the rainfall conditions in a time window of about 5–10 days before the flood event is the Antecedent Precipitation Index (API) [40]. API, which gives a measure of soil moisture index, is calculated at the beginning of the storm for each one of the available rain gauges using the Kohler and Linsley [41] equation:

$${\mathrm{API}}_{\mathrm{i}}={\mathrm{P}}_{\mathrm{i}}+\left({\mathrm{API}}_{\mathrm{i}-1}{\mathrm{K}}_{\mathrm{i}}\right)$$

(2)

where API_i is the Antecedent Precipitation Index for day i, P is for day i, and K is an empirical decay factor less than one. The decay parameter k must be less than one and is usually between 0.85 and 0.98. Cordery [42] recommended an average value of 0.92 and found that k varied from 0.98 in winter to 0.86 in summer. An annually constant value of 0.95 was recommended by Hill et al. [43] and used in the frame of the current approach.

All rainfall-related factors calculated at the stations’ locations were converted into raster datasets using spatial interpolation in GIS to estimate the spatial distributed values for each dynamic factor. Thus, eight additional layers were created to represent: the API index, the total precipitation depth and the maximum rainfall intensities for six relatively low durations (10 min, 20 min, 30 min, 1 h, 2 h, and 3 h) due to the urban character and the small time of concentration of the area; in the outlet of the main Athens basin (Kifissos river entire catchment area, downstream from the study area) the time of concentration is less than 6 h [44]. Regarding spatial interpolation of rainfall data in GIS, geostatistical methods and IDW have proved to be comparable approaches, in particular, for hydrological modeling [45]. For the factors PCP and API, the IDW method is adopted, while for the factors related to the maximum rainfall intensities the ordinary kriging method using an exponential semivariogram was used, based on Earls and Dixon [46].

3. Results and Discussion

3.1. Analysis in GIS; Assumptions and Uncertainties

The analysis in GIS includes the formulation of the two grids that are used for the regression model development (Fishnet400, Fishnet1000), the determination of the two dependent variables (FI and FO) based on the processing on the dataset that includes the flood incidents, and the creation of raster datasets for both static and dynamic factors (i.e., the independent variables). Regarding dynamic factors, freely available 10 min rainfall data from the available rain gauges, which operate across Athens basin and belong to the Hydrological Observatory of Athens, were collected and analyzed for the calculation of the API index, the total accumulated rainfall during each event for each rain gauge, along with maximum rainfall intensity for various durations, including 10, 20, 30 min, and 1, 2, 3 h. Then, using spatial interpolation in the GIS environment, these values were converted into eight raster datasets, corresponding to the factors API, PCP, and RIi as different layers per event. Static factors are created based on the methodology described in the previous paragraph and the main statistics are summarized in

Table 1. Static factor statistics.

Statistics Per Factor	DEM	SSL	FLA	SSO	CNcomp	STR	HOFD	VDCN	Case
Mean	102.14	3.50	11,530.84	70.70	96.25	163.68	142.04	2.77	Not Flooded Cells, 400 m grid
St.dev.	65.96	3.47	30,815.04	26.78	4.12	103.72	57.49	5.32
CV	0.65	0.99	2.67	0.38	0.04	0.63	0.40	1.92
Mean	76.14	3.29	16,950.85	80.92	97.66	159.41	146.23	2.23	Flooded Cells, 400 m grid
St.dev.	49.97	2.61	48,827.91	19.00	2.49	106.60	57.05	4.24
CV	0.66	0.79	2.88	0.23	0.03	0.67	0.39	1.89
Mean	73.24	2.63	19,567.97	78.94	97.52	154.35	133.41	1.69	Not Flooded Cells, 1000 m grid
St.dev.	56.44	1.67	26,576.03	16.68	1.69	60.63	37.59	1.68
CV	0.77	0.63	1.36	0.21	0.02	0.39	0.28	1.00
Mean	125.51	3.82	5962.58	67.19	95.79	154.64	149.23	3.00	Flooded Cells, 1000 m grid
St.dev.	53.27	2.93	13,050.30	21.27	3.35	46.36	40.89	3.74
CV	0.42	0.77	2.19	0.32	0.03	0.30	0.27	1.25

. The main characteristic is the high coefficient of variation in geomorphology-related factors and the fact that the two categories of cells (flooded and not flooded) are not generally related with different values regarding static factors. This is quite reasonable due to the urban character of the area and due to the fact that we investigated a spatially limited built-up area in which common flooding is attributed mainly to heavy precipitation.

There are also a couple of assumptions on this study that may influence the accuracy of the results. First of all, as it can be observed through Table 1, the grid size controls the average values calculated for each factor per cell. This is also the spatial unit used to describe the presence of flooded structures and, thus, the spatial distribution of flooding across the study area. However, this grid is the main tool that can be used to explore statistical correlations in space and the size that is proposed should be oriented to the special spatial characteristics of each case study. Another source of uncertainty derives from inaccuracies of the digital elevation model that influences five additional factors. Finally, regarding dynamic factors, the station network density, as well as the interpolation methods used control the resulting datasets. Figure 3 presents an indicative result regarding the GIS analysis. After the determination of the dependent and independent variables values, the regression analysis follows.

3.2. Regression Analysis

The application of the binary logistic regression (BLR) on a grid basis for three flood events showed that each flood event, as well as the approach of 400 and 1000 m grids lead to different statistics. An analysis that incorporates the total of 18 independent variables and combinations among factors were tested and it was found that certain factors are indeed statistically significant for flood occurrence, meaning that the probability of the Wald statistic is less than 0.05 (sig. < 0.05). This implies that the null hypothesis (i.e., that is the b coefficient equals zero for a variable) was rejected for these particular factors [47]. The b coefficient, which is a regression coefficient, indicates the relative effect of a particular independent variable on the outcome (FO). For both cases examined (400 and 1000 m) and for all events, the model leads to a different hierarchy and different factors.

In the first case corresponding to the event that occurred in 2013, when using a grid of 400 × 400 m² cell size, a model incorporating seven factors seems to have the ability to represent FO (overall predicting ability equal to 88.4). In this combination, the factors RI2h, DEM, SSO, and CNcomp show statistical significance. On the contrary, the overall predicting ability of the BLR applied in a 1000 m grid is 76.5, and this score is not improved in any other combination examined. In the 1000 m grid approach, CNcomp also shows statistical significance, as well as factors RI3h and API, among the dynamic. In the case of the event of October 2014 and with the use of a grid with a cell size 400 × 400 m², there is a model that achieves 77.5% predicting ability using only the DEM, STR, and PCP factors. In this combination, the b coefficient is −0.011, −0003, and 0.204, respectively, while the constant is equal to −8.23 and each factor’s sig. is <0.05. The SOS factor shows statistical significance in the 1000 m grid analysis. Finally, in the case of the event of October 2015, the majority of combinations examined lead to models with predicting ability higher than 77%. In the best combination (87.2% overall predicting ability), the b coefficient is −0.012, −0.005, −0.002, and −0.756, for DEM, STR, HOFD, and RI30min, respectively. All factors apart from HOFD appear sig. < 0.05. The performance of the BLR was evaluated by Nagelkerke R², which was found to vary between 0.274 and 0.653 that is in agreement with previous research (0.468 in [12]). Finally, in the 1000 m grid approach the predicting ability of flooded cells is satisfactory (91.8%).

Amongst the factors examined, maximum rainfall intensities for various durations did not show a systematic statistical significance. On the other hand, CNcomp in most cases shows high absolute b coefficient and sig. < 0.05, although the area is urban. However, due to the fact that as a factor has very similar values in spatial terms, all models that incorporate CNcomp also include at least one dynamic factor (e.g., PCP and API). The finding that each flood event leads to different statistics and to models that consider different variables indicates the necessity for further investigation and highlights that this approach cannot contribute to the determination of a general model for FO assessment.

The Linear Regression Model performed between the 16 factors and the variable FI, regarding the number of flooded locations per cell, also verified the statistical significance of several factors. The importance of the examined variables is assessed through the absolute value of b coefficient, while a positive b coefficient depicts a positive influence and a negative value of b coefficient depicts that the predicted number of flooded locations is related to a low factor value.

The GLM generally verified the results of the BLR, as in most cases the same predictors are considered and the same factors showed statistical significance. However, GLM cases include also rainfall-related factors. For instance, in the event that occurred in 2013 and after using the small grid (400 m), it was observed that the RI2h is the most significant factor (sig. < 0.05). When using the 1000 m grid, the API index is also significant for the model, but RI30min, RI1h, and RI2h are excluded. The factor RI2h is also excluded in the 2014 event (both cases regarding the analysis grid). In the latter event, the variables RI10min and RI3h are the most decisive for a model describing the FI; two factors that are also significant in the 1000 m grid approach. The BLR for the same event emphasizes the contribution of static factors to the model. Finally, in the third event and based on the results obtained for the 400 m- grid, again rainfall-related factors appear to be more significant, as well as the factor regarding surface slope, while RI20min and RI1h are excluded from the model. Similar results regarding the dynamic factors are found in the 1000 m grid approach, while factors related to the distances from the hydrographic network are also significant.

A comparison with other studies [12,48], even though other parameters are incorporated in their models, shows a higher influence of distance to rivers (−0.013), but a smaller influence of soil-related factors (0.028); a factor similar to the soil sealing factor used in the current approach. Among the factors examined, maximum rainfall intensities for various durations, factors related to the distance from the hydrographic network and the surface slope show in general statistical significance. The latter factor is also found significant in the analysis of Diakakis et al. [12]. On the other hand, the contributing area and the factor regarding soil moisture (API) were found to be statistically less significant. Overall, the goodness of fit for GLM results match to a low R².

4. Conclusions and Future Work

Urban development increases flood risk in cities due to local changes in hydrological and hydrometeorological conditions that increase flood hazard, as well as to urban concentrations that increase the vulnerability. For this reason, further research on factors and mechanisms that describe flood phenomena is essential to achieve an integrated and improved flood risk management.

This study analyzed the spatial distribution of flash flooding in the western part of Athens basin during three events, in combination with potentially influencing factors, and developed a methodology to examine which of them and to what degree, could affect the distribution of flooding within the urban area, in terms of flood occurrence and impact. A significant number of flood-influencing factors has been investigated and they are mainly related to the geomorphological and hydrological characteristics of the study area (static factors) and on the rainfall characteristics of the events (dynamic factors). Both binary logistic regression (BLR) and generalized linear model (GLM) tests that include the aforementioned parameters were examined for two analysis grids. Interpretation of the results indicates that simple geomorphological factors inherent in the catchment have an important influence on flooding distribution in the urban environment; however, all combinations include at least one dynamic factor, which leads to a diversity among results regarding the three events. Different grid cell sizes show generally similar but not identical b coefficients for the flood-influencing factors. However, this diversity in the results is the reason why GIS-based techniques are characterized by limitation regarding flood hazard mapping.

Further modifications, such as aggregation of the stormwater network and other relevant information when available, are required in order to incorporate factors that describe pluvial flooding.

The improved reliable results presented in this study suggest that the method can be applied in other regions, to obtain a more holistic and accurate quantification of the influence of the flood-related factors, under different local characteristics.