A GIS-Cellular Automata-Based Model for Coupling Urban Sprawl and Flood Susceptibility Assessment

Abstract

In Urban Planning (UP), it is necessary to take under serious consideration the inhibitors of the spread of a settlement in a specific direction. This means that all those parameters for which serious problems may arise in the future should be considered. Among these parameters are geo-hazards, such as floods, landslides, mud movement, etc. This study deals with UP taking into account the possibility of widespread flooding in settlement expansion areas. There is a large flooding history in Greece, which is accompanied by a significant number of disasters in different types of land use/land cover, with a large financial cost of compensation and/or rehabilitation. The study area is the drainage basin of Erasinos River in the Attica Region, where many and frequent flood events have been recorded. The main goal of this study is to determine the flood susceptibility of the study area, taking into account possible factors that are decisive in flood occurrence. Furthermore, the flood susceptibility is also determined, taking into account the scenarios of precipitation and the urban sprawl scenario in the area of reference. The study of flood events uses the Analytic Hierarchy Process (AHP) model and the urban sprawl model SLEUTH, which calibrates historical urban growth, using open and cost-free data and software. Eventually, flood susceptibility maps were overlaid with future urban areas to find the vulnerable areas. Following, three scenarios of flood susceptibility with the corresponding susceptibility maps and vulnerability maps, which measure the flood susceptibility of the current and future urban space of the study area, are presented. The results have shown significant peaks in the moderate class of flood susceptibility, while, in the third scenario, high values of flood susceptibility seem to appear. The proposed methodology and specifically the output maps can serve as a decision support tool to assist urban planners and hazard managers in making informed decisions towards sustainable urban planning.

1. Introduction

Among the most severe natural disasters are floods. Their effects have a direct impact on people, infrastructure, and the natural environment in general. It is estimated that the flood risk will intensify in the future, as due to climate change, there is a steady increase in the intensity and frequency of floods [1,2,3]. According to EU Directive 2007/60/EC, “flood risk” is the likelihood of a flood event together with the actual damage to human health and life, and the environmental and economic activity associated with a flood event. The quantification of flood risk results either in monetary units or loss of life units, if the losses are measurable, or in qualitative terms in the case of intangible damages (social, environmental, cultural) to the affected areas.

Floods occur due to the rapid accumulation and release of runoff from the upstream to the downstream, which is caused by heavy rainfall [4]. The likelihood of flooding has increased in many areas due to social and morphological changes, such as land-use change due to urban sprawl. In addition, flood risk is increased due to climate change, severe weather conditions (heavy rainfall), and river discharge conditions [5].

Intense and continuous rainfall usually affects urban areas resulting in the occurrence of floods [6]. Urbanization and urban sprawl are causing agricultural land and vegetation change and population growth [7] in pre-existing urban areas close to river flows. A direct relationship has been observed between urban areas and hydrological characteristics, such as increased runoff, increased frequency, and height of flooding [8].

In general, natural disasters have substantial influences in numerous countries worldwide, with many deaths, infrastructure damage, and population relocations. In addition, due to the significant impacts of climate change, these impacts are expected to increase in the coming years in many countries. Although technology and scientific knowledge have grown significantly, natural disasters continue to have disastrous economic and environmental consequences, as well as many human casualties worldwide. The study and monitoring of these phenomena have been consistent in recent decades and have shown increasing trends [9]. In order to successfully monitor these phenomena and manage their impacts, it is very important to create hazard and risk maps for both the natural and the artificial environment [10,11]. It is a reality that half of the population worldwide is living within cities [12]. Flood hazard maps and flood risk maps refer to areas that, based on historical data and/or other available data, have been identified as potentially serious flood risk areas. Most of these maps are based on specific elements, which are the extent of the flood (applied to this work), and the depth of the water. Furthermore, it is understandable that in data-poor places, for the production of such maps, innovative techniques and approaches must be adopted [13]. Therefore, there are a plethora of methodologies concerning the production of flood hazard and risk maps.

Best practices, in terms of urban studies, presupposes consideration of any natural disasters. Thus, urban sprawl is one of the most vital factors that affect flood risk [14] in cities, mainly due to native changes in hydrological/hydrometeorological conditions. Huong and Pathirana [15] noted that the relationship between the growing urban runoff and flooding is better supposed than that between the cyclic impact of urban growth and the urban rainfall via microclimatic changes. The studies of flooding within an urban environment nowadays include Geoinformatics (GIS, RS, and coupling methods) [3,16,17,18,19,20,21]. The scientific growth in terms of technology during the last decades, mainly in the field of Geoinformatics, has presented new opportunities in floods assessment. Nowadays, the present efforts are targeted at improving and evaluating current models, coupling and linking them, and most significantly endorsing new ones founded on novel ideas and tools coming from emerging technologies, techniques, and sciences [22].

