Flood risk assessment and flow modeling of the Stalos stream area

The rise of average temperature is a globally realized and confirmed phenomenon, and there is a general consensus on the scientific community that this is due to anthropogenic stresses, and it leads to the increase of the precipitation events’ severity, as well as duration (Lyu et al. 2016). The frequency of severe flood incidents has increased, necessitating an integrated approach in managing flood risk rather than mitigating the flood effects, utilizing hydraulic and hydrological modeling (Mehta et al. 2014). The term ‘flood’ is most notably used to describe the uncontrolled inundation of an area with water, as a result of the overflow of water from a nearby stream or river over the natural or artificial banks. The basic mechanisms leading to the flooding of an area adjacent to a river or stream are (a) the increase of water flow or (b) the reduction of the cross-section area of the river/stream bed (Sakkas 2004).

Floods that emanate from overflow river systems vary considerably regarding their magnitude and duration. In the case of large rivers, floods can occur long after the rainfall and last for days, weeks, or even months. In small river basins, which exhibit intense topography with thin and impenetrable soil, flash floods can also occur (Georgakakos 1986). A flash flood is defined as a flood which follows shortly – within a few hours – after a heavy or excessive rainfall event. Flash floods can also occur due to rapid snow melting. The karst flash flood is a special kind of flash flood related to the structure and hydraulic properties of a karstic system and has been identified as one of the hazards in karstic terrains. The majority of the karst flash floods appear in arid and semi-arid areas (Marechal et al. 2008).

Several models presented in the past can estimate potential flood damages for different scenarios using flood parameters obtained from historical data, e.g. the Flood Damages Analysis Package developed at the Hydrological Engineering Center (HEC) (USACE1 2002). However, these models cannot be used for real-time flood loss estimation or forecasting as there is no well-established mechanism available for simulating the flood parameters of an actual flood event from the physical characteristics within the flooded area (Dutta et al. 2003). The importance of modeling open channels under varying physical characteristics is unmistakable, especially in examining the constraints of each stream, such as velocity, composite roughness, Froude number, and channel stability, leading to optimal cross-section profiling of the stream (Roushangar et al. 2021). Recent models of real-time flood loss estimation use Geographic Information Systems (GIS) to delineate flood inundation (Consuegra et al. 1995; Jonge et al. 1996). However, these models do not include a physically based model for the flood inundation simulation component; instead, they use GIS to spread the floodwater on the basis of a river model simulation. The main difficulty associated with the flood inundation simulation approach is to obtain an adequate number of flood inundation parameters, which are the most important inputs to a flood loss inundation model. Recently, works that combine a physically based flood inundation simulation model and a loss estimation model have been presented (Dutta et al. 2003, 2006; Jonkman & Vrijling 2008; Kourgialas & Karatzas 2013; Vozinaki et al. 2015). In this paper, a spatial flood risk assessment is combined with a flood inundation model to determine the overall flood risk of the study area and the viability of the GIS spatial flood risk method.

The study area was selected as a suitable area for a pilot level implementation of the flood risk assessment methodology, which was adjacent to historically reported flood events. The study area contains a stream, which had been modified in the past to divert its bed, without the proper engineering studies and measures taken. The stream runs from south to north, originating from the village of Stalos in the form of two separate streams merging into one, before the village of Kato Stalos and draining into the Chania bay.

The watershed has a total area of 3 km², a perimeter of 31.9 km, and spreads at a direction from north to south from the Northern coast of the Chania prefecture to the southern limits of the Stalos village, while at a direction of west to east from Platanias area to Kladissos stream (Figure 1). The studied stream of this area is in the middle of the watershed, originating from two streams in the southern domain of the watershed, ending at the Kato Stalos coast.

Utilizing AcrGIS's Hydrology Toolbox, with additional data from creating a DEM, and rasters for flow accumulation and flow direction, the stream was found to consist of five separate branches, with each having its own sub-watershed. The final modeled portion in the 1-D–2-D coupled model had a length of 814 m.

