High-resolution susceptibility mapping of seismically induced landslides on Ischia island: the 2017 earthquake case study

Abbreviations: CLC: Corine land cover; DEM: digital elevation model; ERT: electrical resistivity tomography; PGA: peak ground acceleration; PVD: Phlegraean Volcanic District; SBAS: small baseline subset; SNAP: seismically-induced landslides susceptibility

1.   Introduction

Ischia island is an active volcano located on the Gulf of Naples (Campania Region, Italy). The area has been affected by intense local seismic activity since prehistoric times ^[1,2,3,4,5]. The seismicity of Ischia is caused by regional tectonics and magmatic dynamics ^[6,7,8] and is characterized by shallow hypocenters (<2 km depth) located mainly in the NW portion of the island. Detailed documentation about earthquakes and related induced effects is available starting from the 13th century ^[9,10]. Considering the seismicity from 1800 to nowadays, the most catastrophic events occurred during the 19th century, on March 4, 1881, and July 28, 1883, with a Mercalli-Cancani-Sieberg (MCS) macroseismic scale of IX and XI, respectively. Those events were particularly destructive as they wiped out the Casamicciola Terme village and its surroundings, causing 120 deaths and 140 injuries in 1881, and 2343 victims in 1883. During the 20th century, only low-magnitude seismicity has occurred. The most recent events were on August 31, 2022, with an M_d of 2.1, April 16, 2020, with an M_d of 2.3, and August 8, 2018, with an M_d of 2.5; the most intense earthquake was on August 21, 2017, with and M_d of 4.0. The epicenter was located once again in the north portion of the island, near the Casamicciola Terme village, with a shallow hypocenter of 1.2 km depth ^[11]. The mainshock was followed by approximately 30 small aftershocks with the most detectable ranging from M_d 0.7 to 1.9. Two victims and significant damage to buildings were caused by the mainshock ^[12,13]. This event also generated several coseismic geological surface effects such as at least 6 landslides and 20 fractures ^[9].

The 2017 coseismic effects, although having been the first observed in recent times, are not the only documented. In the past centuries, the identification of seismically induced landslides was based on the analysis of historical reports, paintings, and archaeological studies that allowed the reconstruction of their chronologies from the 4th century BC to the beginning of the 21st century ^[14]. In this context, the most significant landslides that have been triggered by earthquakes occurred during the 2nd and the 3rd centuries and those dated 1228, 1828, 1863, 1881, and 1883 ^[14].

However, the landslides that affect Ischia are not only induced by earthquakes but especially by rainfall. They can be triggered by very intense or prolonged precipitation, as occurred on November 26, 2022, when a devastating landslide was triggered by 126 mm in 6 hours of rainfall ^[15]. This landslide originated from Mt. Epomeo, crossed the northern side of Casamicciola Terme village, and reached the seashore causing the loss of 12 human lives, the destruction of numerous buildings, the interruption of roads, and the evacuation of over 30 families. This more recent disaster further demonstrates the necessity to develop or update landslide susceptibility maps in order to prevent and mitigate the risks linked to such calamities, independently of the triggering factor.

In the literature, most works have focused on landslides triggered by intense or prolonged rainfall ^{[16,17,18,19,20]}, whereas only a few address seismically induced landslides. In ^[8], based on seismically induced landslides that occurred during the 1883 earthquake, a first susceptibility map was produced covering the entire island with a spatial resolution of 40 m. The map was obtained considering three predisposing factors (the type of outcropping rocks/soils, the terrain slope angle, and the MCS intensity) stated as a significance percentage. In ^[21], a probabilistic approach was proposed integrating topographic and geo-lithological data with seismic input characterized by different annual exceedance frequencies. Finally, a more recent study ^[22] investigated the possibility of identifying areas potentially affected by earthquake-induced landslides by analyzing the coseismic InSAR data of the earthquake of August 21, 2017.

In this work, a high-resolution seismically induced landslide susceptibility map is produced. In particular, considering as the triggering event the earthquake of August 21, 2017, the map is the result of a multidisciplinary approach that integrates and combines, in GIS platform, datasets of different typologies recognized in the literature as predisposing factors to terrain instability. The 2017 earthquake can be considered a good case study as most of the recent and historical seismic events epicenters are localized in the same area as the 2017 one ^[23] and its intensity is representative of the seismicity of the studied area. The susceptibility map is discussed and compared with recent and historical seismically induced ground effects consisting of coseismic fractures, ruptures and cracks, and landslides phenomena ^[12,24,25]. In addition, here for the first time, a test area has been investigated with an electrical resistivity tomography field survey. This geophysical prospecting confirms the high susceptibility to mass movements as a result of our map. Finally, the susceptibility map is combined with 7 years of InSAR temporal series to identify the most predisposed areas of the island.

