Multi-Hazard Susceptibility Assessment Using the Analytical Hierarchy Process in Coastal Regions of South Aegean Volcanic Arc Islands

: Coastal environments are highly recognized for their spectacular morphological features and economic activities, such as agriculture, maritime trafﬁc, ﬁshing


Introduction
Southern European countries, and their low-lying coastal areas, offer Sun, Sea, and Sand (3S), responsible for a significant portion on their total revenues [1]. Coastal tourism is one of the most important economic sectors in Greece, where it accounts for 13% of the country's Gross Value Added (GVA) and 3.8% of employment [2]. Regarding the Greek Population-Housing Census 2021, the South Aegean territory was the only region in Greece that recorded a population increase by 5% [3] from 2011 to 2021. Therefore, the SAVA regions in the context of climate change adaptation, mitigation strategies, and multi-hazard assessment. In terms of geospatial data, this work incorporates for the first-time public road data from Microsoft Bing Maps, combined with the Open Street Map (OSM) data, in order to enhance the accuracy of the transportation network.

Study Areas
The South Aegean Region extends from 24 •  The South Aegean Volcanic Arc islands (SAVA) are distinguished by both subaerial and submerged volcanism and starts from the Saronic Gulf in the west to the Nisyros volcanic field (Figure 1) in the east [40]. The SAVA has been developed over the past 5 million years in the Hellenic Subduction System's (HSS) continental crust due to the northward subduction of the last remnant of the African plate's oceanic crust beneath the active margin of the European plate [41]. Milos and Thira's geological context is dominated by volcanic or volcano-sedimentary rocks, including lavas, dykes, pumice, tuffs, and scoria.
Taking into consideration all the above mentioned hazards, this work introduces an innovative approach that could act as a baseline for multi-hazard approaches and early warning systems in the wider area of the Aegean Sea. This manuscript presents a methodology for integrating and analyzing multiple hazards, in the same geographic area with similar morphological and geological characteristics. Ultimately, the findings of this research could contribute to the effective and holistic management of low-lying coastal regions in the context of climate change adaptation, mitigation strategies, and multi-hazard assessment. In terms of geospatial data, this work incorporates for the first-time public road data from Microsoft Bing Maps, combined with the Open Street Map (OSM) data, in order to enhance the accuracy of the transportation network.
The South Aegean Volcanic Arc islands (SAVA) are distinguished by both subaerial and submerged volcanism and starts from the Saronic Gulf in the west to the Nisyros volcanic field (Figure 1) in the east [40]. The SAVA has been developed over the past 5 million years in the Hellenic Subduction System's (HSS) continental crust due to the northward subduction of the last remnant of the African plate's oceanic crust beneath the active margin of the European plate [41]. Milos and Thira's geological context is dominated by volcanic or volcano-sedimentary rocks, including lavas, dykes, pumice, tuffs, and scoria.

Data Collection
The datasets used for the processing of this study were based on open access products from the Copernicus program and other open-source databases, which are described in the Table 1. Specifically, geological information was obtained from the Hellenic Survey of Geology and Mineral Exploration (HSGME). Additionally, the land use data were processed from the Corine 2018 dataset. The Digital Elevation Model (DEM) was acquired from the National Cadastre in order to extract the required morphological factors for the multi-hazard analysis (Figures 2 and 3). The urban fabric was obtained from previous study [7] and the transportation network from the OpenStreetMap (OSM) and Microsoft [42]. Moreover, Sentinel 2 satellite images were utilized for the calculation of the Normalized Difference Vegetation Index (NDVI) and Bare Soil Index (BSI).

Methodology
According to the following workflow ( Figure 4), a fourfold process was incorporated: (1) the development of a geodatabase; (2) the generation of the multiple hazards individually in terms of landslide, torrential flood, soil erosion and tsunami; (3) the visualization and classification of the selected criteria for each hazard; and (4) the multi-criteria ranking of the derived hazard models utilizing, the Analytical Hierarchy Process (AHP) methodology. The first step was heads-up digitization utilizing high resolution ESRI imagery in an attempt to delineate the bounds of important infrastructure and beaches surrounding the two islands. After that, multispectral Sentinel-2 pictures from 2016 and 2021 were obtained [7] in order to generate the Normalized Difference Vegetation Index (NDVI) and Bare Soil Index (BSI), which are required for soil erosion and torrential flooding approaches.
The components depicted in the right portion of the workflow were incorporated into the relational database, and converted into 10 × 10 m raster datasets. Thus, the findings from each individual hazard model were estimated over the total extent for both islands and then subtracted for the low-lying coastline zone.
In accordance with the applied Multi Criteria Decision Making (MCDM) [53] approach, the AHP was adopted, considering the frequency of the potential natural hazards in Milos and Thira Islands. For the implementation of the AHP approach, all the criteria were reclassified into a single tactical scale from 1 to 5, in order to be comparable to each other and be homogenized. The classification of each criterion was implemented using the natural breaks algorithm (Jenks) method [54], which is one of the most appropriate methods for classification [55].
from the National Cadastre in order to extract the required morphological factors fo multi-hazard analysis (Figures 2 and 3). The urban fabric was obtained from prev study [7] and the transportation network from the OpenStreetMap (OSM) and Micro [42]. Moreover, Sentinel 2 satellite images were utilized for the calculation of the Nor ized Difference Vegetation Index (NDVI) and Bare Soil Index (BSI).