The area displays a varying topography, with the maximum altitude being at 265 m at the south of the watershed, whereas the lowest being 0 m at the northern coastal limits of the watershed. The varying topography leads to multiple seasonal and permanent streams in the area, as denoted in the Digital Elevation Model (DEM) of the study area.

In addition, regarding the surface geological formations, the area exhibits mostly porous formations of varying hydraulic permeability (P1–P3), with additional areas with impermeable formations (A1), as well as karst formations (K1) – mostly limestone, as shown below in Figure 2.

The first part of flood risk assessment at the study area was to determine the factors that comprise flood risk for the area. According to Kourgialas & Karatzas (2011), the significance for the estimation of hazardous areas is a fundamental task in environmental management, started for the assessment of groundwater contamination, and ended up for the study of natural disasters such as flooding, earthquakes, and soil erosion.

The main goal of mapping flood risk areas for a specific domain is to classify it into categories based on the degree of flood risk, with the categories of risk being very high, high, moderate, low, and very low.

To estimate the specific areas of a river basin with increased flood risk, six thematic map layers are needed. Each thematic map is created in the ArcGIS environment, in raster format. Specifically, the six factors transformed in spatial raster data are: flow accumulation, slope, elevation, rainfall intensity, land use, and geology. The produced maps were georeferenced in the GGRS 1987 Greek Grid coordinate system, and had a cell size of 5 m × 5 m. In order to construct the six thematic maps, the following methods were followed, according to Kourgialas (2010) methodology. The factors were selected according to their general relevance to flood hazards, according to the prevailing literature (Eimers et al. 2000; Yahaya et al. 2010).

First, a DEM of the area was constructed, as the basis of creating multiple of the raster layers. The DEM was created using the contour shapefile of the area, and converting it to an elevation raster, by using the Topo to Raster tool, in ArcGIS.

Following the above, a flow accumulation raster map of the area was created, using the ArcGIS Spatial Analyst toolbox. When performing the calculation, the software deposits a certain amount of water into each raster cell and by factoring in slope, flow direction and elevation, determines the spatial accumulation of the water. The water accumulation points were therefore identified, thus comprising the first thematic map. In conjunction with flow accumulation, Kwak & Kondoh (2008) stated that slope and elevation contribute to flood events and their frequency, specifically they are both inversely proportional to the appearance of flood events.

Continuing with the factors contributing to flash flood risk in small river basins, Belmonte & Beltran (2001) have shown that rainfall intensity and soil moisture conditions are the most important factors. Soil moisture in particular is directly correlated with the land use, as different land uses diversify the soil moisture content and retention, as well as the amount and time of precipitation reaching the soil's surface. Finally, geology of the surface layers of soil contributes to flood risk, as it affects runoff generation. A highly permeable top layer of soil leads to reduced runoff, as more moisture is absorbed vertically through the ground, in contrast to a soil layer with low to no permeability, where the majority of the rainfall is retained on the surface, creating runoff.

The above six factors (elevation, slope, flow accumulation, geology, land use, and rainfall intensity) were the factors selected to be combined in order to create the Total Flood risk map, using the ArcGIS software suite. The reason for choosing these six factors is that they contain the appropriate information with regards to flood risk, while parting with additional, often unnecessary information that may increase complexity and utility of the spatial model (Kourgialas & Karatzas 2011).

To construct the six layers, the original data used included contour maps of the area of Stalos, at a 20 m interval (Figure 3). Utilizing the Spatial Analyst toolbox, a DEM of the area was created, in the form of a raster file, with a cell size of 5 m × 5 m (Figure 4).

Using the DEM, a slope raster was created with the Spatial Analyst Toolbox and the Slope tool. The raster cells contain slope values, ranging from 0 (flat) to 90° (vertical) (Figure 5).

The Flow Accumulation raster map was created by utilizing the DEM, and the slope rasters, in addition to implementing the Flow Direction tool and finally the Flow Accumulation tool. The final raster depicts the areas where the runoff accumulates due to slope and elevation (Figure 6).