2.   Geological and volcanological framework

Ischia island represents the westernmost portion of the Phlegraean Volcanic District (PVD) constituted by Campi Flegrei caldera, Procida, and Vivara islands (Figure 1). The island covers an area of approximately 42 km² and its highest point is the Mt. Epomeo (786 m a.s.l.), a volcano structural block that, after the formation of the "Green Tuff" caldera depression, was affected by an uplift of 800–1100 m a.s.l. in the last 33 ka ^[26]. This phenomenon, defined as "resurgence" ^{[26,27,28,29,30]}, represents an almost unique example in volcanic areas ^[10]. The resurgence block boundary assumes a polygonal shape ^[30,31,32] defined by the intersection of NW–SE and N–S faults in the western sector, and NE–SW, E–W faults in the northern sector (Figure 1). The reactivations of regional faults and the generation of new faults due to volcano-tectonic activity tend to reshape the resurgent area ^[33,34]. In fact, the structural evolution of Ischia was defined in ^[35] as a volcano-tectonic horst formed by shallow crustal magmatic intrusion causing the uplifting of Mt. Epomeo ^{[28,29,33,36,37]}. Also, in the eastern sector, N–S and NE–SW structural lineaments are present, corresponding to a series of eruptive vent alignments ^[31].

Figure 1.  Sketch map of the Phlegraean Volcanic District. The historical landslides affect areas represented in yellow (from IFFI catalog, ISPRA 2006); the faults are represented in red lines ^[27,29].

DownLoad: Full-Size Img PowerPoint

From a volcanological point of view, the eruptive activity of Ischia island started more than 150 ka ago ^[10]. The oldest volcanic products, dated 75–150 ka, are represented by lava flows, lava domes, and a sequence of pyroclastic deposits with intercalated paleo soils, mainly exposed along the coastline in the southern part of the island ^[36]. Approximately 55 ka ago, the trachytic ignimbrite eruption of "Tufo Verde" created the calderic depression of Mt. Epomeo ^[26]. Subsequently, the caldera started collapsing with the deposition of "Tuffite" and "Colle Jetto" marine sequence in the central part of the island ^[27]. Between 44 ka and 33 ka ago, the western coast was affected by the "Citara Tuff" eruption, while the central part of the island experienced the resurgence of Mt. Epomeo. Subsequently, periods of quiescence alternated with periods of activity up to 2.9 ka, when the recent volcanism of Ischia began with effusive, extrusive, and explosive activities ^[38]. The most recent volcanic activity was the Arso eruption, dated 1302, and mainly featured by lava flows affecting the eastern sector of the island ^[36,39]. The geomorphological and volcanological aforementioned characteristics make Ischia an island particularly prone to slope instability. In particular, fundamental roles are played by the steep morphologies characterizing most of the island surfaces and the different typologies of volcanic outcropping terrains such as pyroclastites, scoriae, tuffs, and lava flows ^{[7,40,41,42,43]}. In detail, the effusive and explosive eruptions ^[7,30] generated volcanic deposits (either soils and rocks) with various degrees of permeability and erodibility that, in case of earthquakes, intense precipitation, marine erosion, or human actions, can be easily remobilized triggering different typologies of landslides, such as: ⅰ) earthflows/flow-slides when incoherent pyroclastic deposits are involved; ⅱ) rock fall that generally affects compacted pyroclastic deposits, tuff, and lavas; ⅲ) rock and debris slides ^[44,45].

3.   Dataset and method

The multidisciplinary approach applied here analyzes and connects the dataset of parameters recognized in the literature as significant predisposing factors to landslide formation ^{[8,18,21,22,46,47,48,49,50,51,52,53,54,55,56,57]}. The parameters, listed in Figure 2, are represented by a geo-spatial (WGS 84 UTM N33 system) raster dataset covering the entire extent of the island with a spatial resolution of 5 m. The dataset, derived from different sources, was first elaborated and then combined in a GIS environment through a defined equation that calculates the index of susceptibility in each area of 25 m².

Figure 2.  Workflow of the seismically induced landslide susceptibility map.

DownLoad: Full-Size Img PowerPoint

3.1.   Datasets

The datasets of each parameter were derived from the following sources: digital elevation model (DEM), geo-portals, and the literature.

3.1.1.   DEM-derived datasets

Three datasets were derived from a high spatial resolution DEM of the island obtained by the Airborne LiDAR survey acquired in 2012 (http://www.pcn.minambiente.it/mattm). The original spatial resolution of this DEM (2 m) was resampled, through the nearest neighbor algorithm, to 5 m grid cell size to be consistent with the spatial resolution of all datasets. From this data, we derived the slope, relative relief, and aspect parameters:

- The slope was calculated using the algorithm from ^[58]. The five classes of slope (Figure 3a) were defined by applying statistical analysis to Slope values characterizing the source areas of 62 historical seismically induced landslides ^[8,21].

Figure 3.  Parameters used to calculate the index of susceptibility SNAP: a) Slope; b) Relative relief; c) Aspect (red and green colors indicate more and less predisposition to terrain instability, respectively ^[44,48,52]); d) Lithology (codes reported and defined in Table 1).

