The Relationship between InSAR Coseismic Deformation and Earthquake-Induced Landslides Associated with the 2017 M w 3 . 9 Ischia ( Italy ) Earthquake

We investigated the contribution of earthquake-induced surface movements to the ground displacements detected through Interferometric Synthetic Aperture Radar (InSAR) data, after the Mw 3.9 Ischia earthquake on 21 August 2017. A permanent displacement approach, based on the limit equilibrium method, allowed estimation of the spatial extent of the earthquake-induced landslides and the associated probability of failure. The proposed procedure identified critical areas potentially affected by slope movements partially overlapping the coseismic ground displacement retrieved by InSAR data. Therefore, the observed ground displacement field is the combination of both fault slip and surficial sliding caused by the seismic shaking. These findings highlight the need to perform preliminary calculations to account for the non-tectonic contributions to ground displacements before any estimation of the earthquake source geometry and kinematics. Such information is fundamental to avoid both the incorrect definition of the source geometry and the possible overestimation of the coseismic slip over the causative fault. Moreover, knowledge of the areas potentially affected by slope movements could contribute to better management of a seismic emergency, especially in areas exposed to high seismic and hydrogeological risks.


Introduction
On 21 August 2017, an M w 3.9 earthquake (http://cnt.rm.ingv.it/)struck the Ischia Island (southern Italy) (Figure 1).Despite the low magnitude, the event (hereinafter the Ischia earthquake) caused damage and two casualties in Casamicciola village [1].The concurrence of a shallow seismic source (approximately 1 km below the ground), local seismic amplification, and the high vulnerability of buildings explained such unexpected severity.
The epicentral area had already experienced moderate magnitude, shallow seismic events in the past.Major earthquakes occurred at Ischia from AD 1000 onwards, with epicentral intensities varying from MCS (Mercalli-Cancani-Sieberg) V to X.Most of them were located along the northern periphery of the island, close the Casamicciola village [2] (Figure 1).The last and strongest earthquake (MCS intensity = X) occurred in 1883, in the same area of the 21 August 2017 event, causing significant damage and triggering several landslides [3,4].
(MCS intensity = X) occurred in 1883, in the same area of the 21 August 2017 event, causing significant damage and triggering several landslides [3,4].
Space-borne Interferometric Synthetic Aperture Radar (InSAR) and Global Positioning System (GPS) data allowed a complete overview of the coseismic displacement field due to the mainshock to be provided, thus identifying the areas affected by significant deformations [9].
Soon after the Ischia earthquake, several actions were deployed to detect the effects induced by the seismic event on the environment.More than 100 observations across the epicentral area and its surroundings were collected and classified.Coseismic effects consist in ground fractures (small open cracks with vertical offset ≤ 1 cm), ground ruptures (open breaks with vertical offset ≥ 1 cm), gravitational phenomena (small size collapses and small landslides in volcanoclastic deposits), and shaking effects (e.g., failure of drywalls) (Figure 1) [5].The observed fractures, with lengths up to some tens of meters, were attributed to the propagation of the seismogenic motion up to the surface.However, gravitational phenomena cannot be excluded since some ruptures cut across old landslides and reworked soils.
Space-borne Interferometric Synthetic Aperture Radar (InSAR) and Global Positioning System (GPS) data allowed a complete overview of the coseismic displacement field due to the mainshock to be provided, thus identifying the areas affected by significant deformations [9].
The InSAR-derived coseismic displacement field, together with GPS and seismological data, were also exploited to estimate the geometry and kinematics of the seismic source using a finite-fault dislocation model embedded in an elastic and homogeneous half-space [9].The best-fit solution of the source inversion consisted of an E-W striking, south dipping normal fault, with its center located at a depth of approximately 800 m.The solution obtained only considered tectonic phenomena, i.e., the slip along fault planes, neglecting possible non-tectonic causes of the observed coseismic displacement.
Many phenomena can contribute to ground deformation after strong earthquakes, including the fault slip [10], the compaction induced by liquefaction [11,12], and the earthquake-triggered slope instabilities [13][14][15][16][17].All these phenomena play different roles depending on the geological, geotechnical, and hydrogeological features of the epicentral area.Generally, coseismic ground displacements are primarily caused by the slip on the fault, while other contributions are limited both in space and amplitude [13,18].For this reason, the coseismic surface deformation is modeled with linear and non-linear inversion models, which neglect non-tectonic phenomena in favor of fault slip.However, ignoring non-tectonic events is not always acceptable, especially in the case of shallow earthquakes with minor magnitude, as in the Ischia earthquake.
In such a case, the observed ground displacements have a low amplitude (less than 4 cm) and a moderate surface extension (approximately 1 km 2 ) [9] (Supplementary Figure S1).Furthermore, coseismic displacements are located on an area along the northern flank of Mt.Epomeo which is characterized by high earthquake-induced landslide susceptibility [6,19] and has been already affected by catastrophic earthquake-induced landslides in the past [3] (Figure 1 and Supplementary Figure S1).Such observations suggest that non-tectonic displacements caused by earthquake-induced shallow landslides are not negligible and must be opportunely identified, modeled and analyzed to improve any source inversion procedure.
In this work, we applied a procedure to identify areas affected by slope instabilities triggered by the 2017 Ischia earthquake.The method is based on previously published works using a simplified Newmark analysis [19,20] which allows estimation of the landslide spatial extent and the associated probability of failure at a regional scale.The identified areas are then compared with coseismic ground displacements, retrieved from the processing of Sentinel-1 SAR data, to analyze the contribution of tectonic and non-tectonic phenomena to the observed ground displacements.

Geological Setting
Ischia Island constitutes the westernmost active field of a volcanic district comprising the Phlegrean Fields, the Procida Island, and the Mt.Somma-Vesuvius complex.The primary structure of the island is Mt.Epomeo (Figure 1), which constitutes a volcano-tectonic horst with an altitude of 787 m a.s.l.[21].The island's morphology is continuously formed and reshaped by multiple tectonic processes such as eruptions, earthquakes, resurgences, collapses, and landslides [2,7,[22][23][24].As a result, Ischia is composed of various formations and lithologies including volcanic rocks and soils (Figure 1) produced by effusive and explosive volcanic activities occurred since its creation [7,25].The current resurgent phase raised the Mt.Epomeo summit up to its present altitude [26] and was accompanied by the activation of faults and renewal of volcanism that produced over-steepening of slopes [25].The steep angles and the poor geotechnical properties of the outcropping lithologies have favored slope instabilities, resulting in shallow mass movements and large rock and debris avalanches.Indeed, slope instability-related deposits, mainly triggered by local seismic activity, are identified in and around the island [4,25,[27][28][29].In particular, volcaniclastic deposits with variable thickness and strength cover the north-western flanks of Mt.Epomeo [6,29,30] (Figure 1).

Interferometric Synthetic Aperture Radar (InSAR) Data and Analysis
The dataset used for the InSAR analysis of the coseismic crustal deformations occurred after the Ischia event consists of Sentinel-1 (S1), C-Band, Single Look Complex (SLC) SAR images, collected before and after the earthquake (i.e., 16-22 August 2017).We performed two interferograms from the ascending and descending track, respectively, using the SARScape ® software (Sarmap SA, Purasca, Switzerland).The 30 m Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) was used to remove the topography contribution, and Goldstein filtering [31] was adopted to reduce phase noise of InSAR data.Finally, the interferograms were unwrapped by using the minimum cost flow (MCF) phase unwrapping algorithm [32].The retrieved S1 InSAR displacement maps were strongly affected by atmospheric residuals, so we applied a technique to reduce that disturbance.We considered the Generic Atmospheric Correction Online Service (GACOS) (http://ceg-research.ncl.ac.uk/v2/gacos/) [33,34] approach that uses the iterative tropospheric decomposition (ITD) model to separate stratified and turbulent signals from total tropospheric delays, and generate high spatial resolution zenith total delay maps to be used for correcting the InSAR measurements.

Landslide Modeling Approach
The discrimination between the tectonic and non-tectonic contribution to coseismic ground displacements requires accurate modeling of each phenomenon affecting the surface deformations.Hence, analyzing non-tectonic processes calls for methods that differ from those commonly used in the coseismic studies.Manifold approaches are currently available to model the dynamic performance of slopes, i.e., pseudostatic stability analysis, permanent displacement analysis, and advanced numerical simulations [35].Pseudostatic stability analysis is simple, requires few data for calibration, and provides a simple, scalar index of stability.This simplicity stems from a rather crude characterization of the physical processes of dynamic slope behavior that yields results affected by strong uncertainties.Permanent displacement analysis [36] improves simplistic pseudostatic analysis and permits the calculation of earthquake-induced ground movements.Finally, numerical simulations allow introducing geological and geometrical complexities and may account for the correct cyclic behavior of soils affected by seismic loads.However, such simulations are more complicated, require a significant amount of data for calibration, and are hardly applicable at medium to small scale.
The more sophisticated methods do a better job of modeling the dynamic response of the landslide material, thus yielding more accurate displacement estimates, but there is a trade-off in the complexity of the analysis and the difficulty in acquiring the needed input parameters.Therefore, the selection of the appropriate modeling approach is application-dependent, i.e., it depends on the problem scale and the requirement of the analysis.For the Ischia earthquake case study, the target is to identify areas potentially affected by earthquake-induced landslides at a regional level, i.e., over the whole island.According to this rationale, the permanent displacement analysis developed by Newmark [36] represents a good compromise between model complexity and the reliability of results.Such an approach has been already applied successfully for the evaluation of earthquake-induced landslide susceptibility over Ischia Island [19].
In its simplest form, Newmark's method models a potential landslide as a rigid block that slides on an inclined plane (Figure 2a).The block has a known critical acceleration, i.e., the threshold acceleration required to overcome the frictional resistance at the interface between the block and the plane and then initiate sliding.This value generally depends on the geometry of the slope and the frictional resistance at the interface between the slope and the inclined plane.The critical acceleration of a potential landslide block can be estimated as a simple function of the static factor of safety and the landslide geometry according to Equation (1) [36]: where a c is the critical acceleration, g is the gravitational acceleration, FS is the static factor of safety, and α is the angle from the horizontal of the potential landslide block, which is approximated with the slope angle.
is likely to perform during seismic shaking.Indeed, this method yields reasonable results for the dynamic analysis of natural hills and has been widely adopted for the regional assessment of seismic slope stability [20,35,38,39] and for the evaluation of earthquake-induced landslide susceptibility even on Ischia island [19].Therefore, such a modeling approach is acceptable for the regional assessment of earthquake-induced landslides on Ischia Island.

Modeling Procedure
The sequence of steps of the modeling procedure is represented in the flowchart in Figure 3.The workflow is similar to the method proposed in other studies [19,20,39] for the assessment of landslide susceptibility.The whole island is subdivided into a regular grid, and Newmark displacements are estimated for each grid cell.Then, a probability of landslide failure is associated with each computed movement, thus obtaining a map with the landslide failure probability associated with the earthquake of interest.
The procedure requires the sequential calculation of the following parameters: (i) the static FS values; (ii) the ac values; (iii) the Newmark displacement (Dn); and (iv) the associated probability of landslide failure (P(f)).
(i) The FS value in static conditions (i.e., without seismic forces) is calculated as the ratio between the resting and driving forces acting on the slope.A limit-equilibrium approach is adopted, which assumes a simple model of an infinite slope, made with a homogeneous material having both frictional and cohesive strength (Figure 2c).FS is calculated as follows: where the first term of the sum is the cohesive strength, the second term is the frictional strength, and the third term represents the reduction of the frictional strength because of the presence of pore water.Solving Equation (2) requires apriori information about the slope inclination (α) and thickness (t), the effective soil cohesion (c') and friction (φ'), the unit weight of water (γw) and soil (γ), and the portion of the slope that is saturated by water (m).An acceleration time history is compared to the critical acceleration; if the recorded values are below the a c threshold the block is stable; otherwise, the block undergoes permanent displacements.The accelerations exceeding the a c value are integrated twice to obtain the velocity-time history and the cumulative movement of the landslide block, respectively.
Newmark's method considers only rigid-body movements and plastic deformation at the contact between the block and the inclined plane.Such an approach is appropriate for the dynamic performance assessment of shallow, disrupted landslide masses during seismic events (see Jibson [35,37] and references therein).However, it does not provide accurate predictions of landslide displacements measured in the field.Instead, Newmark's method is used to determine how a slope is likely to perform during seismic shaking.Indeed, this method yields reasonable results for the dynamic analysis of natural hills and has been widely adopted for the regional assessment of seismic slope stability [20,35,38,39] and for the evaluation of earthquake-induced landslide susceptibility even on Ischia island [19].Therefore, such a modeling approach is acceptable for the regional assessment of earthquake-induced landslides on Ischia Island.

Modeling Procedure
The sequence of steps of the modeling procedure is represented in the flowchart in Figure 3.The workflow is similar to the method proposed in other studies [19,20,39] for the assessment of landslide susceptibility.The whole island is subdivided into a regular grid, and Newmark displacements are estimated for each grid cell.Then, a probability of landslide failure is associated with each computed movement, thus obtaining a map with the landslide failure probability associated with the earthquake of interest.
The procedure requires the sequential calculation of the following parameters: (i) the static FS values; (ii) the a c values; (iii) the Newmark displacement (D n ); and (iv) the associated probability of landslide failure (P(f )).
(i) The FS value in static conditions (i.e., without seismic forces) is calculated as the ratio between the resting and driving forces acting on the slope.A limit-equilibrium approach is adopted, which assumes a simple model of an infinite slope, made with a homogeneous material having both frictional and cohesive strength (Figure 2c).FS is calculated as follows: where the first term of the sum is the cohesive strength, the second term is the frictional strength, and the third term represents the reduction of the frictional strength because of the presence of pore water.Solving Equation ( 2) requires apriori information about the slope inclination (α) and thickness (t), the effective soil cohesion (c') and friction (ϕ'), the unit weight of water (γ w ) and soil (γ), and the portion of the slope that is saturated by water (m).
strength, may result in statically unstable grid cells (i.e., FS ≤ 1).Some studies propose to iteratively increase the absolute strength values of all lithological units until all slope angles less than 60° are statically stable [20].Other authors suggest avoiding the occurrence of unrealistically high strength parameters, assigning a fixed value FS = 1.01, i.e., barely above equilibrium, to the unstable cells having slope values lower than one [19].Both approaches are equally effective to prevent negative ac values.
(iii) Dn values are then calculated with the Newmark's approach.The latter requires the identification of a strong-motion record for each grid cell to compare with the ac value.For a regional scale analysis, the definition of an acceleration time history for each grid cell is unpractical and inappropriate, since earthquake strong-motion records exhibit a significant spatial variability [20].Therefore, Dn values are calculated through a simplified Newmark analysis that uses an empirical regression equation (Equation ( 3)).The latter combines the ac value of the slope with one or more parameters representing the shacking intensity of the target earthquake [37,40]: where PGA is the peak ground acceleration of the earthquake and σ (=0.51) is the standard deviation of the regression.The estimation of the PGA at each grid cell is not straightforward.Indeed, only one acceleration record of the Ischia earthquake is available at the IOCA station, located at approximately 600 m from the earthquake epicenter (Figure 1).The recorded horizontal PGA, i.e., 0.28 g and 0.19 g in the E-W and N-S directions, respectively [41], gives a picture of the seismic shaking in the near field only.Therefore, the PGA distribution over the whole island has been estimated using the macroseismic earthquake intensity.The latter quantifies the effects of an earthquake on the Earth's surface, anthropic structures, and the environment.An empirical PGA-intensity regression law, specially developed for the Ischia Island by Paparo and Tinti [30], has been used to estimate the PGA distribution over the whole island (Equation ( 4)): (ii) a c value is computed with Equation (1) by combining the FS value for each grid cell with the slope inclination obtained from the available DEM. Equation (1) requires that the slope is statically stable (i.e., FS > 1) to prevent the occurrence of unrealistic negative a c values.Predisposing factors, such as the simultaneous occurrence of moderate-to-steep slopes and lithological units with low strength, may result in statically unstable grid cells (i.e., FS ≤ 1).Some studies propose to iteratively increase the absolute strength values of all lithological units until all slope angles less than 60 • are statically stable [20].Other authors suggest avoiding the occurrence of unrealistically high strength parameters, assigning a fixed value FS = 1.01, i.e., barely above equilibrium, to the unstable cells having slope values lower than one [19].Both approaches are equally effective to prevent negative a c values.
(iii) D n values are then calculated with the Newmark's approach.The latter requires the identification of a strong-motion record for each grid cell to compare with the a c value.For a regional scale analysis, the definition of an acceleration time history for each grid cell is unpractical and inappropriate, since earthquake strong-motion records exhibit a significant spatial variability [20].Therefore, D n values are calculated through a simplified Newmark analysis that uses an empirical regression equation (Equation ( 3)).The latter combines the a c value of the slope with one or more parameters representing the shacking intensity of the target earthquake [37,40]: where PGA is the peak ground acceleration of the earthquake and σ (=0.51) is the standard deviation of the regression.The estimation of the PGA at each grid cell is not straightforward.Indeed, only one acceleration record of the Ischia earthquake is available at the IOCA station, located at approximately 600 m from the earthquake epicenter (Figure 1).The recorded horizontal PGA, i.e., 0.28 g and 0.19 g in the E-W and N-S directions, respectively [41], gives a picture of the seismic shaking in the near field only.Therefore, the PGA distribution over the whole island has been estimated using the macroseismic earthquake intensity.The latter quantifies the effects of an earthquake on the Earth's surface, anthropic structures, and the environment.An empirical PGA-intensity regression law, specially developed for the Ischia Island by Paparo and Tinti [30], has been used to estimate the PGA distribution over the whole island (Equation ( 4)): log(PGA) = 0.764 + 0.176 where I MCS is the macroseismic intensity defined according to the MCS scale [42], and σ (= 0.222) is the standard deviation of the regression.Finally, the PGA values were interpolated to obtain a continuous PGA distribution.We used a natural neighbor (NN) interpolation technique [43], which performs well with scattered and irregularly spaced data [44].
Given the PGA distribution, Equation (3) allows estimating the downslope-oriented D n values caused by the earthquake at each grid cell.D n predictions do not necessarily correspond to measurable slope movements in the field since the adopted approach is a simplification of the real dynamic behavior of soils under seismic loads.Commonly, the effectiveness of modeled displacements is evaluated through their possible effect on a potential landslide.In particular, movements in the range of 15-100 cm are expected to induce ground cracking for deep, massive landslides, while shallow, small landslides are triggered for lower movements, i.e., 2-15 cm (see Jibson [35] and references therein).
In this work, we avoided fixing a D n threshold value, but we interpreted the modeled displacements as an index of potential instability to be correlated with field observations.It is known that landslides preferentially happen in areas where landslides have already occurred in the past [45].Therefore, (iv) we exploit such a finding by spatially comparing the modeled displacements with the digital inventory of landslides triggered by past earthquakes over the study area [20] to construct a probability function relating the predicted D n values to a probability of failure (P(f )).
The procedure consists of [20]: distribution curve [20,46] (Equation ( 5)): where P(f ) is the probability of failure associated to each mean D n /D n_max value and m, a and b are regression constants to be determined through a best-fit solution between the model and the data.Such an approach allows providing a continuous curve function that associates a probability of failure P(f ) to each normalized Newmark displacement [20]; (e) associate the likelihood of landslide failure P(f ) to each grid cell using Equation (5).
It is expected that the P(f ) will increase for increasing D n /D n_max values, i.e., the number of grid cells falling on landslide areas increases with increasing Newmark displacements.
The modeling procedure requires a set of eight parameters, i.e., two state parameters (γ and γ w ); three geometrical parameters (t, α, and m); two strength parameters (c' and ϕ'); and one shaking parameter (PGA).

Coseismic Ground Displacements from InSAR
The coseismic displacement maps, obtained along ascending and descending orbits (Figure 4) present a negative displacement lobe (i.e., ground movements away from the satellite sensor) extending over an area of approximately 1 km 2 , south of Casamicciola and along the northern flank of Mt.Epomeo.The maximum displacements, ranging between −2 and −4 cm along the sensor line of sight (LoS), fairly overlap the area where surface movements and cracks were detected after the earthquake (Figure 1).The atmospheric correction performed with the GACOS tool was particularly useful to correct the displacement map for the descending interferogram.Indeed, the uncorrected wrapped interferogram (Supplementary Figure S2a) shows the presence of some interferometric fringes correlated with the topography, i.e., only stratified atmospheric delay is visible.The corrected wrapped interferogram (Supplementary Figure S2b) shows fringes associated with the coseismic ground deformations in the vicinity of Casamicciola, accordingly both with the ascending result and literature knowledge [9].wrapped interferogram (Supplementary Figure S2a) shows the presence of some interferometric fringes correlated with the topography, i.e., only stratified atmospheric delay is visible.The corrected wrapped interferogram (Supplementary Figure S2b) shows fringes associated with the coseismic ground deformations in the vicinity of Casamicciola, accordingly both with the ascending result and literature knowledge [9].The LoS displacement map along the ascending orbit (Figure 4a) shows positive displacements (i.e., ground movements towards the satellite sensor) diffused over the south-east of the island.Further residual atmospheric phase delay causes this signal since it is not present along the descending orbit (Figure 4b).Such disturbance does not affect the coseismic ground displacements since both the generated displacement maps present a consistent range increase pattern localized in the same area south of Casamicciola Terme, independently from the orbit pass, look angle, and interferometric pair, as confirmed by the detailed analysis on InSAR coseismic displacements performed by De Novellis et al. [9].

Landslide Failure Probability Estimation
The procedure described in Section 2 was applied to the Ischia earthquake in order to locate areas potentially affected by earthquake-induced landslides.

Available Dataset
The parameters required for the analysis were selected from up-to-date literature data available for Ischia Island.Geotechnical parameters, i.e., the state (γ) and strength (c' and φ') properties of the different lithologies covering the island, were taken from a detailed lithological map, 1:5000 scale, of the island (Figure 5a) [19].Each lithological unit consists of similar geotechnical properties obtained by merging available geological maps and fieldwork data, whose values are listed in Table 1.The range of the strength parameters is relatively broad; higher strength belongs to rocky lithologies, while lower strength applies to slope debris that covers most of the epicentral area (13, 14 and 15 in Figure 5a).Even though somewhat subjective, these parameters have been already successfully used as a first approach for the assessment of earthquake-induced landslide susceptibility on Ischia Island at a regional scale [19].The LoS displacement map along the ascending orbit (Figure 4a) shows positive displacements (i.e., ground movements towards the satellite sensor) diffused over the south-east of the island.Further residual atmospheric phase delay causes this signal since it is not present along the descending orbit (Figure 4b).Such disturbance does not affect the coseismic ground displacements since both the generated displacement maps present a consistent range increase pattern localized in the same area south of Casamicciola Terme, independently from the orbit pass, look angle, and interferometric pair, as confirmed by the detailed analysis on InSAR coseismic displacements performed by De Novellis et al. [9].

Landslide Failure Probability Estimation
The procedure described in Section 2 was applied to the Ischia earthquake in order to locate areas potentially affected by earthquake-induced landslides.

Available Dataset
The parameters required for the analysis were selected from up-to-date literature data available for Ischia Island.Geotechnical parameters, i.e., the state (γ) and strength (c' and ϕ') properties of the different lithologies covering the island, were taken from a detailed lithological map, 1:5000 scale, of the island (Figure 5a) [19].Each lithological unit consists of similar geotechnical properties obtained by merging available geological maps and fieldwork data, whose values are listed in Table 1.The range of the strength parameters is relatively broad; higher strength belongs to rocky lithologies, while lower strength applies to slope debris that covers most of the epicentral area (13, 14 and 15 in Figure 5a).Even though somewhat subjective, these parameters have been already successfully used as a first approach for the assessment of earthquake-induced landslide susceptibility on Ischia Island at a regional scale [19].[6] and the location of the landslide source areas [19].The yellow star identifies the Ischia earthquake epicenter.

Results of the Procedure Applied to Ischia Island
Ischia Island was discretized with a 10 m × 10 m grid (463,544 cells), according to the available DEM resolution.Attributes about the lithology type, unit weight, cohesion, friction angle, slope inclination, thickness, and PGA were assigned to each grid cell.
Figure 6a shows the resulting PGA distribution map of the Ischia Island, obtained by using  [6] and the location of the landslide source areas [19].The yellow star identifies the Ischia earthquake epicenter.
Table 1.Lithological units (for the reference number see Figure 5a) and main geotechnical parameters [19].The slope map (α) of the island (Figure 5b) was obtained from the 10-m resolution TINITALY DEM [47,48].The thickness of the slab involved by failure (t) was set to 3 m, according to the average depth of earthquake-induced landslides that have occurred at Ischia in the past [6,29,49] and to the thickness values adopted by Caccavale et al. [19] for the assessment of earthquake-induced landslide susceptibility on Ischia Island.We assume that the water table is at or below the failure surface (i.e., m = 0) since the earthquake occurred during the dry season (August 2017).Consequently, the third term of Equation ( 2) is zero.
The macroseismic intensity for the Ischia earthquake was estimated soon after the event (Figure 5c) [1] according to the European Macroseismic Scale (EMS).Evaluating the PGA with Equation (4) requires the macroseismic intensity defined according to the MCS scale.However, both EMS and MSC scales have 12 degrees of damage, thus, the equivalence between them is roughly one to one.Differences derived from the use of both intensity scales can substantially outweigh any actual differences in the scales themselves.Therefore, Equation (4) can be used to estimate PGA values according to both EMS and MCS intensity scales [50].
Finally, Figure 5d shows the inventory of the Ischia earthquake-induced landslides happened in the past.It mainly refers to the shallow landslide events that were inferred by geomorphological and archaeological studies, as well as from historical chronicles (see Rapolla et al. [6] and references therein).The mapped landslides consider both source and runout areas.However, Newmark displacements must be compared with landslide source areas only [20].Therefore, we assumed as landslide source areas those provided by Caccavale et al. [19] (Figure 5d), which are defined on a morphological basis.

Results of the Procedure Applied to Ischia Island
Ischia Island was discretized with a 10 m × 10 m grid (463,544 cells), according to the available DEM resolution.Attributes about the lithology type, unit weight, cohesion, friction angle, slope inclination, thickness, and PGA were assigned to each grid cell.
Figure 6a shows the resulting PGA distribution map of the Ischia Island, obtained by using Equation (4), plus its standard deviation (σ).We observe that the PGA estimate at the location of the IOCA seismic station (i.e., PGA ≈ 0.21 g) is in agreement with the horizontal PGA recorded during the Ischia earthquake (approximately 0.2 g, see Section 3.3) [41].Preliminary calculations of the FS map highlight that some areas of the island (approximately 8% of the grid cells) are statically unstable (i.e., FS ≤ 1).Such instability is due to the combination of steep slopes and low strength properties.However, most of these areas are affected or have been affected by active landslide phenomena [51].Therefore, we preferred to assign them a minimal fixed value of FS = 1.01, i.e., barely above equilibrium [19].The resulting static FS map (Figure 6b) shows that lower FS values correspond to grid cells with steep slopes and weak material properties.In detail, Preliminary calculations of the FS map highlight that some areas of the island (approximately 8% of the grid cells) are statically unstable (i.e., FS ≤ 1).Such instability is due to the combination of steep slopes and low strength properties.However, most of these areas are affected or have been affected by active landslide phenomena [51].Therefore, we preferred to assign them a minimal fixed value of FS = 1.01, i.e., barely above equilibrium [19].The resulting static FS map (Figure 6b) shows that lower FS values correspond to grid cells with steep slopes and weak material properties.In detail, slope areas close to the earthquake epicenter show low FS values that increase by decreasing the slope angle.
The insertion of the FS and α values (Figures 5b and 6b) into Equation ( 1) allowed calculating the a c values for each grid cell.The resulting a c map (Figure 6c) can be considered as a seismic landslide susceptibility map since it portrays a measure of intrinsic slope properties independent of any ground-shaking scenario.Lower a c values belong to areas with low FS and high slope angle.
Finally, D n values were computed for each grid cell by combining PGA and a c maps (Figure 6a,c, respectively) into Equation (3).Predicted D n values were divided by the maximum displacement estimate (D n_max ).The resulting D n /D n_max map is shown in Figure 6d, where surficial displacements are predicted over 12% of the whole island.Most of the predicted values are low (D n /D n_max < 0.1) and are located in the far field of the earthquake area, close to the Fontana, Serrara Fontana and Barano d'Ischia villages (Figure 6d).On the contrary, higher D n /D n_max values are located in the near field, on the north-western flank of Mt.Epomeo, close to the Casamicciola and Lacco Ameno villages.
The normalized displacements were compared with the landslide source areas of Figure 5d.In detail, D n /D n_max values were grouped into 50 classes (∆D n /D n_max = 0.02).Then, the proportion of the cells falling on landslide source areas, i.e., P(f ), was calculated for each class.Figure 6e plots the P(f ) distribution as a function of the computed D n /D n_max values.The data (red dots) show that P(f ) values increase almost monotonically for increasing D n /D n_max .This increase is faster for low D n /D n_max values and then levels off at approximately D n /D n_max = 0.5 and P(f ) ≈ 90%.
The best-fit curve according to Equation ( 5) is plotted as a continuous black line in Figure 6e, showing a good fit with the data (R 2 = 0.98, RMSE = 0.04).This curve was used to assign the P(f ) value to each D n /D n_max estimation for each grid cell, thus obtaining the P(f) map of Figure 6f.This map shows that the areas with the higher P(f ) values are located close to the epicenter of the Ischia earthquake.

Discussion
The proposed procedure allowed areas potentially affected by earthquake-induced landslides after the Ischia earthquake to be identified.The chart in Figure 6e demonstrates the effectiveness of the adopted method.Indeed, the higher the predicted displacement, the more significant the number of grid cells falling on landslide source areas (P(f )).This result agrees with the legacy effect of overlapping landslides [45].
Results highlight that, for the assumed slope geometries, strengths, and PGA distribution, the dynamic mobilization of the slopes in the near field of the Ischia earthquake epicenter is very likely (Figure 6f).The comparison with coseismic displacement maps of Figure 4a,b shows that areas with a P(f ) higher than 0.7 overlap the maximum displacement lobe inferred from InSAR measurements (Figure 7a,b).This analysis sets the high probability of a slope instability occurrence, though it does not allow for the provision of an absolute quantification of this phenomenon.It must be considered that InSAR data show deformations also where the P(f ) is close or equal to zero (Figure 7b), i.e., where only tectonic phenomena could be invoked to justify the observed deformation.It is, therefore, straightforward to assume that the observed ground displacement, despite its low intensity, could be attributed to the joint effect of tectonics and slope instabilities.Therefore, the movements detected by InSAR data can be interpreted as the combination of two different phenomena, i.e., the fault slip and the dynamic mobilization of slopes.
The performed analysis has a twofold implication.First, the proposed procedure allowed to identify critical areas potentially affected by slope movements after earthquakes.This information could contribute to better management of the seismic emergency, especially in areas exposed to high seismic and hydrogeological risks.Secondly, the ground displacements produced by the Ischia earthquake are the combination of both deep and shallow phenomena.Such information has to be adequately considered as far as the InSAR displacements are exploited for the research of the fault parameters with source inversion techniques.Indeed, non-tectonic shallow movements could introduce unknown biases in the estimation of the retrieved fault parameters.In general, when coseismic ground displacements are used to estimate the earthquake causative fault, a preliminary analysis must be included to locate areas affected by non-tectonic deformations, and eventually proceed through one of the following steps: (i) exclude displacements of those areas from the inversion procedure; (ii) weight the uncertainties associated to each displacement observation; (iii) model both contributions separately.The performed analysis has a twofold implication.First, the proposed procedure allowed to identify critical areas potentially affected by slope movements after earthquakes.This information could contribute to better management of the seismic emergency, especially in areas exposed to high seismic and hydrogeological risks.Secondly, the ground displacements produced by the Ischia earthquake are the combination of both deep and shallow phenomena.Such information has to be adequately considered as far as the InSAR displacements are exploited for the research of the fault parameters with source inversion techniques.Indeed, non-tectonic shallow movements could introduce unknown biases in the estimation of the retrieved fault parameters.In general, when coseismic ground displacements are used to estimate the earthquake causative fault, a preliminary analysis must be included to locate areas affected by non-tectonic deformations, and eventually proceed through one of the following steps: (i) exclude displacements of those areas from the inversion procedure; (ii) weight the uncertainties associated to each displacement observation; (iii) model both contributions separately.4b); The gridded pattern identifies areas with a landslide probability P(f ) ≥ 70% (black grid) and regions with a landslide probability P(f ) ≤ 70% (grey grid).Areas with P(f ) ≤ 10% are not shown.The yellow star identifies the Ischia earthquake epicenter.
The grid cells affected by slope deformations (Figures 6f and 7b) overlap fairly with the location and geometry of the observed coseismic fractures and collapse close to Casamicciola and Lacco Ameno villages (Figure 1) [5].However, areas with greater P(f ) values are located along Mt.Epomeo's flank, where no significant coseismic effects were observed.Such a discrepancy could be caused by the unmodeled lithological variability, which could lead to an overestimation of the areas affected by slope displacements; as well as to the presence of coseismic effects that have not been identified because of the impervious area.
Finally, it is worth discussing the significance of the results regarding the modeling assumptions and the parameters adopted in the analysis.
The assumed simplified rigid-block analysis does not allow quantifying accurately the amplitude of the earthquake-induced landslide movements; this is the reason why we normalized the computed displacement values.Indeed, estimating earthquake-induced displacements requires: (i) a detailed geological and geotechnical characterization of the outcropping lithologies; (ii) a rheological model accounting for the complex elastic-plastic behavior of soils subjected to cyclic loads [52,53]; (iii) a proper characterization of the seismic load.The performed analysis is adequate for the required target, i.e., the identification of areas potentially affected by earthquake-induced landslides at a regional scale.Different approaches at smaller scales must be used if an accurate estimation of the amplitude of earthquake-induced slope displacements is required [35,54,55].
The parameters adopted in the analysis, i.e., the slope angle, the depth of the failure surface, the PGA distribution, and the soil strength properties, present some uncertainties that must be investigated using a sensitivity analysis.The accuracy of the slope angle dramatically affects the outputs since the slope value is used twice in the proposed procedure (Equations ( 1) and ( 2)).Generally, slope maps from DEMs tend to underestimate steep slopes when the DEM resolution decreases.The relatively high resolution (10 m) of the adopted DEM allows for a reasonable estimation of local morphological variations without excessive smoothing of the slope values.Regardless, according to Equations ( 1) and ( 2), increasing the DEM resolution would increase the amplitude and spatial extent of the earthquake-induced landslide displacements, thus supporting the results obtained.
The sensitivity of the Newmark displacement spatial distribution to the other parameters, i.e., the slope thickness, the soil strength properties, and the PGA distribution, has been evaluated parametrically.Regarding the slope thickness (t), Equation (2) shows that if the slope thickness decreases, the FS values increase, and the spatial extent of Newmark displacements reduces.
Therefore, the calculation of D n /D n_max values is performed by halving the assumed t value (i.e., from 3 m to 1.5 m).The results (Figure 8a) show that the slope thickness has a negligible effect on the computed spatial distribution of surficial ground movements.Indeed, the grid cells affected by Newmark displacements decrease of approximately 25% in the far field of the earthquake epicenter, while the near-field movements are still predicted.Regarding PGA, low acceleration values lead to small ground movements (Equation ( 4)).Thus, we decreased the amplitude of the PGA pattern in Figure 6a by subtracting the standard deviation (σ = 0.222) in Equation (4).In this way, the PGA reduces by approximately 30%.Even in this case, the PGA has a negligible effect on the computed surficial movements since the grid cells affected by Newmark displacements decrease by roughly 16% in the far field areas only (Figure 8b).Finally, we increased the strength parameters (i.e., c' and ϕ) by 50% concerning the values listed in Table 1.Such an increase allows statically unstable areas (i.e., FS ≤ 1) for slopes less than 60 • to be avoided [20].Increasing the strength of lithologies has a substantial effect on the displacement prediction.Indeed, the pixels affected by Newmark displacements (Figure 8c) decreased by approximately 70% concerning those shown in Figure 6d.However, even assuming unrealistically high strength values, Newmark displacements are still predicted close to the epicentral area of the Ischia earthquake, thus confirming the presence of non-tectonic ground deformations.

Figure 2 .
Figure 2. Analogy between the block on the inclined plane and the landslide mass.(a) Representation of the Newmark [36] sliding block model; (b) infinite slope scheme.

Figure 2 .
Figure 2. Analogy between the block on the inclined plane and the landslide mass.(a) Representation of the Newmark [36] sliding block model; (b) infinite slope scheme.

Figure 3 .
Figure 3. Flowchart showing the required steps to produce the landslide failure probability map.

Figure 3 .
Figure 3. Flowchart showing the required steps to produce the landslide failure probability map.
(a) normalize the D n values at each grid cell by calculating the ratio between D n and the computed maximum displacement (D n_max ) from Equation (3) for the Ischia earthquake, i.e., D n /D n_max ; (b) group the D n /D n_max values into sequential intervals, e.g., the first interval contains D n /D n_max values ranging between 0 and 0.1, the second interval between 0.1 and 0.2, and so on; (c) compute for each interval the proportion of grid cells falling on areas affected by earthquake-induced landslides in the past.Such values correspond to the probability of landslide failure associated with each D n /D n_max interval, i.e., P(f ); (d) plot the obtained P(f) values versus the mean D n /D n_max of each interval and fit with a Weibull

Figure 4 .
Figure 4. Coseismic displacement maps from S1 SAR data of the Ischia earthquake.(a) Line of sight (LoS) displacements along the ascending orbit; (b) LoS displacements along the descending orbit.The yellow star identifies the Ischia earthquake epicenter.

Figure 4 .
Figure 4. Coseismic displacement maps from S1 SAR data of the Ischia earthquake.(a) Line of sight (LoS) displacements along the ascending orbit; (b) LoS displacements along the descending orbit.The yellow star identifies the Ischia earthquake epicenter.

Figure 5 .
Figure 5.The available dataset for the Ischia Island.(a) Geolithological map, scale 1:5000 form Caccavale et al. [19]; (b) Slope map derived from the 10-m resolution TINITALY digital elevation model (DEM); (c) Macroseismic intensity estimated according to the European Macroseismic Scale (EMS) after the 21 August 2017 earthquake [1]; (d) Map of areas affected by past earthquake-triggered landslides[6] and the location of the landslide source areas[19].The yellow star identifies the Ischia earthquake epicenter.

Figure 5 .
Figure 5.The available dataset for the Ischia Island.(a) Geolithological map, scale 1:5000 form Caccavale et al. [19]; (b) Slope map derived from the 10-m resolution TINITALY digital elevation model (DEM); (c) Macroseismic intensity estimated according to the European Macroseismic Scale (EMS) after the 21 August 2017 earthquake [1]; (d) Map of areas affected by past earthquake-triggered landslides[6] and the location of the landslide source areas[19].The yellow star identifies the Ischia earthquake epicenter.

Figure 6 .
Figure 6.Results for Ischia Island.(a) Estimated horizontal peak ground acceleration (PGA) distribution according to Equation (4); (b) static factor of safety (FS) distribution according to Equation (2); (c) ac distribution according to Equation (1); (d) Dn/Dn_max distribution according to Equation (3); (e) P(f) curve (black line) obtained by best fitting Newmark displacement data (red dots) with Equation (5); (f) P(f) map associated with the calculated Newmark displacements.The yellow star identifies the Ischia earthquake epicenter.

Figure 6 .
Figure 6.Results for Ischia Island.(a) Estimated horizontal peak ground acceleration (PGA) distribution according to Equation (4); (b) static factor of safety (FS) distribution according to Equation (2); (c) a c distribution according to Equation (1); (d) D n /D n_max distribution according to Equation (3); (e) P(f) curve (black line) obtained by best fitting Newmark displacement data (red dots) with Equation (5); (f) P(f) map associated with the calculated Newmark displacements.The yellow star identifies the Ischia earthquake epicenter.

Geosciences 2018, 8 , 19 Figure 7 .
Figure 7.Comparison of spatial patterns of InSAR coseismic ground displacements and landslide failure probability.(a) 3D view of Ischia Island; (b) detail of the area bounded by the dashed black rectangle in panel a.The color map shows the coseismic LoS displacements from S1 descending data (Figure4b); The gridded pattern identifies areas with a landslide probability P(f) ≥ 70% (black grid) and regions with a landslide probability P(f) ≤ 70% (grey grid).Areas with P(f) ≤ 10% are not shown.The yellow star identifies the Ischia earthquake epicenter.

Figure 7 .
Figure 7.Comparison of spatial patterns of InSAR coseismic ground displacements and landslide failure probability.(a) 3D view of Ischia Island; (b) detail of the area bounded by the dashed black rectangle in panel a.The color map shows the coseismic LoS displacements from S1 descending data (Figure4b); The gridded pattern identifies areas with a landslide probability P(f ) ≥ 70% (black grid) and regions with a landslide probability P(f ) ≤ 70% (grey grid).Areas with P(f ) ≤ 10% are not shown.The yellow star identifies the Ischia earthquake epicenter.

Figure 8 .
Figure 8.Estimated normalized Newmark displacements after being (a) halved the soil thickness; (b) decreased the PGA values; and (c) increased the strength parameters.The yellow star identifies the Ischia earthquake epicenter.

Figure 8 .
Figure 8.Estimated normalized Newmark displacements after being (a) halved the soil thickness; (b) decreased the PGA values; and (c) increased the strength parameters.The yellow star identifies the Ischia earthquake epicenter.