For the land use raster, the Corine Land Use survey was the source for data (Figure 7). Since there were multiple types of land use, they were subsequently categorized according to their effect on water permeability to the ground.

Finally, the geological formation data of the area were procured by the Geoenvironmental Engineering Lab of the Environmental Engineering School of the Technical University of Crete. The area mostly comprises of porous formations of varying hydraulic conductivity (P1: medium to high permeability, P2: medium permeability, P3: medium to low permeability), as well as various karst (K1) and impermeable (A2) formations, albeit at a significantly lower percentage (Figure 9).

In addition to the spatial flood risk assessment of the area, the simulation of a stream located in the study area was performed, using hydrologic modeling techniques. Flood events in general are considered as the overflow of water in rivers or streams, due to an increased flow rate, which in turn is caused by heavy precipitation events. Floods are characterized as a violent phenomenon of short duration and usually confined to a specific domain (Goumas 2019). To predict and manage flood events, a statistical analysis of the hyrometeorological data of the study area is necessary, followed by the use of hydrologic and hydraulic simulation models (Bellos et al. 1992). Thus, the first step in understanding flood phenomena and the peak flow is the analysis of spatial and temporal distribution of rainfall. The resulting rainfall time series are transformed to flood peaks or hydrographs.

In a watershed area simulation, especially in areas that experience frequent flood event occurrence, models are mainly used for the determination of vulnerable zones, as well as the determination of rainfall threshold which leads to an extreme event, enabling authorities to design, plan, and construct the appropriate protective structures (Stathis 2004).

For the purposes of this study, the hydraulic model MIKE HYDRO River suite was chosen, due to the options and customizability of its parameters. The suite offers a complete solution for flow modeling in rivers, as the main function for the River module is the determination, construction, and execution of 1-D hydraulic river models, in various scenarios where hydraulic modeling is necessary. For this study, the stream of the study area was modeled, utilizing the data from the water level data logger, installed downstream, as calibration data. The input data were the rainfall time series from the installed meteorological station upstream in the watershed.

In addition to the rainfall and water level data, the data for the spatial characterization of the study stream were also introduced into the model, as depicted in Figure 10, in order to accurately model the riverbed in the MIKE HYDRO River software, after the necessary digitization and processing. Fifteen cross-sections were logged, especially for the downstream part of the stream. For example, Figure 11 shows the analytical cross-section data as it was input into the MIKE FLOOD model interface. Individual markers are placed for the center and the bank points, as well as additional ones for which data were recorded.

The rainfall was inserted as a time series, aggregated to suit the model's time step, in order to match the calibration data of the downstream water level meter (Figure 12).

After successfully implementing the 1-D model of the study area, it was coupled with a 2-D model in order to determine the flood zones of the stream. The input data included the above 1-D model, a 2-D bathymetry input, and a coupling element of the MIKE Flood package. This enables the modeling of lateral flow of water during flood events of excessive rainfall.

Initially, the bathymetry was inserted in the model by converting the DEM into a compatible .dfs2 format, followed by the creation of a bathymetry boundary around the stream indicating the desired simulation area. The selected elevation was 80 m in order to exclude simulation artificially outside that area (Figure 13).

A simple MIKE21 model was implemented, without any additional parameters. Only simulation time and time step were chosen, as were the data output types. The simulation of the flood area takes place during the MIKE Flood simulation, negating the need for a full MIKE21 model setup. The Mike Hydro (1-D) element includes all the simulation data and parameters, with the MIKE21 setup only contributing the necessary 2-D modules for the coupling MIKE Flood model setup. Finally, the total simulation time was a 3-month wet period from 1 October 2019 to 10 January 2020, due to the plethora of rain data. It must be noted here that Western Crete (the greater domain of the study area) being on the West end of the island of Crete, has no defined wet and dry periods. Usually, the majority of rain events occur between October and April of each year.

The coupling, as shown in Figure 14 connects the 1-D simulation with 2-D data, in order to simulate flood inundation in a selected study area, via a 1-D–2-D coupled simulation. In essence, the 1-D simulation is executed and calibrated.