DownLoad: Full-Size Img PowerPoint

- The relative relief was calculated as the difference between the maximum and minimum digital elevation matrix, both obtained in a 5 × 5 window. The four classes of relative relief (Figure 3b) were defined through the statistical analysis of the matrix values characterizing the previous 62 source areas.

- The aspect was calculated considering the direction of the maximum terrain slope. Each pixel of the obtained raster matrix stores the numerical value of the angle between the maximum slope direction and the geographic north. Figure 3c illustrates the aspect through the eight main cardinal directions (N, NE, E, SE, S, SW, W, NW).

3.1.2.   Geo-portals-derived datasets

Lithology and land use datasets, available in polygon vector format through the geo-portals, were converted in raster format with a spatial resolution of 5 m, following the topological criteria on which the vector-raster transformation is based ^[58]:

- The lithology dataset comes from the "Campania" geo-portal (https://www.distrettoappenninomeridionale.it). This layer, consisting of 146 polygonal features, was grouped in 15 different lithology types (Table 1 and Figure 3d) according to the values of soil cohesion parameter reported in ^[21].

Table 1.  Lithology listed for the crescent value of soil cohesion parameter.

Lithology/Code (LC)	Lithology description	Cohesion (kPa)
1	Sandy and silty beach deposits	0
2	Dump materials and reworked soils	0
3	Slope debris	2
4	Eluvial-colluvial pyroclastic deposits	2
5	Landslide, debris, and colluvial deposits	2
6	Loose pyroclastic deposits	5
7	Pumice and ash layers	10
8	Marine silty sandstones with tuffites	30
9	Breccia	100
10	Lithified yellow tuff	325
11	Tuff with pyroclastic deposits	370
12	Welded tuff	425
13	Lithified pyroclastic deposits	600
14	Lava with welded scoriae	800
15	Lava	1000

---

 | Show Table

DownLoad: CSV

- The land use dataset was obtained by updating the 2018 Corine land cover (CLC) data, derived from Sentinel-2 satellite images at a spatial resolution of 10 m (https://land.copernicus.eu), with a vector layer digitized from 2020 ortho-photos acquired by CGR Spa with a spatial resolution of 20 cm. The updated layer consists of 322 polygon features described through 12 different typologies of land use (Figure 4a).

Figure 4.  Parameters used to calculate the index of susceptibility SNAP: a) Land use, grouped into six main classes according to their different influence on terrain instability ^[53,55,60]: (WD), high scrub (HS), low scrub (LS), cultivated soil (CS), barren soil (BS), and anthropic (AN); b) PGA; c) hydrothermal activity, divided into fumaroles fields (FM), hot springs (HS), and cold springs (CS); d) Faults, in which PF stands for proximal faults and NP for non-proximal faults.

DownLoad: Full-Size Img PowerPoint

3.1.3.   Literature-derived datasets

Datasets coming from the literature are the peak ground acceleration (PGA), the hydrothermal activity, and faults distribution:

-The PGA dataset, represented by a raster map (Figure 4b), was obtained by elaborating the macroseismic intensity data I_MCS of the earthquake of August 21, 2017 ^[22]. Data, represented by geospatial points, were interpolated using the inverse distance weighted (IDW) algorithm to obtain a continuous spatial surface in raster format of I_MCS with a 5 × 5 m spatial resolution. Subsequently, the following equation (1) from ^[59] was applied to I_MCS raster to derive the PGA values:

where I_MCS indicates the macroseismic intensity and σ is the standard deviation (σ = 0.222).

- The hydrothermal activity dataset was obtained as a vector layer where each polygon represents the area where fumaroles and springs (hot and cold) are located (Figure 4c). In detail, starting from the volcanological map of Ischia ^[37], the points representative of the fumarolic fields (FM), hot springs (HS), and cold springs (CS) were digitalized and then buffered with a radius of 50 m in order to account not only for the possible uncertainties on their exact location but also the dispersion of weathering around the springs.

- The faults dataset (Figure 4d) was created as a polygon vector layer by applying buffers of 25 m to the fault lines derived from the maps in ^[7] and ^[32]. This buffer accounts for possible uncertainties regarding faults' location. The buffer areas represent zones potentially affected by two typologies of faults: proximal to the 2017 earthquake epicenter (PF) and non-proximal (NP).

3.2.   SNAP index and weight assignment