Landslide Susceptibility
Mapping areas susceptible to landslides is an important tool for disaster management and development activities, such as the planning of infrastructure, settlement, and agricultural activities. In this study, the Rock Engineering System (RES) method was implemented [56], which utilizes an Interaction matrix to determine the interactions between different factors, based on their importance. The RES method is a semi-quantitative method usually applied to problems concerning rock engineering by determining the process through which a given factor causes the problem that is being researched [57][58][59][60][61]. Specifically, this is achieved by studying the way any given causal factor interacts with another in a clockwise manner ( Figure 5). Furthermore, this method utilizes an Interaction Matrix (see Section 4.1), to determine the degree of interaction between pairs of causal factors. An interaction value of 1 represents no causal interaction between the two factors, while the higher the value, the easier that one of the factors can affect the other.

Conditioning Factors
The following factors were selected for Milos and Thira: the lithology, the elevation, the slope gradient and curvature, the annual solar radiation, the proximities to tectonic structures, drainage and transportation networks, seismic parameters, and the mean annual precipitation.

• Lithology
Lithological information is one of the most critical factors influencing an increase in the likelihood of a landslide occurrence. Depending on geotechnical, stratigraphic, mechanical or hydrogeological characteristics, it is directly related to a slope's resistance to landslide events [60,62]. Consequently, formations such as flysch are mechanically more unstable and prone to trigger landslide phenomena. Likewise, carbonate rocks, such as limestones, are capable of containing a great amount of water and are very susceptible to dissolution by water and by extensive landslide events. In our case, lavas, dykes and pumice were considered as volcanic rocks, while tuffs and scoriae were grouped as volcano-sedimentary rocks. The Quaternary and more recent sediments for both islands were generally identified as "loose sediments", and scree and mudflows as "Scree & debris flows" ( Table 2).