Finally, for the flood inundation simulation, an appropriate time period was chosen. During the study period, the stream did not overflow, a fact confirmed by recorded data as well as observations. To properly estimate the flood zone in a potential/possible flood event, the rainfall data were overestimated and used as an extreme case input for the Flood model.

Regarding the computational parameters of the 1-D model, Manning's n formula was implemented, with a uniform distribution value of 0.05, which was derived after consulting empirical evidence and previous simulations of the region, in addition to fine tuning during calibration. The time step was fixed, at 10 s, instead of using a dynamic option, since the dynamic option yielded out of bounds errors, with the high volume of water that resulted from rainfall spikes.

Calibration was performed manually, by extracting the water level time series from the software, performing aggregation and editing, in order to directly compare them to observed levels from the installed water level station in the stream. The data were compared to the respective time moment, for as much direct comparison as possible. The time series aggregation was performed using the Hydrognomon software, which is a useful hydrological time series editing tool, developed by the National Technical University of Athens, and distributed openly.

Adequate calibration was achieved by observing the R² value of the two time series, and adjustments were made to the scale of the time series, to adjust the percent of surface runoff in comparison to the total rain volume. In addition, small adjustments were performed to the Manning n value to better approximate the true value, when flow seemed inadequate between control points.

The six factors chosen (elevation, slope, flow accumulation, rainfall intensity, land use, and geology) were subsequently reclassified and given five discrete levels of flood risk: very low, low, moderate, high, and very high. However, four of the factors have numerical values (elevation, slope, flow accumulation, and rainfall intensity), and the remaining two (land use and geology) have descriptive values.

For the numerical factors, the values of each one were classified using the Natural Breaks (Jenks) method, which – according to Smith (1986) – determines the classes by finding the adjacent feature pairs between which there is a relatively large difference in data values.

For the descriptive factors of geology and land use, the classification was based on each descriptive values relative contribution to flood risk. In the cases of land use, urban areas posed a high flood risk value, whereas woodlands and areas with dense flora posed a low flood risk value. The geology classification was performed by taking into consideration the permeability of each formation; specifically, there is an inverse correlation between permeability – hydraulic conductivity and flood risk levels.

The classification of each factor's flood risk required calculations, which are presented in Table 1, with the rating of the flood risk levels in points being: very high = 10, high = 8, moderate = 5, low = 2, and very low = 1.

Having completed the calculations, as presented in Table 1, it is evident that the elevation factor contributes the most to the total flood risk, at a percentage of 29.03%, followed by land use and geology at 19.35%. Slope contributed lower to flood risk with a percentage of 12.9% and the factors of rainfall intensity and flow accumulation scored the lowest scores, with a 9.68% for both factors. These percentages represent each factor's contribution to the total flood risk. Each raster map is subsequently reclassified to a scale of 1, 2, 5, 8, and 10, representing the five levels of flood risk, and the result is presented in Figure 16(a)–(f).

To fully evaluate the flood risk of the study area of Stalos, all of the six factors must be combined into one overall flood risk map. Each factor contributes to the final map by the amount of its weight calculated in Table 1, and all the raster maps are combined into the final map, with each of their respective percentages in the last column of Table 1.

The simulation period for the hydraulic model was from October 2018 to January 2019. The rainfall was transformed to surface runoff and used as an input boundary condition, while the validation and calibration of the model was performed by comparing the output at the stream point at which the water level metering station was installed. The extracted water level data closely approximates the measured water level, as is evident in Figure 18, as well as the simulated vs. observed graph in Figure 19.

As is evident, the maximum water level both observed as well as recorded was 0.4 m, well below the 1.6 m bank height of the stream at the narrowest point, where the water level meter was installed. This, coupled with the high R² value, indicates a very low probability of a flood event caused by extreme rainfall. The 1-D modeling results above were after extensive calibration, comparing the simulated water level and the actual water level recorded by the stream water level monitoring station installed upstream. The point at which the simulated water level was extracted from the model is the exact same point where the stream monitoring station is installed, at which chainage point an accurate cross-section was measured in situ, eliminating additional variables and ensuring consistency. Therefore, a significant correlation between the two parameters with minimal variance ensured sufficient calibration for the purposes of this study.