The SNAP index (Seismically iNduced lAndslides suscePtibility) was computed by applying equation (2) and performed on a GIS environment. It combines the eight previously described parameters as raster by using the map algebra ^[58]. For each parameter and respective class of values, a weight (primary and relative, respectively), which reflects the inferred influence on generating terrain instability, is assigned. The SNAP equation is as follows:

where i represents the parameter and j its class; p_i is the primary weight (the weight assigned to each parameter), and w_ij is the relative class weight (the weight assigned to each class defined for each parameter). The weights values are ranked on a scale from 1 to 10 for both p_i and w_ij, where 1 stands for the lowest influence and 10 stands for the highest one. Weights are assigned according to both statistical analysis and available literature on Ischia island and similar geo-volcanological contexts. Regarding the primary weights (p_i), we assigned the highest value (10) to the PGA since it represents one of the most critical ground motion factors for triggering landslides during an earthquake ^{[8,61,62,63,64,65,66]}. Several works on Ischia ^[18,21] suggest that seismically induced landslides tend to occur along steep slopes mantled by loose volcaniclastic material having very poor geotechnical properties (with extremely low cohesion factor, to be considered null). According to such studies, the slope, relative relief, and lithology are considered the most relevant parameters to generate terrain instability. Slope is recognized as the first stability parameter since it directly controls the balance between resistant and gravitational forces causing any type of landslide ^[48,54]. For this reason, the assigned weight (9) is slightly higher than the weight assigned to both relative relief and lithology (8), followed by hydrothermal activity with a weight equal to 7. This parameter is the expression of vertical pathways of hydrothermal fluids and creates local rock weaknesses that, subjected to a seismic event, could make the terrain more prone to slide ^[67,68]. The weights assigned to land use and fault parameters are equal to 6. In one case, the value is given taking into account the relation between the land use and the slope stability according to ^[69]. Regarding the fault parameter, the assigned weight justifies that the landslides may occur in proximity to active faults ^[57] as considered potential weakness zones. Finally, for the aspect parameter, we defined a weight equal to 3 because this parameter is relevant only for some slope and land-use combinations.

Table 2 shows, for each parameter (i) and its class (j), the primary weights (p_i) and the relative class weights (w_ij) assigned according to literature data and statistical approach. The five classes of the slope parameter are defined by applying a statistical analysis on 62 historical seismically induced landslide source areas ^[8,21]. Weights have been assigned considering the slope thresholds recognized for both earthflows/flow-slides and rock fall generation ^[49,50,52]. The relative relief parameter is strongly linked to the slope and connected to erosion activity ^[70], which influences landslide triggering. The setup classes are 4, and the respective weights are defined with the same approach used for the slope parameter. The aspect parameter is represented by 8 standard classes that are grouped into 2 classes related to the northern and southern portions of the island. The south areas are more erodible than those facing north ^[47,51] because they tend to have more solar radiation, a greater thermal excursion, and less vegetation. For the lithology parameter, we assigned the weight to each class considering that the poorer the geotechnical properties, the higher the probability for a landslide to occur ^{[18,21,71,72]}. In particular, we consider cohesion as a reference geotechnical property (Table 1). For the land use parameter, the weights assigned to the classes take into account their influence on instability. In our case, we considered ⅰ) the stabilizing effect of woody roots ^[53], ⅱ) that urbanization can interfere with loose deposits equilibrium ^[57], and ⅲ) that sparse vegetation soils are more subjected to erosion ^[60]. For the PGA parameter, four equi-interval classes were defined and weighted considering the linear trend between MCS intensity and susceptibility proposed in ^[8]. Regarding the hydrothermal activity parameter, the classes distinguish cold springs from hot springs and fumarolic fields considering that high temperatures weather deposits and rocks more efficiently than cold ones ^[73,74]. The maximum weight assigned to fumarolic fields is related to the typical chemical weathering associated with this phenomenon. Finally, the faults parameter consists of two classes defined according to faults proximity to the epicenter of the earthquake of August 21, 2017: faults proximal to the epicenter (PF), which have a high possibility to be reactivated in case of an earthquake, and non-proximal faults (NP).

Table 2.  Susceptibility parameters, the assigned weight (p_i), and the classes (j) and their respective weights (w_ij). LC is the lithology code from Table 1.

Parameters	pi	j	wij
Slope (SL)	9	SL < 10°	1
		10° < SL < 20°	3
		20° < SL < 36°	8
		36° < SL < 45°	9
		SL > 45°	10
Relative relief (RR)	8	RR < 3m	2
		3 m < RR < 9 m	10
		9 m < RR < 13 m	4
		RR > 13m	1
Aspect	3	N, NE, NW, W	4
		S, SE, SW, E	8
Lithology	8	LC 14–15	1
		LC 12–13	3
		LC 9–11	5
		LC 1–8	10