In this context, urban sprawl models which use semi-physical systems, such as the cellular automata SLEUTH, have been used in recent years to produce scenarios in terms of flood risk assessment within an urban environment. This kind of model is based on the computation of a discretized space that encompasses a cell state, a time step, and a set of transition rules to other neighboring cells within a study area [23]. Cellular automata have been used for Flood Assessment [24,25,26,27,28,29], while SLEUTH has been used to examine the impact of urban sprawl and land-use changes on floods and vice versa [12,14,16,30,31,32,33,34]. Another parameter related to flooding within urban environments near coasts is sea-level rise. Sea-level rise, together with the growing population in coastal areas, increases the challenges in terms of urban planning [12,30,35,36].

The main objective of the work is to investigate the susceptibility to flooding of the catchment area of the Erasinos River, located in eastern Attica, taking into account the dynamics of rainfall and land use, as well as the morphological characteristics of the area. This is achieved by selecting rainfall and land-use scenarios in the study area. The study area has a high probability for intense future urbanization, while it is also characterized by frequent river overflows. The proposed analysis attempts a holistic approach to urbanization issues in flood-prone areas, which can contribute to more efficient spatial planning.

For the evaluation of the flood susceptibility of the study area, the key is the selection of the criteria related to the phenomenon as well as the determination of the gravity of each criterion. The creation of maps that provide information about the geographical/geomorphological conditions of the area under study, in combination with the characteristics of the flood phenomenon, helps in proper forecasting and management. It is also very important to communicate the final maps to groups of scientists, as well as to the general public in an interactive way so that as much data as possible are integrated into the modeling process, to better interpret the situation and then better deal with it.

The present work uses the Analytical Hierarchy Process (AHP) with weighted overlay, Geographic Information Systems (accounting QGiS), and the urban deployment model with cellular automata (SLEUTH software). All the used data is free, and the software is not commercial. Therefore, this methodology can be reproduced and transferred to other study areas without any economic cost.

2. Materials and Methods

2.1. Study Area

The present work focuses on the Attica region of Greece (a Mediterranean area where the country’s capital is located) and specifically on a sub-basin of the Erasinos River catchment area (Figure 1). Most Municipalities from Eastern Attica are included in the study area, such as Markopoulo—Mesogeia, Spata—Artemis, Koropi, and Paiania.

The residential areas of eastern Attica are being transformed into main residential areas, while the countryside and traditional crops in the Mesogeia are constantly changing due to urbanization. In this area (which is touristy and with a seaside section), for years, there has been a strong urban spread [35].