As an output of the 2-D couple model, the floodplane is shown in Figure 20.

The output of the 2-D coupled model was visualized in the GIS software, where the water depth points were displayed in ranges, depicting the potential floodplane. It must be noted that the floodplane depicts a beyond-extreme event, such that is not supported by current rain conditions; nevertheless, it was decided to identify and model an extreme event to depict a potential floodplane and which areas near the stream would be affected. The scale bar in Figure 20 is in meters, as well as the water depth in the legend.

The stream level which was used as input for the 2-D coupled model, was the output of the initial 1-D model, properly validated and subsequently calibrated to real collected data of 2 years, picking the wet periods of the season. The validity of the 1-D model carried over to the 2-D model, which was properly validated. However, there are no historical floodplane measurements and observations in order to properly assess the extent of the simulated floodplane. However, it is noted that in recent recorded flood events of the area, in none of them did the river overflow, therefore validating the assumption that the area is not susceptible to flood events due to river overflow, rather to poor infrastructure.

As mentioned, the final 814 m of the stream were modeled utilizing the coupled model. This is due to inaccessibility to the upper stream chainage and physical blockages that prevent accurate measuring of the stream cross-sections upstream. The stream is poorly if at all maintained, therefore accessibility is near impossible upstream.

The use of GIS software in the presented methodology allows users to determine the spatial extent of flood risk zones, by using readily available data. This enables numerous public organizations and private individuals to assess – even at a preliminary level – the flood risk zones of an area of interest, and determine the factor or factors that mostly affect the resulting flood risk. In addition, the use of GIS software is becoming increasingly common, with open-source GIS platforms being readily available and fully compatible with this methodology, with no cost allocated to software procurement and with extended support, facilitating migration to said platforms. The necessary data, such as altitude maps, land use maps, rainfall and geology are also becoming easier to access, with the implementation of online repositories free to access, enabling stakeholders to follow and implement the presented methodology of this study.

In this study, the spatial flood risk assessment was implemented for the area of Stalos in Nea Kydonia, Chania, Greece, with a specific focus on the Stalos stream and adjacent to it areas. Despite being coastal and to the north, the study area contains mountain ridges to the south, subsequently with lower flood risk. Indeed, the spatial flood risk methodology confirmed a higher flood risk of the coastal areas of the domain compared to the southern mountain areas, due to higher elevation, higher slope and less land uses that contribute to flood risk (e.g. urban areas).

The purpose of modeling the stream with the MIKE HYDRO River software was to accurately examine the rainfall effect to the stream's water level. The studied stream was in the center of the study area, and chosen due to a number of recent flood events (2015 and 2017). Despite the increase of water level during severe rainfall events with long duration, the water level never surpassed 25% of the maximum flow depth, before overflowing, which is an observation confirmed by the observed water level data. This supports the initial hypothesis that the flood reports for the area were not due to overflowing of the stream, rather than inefficient surface runoff management, with the main issue being over urban areas without the appropriate rainwater management infrastructure implementation. Additionally, a 2-D coupled element was introduced, in order to determine the flood plain of the area in cases of extreme rainfall. As a result, the implemented methodology has determined that for the study area there is no flood risk as posed by the current rain fall data, and the flood risk will become apparent for a significant increase in the rainfall amount in the area. However, the study still produces valuable data in the dangers imposed by excessive rainfall and improper implementation of rainfall runoff management infrastructure, as well as the changes proposed in land use, where a more significant portion of land is covered by impermeable structures.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The work presented in this paper was performed and supported by the ERMIS-F Interreg V-A Greece-Cyprus research program (2018–2020), and by the contributing team of the Geoenvironmental Engineering Laboratory of the School of Environmental Engineering - Technical University of Crete, under the guidance of professor Dr. G.P. Karatzas.