Land use	6	Woods	3
		High scrub	4
		Cultivated soil	7
		Anthropic	8
		Low scrub	7
		Barren soil	8
Peak ground acceleration (PGA)	10	PGA < 0.1 g	1
		0.1 g < PGA < 0.2 g	2
		0.2 g < PGA < 0. 3 g	5
		PGA > 0. 3g	10
Hydrothermal activity	7	Cold springs	1
		Hot springs	4
		Fumarolic fields	10
Faults	6	Proximal to epicenter	8
		Non-proximal to epicenter	4

---

 | Show Table

DownLoad: CSV

4.   Results

The resulting SNAP index map classifies the entire Ischia territory using five classes relating to seismically induced landslide susceptibility. Figure 5 shows the classes obtained by using the Jenks Natural Breaks algorithm: very low, low, moderate, high, and very high SNAP. Areas with very low (75 ≤ SNAP ≤ 165) and low (165 < SNAP ≤ 218) susceptibility are poorly exposed to seismically induced landslides. They cover ~34% of the entire territory and are mainly localized in the west and east portion of the island. Areas with moderate susceptibility (218 < SNAP ≤ 273) are located mainly in the southeastern coastal zones where the generation of landslides induced by an earthquake is due to the presence of significant slopes (>25°). These areas affect 33% of the island. Finally, the zones with high (273 < SNAP ≤ 315) and very high (SNAP > 315) susceptibility cover the remaining 33% of the island, involving mainly its central portion (Mt. Epomeo) and a few coastal areas located south and east.

Figure 5.  SNAP susceptibility map for seismically induced landslides in the Ischia island.

DownLoad: Full-Size Img PowerPoint

Further analyzing the map, it is observed that the areas with very high susceptibility are located along the northern, western, and eastern flanks of Mt. Epomeo, at elevations above 50–100 m a.s.l.. These areas show steep slopes, high relative relief, loose material (in particular debris and colluvial deposits in the west and north zones, pyroclastic deposits in the eastern zones), fumarolic fields or high-temperature springs, and faults bordering the resurgence block of Mt. Epomeo. The maximum values of susceptibility are found in the northern zones, inside Casamicciola Terme municipality, where the highest values of PGA was recorded due to the proximity to the epicenter of the August 21, 2017 seismic event. Focusing on Mt. Epomeo, the susceptibility of the southern flank is high rather than very high, as in the northern flank, despite the presence of very steep slopes, faults, and high values of relative relief due to gullies and drainage erosional channels. This difference can be explained by the southern flank being affected by moderate PGA values and showing a greater distribution of areas with slopes < 25° than the north flank. Moreover, the areas from moderate to high SNAP involve the entire SE sector of the island, in particular the cliffed coasts that are featured by slopes > 36° and relative relief > 9 m. In addition, the areas of the SE sector are characterized by low scrub and cultivated soils. Finally, the low and very low SNAP zones are distributed in the remaining territory of the island, with more gentle slopes, such as the base of Mt. Epomeo, or in several areas very close to the coast, such as the north and west ones. In general, the low and very low SNAP areas show low PGA values and lithology with cohesion values > 100 kPa.

In addition, the SNAP distribution of each municipality (Figure 6) shows that Casamicciola Terme presents the largest area with very high SNAP (1.4 km² corresponding to approximately 25% of the entire municipality), followed by Lacco Ameno with 0.3 km², corresponding to 12% of its territory. Combining the two most relevant classes (high and very high) into a new class, here defined as "critical class", approximately 60% of Casamicciola Terme seems to be vulnerable to seismically induced landslides. The second most exposed municipality is Serrara Fontana, with 51% of the territory involved, followed by Barano d'Ischia (36%), Lacco Ameno (35%), Forio (23%), and Ischia (13%).

Figure 6.  Distribution of the SNAP classes in Ischia municipalities.

DownLoad: Full-Size Img PowerPoint

5.   Discussion

In this section, we discuss the obtained SNAP map, comparing our results with recent and historical seismically induced ground effects, the analyses obtained by a dedicated ERT field survey, and the deformation data derived from 2015–2022 InSAR analysis.

5.1.   Seismically induced ground effects comparison

The comparison between the SNAP map and the seismically induced ground effects was made through geospatial datasets coming from available databases at different scales, from national to local. The datasets were overlapped with our susceptibility map (Figure 7) and highlight that:

Figure 7.  Ground seismic effects ^[12,24,25] and seismically induced landslide source areas ^[21] overlapped on the SNAP map. The symbols inside the square of the CFTIlandslides are related to the degree of location accuracy (cross: poor; point: intermediate; no symbol: excellent/good).

DownLoad: Full-Size Img PowerPoint