Methodology
According to the following workflow ( Figure 4), a fourfold process was incor (1) the development of a geodatabase; (2) the generation of the multiple hazards ually in terms of landslide, torrential flood, soil erosion and tsunami; (3) the visu and classification of the selected criteria for each hazard; and (4) the multi-criteria of the derived hazard models utilizing, the Analytical Hierarchy Process (AHP) ology. The first step was heads-up digitization utilizing high resolution ESRI im an attempt to delineate the bounds of important infrastructure and beaches surr the two islands. After that, multispectral Sentinel-2 pictures from 2016 and 2021 tained [7] in order to generate the Normalized Difference Vegetation Index (ND Bare Soil Index (BSI), which are required for soil erosion and torrential floo proaches.    agricultural activities. In this study, the Rock Engineering System (RES) method was implemented [56], which utilizes an Interaction matrix to determine the interactions between different factors, based on their importance. The RES method is a semi-quantitative method usually applied to problems concerning rock engineering by determining the process through which a given factor causes the problem that is being researched [57][58][59][60][61]. Specifically, this is achieved by studying the way any given causal factor interacts with another in a clockwise manner ( Figure 5). Furthermore, this method utilizes an Interaction Matrix (see Section 4.1), to determine the degree of interaction between pairs of causal factors. An interaction value of 1 represents no causal interaction between the two factors, while the higher the value, the easier that one of the factors can affect the other.

Conditioning Factors
The following factors were selected for Milos and Thira: the lithology, the elevation, the slope gradient and curvature, the annual solar radiation, the proximities to tectonic structures, drainage and transportation networks, seismic parameters, and the mean annual precipitation. •

Slope gradient
Numerous studies have shown that areas with steeper slopes and a more distinct topography are significantly more susceptible to landslides. Therefore, the steeper the slope, the greater the likelihood of a landslide [63][64][65]. The slope gradient layer was derived from a resampled DEM with a 10 × 10 m 2 cell size, categorized into five classes ( Table 2).

Slope Curvature
The shape of an area's slopes can control the water accumulation, thus affecting their landslide susceptibility. Shape can be classified into concave or convex, where concave slopes tend to hold more water, thus being more saturated and more unstable than convex slopes. [63,66]. This makes concave slopes more likely to cause a landslide. In general, the more negative a value, the steeper the slope, whereas the closer a value is to 0, the flatter the surface (Table 2).

•
Annual Solar Radiation (WH/m2) Areas exposed to solar radiation for extended periods of time annually tend to lose much of their moisture quickly. This creates porous gaps between the rock's grains which contain no water, which would normally widen them, thus decreasing the slopes stability. Therefore, slopes with high annual solar radiation exposure are less susceptible to landslides [67]. The solar radiation layer was categorized into five classes ( Table 2).

• Proximity to tectonic structures (m)
In most instances, tectonic fabric is directly related to slope failure. Specifically, slopes close to faults are more likely to be affected and mechanically weakened, resulting in a landslide. Landslides are more likely to affect regions close to tectonic structures; therefore, it is essential to analyze this factor. Additionally, a raster of Euclidean distances to faults was calculated, which was then classified into five classes ( Table 2).

•
Distance from the main road network For the creation of distance from the road network, a shapefile of the main road network, combing OSM and Microsoft datasets for the two islands, was utilized, using Euclidean distance. The resultant distance raster was reclassified into five classes ( Table 2). It is important to underline that, when a road is blocked by a landslide event, the entire network is affected, as transportation becomes more difficult, especially in islands, where the road network is limited.

• Mean annual precipitation
From 1971 to 2022, precipitation data were collected from the Climate Atlas of Greece provided by the National Meteorological Service, in order to create a rainfall map. It should be noted that, although precipitation heights are low, a small amount of rain can be sufficient to trigger landslides if the duration of rainfall exceeds a certain threshold [68]. The generated dataset was also divided also into five classes (Table 2).

• Seismic parameters
During an earthquake, a massive amount of stress is caused to the geological formations of the area, thus reducing their stability, leading to landslide events. Since an earthquake can make an area highly unstable [69], a high earthquake density is going to increase the susceptibility of the area drastically. To calculate the earthquake density, an earthquake epicenter inventory for magnitudes greater than 3.5 was collected from the National and Kapodistrian University of Athens from 1900 to 2022, and classified into five categories ( Table 2).
The second seismic parameter considered for landslide susceptibility was earthquake depth [69]. Lower earthquake depths tend to cause a greater negative effect on slope stability. As a result, an earthquake depth layer was created, separated into five classes ( Table 2).

FFPI Methodology
In accordance with the study conducted by Durlević et al. [70], the Flash Flood Potential Index (FFPI) was adopted in an attempt to assess the terrain susceptibility to torrential flooding. This index represents the probability that a flash flood may occur in a specific research region, illustrating its susceptibility and vulnerability to this hazard. This methodology is based on four conditioning factors, in terms of their effect on surface runoff: the terrain slope, the land cover, the soil type and the vegetation density [71,72]. These datasets were classified based on runoff and flash flood potential and assigned relative FFPI ratings using the following Equation (1): where S is the Slope, LC is the Land Cover, ST is the Soil Type, V is the Vegetation density and the variations of n are the respective weights of each factor.

Torrential Flood Factors
Torrential floods and flash floods are natural hazards of great complexity. In a floodbased analysis, the parameters regarding the topography of the study area are considered to be the most crucial.

• Slope gradient
Among all of the potential causes of flash floods, topographical gradient is a key parameter, due to the fact that a steep slope, depending on the angle, can create surface water runoff at higher speeds, leading to high water accumulation at the base of those slopes. Consequently, this increases the likelihood of a flash flood occurring significantly, yet the slopes themselves do not typically become flooded. Instead, the low-lying areas that are located at the bottom of the slope are those that are confronted with the greater repercussions of a potential flash flood. As a result, the locations that have a low slope gradient were assigned with a greater relative FFPI in this study.
When determining the FFPI of the slope gradient, the method of Minea et al. [73] was considered, where slope values were calculated in degrees and then reclassified in five classes, as presented in Table 3.

87
The type of an area's soil is considered as responsible for its ability to either drain the surface water, or allow it to move as runoff. Soils such as clays have extremely low permeability and cause increased runoff, resulting in flash flooding. However, soils such as sands are highly permeable, due to larger gaps between the grains, thus being more easily infiltrated by surface water, and preventing flash flood events [74].
The mechanical and hydrological properties of the soil type of a study area control its susceptibility to torrential floods in different ways. Due to the lack of soil maps, lithological datasets were used in an attempt to convert the lithology into soil categories. The correlation between the simplified lithologies and the different soil types was based on research. For example, volcanic or volcano-sedimentary rocks were considered as Andosols, loose sediments and deposit varieties were classified as Alluvials, carbonate rocks were connected to Regosols, and sedimentary rocks, such as sandstones, were considered as Cambisols. Specifically, in this work, the Alluvials are considered the most susceptible soil type, while andosols present the lowest susceptibility, due to increased runoff. As a result, the soil type was also divided into five distinct groups (Table 3).

•
Land cover/Land use The use of land can affect the ability of the ground surface to contain water or cause runoff, similar to the effect of soil type [75][76][77]. Anthropogenic artificial surfaces such as roads, ports or airfields can increase the impermeability of the ground significantly, thus causing more runoff and allowing easier flooding [70,74]. Likewise, areas with little to no vegetation have reduced ability to absorb water, as runoff travels at extreme speeds, thus causing flash floods before it can infiltrate the ground. Vegetated regions, on the other hand, have a natural slowing effect on surface water flow, due to the fact that they allow water ground infiltration, preventing flooding. The necessary data were collected from Corine (2018) ( Table 3).

• Vegetation density
One of the most important factors that may prevent the occurrence of a flash flood is the existence of vegetation [78]. Particularly, trees and forests, which are known for their stability, can decrease the rate of surface water runoff when they come into contact with it. This becomes even more apparent in areas where there is a higher density of plants and trees [76]. The vegetation density layer was created using the method suggested by [70], in which the Bare Soil Index (BSI) of the study area is calculated in order to determine the areas covered by vegetation by separating the bare areas [79]. This was accomplished using satellite images derived from Sentinel 2, from 2016 and 2021, using the following Equation (2): where B11 is the Short-wave infrared band, B4 is the Red band, B8 is the Near Infrared band and B2 is the Green band of the Sentinel-2 satellite system. The following Equation, suggested by [80], was utilized to calculate the vegetation density (V) based on the BSI index.
The results were categorized into five classes, as represented in the following Table.  Table 3. Torrential flood causal factor classification and ranks in Milos and Thira islands.

Soil Loss Susceptibility Assessment-Revised Universal Soil Loss Equation (RUSLE)
A great variety of soil erosion models have been developed over the years, used in the estimation of soil loss and its impact, along with the increase in land conservation efficiency. One of the most commonly utilized methods of soil erosion rate prediction is the Revised Universal Soil Loss Equation (RUSLE), first published in 1991, a modified version of the original Universal Soil Loss Equation (USLE). The main difference and advantage of RUSLE is the ability to estimate the soil loss caused by both surface and rill erosion [81][82][83][84][85][86]. The RUSLE model estimates the (annual) average soil loss (A) in tons per hectare per year (t ha −1 y −1 ), based on six major numerical factors (Equation (4)): where R is the rainfall erosivity, K is the soil erodibility, LS is a factor that usually combines the length (L) and the steepness (S) of slope, C expresses the land cover and P is the erosion control practice factor [85].

RUSLE Factors
Soil erosion is a natural hazard like many others, but unlike other risks, it is a dependent variable, making investigation challenging. To evaluate the vulnerability of Milos and Thira to soil loss, the following factors were investigated. •

Rainfall (R) factor
The R factor in the RUSLE equation is used to determine the kinetic energy of the rainfall intensity [87,88]. An empirical Equation was implemented due to the lack of intensity data availability, in order to calculate the R factor.
where P is the mean annual rainfall in mm of the two study areas. This equation requires long-term continuous rainfall data from meteorological stations within the study area. In addition, this study adopts the approach of [89] by using the mean annual rainfall from monthly data to reduce the rainfall distribution variation.

• Soil Erodibility (K) factor
One of the main issues in estimating soil erosion is the lack of data on soil features. K-factor (MJ-1 mm −1 ), as used in the soil erosion model, the Universal Soil Loss Equation (USLE) and its revised version (RUSLE), can be used as a key factor to model soil, based on soil erodibility. More particularly, this factor, which expresses the susceptibility of a soil to erosion, is related to soil properties such as organic matter content, soil texture, soil structure and permeability. The K factor calculation based on LUCAS datasets is provided from JRC [52].

• Topographic Slope Length & Steepness (LS) factor
The LS factor indicates the effect of topography on soil erosion and describes the slope length (L) and the steepness (S) parameters. The former (L) is defined as the point of departure of the surface runoff to the point, while the latter (S) describes the behavior of soil erosion with an inclination angle [83].
Ref. [90] developed an approach that combines both parameters to estimate the LS factor. In this way, input data is distinguished in the upslope contributing area per unit width, determined by the flow accumulation, the pixel size, and the slope, while the outcome is unitless (Equation (6)).
where U is the flow accumulation multiplied by the pixel size, L0 is the slope length (22.13 m), β is the slope in degrees, S0 is the slope percentage (9%), m is sheet erosion ranging from 0.4 to 0.6, accordingly to the slope intensity, and n is rill erosion ranging from 1 to 1.3. The values of the indicators used are m = 0.5 and n = 1.2 [87].
• Cropping and Land-Cover (C) factor The C factor describes the influence of soil cover, crop and management territorial loss in relation to the territorial loss in bare fallow land areas [86]. It is a dimensionless factor ranging from 0, in case of intense coverage, to 1 when there is significant lack of coverage in vegetation.
Due to the important relationship between soil and vegetation coverage, Durigon et al. [91] provided the following Equation (7) to describe the C factor based on the Normalized Difference Vegetation Index (NDVI): • Conservation Practices (P) factor The support Practice (P) factor is defined as the ratio of soil loss in a specific soil conservation practice to a rough field. Thus, factor P is important to take into consideration, as it can provide information about what practices are more advantageous for soil conservation [84].
Concerning the values of the P factor, these range from 0.2 to 1, with low values corresponding to greater control of soil erosion. P values can be extracted from the literature, while in other cases P is empirically considered. For example, a value of 1 for the P factor means that either there are no support practices or there are conventional techniques. On the other hand, a value of 0.25 shows the potential for this management factor to reduce soil by a 75% loss [84]. P factor can also be estimated based on the Corine Land Cover (CLC) data. According to Yang et al. [92] for Corine land categories (211, 212, 221-223, 231, 241-243), P is estimated as 0.5, which according to David [93] corresponds to "minimum plowing", while all other classes were given the value of 1.

Tsunami Run-Up Scenario
Tsunami events are extremely destructive natural hazards, which may occur rarely, but cause extreme loss of human life, as well as of structural properties. This makes the estimation of the parts of a coastal zone more susceptible to tsunami events a necessity when studying the hazard susceptibility of an area. It is worth mentioning that tsunamis can be generated from underwater earthquakes, volcanic eruptions, landslides, or extraterrestrial impacts such as asteroids. In our work, the tsunami run-up model was based on seismogenic tsunamis in an attempt to visualize the worse-case-scenario, as suggested by the research of Batzakis et al. [38]. This model was created using tectonic plate seismic data based on the earthquake magnitude most associated with high tsunami run-up. This specific threshold was suggested by Iida (1963) as approximately Mw~8.5 [94]. These magnitude values were associated with a wave of 10m mean maximum run-up, which has been confirmed by several studies in the coastal zones of the Cyclades islands [38]. Through the application of a GIS, the elevation pixel values of the utilized DEM were transformed into inundation depths, divided into five categories using natural breaks (Jenks).

Analytical Hierarchy Process
The total susceptibility of coastal regions to natural hazards was estimated by using the Analytical Hierarchy Process (AHP) method. The AHP is a decision-making analytical methodology, which excels in solving multi-criteria problems [95,96]. In our work. the four hazards represent the problem variables.
As a second stage, it was necessary to compare these variables in pairs (pair-wise comparison) based on their impact on the overall problem. This pair-wise comparison approach places the separate hazards in a hierarchy, via a comparison matrix [96]. Once the comparison matrix is completed, the weights of importance for each hazard value were calculated, which represent their influence. This was performed by dividing the geometric mean of each line (u i ) by the sum of the geometric means of all rows (u k ) of the matrix, via the following Equation (8): The consistency of the pairwise comparisons of the square matrix can be evaluated via the Consistency Index (CI) and the Consistency Ratio (CR). The CI was calculated via Equation (9): where λ max is the largest eigenvalue of the matrix and n is the order of the matrix. Meanwhile, the CR is calculated through the following Equation (10): where the RI is the randomness index value that depends on the order of the matrix published by Saaty [95], and CI is the Consistency Index. The CR needs to be lower than 0.1 and realistically higher than 0 in order for the comparison to be considered consistent, and thus successful. In this work, CR value was calculated with an acceptable consistency at 0.01. The calculated weights were used as multiplying factors for each classified criterion, in order to create the final susceptibility map.

Results
In this study, four criteria corresponding to the most common hazards on both islands were selected for homogenization and incorporation into the presented AHP implementation.

Susceptibility to Landslides
Applying the RES to determine the susceptibility of Milos and Thira islands to landslides, a 10 × 10 interaction matrix (Table 4) was implemented. It is important to note that the same comparison ranks were used for the causal factors of both islands due to their similarities, regarding their interaction intensity (C+E). The factors were categorized as P1 lithology, P2 mean annual precipitation, P3 slope gradient, P4 solar radiation, P5 earthquake density, P6 earthquake depth, P7 distance from the road network, P8 proximity to tectonic structures, P9 proximity to the drainage network and P10 slope curvature.  After determining the interactions between the 10 causal factors, a sample of random landslides equal to 70% of the landslide inventory was used to calculate the weights of importance for each factor. The landslide inventory was mainly obtained using optical observations in Google Earth Pro and bibliographic data [63,97]. A total of 162 landslides were observed in Milos and Thira, respectively, within the boundaries of the low-lying coastal areas.
The remaining 30% were used to validate the results of the landslide susceptibility model. The C+E% values were divided by the maximum susceptibility value (Max Pij) of each factor in order to calculate the necessary weights (Wi) ( Table 5).
Once the weights of importance were calculated, the factor layers were multiplied by their respective weights and then added together to generate the landslide susceptibility models.
Based on the RES analysis and processing of the existing landslides inventory in GIS, the landslide susceptibility maps of Milos and Thira Islands ( Figure 6) were generated as illustrated in the following figures.
Due to the similar morphology of both islands, which is characterized by relatively steep terrain slopes, landslides are the most common threat. However, only a small percentage of the entire low-lying coastal zones fall within the high and very high susceptibility classes. In particular, on the island of Milos, the high (13.79%) and extremely high (6.25%) susceptibility classes cover 19%, while on the island of Thira they cover only 13% (Table 6).  Due to the similar morphology of both islands, which is characterized by relatively steep terrain slopes, landslides are the most common threat. However, only a small percentage of the entire low-lying coastal zones fall within the high and very high susceptibility classes. In particular, on the island of Milos, the high (13.79%) and extremely high (6.25%) susceptibility classes cover 19%, while on the island of Thira they cover only 13% (Table 6).
Specifically, the beaches of Agia Kyriaki (Ag. K), Filotas (Fi), and Palaiochori (Pal) on the southeast side of the island of Milos were the most vulnerable. Concerning the island of Thira, areas of high threat were identified in the southern parts of the coastal region in the Perissa (Pe) and Kamari (Ka) settlements, as well as Red beach (Rb) (Figures 7 and 8).  Specifically, the beaches of Agia Kyriaki (Ag. K), Filotas (Fi), and Palaiochori (Pal) on the southeast side of the island of Milos were the most vulnerable. Concerning the island of Thira, areas of high threat were identified in the southern parts of the coastal region in the Perissa (Pe) and Kamari (Ka) settlements, as well as Red beach (Rb) (Figures 7 and 8).

Susceptibility to Torrential Floods
Once all the conditioning factors were processed into a common scale, the four thematic layers were combined through map algebra, using the following Equation (11): (Slope + Land cover + Soil type + Vegetation density ) 4 (11) Due to their significant influence on the coastal area, the causal factors studied were considered equal when calculating the FFPI values of the islands. The resultant values were on a scale of 1 to 5, representing the increasing flash flood potential and overall torrential flood susceptibility of the study areas (Figure 9).
A large part of the coastal region of Milos Island has a high (24.04%) and very high

Susceptibility to Torrential Floods
Once all the conditioning factors were processed into a common scale, the four thematic layers were combined through map algebra, using the following Equation (11): FFPI = (Slope + Landcover + Soiltype + Vegetationdensity) 4 (11) Due to their significant influence on the coastal area, the causal factors studied were considered equal when calculating the FFPI values of the islands. The resultant values were on a scale of 1 to 5, representing the increasing flash flood potential and overall torrential flood susceptibility of the study areas (Figure 9). GeoHazards 2023, 4, x FOR PEER REVIEW 19 of 33 Figure 9. Torrential flood susceptibility maps in Milos (a) and Thira (b) islands. Figure 9. Torrential flood susceptibility maps in Milos (a) and Thira (b) islands.
A large part of the coastal region of Milos Island has a high (24.04%) and very high (50%) susceptibility to the occurrence of torrential floods (Table 7). In the north and west of the island, very low, low, and moderate susceptibility are predominant. In terms of the onshore regions of Thira Island, the eastern and southern regions have high (10.78%) and extremely high (19.50%) susceptibility values. The medium susceptibility class has the highest coastal coverage at 59.89%, followed by the very low and low classes at 5.11% and 4.71%, respectively.

Soil Loss
The soil erosion model was created by using the RUSLE equation to combine the soil loss causal factors. After generating the required layers, the result was a soil loss map of Milos and Thira's coastal zones ( Figure 10).
The erosion susceptibility maps ( Figure 10) depict low potentiality in both islands. Although, most of the coastal areas are characterized by low erosion that ranges between 94.94-95.28% (Table 8), the high and very high susceptibility values were identified near to the drainage network [89].

Susceptibility to Tsunami
As stated by Batzakis et al. [38], the inundation depth of the zones required to be estimated in order to develop the tsunami susceptibility model. The results were limited to a maximum elevation of 10 m, because this is considered the maximum possible tsunami run-up, and classified into five categories. The susceptibility maps (Figure 11) illustrate the coastal regions that are more vulnerable to tsunami. The high and very high classes, in particular, cover up to 20% of Milos Island and approximately 26% of Thira Island (Table 9). Specifically, in Milos Island, the high susceptibility classes cover the central onshore areas at Alikes (Al), Adamas (Ad) settlement and Rivari Lagoon (Ri). as well as on the northern side in Papafragas (Pap) and Polonia (Po) settlements. In Thira Island, the most susceptible areas are indicated on the southeastern sides in Exomitis (Ex), Perissa (Pe) and Kamari (Ka) settlements.

Soil Loss
The soil erosion model was created by using the RUSLE equation to combine the soil loss causal factors. After generating the required layers, the result was a soil loss map of Milos and Thira's coastal zones ( Figure 10). The erosion susceptibility maps ( Figure 10) depict low potentiality in both islands. Although, most of the coastal areas are characterized by low erosion that ranges between 94.94-95.28% (Table 8), the high and very high susceptibility values were identified near to the drainage network [89].

Susceptibility to Tsunami
As stated by Batzakis et al. [38], the inundation depth of the zones required to be estimated in order to develop the tsunami susceptibility model. The results were limited to a maximum elevation of 10 m, because this is considered the maximum possible tsunami run-up, and classified into five categories. The susceptibility maps (Figure 11) illustrate the coastal regions that are more vulnerable to tsunami. The high and very high classes, in particular, cover up to 20% of Milos Island and approximately 26% of Thira Island (Table 9). Specifically, in Milos Island, the high susceptibility classes cover the central onshore areas at Alikes (Al), Adamas (Ad) settlement and Rivari Lagoon (Ri). as well as on the northern side in Papafragas (Pap) and Polonia (Po) settlements. In Thira Island, the most susceptible areas are indicated on the southeastern sides in Exomitis (Ex), Perissa (Pe) and Kamari (Ka) settlements.   Figure 11. Tsunami susceptibility maps in Milos (a) and Thira (b) islands. Figure 11. Tsunami susceptibility maps in Milos (a) and Thira (b) islands.

Total Susceptibility to Hazards
The AHP methodology was utilized to determine the overall susceptibility to multiple natural hazards using a 4 × 4 comparison matrix (Table 10), in order to calculate the weights of importance for each criterion (hazard) using Equation (12). According to the applied scenario, landslides and torrential floods are the most prevalent threats, with landslides being the most prevalent. This is because both islands are frequently affected by landslides, as evidenced by the large number of past landslide events in the two islands, as well as by the numerous debris and mud flows mapped on the geological maps. In addition, although less frequent than landslides, flood events have been recorded, especially in Thira in recent years, making this a critical hazard. As the coastal environment is subject to a constant barrage of erosion processes, soil loss is a greater threat than severe seismogenic tsunamis, due to their high return period. According to [98], the return period values for shallow earthquakes of magnitude M ≥ 6.0 in the South Aegean is more than 200 and/or 400 years. Particularly, landslides were assigned with the greatest weight of importance value (0.4669). Intense precipitation is rare in the Cyclades, but in the context of climate change, the islands are prone to intensive torrential flood events that can cause damage, especially in low-lying areas. This hazard was assigned a weight value of 0.2776. Moreover, soil erosion was assigned with 0.1603 and tsunami with 0.0953. Based on the obtained weight values (Equation (12)), total susceptibility to the assessed hazards was calculated, utilizing the following Equation: where L represents landslides, F torrential floods, S soil erosion and T tsunami. The results (Figures 12 and 13) were reclassified into five classes using the natural breaks (Jenks) method. Green color illustrates the areas with the lower susceptibility score and red color the higher, respectively. GeoHazards 2023, 4, x FOR PEER REVIEW 23 of 33   The results indicated that the highly scored areas were located in the Adamas (Ad) settlement, Sarakiniko (Sa) beach and the Rivari Lagoon (Ri) (Figure 14).  Particularly, the analysis of the results for Milos illustrated ( Table 11) that 25.85% (0.24 km 2 ) of the built-up areas, 56.46% (0.05 km 2 ) coverage of the ports, 3.77% of the airport, 35.48% (0.35 km 2 ) of beaches and 20.44% (0.05 km 2 ) of the main transportation network are characterized by high susceptibility (Figure 14). Regarding the results for Thira Island, the highly scored areas were identified on the west side, where geomorphology is characterized by high relief slope. Furthermore, high values were also detected in the Perissa (Pe) settlement, which is located on the northeast side of Thira ( Figure 15). Specifically

Result Validation
The validation of the total susceptibility is a complex problem; therefore, the models' validation was accomplished by validating each hazard individually. Based on 30% of the landslide inventory, the percentage success of RES methodology was 97.37%.
In the case of susceptibility to torrential floods, areas with high and very high zones were confirmed from recent flood events near to the airport area and Kamari settlement. There was no data for validation on Milos. However, recent news from online local newspapers have announced the establishment of new anti-flooding works in Milos for the upcoming period.
Despite the fact that the soil loss values are low, the higher values were identified in specific areas and especially along the existing drainage network, following the pattern of the LS factor. Particularly, high erosion values were detected between 2016 and 2021 in Kamari (Figure 16c,d) settlement (eastern Thira) and adjacent to the southeastern part of the airport (Figure 16a,b).

Result Validation
The validation of the total susceptibility is a complex problem; therefore, the models' validation was accomplished by validating each hazard individually. Based on 30% of the landslide inventory, the percentage success of RES methodology was 97.37%.
In the case of susceptibility to torrential floods, areas with high and very high zones were confirmed from recent flood events near to the airport area and Kamari settlement. There was no data for validation on Milos. However, recent news from online local newspapers have announced the establishment of new anti-flooding works in Milos for the upcoming period.
Despite the fact that the soil loss values are low, the higher values were identified in specific areas and especially along the existing drainage network, following the pattern of the LS factor. Particularly, high erosion values were detected between 2016 and 2021 in Kamari (Figure 16c,d) settlement (eastern Thira) and adjacent to the southeastern part of the airport (Figure 16a Finally, the tsunami susceptibility model depicted that low-lying coastal areas on both islands are vulnerable to tsunami, especially the areas adjacent to the coastline and with lower elevation values. Specifically, Kamari, Perissa and the eastern part of Thira are more susceptible to tsunami hazard. It is worth mentioning that Kamari settlement was developed after the 1956 earthquake, due to the relocation of the population from other settlements of the island that were destroyed by this very strong earthquake, following the initial Hellenic Anti-Seismic Regulation in 1959 [99]. Moreover, terrain vulnerability to tsunami in Milos is higher in Adamas (central port), Alikes (central Milos) and the southeastern part, where high tourism activity and marine activities are concentrated annually.

Discussion and Conclusions
The SAVA region is one of Europe's most seismic and volcanic landscapes, making it a one-of-a-kind physical laboratory for researchers and visitors. Additionally, reports have shown that the Aegean Sea will be one of the most vulnerable areas in the Mediterranean basin to global warming and climate change scenarios, with more extreme meteorological events, heat waves and droughts. Apart from that, the Cyclades complex and the South Aegean region was the only area in Greece with a population rise of more than 5% from 2011 to 2021. Thus, the study areas are a hub of human-caused activity, especially during the tourist season, which starts from April and ends in October, annually.
This work introduces a workflow to determine the most susceptible regions, in terms of multi-hazard susceptibility, for low-lying coastal areas in Thira and Milos islands. A multi-hazard methodology is a complex work, due to the fact that single hazards with different approaches must be compared in a specific period and specific geographic  Figures  (a,b) correspond to the airport area, while (c,d) depict the high urbanization rate in the Kamari settlement.
Finally, the tsunami susceptibility model depicted that low-lying coastal areas on both islands are vulnerable to tsunami, especially the areas adjacent to the coastline and with lower elevation values. Specifically, Kamari, Perissa and the eastern part of Thira are more susceptible to tsunami hazard. It is worth mentioning that Kamari settlement was developed after the 1956 earthquake, due to the relocation of the population from other settlements of the island that were destroyed by this very strong earthquake, following the initial Hellenic Anti-Seismic Regulation in 1959 [99]. Moreover, terrain vulnerability to tsunami in Milos is higher in Adamas (central port), Alikes (central Milos) and the southeastern part, where high tourism activity and marine activities are concentrated annually.

Discussion and Conclusions
The SAVA region is one of Europe's most seismic and volcanic landscapes, making it a one-of-a-kind physical laboratory for researchers and visitors. Additionally, reports have shown that the Aegean Sea will be one of the most vulnerable areas in the Mediterranean basin to global warming and climate change scenarios, with more extreme meteorological events, heat waves and droughts. Apart from that, the Cyclades complex and the South Aegean region was the only area in Greece with a population rise of more than 5% from 2011 to 2021. Thus, the study areas are a hub of human-caused activity, especially during the tourist season, which starts from April and ends in October, annually.
This work introduces a workflow to determine the most susceptible regions, in terms of multi-hazard susceptibility, for low-lying coastal areas in Thira and Milos islands. A multihazard methodology is a complex work, due to the fact that single hazards with different approaches must be compared in a specific period and specific geographic boundaries. Our pair-wise comparison and applied methodologies for terrain susceptibility to landslide, torrential flood, soil erosion, and tsunami yielded the following results.
Initially, the landslide susceptibility map was derived by implementing and integrating the RES approach. Especially, for this hazard a landslide inventory was created in order to validate the model on both volcanic islands. The calculated results illustrated that more than 55% of the terrain was found within the very low and low susceptibility zones, and more than 10% was found within the high and very high susceptibility zones (Table 6). In particularly, the visualized outcomes highlighted those southeastern areas of Milos, such as Ag. Kiriaki, Filotas, Palaiochori and the western part of Adamas settlement, as the most landslide-affected areas ( Figure 6). Additionally, the areas most prone to landslides are located in southeastern regions of Thira, such as Perissa, Kamari, Red beach, Athinos port (central region), and Oia's desalination plant (northwest side). These findings were driven due to the interaction and weighting analyses of 10 parameters, such as seismicity, lithological conditions, terrain curvature, soil radiation, etc.
According to Figure 9, the high and very high flood susceptibility zones are located in low-lying regions. The results illustrated that more than 10% of the terrain was found within the very low and low susceptibility zones, and more than 70% (Milos) and 30% (Thira) were found within the high and very high susceptibility zones ( Table 7). The northern areas of Milos, such as Sarakiniko, Papafragas and Pollonia settlements, central regions, e.g., Adamas and Alikes, and the western area of Firiplaka settlement are occupied by agricultural activities (Figure 9). Additionally, the most flood-affected areas are located in southeastern regions of Thira such as Perissa, Kamari and Exomitis settlements (Figure 9).
The calculated values were quite low on both islands in terms of soil loss assessment. Specifically, more than 90% was mapped within the low and very low susceptibility zones. In contrast to these results, regions with high soil loss values were identified in the vicinity of the existing drainage network following an LS pattern. According to Figure 10, prone areas were identified in southeastern regions such as Agia Kiriaki, Paliochori, Spathi and Thiorichio. Moreover, in Thira, eastern and southeastern areas were found to be more susceptible to flood, such as Exomitis, Perissa, Kamari, and areas near to the southeast side of the airport.
The tsunami susceptibility map was generated based on the conversion of low elevation values up to 10 m. Coastal regions of the Cyclades islands complex have been affected historically by a run-up of 10 m. According to the worse-case run-up scenario of 10, the results illustrated that more than 55% of the terrain was found within the very low and low susceptibility zones in both case study areas (Table 9). In particular, the visualized outcomes highlighted the central area of Milos, such as Adamas and Alikes, and the northeastern side, such as Pollonia settlement ( Figure 11). Additionally, the tsunami-affected areas are located in eastern and southeastern regions of Thira, such as Perissa, Kamari and Exomitis, and adjacent to the airport facility ( Figure 11).
According to Figures 12 and 13, two multi-hazard susceptibility maps were produced within the low-lying coastal limits. Although the percentage of the most vulnerable regions on both islands is less than 25%, the majority of ports and transportation networks are located in the high and very high susceptibility zones. According to Table 11 and the adopted scenario, high and very high zones include more than 30% of the built-up area, 20% of the transportation network, and 50% of the seaports. Notable is the fact that 33% of Thira's airport and more than 30% of Milos's beaches are vulnerable to severe hazards. Further to that, the aforementioned findings highlight the necessity for scenarios and potential threats modelling, during a potential hazard event in low-lying areas.
This work has some limitations due to the fact that it had to compare four specific hazards with selected criteria. Furthermore, the ranking of the criteria in AHP analysis was based on the appearance frequency of the hazards. Finally, this study presents only one adopted scenario, considering landslide susceptibility as the most important hazard.
Summarizing, this study highlights the importance of using earth observation and geo-information methods to investigate hazard assessment in low-lying coastal areas, on a regional scale. Future research could be based on multi-hazard integration with other hazards, such as volcanic activity and hydrothermal investigation in southern Milos, where the EU Pathfinder project RAMONES will research the offshore area in 2023 using innovative prototype instruments. Therefore, advanced remote sensing methods in terms of ground subsidence (e.g., InSAR methodology) for shoreline changes utilizing the Digital Shoreline Analysis System (DSAS) will be a part of future work, in order to investigate potential correlation. The identification of more vulnerable areas to hazards using multihazard approaches could be used as a baseline for civil protection, regional planning and decision making. Thus, the results of this work could be combined with data from innovative technological instruments such as High-frequency (HF) radar that monitors the sea surface velocity and provides information for tsunami early warning and SLR [100].