In the Attica region, as mentioned in the National Flood Management study, (https://floods.ypeka.gr/egyFloods/gr06/report/I_3_P05_EL06.pdf—accessed date: 5 October 2021) the Potentially High Flood Risk Zones occupy 162.51 km². The study area includes most of the aforementioned Flood Risk Zones.

The Erasinos River is the recipient of the runoff of the central region of the Mesogeia (wider region of Paiania, Koropi, and southern zone of Spata) and the southern region of Mesogeia (Markopoulo, Kouvaras). Specifically, it gathers the outflow of three main streams: Agios Georgios, Markopoulo, and Erasinos itself, which is the final recipient. The catchment area occupies an area of 204 km² and is bounded by the watershed line of the Rafina stream (north), the mountain range of Ymittos (west), and the ridges between the hills Profitis Ilias, Mavrovouni, Stavraetos, Paneion Mt, and Merenda Mt.

In the study area, within the Potentially High Flood Risk Zones, the largest area is occupied mainly by orchards and arboriculture at a rate of 47.48%, followed by broad linear crops at a rate of 29.49%, areas with bare soil at a rate of 11.94% and finally the urban areas with dense construction with a percentage of 7.33%.

More specifically, in the study area, the largest area is occupied by olive groves, vineyards, and forests, and the smallest area is urban areas.

In addition, 5 areas with industrial facilities are located in this area, while there is also the archeological site of Vravrona, where the Temple of Artemis is housed. It is basically a small valley, between small hills and the sea. The Erasinos River, which flows through this valley, floods the temple when it overflows.

The climate is characterized as Mediterranean. Only in the highlands is the climate considered mountainous. Specifically, the average annual rainfall ranges from 350 mm in the Attica basin to 1000 mm in the mountainous parts (Parnitha) [36], rainy days range from 50 to 100 per year, snowfall is rare in the coastal areas and increases inland, depending on the altitude and distance from the sea the average annual temperature ranges from 16 to 18 °C, the annual range is about 16 °C, and the average relative humidity is between 56–64% (https://floods.ypeka.gr/egyFloods/gr06/report/I_3_P05_EL06.pdf—accessed date: 5 October 2021).

The area has been included in the Natura 2000 network, as well as in the network of the conservation of natural habitats and of wild fauna and flora based on the Council Directive 92/43/EEC of 21 May 1992.

Several floods have been recorded in the Erasinos river basin, which presents fewer and milder floods compared to other rivers and streams in Eastern Attica. Nevertheless, significant flood phenomena have been recorded that destroyed the ancient temple of Vravrona Artemis, which was covered by the sediments of the flooded field, etc. Specifically, in the last 3 decades, at least 12 floods have been recorded (1991, 1993, 1998, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2013, 2014), some of limited and some of greater extent. The damages that have occurred and their frequency show that the problem is real (https://mthymettosgreece.com/topothesies/erasinos-potamos/—accessed date: 5 October 2021).

The studied sub-basin through the prediction of the SLEUTH model presents the highest rates of population change in the period 2011–2031 and the potentially greater urban spread compared to the rest of Attica [37].

2.2. Data

The methodology presented in this paper is based on a set of data, the vast majority of which are physical geographical. Floods are directly related to the topography of an area, the geological/hydrolithological background of the area, its hydrographic network, and land use [38,39]. Thus, for example, an area with a strong relief, without vegetation and non-water permeable geological formations, presents a high probability of flooding, especially in the case of a sudden storm of very high intensity. In other words, the flood susceptibility of an area depends to a large extent on the physical geographical/geological conditions that prevail in it. In the present work, the data used are related to topography (Digital Elevation Model (DEM)), hydrolithology, and land use of the study area. Furthermore, for rainfall, which is the main activating factor for the occurrence of flood, values of 6-h rainfall (8:00–14:00) or 12/12/14 from 6 meteorological stations (Glyka Nera, Koropi, Markopoulo, Pallini, Spata, Airport) were used, and the Spline interpolation method was used to create a continuous level of rainfall in the study area. The basic reason for choosing the Spline method was because the study area is exactly inside the meteorological stations, while it is a two-dimensional minimum curvature (plate interpolation) [40,41].

Rainfall data were provided by the Hellenic National Meteorological Service for the 6 meteorological stations.

The following

Table 1. Data and parameters of the methodology.

Data	Type	Spatial Resolution	Source
Altitude	Raster	30 × 30 m	Aster GDEM *
Hydrolithology	Vector (Polygon)	1:50,000	EAGME **
Land Use	Vector (Polygon)	1:5000	OPEKEPE ***
Slope	Raster	30 × 30 m	Aster GDEM
Distance from hydrographic network	Raster	30 × 30 m	Aster GDEM
Precipitation	Raster	30 × 30 m	Hellenic National Meteorological Service—HNMS ****
Flow accumulation	Raster	30 × 30 m	Aster GDEM

* https://asterweb.jpl.nasa.gov/gdem.asp—accessed date: 5 October 2021; ** https://www.eagme.gr/—accessed date: 5 October 2021; *** https://www.opekepe.gr/en/—accessed date: 5 October 2021; **** http://emy.gr/emy/en—accessed date: 5 October 2021.

presents the parameter variables used in this work.

2.3. Methodological Framework

Flood susceptibility maps were created by applying Multi-Criteria Analysis (MCA), which aimed to prioritize the gravity of the parameters. For this reason, the Analytical Hierarchy Model was used in order to estimate the weights for each factor. Moreover, an attempt was made to estimate the vulnerability of the area to a flood event based on both the current land cover situation and the estimated future coverage taking into account the results from the urban sprawl model SLEUTH. Figure 2 shows the flow chart of the proposed methodology.

2.3.1. The AHP Model

The Analytic Hierarchy Process (AHP) model is a method of converting subjective scores of relative weight into a set of weights [42]. One of its advantages is that it handles efficiently and simply many evaluation criteria, both quantitative and qualitative, with the pairwise method comparisons [43,44]. More specifically, it includes a table of binary comparisons and develops one scale of preference between factors. The scale of preference in question is formed based on the assessments of the designer and is taken into account in creating weights in the respective factors.

The Analytic Hierarchy Process (AHP) is one of the Multi-Criteria Decision-Making techniques. Due to its extensive range of application areas, it has been an overwhelming investigation topic for various researchers around the world. AHP takes into account the relative priorities of variables trying to maximize the quality of the results. The issue is that this method cannot be used in some cases (e.g., linear equations). It is a very simple method to use, but when a problem has many variables, then the computational requirement increases [45,46,47]. There are a plethora of models and methodologies that can provide the required results in terms of weights of parameters [48,49,50,51,52,53,54].

A weighted overlay is then applied to produce the flood susceptibility map. At this stage, all the factors are reclassified to a regular scale of 1–5, where 1 corresponds to very little susceptibility (e.g., regions with very steep slopes) and correspondingly, 5 corresponds to a great susceptibility (e.g., flat areas or areas with a very small slope) before superimposition. The reclassification is followed by the weighted overlay using for each factor the weights that emerged in the first stage of the Analytical Hierarchy.

Furthermore, the pairwise matrix with all the selected parameters was sent to experts in the field, and the final pairwise comparison matrix for the AHP method is presented in

Table 2. Data and parameters of the methodology.

	Flow Accumulation	Altitude	Precipitation	Land Use	Distance from Hydrographic Network	Slope	Hydrolithology
Flow Accumulation	1	2	2	2	3	4	6
Altitude	1/2	1	2	2	3	4	6
Precipitation	1/2	1/2	1	1	2	3	3
Land Use	1/2	1/2	1	1	2	3	3
Distance from hydrographic network	1/3	1/3	1/2	1/2	1	4	3
Slope	1/4	1/4	1/3	1/3	1/4	1	3
Hydrolithology	1/6	1/5	1/3	1/3	1/3	1/3	1

.

In particular, the first step of the method concerns the creation of the hierarchical structure (tree), with the top level depicting the overall objective of the analysis. The intermediate levels represent the factors that influence the sought objective, while the bottom of the tree includes the variables that influence these factors. Then, the relative importance of every criterion in each level of the tree is defined via their pairwise comparison. The comparison between each pair of criteria is made with relative grades expressed as absolute numbers. This scale of numbers indicates how many times more important or dominant one element (criterion, variable, etc.) is over another. The weight of each hierarchy element is calculated by the comparative matrices and verified by a consistency check.

$$A=\left[\begin{array}{ccc}{a}_{11}& \cdots & a_{1n}\\ ⋮& \ddots & ⋮\\ a_{n1}& \cdots & a_{nn}\end{array}\right]$$

(1)

where $a_{11}$,…, $a_{nn}$= 1, and α₂₁ = 1/α₁₂, etc.

Taking into consideration the above principles, the number of comparisons required were n(n − 1)/2, where n is the number of elements.

The values of the normalized matrix were determined based on Equation (2) and priority vectors indicating weights of elements from Equation (3).

$$w_{ij}=\frac{a_{ij}}{{{\displaystyle \sum }}_{i=1}^n{a}_{ij}}$$

(2)

$$w_i={\displaystyle \sum }_{i=1}^n{w}_j{a}_{ij}$$

(3)

where

$$w_j=\frac{{{\displaystyle \sum }}_{i=1}^n{w}_{ij}}{n}$$

(4)

The consistency estimation of the produced weights involves the calculation of the consistency ratio (CR). CR indicates the probability that the matrix values have been randomly generated [55].

Saaty [42] argued that when a matrix has a CR value greater than 10%, it should be re-evaluated [2]. Consistent matrices have CR ≤ 0.10. The calculation of the consistency ratio was conducted through the following Equation (5).

$$CR=\frac{CI}{RI}$$

(5)

The Consistency Index (CI) was calculated as shown below (Equation (6)).

$$CI=\frac{{\lambda }_{max}-n}{n-1}$$

(6)

where ${\lambda }_{max}$ is the average value of the consistency vector, and n = the number of criteria. The random index (RI) represents the CI value of a pairwise comparison matrix, which is randomly generated and depends on the number of the compared elements.

After having calculated the weight of each element in each corresponding level of the hierarchy (criteria, sub-criteria, variables) and checked its consistency, the final step of the method follows, which is the synthesis of the main goal of the analysis (the top of the tree structure). This can be accomplished via the following formula Equation (7).

$$G={\displaystyle \sum }_{i=1}^n{w}_i{c}_i$$

(7)

where G is the main goal of the analysis, $w_i$ is the weight of a criterion, $w{c}_i$ is the value of a criterion and i = 1, 2, 3…., n is the number of a criterion.

At this point, it should be highlighted that AHP was applied using GIS and refers to spatial data/criteria. The final aim of this methodology part is the formation of a spatial index map via Weighted Linear Combination (WLC), where the weights of the spatial criteria are calculated by AHP. Therefore, the factors of Equation (7) refer to the pixel values of the spatial layers (goal and criteria) and were transformed as presented below (Equation (8)).

$$G^p={\displaystyle \sum }_{i=1}^n{w}_i{c}_i^p$$

(8)

where $G^p$ is the value of each pixel of the main goals is spatial layer and $c_i^p$ is the value of each pixel of each criterion’s spatial layer [55].

The final formula for this current methodology is:

$$\begin{array}{ll}{\displaystyle {\displaystyle \sum }_{i=1}^7}{w}_i{c}_i=& Flow accumulation\times 0.2857+Altitude\times 0.2279+Precipitation\\ & \times 0.1415+Distance from the hydrographic network\\ & \times 0.1043+Slope\times 0.0601+Hydrolithology\times 0.0391\end{array}$$

(9)

Flow accumulation is considered the most important factor [39] and especially for the study area, as several streams are collected in the main bed of the Erasinos River. Next is the altitude, which does not fluctuate much, as it is low in most of the study area. The following are the land uses and the intensity of the rainfall that were considered to be of equal importance for the occurrence of the phenomenon [56] but also with the historical flood events of the area. Most of the flood events recorded in the study area were located in urban areas and under conditions of heavy rainfall. The distance from the hydrographic network is defined with less comparative importance in the event of a flood in the specific area. Immediately after are the slopes, which are in general mild in the area and are therefore of lesser importance. Finally, hydrolithology is classified last, as most of the land is covered by forests, crops and settlements, without sparse vegetation.

2.3.2. The SLEUTH Model

Based on the mechanism of cellular automata, a diverse range of models have been developed and improved to simulate and predict future urban growth over the past few decades [57]. Among these, the SLEUTH model has been widely applied due to its comprehensive evaluation and high accuracy [58,59] and, as a simple single-land-use model, is unrivaled for rapid exploration and simulation of urban growth, pattern and form, as their enduring popularity attests [57,60,61,62,63,64,65].

SLEUTH model was developed by Keith C. Clarke [54], in research of the University of California’s Department of Geography, which started initially under the supervision of the USGS (United States Geological Survey) and then of the NSF (National Science Foundation). The SLEUTH model is essentially a combination of the UGM (Urban Growth Model) and the LCD (Land Cover Deltatron) model. UGM is a C programming program, which works with the operating system UNIX (Linux environment), while the LCD is a land cover model. The LCD has built-in code, which can be read by UGM. So, rgw SLEUTH model was created, which is equipped with everything is needed to capture and predict a dynamic phenomenon, such as urban sprawl.

SLEUTH is an acronym made up of its input datasets:

• S lope;
• L and uses;
• E xclusion (Building Prohibition Areas);
• U rbanization;
• T ransportation;
• H illshade.

Specifically, using the QGIS software, the corresponding thematic layers are created. Utilizing the input datasets, the SLEUTH model sets out five growth parameters:

• Diffusion (determines the minimum and automatic probability of urbanization);
• Spread (determines the part of the development that derives from existing urban centers);
• Breed (defines each new urban place to develop into an urban center);
• Slope resistance (determines the reduction in urbanization due to ground slopes);
• Road gravity (determines the urbanization, which follows the road network).

Among the models based on cellular automata, the SLEUTH model is the most accepted. There were many reasons for choosing SLEUTH in this project: it is available free online and with technical support, it includes a fairly powerful mechanism for historical calibration, it can simulate future development scenarios, its results are in the form of GIF maps, which are effective visualization tools, and finally it has been successfully applied in many recent studies on urban development [66,67,68].

It is required to simulate urban growth as a key asset to better inform planners and authorities in landscape design and urban planning. Figure 3 summarizes the process followed by the SLEUTH model for final results’ production.

Particular attention should be paid to calibration, as it is the most crucial step of the model. Essentially the model calibrates the historical data given in order to calculate the appropriate parameters, based on which an urban sprawl has been created (Figure 4).

The calibration, which works according to the Lee Salee Metric method, is divided into four phases. Each phase requires specific changes to the files and images data. The function of the SLEUTH model is specific. Initially, it identifies the parameters that characterized the existing urban sprawl and based on them, predicts the upcoming one. Therefore, the trend of urban development, whether reduced or increased, is also reflected in the forecast [37].

In this current work, SLEUTH provided urban sprawl maps for the period 2011–2030.

2.3.3. Precipitation Scenarios

For the needs of this work, 2 rainfall scenarios were used. The first rainfall scenario results from the rainfall values recorded at the 6 meteorological stations that include the study area on 12/12/14, where there was heavy rain, and as a result, many areas in the Erasinos river basin were flooded.

Based on the six-hour rainfall values on the day in question and through the Spline interpolation method, the corresponding distribution of the cumulative 6-h rainfall was obtained.

The second rainfall scenario concerns the worst rainfall scenario in terms of increasing flood susceptibility, out of a wide range of 18 hypothetical rainfall scenarios. The scenarios in question arise by increasing the actual rainfall values of the 6 meteorological stations of the 1st scenario by 10%, 30%, and 50%. Therefore, 18 precipitation maps were produced according to the above scenarios (Figure 5).

In the above Figure, the first image is for an increase of 10% precipitation at the Glyka Nera station, 2 for an increase of 30% precipitation at the Glyka Nera station, 3 for an increase of 50% precipitation at the Glyka Nera station, 4–6 for an increase respectively at the Koropi station, 7–9 at the Markopoulo station, 10–12 at the Pallini station, 13–15 at the Spata Station and 16–18 at National Airport station.

The above maps were categorized on a common susceptibility scale to check which scenario resulted in greater flood susceptibility to the rainfall factor. From the generated maps, the area (in km²) was calculated per susceptibility category (very small, small, medium, large, and very large).

Based on the susceptibility data, we compared the maps that showed the square kilometers in the middle (9th map), in the high (15th map), and very high category (18th map) of flood susceptibility, in terms of the rainfall factor.

It was observed that regarding the important categories of flood susceptibility (moderate, high, very high), only one category with the highest percentage appeared per map.

Specifically, the percentage in the middle category appeared higher on the 9th map compared to the rest (16.3% compared to 14.9% on the 15th map and 9.1% on the 18th).

The percentage in the high category appeared higher on the 15th map compared to the rest (14.7% compared to 12.9% on the 9th map and 5.9% on the 18th).

Finally, the percentage in the very high category appeared higher on the 18th map compared to the rest (10.4% compared to 8.2% on the 15th map and 2.4% on the 18th).

Based on the above and for the selection of the map that will most affect the flood susceptibility of the study area, in terms of the rainfall factor, the percentages of the maps were compared cumulatively for the high and very high categories. In this case, the percentage of the 15th map in the high and very high zone was 22.9%, compared to 16.3% and 15.3%, respectively, for the 18th and 9th maps. Considering that the percentage in the middle category of the 15th map was also large, the susceptibility map in question (an increase in rainfall values by 50% at Markopoulo station) was selected for the worst-case scenario of flood susceptibility to rainfall.

The flood response of the area in the scenarios with the actual rainfall values and an increase of 10% and 30% of the actual rainfall values appeared relatively limited, while the scenario with an increase of 50% of the real values brought a substantial response to the flood susceptibility area.

2.4. Validation

The process of validation of the model concerns the model itself and consists of investigating through the program, which has been verified, the question of whether the model is an acceptable/satisfactory representation of the system being simulated and therefore whether it incorporates logical details and assumptions that the analyst initially wanted to include [70]. That is, the accuracy of the system transfer to the model was investigated, or at least whether it satisfactorily approaches the operation of the system in relation to the entities, events, and activities of interest.

The validation of the results of this methodology was based on actual flood data. In particular, the scenarios compared with a flood event that happened in 12 December 2014. In addition, the historical flood events of the last decades were used to examine the validity of the model.

3. Results

3.1. Main Results of the Methodology

All the produced maps obtained by this methodology were at a spatial resolution of 30 × 30 m, mainly due to the ASTER GDEM and its derivatives.

As already mentioned, SLEUTH provided urban sprawl maps for the period 2011–2030. Figure 6 presents these produced maps.

The final map of the year 2030 was used for the needs of this current methodology.

In the present work, three different scenarios of flood susceptibility and urbanization were examined.

The first scenario took into account the seven factors that can affect flood susceptibility (Figure 7) and used the models of analytical hierarchy and urban sprawl.

More specifically, a flood susceptibility map was created based on the factors of Table 1 and with the values of 6-h rainfall for 12 December 2014 (where due to heavy rain, large areas in Eastern Attica were flooded). In addition, a layer was created that identifies urban areas. To achieve this, the layer of land used was categorized into two categories, urban and non-urban. By isolating the category of urban space, the desired layer was created. Using the same layer of flood susceptibility again, the layer of the urban area was changed, as it was formed with the help of the SLEUTH model. The aim was to capture future vulnerable urban areas from the flood phenomenon, assuming that the susceptibility of the area will not change in the future. The following Figure 8 shows (a) the flood susceptibility map, (b) the current state of the urban fabric on the previous map, and (c) the urban sprawl with the SLEUTH estimation incorporated in the susceptibility map.

With the creation of the above maps of the present scenario, it was possible to quantify the urban area in relation to its susceptibility to flooding. Below is

Table 3. Scenario 1 for the assessment of susceptibility (in km²).

Susceptibility	Current State	Urban Sprawl
Very low	0	0.09
Low	1.657	23.102
Moderate	18.118	73.95
High	0.624	7.99
Very high	0	0.022

with the total area occupied by the urban area by category of flood susceptibility, in the cases of the current situation and the urban spread provided by SLEUTH.

Cumulatively, it was observed that a significantly larger urban area (84.75 km²) fell into the categories of high flood susceptibility in the case of urban sprawl compared to the current situation. Figure 9 presents these results.

Noteworthy is the large increase observed in the case of urban sprawl, mainly in the moderate category of flood susceptibility but also the occurrence of 8 km² versus 0.6 km² in the high category of flood susceptibility of the current situation.

The second scenario also took into account the seven factors related to flood susceptibility and used the models of analytical hierarchy and urban sprawl. The main difference with Scenario 1 is that a different level of rainfall was chosen to create the flood susceptibility map, while the other factors remained the same.

More specifically, different cases of rainfall were examined by increasing the rainfall values of the six meteorological stations for 12 December 2014, and the case of an increase of 50% of the values of Markopoulo station was selected as the most unfavorable for the flood susceptibility of the study area. The layer of flood susceptibility was then used by superimposing the level that defined the existing urban areas in the 1st scenario. This process made it possible to quantify the urban area in relation to its susceptibility to flooding. Using the same layer of flood susceptibility again, the level of the urban area was changed, as it was formed with the help of the SLEUTH model. The aim was to capture the future vulnerable urban areas from the flood phenomenon, assuming that the background conditions of the area will not change significantly in the future. The following Figure 10 shows (a) the flood susceptibility map, (b) the current state of the urban fabric on the previous map, and (c) the urban spread with the SLEUTH estimation on the susceptibility map.

With the creation of the above maps, it became possible to quantify the urban area of the 2nd scenario in relation to its susceptibility to flooding. Below is

Table 4. Scenario 2 for the assessment of susceptibility (in km²).

Susceptibility	Current State	Urban Sprawl
Very low	0	0.138
Low	1.471	21.116
Moderate	16.413	74.326
High	2.513	9.506
Very high	0	0.025

with the area of the urban area per category of flood susceptibility, in the cases of the current situation and the urban sprawl.

Noteworthy, in this scenario as well, the large increase in the case of urban sprawl mainly in the moderate category of flood susceptibility but also the occurrence of 9.5 km² versus 2.5 km² in the high category of flood susceptibility of the current situation was observed.

The third scenario aimed to map (visualize) the possible change in flood susceptibility of the study area, which was calculated and depicted in the relevant map of the first scenario, taking into account the possible change in land use due to future urban sprawl for the year 2030. The response of the flood susceptibility of the study area with the possible change in land use was due to the increase in the urban area. Therefore, a third flood susceptibility map was created, using the seven factors of the first scenario, with the only difference that a different layer of land use was used. Via SLEUTH urban sprawl model, land use is changing, and large areas of crops and other uses are being converted into urban ones. By shaping the specific layer of land use, the necessary reclassification to the common scale was made to create the map of flood susceptibility. Figure 11 shows the flood susceptibility map.

The graph of Figure 12 was produced from the three flood susceptibility maps.

The above graph compares the three flood susceptibility scenarios based on the distribution of their surface in the five flood susceptibility categories. From the produced flood susceptibility maps, it was observed that the third flood susceptibility category (moderate) included the largest area of the study area for all three scenarios. Moreover, the category of very high flood susceptibility included an area of less than 0.5 km² for all three scenarios, while about twice the area was included in the category of high flood susceptibility of the third scenario, compared to the other two.

It was generally observed that the areas with the highest flood susceptibility were the urban areas that received the highest rainfall values.

3.2. Validation of the Methodology

It was necessary to compare the three flood susceptibility maps-scenarios with the flood behavior of the study area to check the validity of the results of the scenarios.

More specifically, the first scenario was checked based on the flood event of 12 December 2014. The flooded areas were recorded near the mouth of the Erasinos River, in the archeological site of Vravrona, as well as in arable lands on either side of its riverbed. The above areas were included in the high and very high flood susceptibility category of the flood susceptibility map of the first scenario, verifying its validity.

Regarding the second flood susceptibility map, which took into account the worst-case scenario of rainfall values, its validity was examined based on the recordings of historical flood events of the last decades.

Finally, the third flood susceptibility map, which concerned the future change in land use due to urban sprawl, confirmed the increase in the flood susceptibility of the urban space in the study area.

Below (Figure 13) is the cartographic comparison of the first flood susceptibility map with its flood area during the flood of 12 December 2014 and the comparison of the second scenario with the recorded flood events of the last decades.

The flood susceptibility map of the study area in the current situation (1st scenario) took into account the actual rainfall values recorded by six meteorological stations on 12 December 2014.

To check the results, the cartographic comparison of the flood susceptibility map with the area that flooded on the day of reference was selected, according to reports from residents and competent public bodies.

With the overlay, a comparison of the first flood susceptibility scenario with the actual flood event could be made (Figure 13).

The flooded areas were concentrated at the mouth of the Erasinos river, in the archeological site of Vravrona, and the surrounding cultivated areas. As can be seen from the comparison map, the flooded areas coincided with a large part of the very high and high flood susceptibility categories of the first scenario. This comparison verified the first susceptibility scenario and its validity.

4. Discussion

4.1. Main Aspects of Discussion

The methodology has provided significant results to be considered. In particular, from the comparison of the urban area of the 1st scenario in each category of flood susceptibility, for the current situation and the case of urban sprawl, it was observed that significantly larger extents of the urban area (84.75 Km² or 415%) fell into categories of flood susceptibility in the case of urban sprawl, compared to the current situation. Remarkable was the large increase observed in the case of urban sprawl, mainly in the moderate category of flood susceptibility but also the appearance of 8 km² versus 0.6 km² in the high flood susceptibility category of the current situation.

Furthermore, from the comparison of the urban area of the 2nd scenario in each category of flood susceptibility, for the current situation and the case of urban sprawl, it was observed that significantly larger urban areas (84.75 km² or 415%) fell into categories of flood susceptibility in the case of urban sprawl, compared to the current situation. Here as well, the large increase observed in the case of urban sprawl mainly in the moderate category of flood susceptibility but also the appearance of 9.5 km² versus 2.5 km² in the high flood susceptibility category of the current situation must be underlined.

In general, the fact that urban models are unable in most cases to capture, incorporate, and co-assess the main driving forces leading to urban growth [71,72] must be highlighted. Socioeconomic factors, such as population growth and economic development that mainly affect and drive urbanization [73], are left out of these models. Thus, geospatial and socioeconomic data amplified with an obligation of forceful urban processes are required in the urban policy assessment [74]. As expected, this is SLEUTH’s limitation as well [75]. So, even though a lot of studies on urban sprawl can be found for various areas throughout the world [76,77,78,79,80], it is hard to follow a generic modeling approach. As mentioned above, this is not only due to the different background characteristics (e.g., land use, slope, etc.) of each area but mostly due to the unique socioeconomic and demographic regime and trend of each study site [81]. Having said that, it can be easily deduced that this is a limitation of the current research work as well.

Furthermore, another limitation of this work, as with more research projects that deal with data analysis, lies under the specific characteristics of the selected and used dataset. In this research work, leaving aside data gathered from national/local sources and repositories (e.g., land use cover obtained by OPEKEPE, hydrolithology map from EAGME, and rainfall data-series from the National Weather Service) that cannot be compared and evaluated in a common and wide (global) frame, satellite data from the ASTER GDEM initiative were exploited. The DEM and all the other spatial layers derived from it have a 30 m pixel size as a base of reference, and thus the produced results of the analysis depend on the accuracy and generalizations of this spatial scale. In general, there are various studies on urbanization/urban sprawl and natural hazards (flooding, landslides, sea-level rise, etc.) leaning on remote sensing and mainly satellite data [14,16,36,66,82] and these approaches will keep multiplying as contemporary satellite projects producing state-of-the art and highly detailed data keep rising. In this light, the present work must be currently evaluated based on ASTER satellite data characteristics and re-evaluated in the future using more detailed satellite products (e.g., SENTINEL data).

In addition, there are a plethora of studies concerning natural hazards and their impact on socioeconomic, spatial aspects of human life [14,30,34]. The need to take under consideration natural phenomena and mainly extreme events and disaster risk in policymaking, design, and management of various aspects of human life are already imprinted in contemporary research trends. It is nowadays of utmost importance to work on holistic and combinatorial approaches to cope with climate change pressures, efficiently adapt to their impact, and design the future world that humanity will safely live and prosper. This paper worked on three different scenarios coupling different rainfall events (one recorded and one extreme), and different land uses (current state and more urban areas) to evaluate the cumulative and both ways effect of the dipole flooding—urbanization. Although the results, as presented in the previous section, give a clear picture of the effect of flooding on urban sprawl and vice versa, the model must be expanded and tested in the three IPCC’s scenarios (RCPs).

4.2. Limitations and Future Work

The specific analysis presented in this study, like all corresponding analyzes, is subject to specific limitations. These restrictions apply to specific parts of the analysis as well as to the used data. In particular, geospatial and socioeconomic data augmented with a requirement of dynamic urban processes are required in the urban policy assessment. Actually, this is a strong limitation of the SLEUTH model itself. This model cannot, at the moment, assign the different background characteristics of each study area but mostly the unique socioeconomic and demographic regime and trend of each of them.

Additionally, another limitation of this work lies in the specific characteristics of the selected and used dataset. Most of the time, the used data cannot be compared and evaluated in a common and global frame. In this work, all the data did not come from a common source with similar characteristics in terms of accuracy.

In this light, future work should take into account these limitations. In addition, it should include a comparison with all existing and available studies in order for the conclusions to be even safer.

Overall, the methodology of this work, despite its limitations and considering that it can be further improved as discussed above, can help decision-makers significantly in urban planning [63], design proper mitigation measures in vulnerable urban areas, and set operationally efficient policies for optimum urban adaptation and management under future climate regime and mainly against extreme natural phenomena and disasters. Directive 2007/60/EC on the assessment and management of flood risks, and a lot of the SDGs set by the UN (e.g., SDG 3 on good health and well-being, SDG 11 sustainable cities and communities, etc.) can be supported by Decision Support Systems (DSS) that integrate interdisciplinary and combinatorial Hazard-Urban Sprawl (HUS) models [12], which can also be fed with climate change scenarios and produce educated predictions [33,82]. This work steps on these pillars.

5. Conclusions

The Erasinos river basin is the largest part of the Eastern sector of Attica Region. Compared to the rest of Attica, it is distinguished by the preservation of traditional land uses the development of animal husbandry and agriculture. But in recent years, the urban space has been growing at a very fast pace, while some crops have been being transformed into urban space.

Various records from the last decades show that the area has been affected by several floods. Most floods usually occurred in the archaeological site of Vravrona, located at the mouth of the Erasinos River, but also in the settlements of the area, with the main cause of their occurrence being heavy rainfall.

The flood susceptibility of the studied area increases along with the increase in rainfall values compared to those already recorded in the area, affecting mainly the areas on either side of the Erasinos main riverbed, the Archaeological site of Vravrona, the Markopoulo Hippodrome, and parts of Kouvaras, Kalivia, and Markopoulo settlements. Furthermore, with the possible land-use changes that the area will undergo over time, taking into account the factor of urban sprawl, the aforementioned areas and parts of the settlements of Koropi and Karella will be affected to a greater extent. In addition, the flood vulnerability of urban space increases significantly in the case of urban sprawl.

This is a methodology to be considered by, mainly, Urban Planners toward suggesting that urban sprawl should not be extended to areas that may be vulnerable to flooding in the future or indicating the need for proper preventive and protective measures and mitigation actions to be adopted if urban sprawl becomes a reality in these vulnerable areas. These maps can be used in conjunction with flood hazard maps wherever possible in order to depict flood-prone areas in their full scale.