- More than 80% of the landslides and ground cracks from the CEDIT database ^[24] fall in high and very high SNAP areas. CEDIT is the Italian catalog of earthquake-induced ground failures, which maps the earthquakes' epicenters and related ground effects for the Italian territory from 1125 to 2022.

- 82% of rockfalls and slow earthflows detected from the CFTIlandslides, the Italian database for historical earthquake-induced landslides ^[25], fall in high and very high SNAP areas. This database localizes (with different degrees of accuracy) the landslides at national, regional, and local scales, combining the existing databases and integrating them with the results of previous studies and research activities.

- 85% of source areas related to 62 historical seismically induced landslides mapped by ^[21] result in the high (48%) and very high (37%) SNAP classes. Moreover, it is noteworthy that most of the source area boundaries delineate the limits of areas classified with high/very high susceptibility. This is particularly evident in the northern flank of Mt. Epomeo.

- Finally, 47 ground effects (fractures, earth slides, landslides, and rockfall) induced by the 2017 earthquake and mapped by ^[12] as geo-coded points were superimposed on the SNAP map. In detail, 46 points (98%) fall inside the areas classified as high and very high susceptibility, despite a part of these points (20%) being located where the PGA does not reach the highest-class values (Figure 4b).

5.2.   Field ERT analysis

In order to directly investigate and characterize the subsoil elements influencing mass movements induced by earthquakes, a dedicated electrical resistivity tomography (ERT) field survey was applied on a specific test site among those identified in the susceptibility map showing significant deformation ^[75] and located in an accessible location. The test site was located in the municipality of Serrara Fontana, in the south sector of Mt. Epomeo, at an elevation of ~200 m from Maronti beach (Figure 8). This area is affected by tectonic lineaments, lies at the edge of an ancient paleo-landslide, and has a strong energy relief. Furthermore, an earthquake event arose close to the Serrara Fontana site (Maronti beach) as one of the higher intensity events for the island (M = 2.6 in 1997, see ^[23]). Finally, in southern Ischia Island, offshore geophysical investigation and marine sampling show that some of the higher mega-landslides occurred in the area ^[76,77]. These authors refer to giant debris avalanches and debris flow events that also occurred in historical time (3 ka B.P. and 2.4 B.P.). Thousands of blocks can still be observed on the sea floor, extending 40–45 km from the island. The age of the debris flow event aligns with a historical account by the Greek historian Strabo (64 BC–21 AD), which describes a sudden collapse of one of the island's flanks accompanied by a tsunami. This event is reported to have occurred at the beginning of the Greek colonization of the island, around 600–700 B.C. ^[76,77].

Figure 8.  Southern flank of the Ischia island geological map 1:25.000 (F.154 sheet, ISPRA Geo-logical Map). The ERT profiles (green lines) overlain on the main outcropping geological units.

DownLoad: Full-Size Img PowerPoint

The Serrara Fontana outcropping geological units have poor mechanical properties such as debris flow deposits, rock avalanche, and hydrothermally altered pyroclastic deposits (Carta Geologica d'Italia 1:10.000 – "Ischia island" F.154 by ISPRA in Figure 8). The above-described characteristics highlight the test site's predisposition to mass movements that, together with the high energy relief, support its high susceptibility level.

The ERT technique aims to reconstruct the subsoil resistivity distribution through multi-electrode resistivity measurements. These resistivity (Ω·m) and their inverse, conductivity (mV/V), are properties of the material and provide information on the ability of the electric current to flow in the subsoil. Injection electrodes send electric currents into the ground generating an electric field, whereas measuring electrodes detect the potential difference at different distances from the injection electrodes ^[78]. During the ERT field survey, resistivity (ρ) and induced polarization (IP) measurements ^[79,80,81] have been acquired along three profiles (green lines in Figure 8). The measurements were carried out using the IRIS SYSCAL Pro Switch 72 resistimeter with a cable setting of 48 electrodes placed every 5 m and extending for approximately 240 m in the NE–SW direction. The maximum reached exploration depth was approximately 50 m. The three profiles are almost consecutive along the NE–SW direction (Figure 8) and allow an exploratory continuity of the western "Serrara Fontana" flank of approximately 420 m along and 100 m across the slope. The real resistivity distributions were generated by inverting the measurements of the apparent resistivity using the RES2DINV software ^[82,83].

The inversion of the ERT measurements reached acceptable RMS misfits for profiles 1 and 2 (the absolute errors are 1.4 and 3.3, respectively), whereas the misfit for profile 3 is higher (absolute error = 7.4). However, the tomography of profile 3 appears to be coherent with those of the other two.

The tomographies are shown in Figure 9 and their interpretation is in Figure 10. The resistivity patterns of the three profiles show remarkable lateral and vertical variations of the electrical properties (intensity range of 1.3 Ω·m < ρ < 430 Ω·m) and noticeable discontinuities. Particularly, the tomographies reveal eight consecutive electro-layers named A, A₁, B, C, D, E, F, and G. The most striking feature is the noteworthy vertical electrical contacts between the A–B and A–C electro-layers (F₁ and F₂ in profile 1 southern sector, Figures 9 and 10). These geophysical discontinuities are most likely related to the fracturing separating geological units with different electrical properties. In fact, an evident tectonic fault can be observed to outcrop a few tens of meters southwestward the profile 1 (Figure 8), whereas the limit between the "TSP" unit (Serrara-Cava Petrella Tuffs) and the "I_b" unit (mixed deposits, debris flow and/or torrential) runs parallel to the profile a few meters westward (Note Illustrative della Carta Geologica d'Italia 1:25000, Foglio 464, ISPRA). Moreover, both the B and C electro-layers (profile 1) have a similar geometry and resistivity range (50 Ω·m < ρ < 80 Ω·m) and can represent the same geological formation or part of it, most likely the TSP unit. The latter is composed of pyroclastites emitted from eruptive centers located close to the profiles. The A unit has a very low resistivity (ρ < 10 Ω·m) and involves 2/3 of profile 1 and one-half of profile 2; it may represent the hydrothermally altered portion of the BSR unit ("Bocca di Serra unit"). This consists of rock avalanches with very coarse clastics, blocks, and mega blocks (metric to decametric) belonging to hydro-thermalized tuffs and epiclastic rocks. They are immersed in a coarse matrix with sands and lapilli and hydrothermalized tuffs. However, a cylindrical-like shape feature, the E electro-layer, occurs in the northeastern half of profile 2 and in the southwestern side of profile 3 with similar thickness and geometry (10 Ω·m < ρ < 27 Ω·m and 10 Ω·m < ρ < 50 Ω·m, respectively, in Figure 9). Therefore, the E unit could represent one of these decametric elements of the BSR unit intercepted at different elevations by the two profiles. In the northeastern side of profile 2 and 3, at the same depth interval of the A unit, the A₁ electro-layer occurs with the same resistivity range as the former. The A₁ unit also corresponds to a higher chargeability value area (30–130 mV/V) of the induced polarization section, following the same geometry and depth interval of the resistivity distribution (Figure 9). The contemporaneous presence of very-low resistivity and higher chargeability indicates a high ion concentration in the subsoil inside the BSR unit, responsible for the battery-like behavior of the A₁ unit. This may be related to a water trap occurring approximately 10 m below the ground and caused by hydrothermal water upwelling. Alternatively, fractures, cavities, and caves are known to host saline concretions several meters thick on Ischia island. These saline concretions can increase the subsoil chargeability. The D and F electro-layers have variable thickness ranging from a few meters to 20 m and represent the most recent lithologies not affected by the F1 and F2 discontinuities. The "Ib" geological unit is largely cropping out along this side of the island (Figure 8) and is associated with debris and mud flows deposits, whereas the uppermost stratigraphic "b2" unit ("eluvial-colluvial blanket") is related to incoherent loamy-sandy soils. The varied composition of these units can account for the different electrical properties of the superficial D and F electro-layers, ranging from very resistive (up to 430 Ω·m) to very conductive (<10 Ω·m).

Figure 9.  ERT profiles together with the main electro-layer units. Profile 1 and 2 are consecutive along the NW-SE direction, profile 3 is shifted approximately 100 m southeastward. Profile 3 also has the IP section.

DownLoad: Full-Size Img PowerPoint

Figure 10.  Interpretative geological sections of the ERT profiles. The figure vertical-to-horizontal ratio is 1.

DownLoad: Full-Size Img PowerPoint

The high resistivity G electro-layer is placed at the northeastern boundary of profile 3, not completely sampled by the measurement.

The discussed elements highlight a paleo-landslide structure reaching at least 50 m depth close to the southwestern Mt. Epomeo boarding fault, together with the pyroclastites. Therefore, the ERT investigations show characteristics such as discontinuities, faults, cracks, and volcanic lithologies with different geometries and bedding in a very short space. These characteristics, together with high energy relief and the existence of a paleo-landslide, represent at a small scale the predisposing factors to mass movements, which can most likely occur on a wider scale.

5.3.   InSAR data

A further investigation of areas with very high SNAP was conducted combining the ground displacement obtained through the InSAR processing of Sentinel-1 satellite data. In detail, the small baseline subset (SBAS) multitemporal method ^[84] has been adopted to process Sentinel-1 data acquired along both the ascending and descending satellites track, from January 2015 to January 2022. Data were processed at 15 m of ground resolution by using the Sarscape software (Sarmap SA, 5.6 version). The interferogram stacks were obtained with orbital separation (spatial baseline) not exceeding 211 m and maximum temporal distance between two satellite passes (temporal baseline) of 48 days. These settings resulted in 581 and 562 interferograms for descending and ascending orbit, respectively. The geophysical model derived from ^[75] and characterizing the August 21, 2017 earthquake was introduced into processing to avoid the influence of the coseismic displacement within the long-term signal. The availability of both ascending and descending acquisition geometries allowed us to calculate the horizontal (east-west) and vertical components of the ground motion ^[85] in the time period common to the two elaborations, i.e., 12/01/2015–10/01/2022. However, the vertical displacement map does not cover entirely the Mt. Epomeo, as the processing does not resolve the summit area. The errors associated with the vertical displacement map are ± 3 mm/yr and were estimated using the formula from ^[86], which considers parameters such as interferometric coherence and wavelength. The map of vertical displacement rates of the Ischia, characterized by only negative values, was then combined with the susceptibility map considering the very high and high class (Figure 11). In such combination, only the areas characterized by vertical displacement rates lower than the derived error (−3 mm/y) have been considered. In detail, such areas show a deformation rate in subsidence from −3 mm/y to −30 mm/y in the period 2015–2022. A new critical class was defined by overlapping the subsidence areas with the very high SNAP. This class identifies the areas most prone to seismically induced landslides (Figure 11). The north flank of Mt. Epomeo provides the most interesting results, involving, especially, the Casamicciola Terme municipality and confirming what arose from Figures 5 and 6. The critical class also covers the west flank, extending to Lacco Ameno village, and the W–NW flank, affecting the Forio municipality. Finally, several sparse critical areas are in the south of Mt. Epomeo, in the gullies zone of Serrara Fontana, and in some hilly parts of the Barano d'Ischia village. The only municipality not affected by critical areas is the Ischia municipality.

Figure 11.  The critical class areas overlapped with the vertical displacement map derived from InSar Sentinel-1data. In the background, the morphology of the island.

DownLoad: Full-Size Img PowerPoint

6.   Conclusions

In this work, a very-high spatial resolution seismically induced landslide susceptibility map of the Ischia island is presented, considering the earthquake of August 21, 2017, as the reference event. The map is the result of a multidisciplinary approach based on the different typologies of data considered relevant for the classification of the areas more prone to generate landslides. Eight parameters (slope, relative relief, aspect, lithology, land use, PGA, hydrothermal activity, and faults) were combined in a GIS environment to build a SNAP index able to quantify seismically induced landslide susceptibility. The SNAP index classifies the entire territory of the island into five degrees of susceptibility. The zones with very low and low SNAP are mainly distributed in the northern and western coastal areas where slopes are gentle, and the PGA values are low. The moderate class covers the areas at S and SE, where gullies associated with drainage system channels and cliffed coasts are present, respectively. Finally, the areas with high and very high susceptibility cover 33% of the island, mostly affecting the NE, N, and NW flanks of Mt. Epomeo. Such results are compared with historical and recent seismically induced ground effects from several databases. These comparisons highlight that the generated susceptibility map is consistent, as more than 81% of the effects stored in the databases fall in the high and very high SNAP classes. Moreover, the ERT investigation, carried out in high SNAP areas of the southern flank of Mt. Epomeo, reveals tectonic and litho-stratigraphic characteristics that confirm a significant predisposition to mass movements. Finally, combining the very high SNAP areas with the subsidence areas obtained from the 2015–2022 Sentinel-1 InSAR data, a further critical class was identified, confirming that the north flank of Mt. Epomeo is the most vulnerable to landslides. The Casamicciola Terme municipality is the most exposed, whereas the Ischia municipality does not contain critical areas. Such results represent an improvement in the knowledge of the areas potentially exposed to seismically induced landslides in Ischia island in case of future earthquakes comparable to the 2017 one.

Author contributions

RG: Conceptualization, Methodology, Validation, Formal analysis, Writing – Original Draft, Visualization. MB: Conceptualization. Methodology, Validation, Formal analysis, Writing – Review & Editing, Supervision. LB: Validation, Formal analysis, Writing – Review & Editing. RDR: Investigation, Validation, Formal analysis, Writing – Review & Editing. VS: Investigation, Writing – Review & Editing. LC: Validation, Writing – Review & Editing. CT: Validation, Formal analysis, Writing – Review & Editing. LC: Investigation, Writing – Review & Editing. CS: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Writing – Review & Editing, Project administration, Supervision.

Acknowledgments

The present work is supported by the Research Project FRASI-Integrated and multi-scale approach for the definition of seismically induced landslides hazard in the Italian territory, funded by the MATTM. Copernicus Sentinel-1 data are freely distributed by the European Space Agency and available at the Copernicus Open Access Hub.

Conflict of interest

Marina Bisson and Claudia Spinetti are guest editors for AIMS Geosciences and were not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.