Multitemporal and Multisensor InSAR Analysis for Ground Displacement Field Assessment at Ischia Volcanic Island (Italy)

: Volcanic islands are often affected by ground displacement such as slope instability, due to their peculiar morphology. This is the case of Ischia Island (Naples, Italy) dominated by the Mt. Epomeo (787 m a.s.l.), a volcano-tectonic horst located in the central portion of the island. This study aims to follow a long temporal evolution of ground deformations on the island through the interferometric analysis of satellite SAR data. Different datasets, acquired during Envisat, COSMO-SkyMed and Sentinel-1 satellite missions, are for the ﬁrst time processed in order to obtain the island ground deformations during a time interval spanning 17 years, from November 2002 to December 2019. In detail, the multitemporal differential interferometry technique, named small baseline subset, is applied to produce the ground displacement maps and the associated displacement time series. The results, validated through the analysis and the comparison with a set of GPS measurements, show that the northwestern side of Mt. Epomeo is the sector of the island characterized by the highest subsidence movements (maximum vertical displacement of 218 mm) with velocities ranging from 10 to 20 mm/yr. Finally, the displacement time series allow us to correlate the measured ground deformations with the seismic swarm started with the Mw 3.9 earthquake that occurred on 21 August 2017. Such correlations highlight an acceleration of the ground, following the mainshock, characterized by a subsidence displacement rate of 0.12 mm/day that returned to pre-earthquake levels (0.03 mm/day) after 6 months from the event.


Introduction
Volcanic islands are often affected by natural hazards, such as eruptions, earthquakes, landslides and mass movements, frequently interrelated to each other and affecting the same zones. For this reason, it is critically important to continuously study these areas, trying to deeply understand natural hazards, leveraging multidisciplinary knowledge and aiming to mitigate the risk of devastating effects risk.
Ischia Island is a quiescent volcano located in Southern Italy within the Gulf of Naples ( Figure 1) and it represents an interesting case study for natural hazards interconnections.
The island is characterized by a long eruptive history, lasting more than 150 ka [1]. The oldest volcanic deposits, individuated along the coastline, represent the remnants of an ancient caldera generated by a massive collapse between 150 and 75 ka. Stratigraphic evidence highlighted that a significant ignimbritic eruption dated 55 ka produced a further caldera collapse and the emplacement of the Green Tuff, the best-known volcanic deposit on evidence highlighted that a significant ignimbritic eruption dated 55 ka produced a further caldera collapse and the emplacement of the Green Tuff, the best-known volcanic deposit on the island. After this event, the morpho-geodynamic evolution of the island was essentially controlled by a volcano-tectonic activity, which uplifted the Mt. Epomeo block inside the pre-existing caldera depression over the past 33 ka [2]. This 900-1100 m uplift is a quite unique case of intra-calderic resurgence, very rarely recognized in other volcanic areas [3,4]. The last eruptive activity that occurred in 1302 AD is represented by the Arso eruption. This event produced a 2.5 km length lava flow in the eastern sector of the island with a volume of 0.03 km 3 [1] and triggered numerous landslides [5]. The seismic activity of Ischia Island is known with a good degree of confidence starting from the 13th century [5] and it is mainly concentrated in the northern area of the island (Casamicciola Terme and Lacco Ameno localities). Analysing the effects of the historical earthquakes, it was observed that even low magnitudes events, e.g., the Mw 2.8 earthquake that occurred on Casamicciola Terme in 1863 with an epicentral intensity MCS (Mercalli-Cancani-Sieberg) equal to "V", were able to activate ground movements [5]. Following the 1883 destructive event [6] and, subsequently, sporadic and low-magnitude earthquakes, the most recent and severe event is dated 21 August 2017. It was characterized by a seismic swarm started with a Mw 3.9 earthquake [1,5]. In this case, the conjunction of a very shallow seismic source, local ground amplification and the poor resilience of the buildings, caused unexpected severity [7]. The occurrence of this earthquake allowed scientists to acquire and integrate different geophysical and geological datasets (e.g., historical seismicity, macroseismic observations, geodetic and topographic instrumental information) as well as a detailed mapping of the coseismic effects [5].
In this work, the areas on Ischia Island affected by ground displacement are investigated for the first time in a 17-year time span by applying an InSAR technique to several datasets acquired during Envisat, COSMO-SkyMed and Sentinel-1 satellite missions. We started with those zones historically manifesting slope instabilities [1], such The seismic activity of Ischia Island is known with a good degree of confidence starting from the 13th century [5] and it is mainly concentrated in the northern area of the island (Casamicciola Terme and Lacco Ameno localities). Analysing the effects of the historical earthquakes, it was observed that even low magnitudes events, e.g., the Mw 2.8 earthquake that occurred on Casamicciola Terme in 1863 with an epicentral intensity MCS (Mercalli-Cancani-Sieberg) equal to "V", were able to activate ground movements [5]. Following the 1883 destructive event [6] and, subsequently, sporadic and low-magnitude earthquakes, the most recent and severe event is dated 21 August 2017. It was characterized by a seismic swarm started with a Mw 3.9 earthquake [1,5]. In this case, the conjunction of a very shallow seismic source, local ground amplification and the poor resilience of the buildings, caused unexpected severity [7]. The occurrence of this earthquake allowed scientists to acquire and integrate different geophysical and geological datasets (e.g., historical seismicity, macroseismic observations, geodetic and topographic instrumental information) as well as a detailed mapping of the coseismic effects [5].
In this work, the areas on Ischia Island affected by ground displacement are investigated for the first time in a 17-year time span by applying an InSAR technique to several datasets acquired during Envisat, COSMO-SkyMed and Sentinel-1 satellite missions. We started with those zones historically manifesting slope instabilities [1], such as the sector located south of the Lacco Ameno municipality (Figure 2), which, together with Casamicciola Terme, was the most damaged area during the Mw 3.9 earthquake from 21 August 2017 [8].
as the sector located south of the Lacco Ameno municipality (Figure 2), which, together with Casamicciola Terme, was the most damaged area during the Mw 3.9 earthquake from 21 August 2017 [8].

Figure 2.
Hillshade relief of Ischia with the epicentres of the 2017 seismic sequence [9] and the area examined through the displacement time series analysis in addition to the municipalities and their limits, contour lines, faults [10] and landslides phenomena mapped within the island (IFFI landslides inventory). Geographic coordinate system is UTM 33N, WGS 84.
In the next sections, the multitemporal and multisensor interferometric processing results are presented, followed by the GPS comparisons. Moreover, for the first time, the ground deformation effects caused by the 2017 seismic swarm are investigated over time in the time series analysis and discussion sections.

Geo-Volcanological Setting
Ischia is a volcanic island belonging to the Phlegraean Volcanic District (PVD), which comprises the Campi Flegrei caldera, Procida and Vivara islands ( Figure 1). PVD is the most widespread active volcanic system of the Mediterranean area, established inside the Campanian Plain half-graben following the opening of the Tyrrhenian Sea [11] started in the late Miocene (9-10 Ma) [12]. The volcanism of Ischia is related to the extensional tectonics that caused the eastward and anticlockwise migration of the Italian peninsula, with the consequent activation of NW-SE and NE-SW normal faults systems, which permitted magmas to reach the surface, feeding volcanism. The island is mainly composed by alkali-trachytic volcanic rocks [13], landslide deposits and sedimentary terrigenous rocks that testify its complex history, characterized by the occurrence of volcanic, tectonic and slope instability phases [1]. The eruptive activity of Ischia lasted from the Late Pleistocene up to historical times [13] and the oldest evidences were dated back to about 150 ka, whereas the most recent outcrops were created by the 1302 eruption (Arso lava flow) [14].
The central part of the island is dominated by Mt. Epomeo, a volcano-tectonic horst (787 m a.s.l.) whose origin is due to the resurgence of the caldera formed after the large explosive eruption (55 ka BP), which emplaced Mt. Epomeo's Green Tuff [14]. The  [9] and the area examined through the displacement time series analysis in addition to the municipalities and their limits, contour lines, faults [10] and landslides phenomena mapped within the island (IFFI landslides inventory). Geographic coordinate system is UTM 33N, WGS 84.
In the next sections, the multitemporal and multisensor interferometric processing results are presented, followed by the GPS comparisons. Moreover, for the first time, the ground deformation effects caused by the 2017 seismic swarm are investigated over time in the time series analysis and discussion sections.

Geo-Volcanological Setting
Ischia is a volcanic island belonging to the Phlegraean Volcanic District (PVD), which comprises the Campi Flegrei caldera, Procida and Vivara islands ( Figure 1). PVD is the most widespread active volcanic system of the Mediterranean area, established inside the Campanian Plain half-graben following the opening of the Tyrrhenian Sea [11] started in the late Miocene (9-10 Ma) [12]. The volcanism of Ischia is related to the extensional tectonics that caused the eastward and anticlockwise migration of the Italian peninsula, with the consequent activation of NW-SE and NE-SW normal faults systems, which permitted magmas to reach the surface, feeding volcanism. The island is mainly composed by alkalitrachytic volcanic rocks [13], landslide deposits and sedimentary terrigenous rocks that testify its complex history, characterized by the occurrence of volcanic, tectonic and slope instability phases [1]. The eruptive activity of Ischia lasted from the Late Pleistocene up to historical times [13] and the oldest evidences were dated back to about 150 ka, whereas the most recent outcrops were created by the 1302 eruption (Arso lava flow) [14].
The central part of the island is dominated by Mt. Epomeo, a volcano-tectonic horst (787 m a.s.l.) whose origin is due to the resurgence of the caldera formed after the large explosive eruption (55 ka BP), which emplaced Mt. Epomeo's Green Tuff [14]. The resurgent block, tilted southward and split in smaller portions, is bordered by a system of faults and fractures trending mostly NW-SE, NE-SW and N-S.
Currently, the island morphology is mainly driven by slopes instability affecting Mt. Epomeo's flanks. Along such slopes, the mobilization of volcanic and sedimentary material can generate landslides phenomena that, according to recent studies [1], have been grouped The numerous historical landslides that occurred in the densely populated island (more than 60,000 inhabitants spread over less than 50 km 2 ) have been recorded and organized in the Inventario dei Fenomeni Franosi in Italia (IFFI) catalogue ( Figure 2). These phenomena generally occur following extreme meteorological events, eruptions or earthquakes, as it happened during the last seismic sequence at the end of August 2017. Indeed, many coseismic effects were detected in conjunction with this event along Mt. Epomeo's flanks, including rockfalls and debris flows in volcanoclastic deposits [8].
Therefore, the detailed study of Ischia's ground deformations, paying attention to Mt. Epomeo's slopes instability, is considered crucial to understand the complex dynamic of the island, often influenced by seismic events.

Datasets and Methods
The numerous synthetic aperture radar (SAR)-dedicated space missions since the launch of ERS-1, have provided a huge amount of data available for interferometric processing. In detail, using a multitemporal series of SAR images acquired from slightly different satellite positions and in distinct time instants, it is possible to apply differential interferometry techniques (DInSAR). The use of a great number of data allows us to obtain parameters not measurable with a single interferogram, granting the description of the temporal evolution of the deformations [15].
The two main algorithms used in the multitemporal DInSAR techniques are known as persistent scatterers [16] and small baseline subset (SBAS) [17]. The first one is based on the identification and displacement measurements of point targets, objects characterized by high temporal stability of the backscattered signal and by a well-defined geometry like corner reflectors, buildings and stable rocks. Therefore, these points are known as the persistent scatterers. Instead, the second algorithm performs the analysis on natural distributed targets and not very geometrically defined objects, allowing the maximum spatial coverage of the outputs within the study area [18]. Detailed information on the SBAS algorithm was provided in [17,19].
The SBAS approach has been adopted here to process several datasets spanning a frame of 17 years from the end of 2002 to the end of 2019 quite continuously. The datasets comprise stacks of Envisat (11/2002-08/2010), COSMO-SkyMed (02/2011-08/2017) and Sentinel-1A (01/2015-12/2019) data, covering the entire area of Ischia Island. With the aim to obtain correct ground deformations, the 30 m Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) has been used to remove the topography contribution from the generated interferograms. All the available data have been processed with SARscape software (ENVI platform). The ascending and descending tracks of both Envisat and Sentinel-1 datasets were processed in order to evaluate either the measurement along the satellite's line of sight (LOS) and the vertical and the horizontal (east-west) components of the displacement. Moreover, the obtained displacement time series allowed us to identify areas affected by ground movements and to measure their relative velocities.
The average velocities resulting from the SAR processing have been validated through the comparison with the displacement rates of the global positioning system (GPS) measurements, acquired in correspondence with the INGV monitoring geodetic network [20,21]. Table 1 summarizes the main characteristics of the processed data: name of the satellite mission, acquisition geometry, track of the satellite, number of images used during each processing, number of interferograms evaluated as reliable for the generation of the outputs, time interval considered in each processing, ground resolution and angle of incidence at the scene centre.

Results and GPS Comparisons
The results and the comparative graphs are shown in Figures 3-7. In particular, the ground displacement maps, results of the InSAR data processing, are shown in Figures 3, 5a and 6 where the ground displacement rate is in millimetres per year (mm/yr). In detail, the negative LOS velocities indicate ground movements away from the satellite sensor, while positive LOS velocities imply ground movements approaching the satellite sensor. All the displacement maps display the highest negative LOS velocities in: (i) the northwestern sector of the island, (ii) along the north side of Mt. Epomeo and (iii) in the area south of the Lacco Ameno and Casamicciola Terme municipalities. It is evident that the displacement maps obtained by Sentinel-1A ( Figure 6) allow the greatest areal coverage. These results also detect significant movements near the southern slope of Mt. Epomeo.
In order to validate such results, the availability of 20 geodetic stations already existing in the island was exploited. The comparisons were made between the average displacement rates derived by using the SBAS technique and those obtained from the GPS measurements reprojected along the relative satellite LOS. This reprojection was based on the incidence line of sight (ILOS) and azimuth line of sight (ALOS) angles. In addition, the error associated with the GPS measurements was appropriately reprojected along the satellite LOS. In order to calculate the interferometric results errors (σInSAR), which are derived from parameters such as coherence and wavelength, the following formula (1), where γ is the interferometric coherence and λ is the wavelength, was used [22,23]: The comparisons are shown in the comparative graphs in Figures 4, 5b, and 7 for Envisat, COSMO-SkyMed and Sentinel-1A, respectively. The accuracies derived from reprojected GPS displacement rates and those derived from the SAR data processing were compared by calculating the root mean square error (RMSE). In detail, the accuracy was calculated as the mean square discrepancy between the values of the GPS data and those measured in the SAR data processing at each measurement station ( Table 2).       . Comparison between GPS displacement rates reprojected on Envisat descending (blue triangles) and ascending (red triangles) LOS and, respectively, displacement rates evaluated through SBAS processing: descending track (light blue diamonds) and ascending track (orange diamonds).    Comparison between GPS displacement rates reprojected on Sentinel-1A descending (blue triangles) and ascending (red triangles) LOS and, respectively, displacement rates evaluated through SBAS processing: descending track (light blue diamonds) and ascending track (orange diamonds).
The availability of both descending and ascending acquisition geometries for Envisat and Sentinel-1A datasets allowed the calculation of the horizontal (east-west) and vertical components of the displacement [24].
The maps of the vertical displacement rate show negative values if the ground motion is downward, while positive values if surface uplift occurs. The maps of the displacement rate in the horizontal direction display negative values if the motion is oriented to the west, while positive values if the movement is oriented towards the east. Figure 8 shows both the horizontal (Figure 8a,c) and vertical (Figure 8b,d) displacement maps generated for the Envisat (Figure 8a,b) and Sentinel-1A (Figure 8c   Comparison between GPS displacement rates reprojected on Sentinel-1A descending (blue triangles) and ascending (red triangles) LOS and, respectively, displacement rates evaluated through SBAS processing: descending track (light blue diamonds) and ascending track (orange diamonds).
The availability of both descending and ascending acquisition geometries for Envisat and Sentinel-1A datasets allowed the calculation of the horizontal (east-west) and vertical components of the displacement [24].
The maps of the vertical displacement rate show negative values if the ground motion is downward, while positive values if surface uplift occurs. The maps of the displacement rate in the horizontal direction display negative values if the motion is oriented to the west, while positive values if the movement is oriented towards the east. Figure 8 shows both the horizontal (Figure 8a,c) and vertical (Figure 8b,d) displacement maps generated for the Envisat (Figure 8a,b) and Sentinel-1A (Figure 8c  Comparison between GPS displacement rates reprojected on Sentinel-1A descending (blue triangles) and ascending (red triangles) LOS and, respectively, displacement rates evaluated through SBAS processing: descending track (light blue diamonds) and ascending track (orange diamonds).
The availability of both descending and ascending acquisition geometries for Envisat and Sentinel-1A datasets allowed the calculation of the horizontal (east-west) and vertical components of the displacement [24].
The maps of the vertical displacement rate show negative values if the ground motion is downward, while positive values if surface uplift occurs. The maps of the displacement rate in the horizontal direction display negative values if the motion is oriented to the west, while positive values if the movement is oriented towards the east. Figure 8 shows both the horizontal (Figure 8a,c) and vertical (Figure 8b,d) displacement maps generated for the Envisat (Figure 8a,b) and Sentinel-1A (Figure 8c These results were also validated by comparing the horizontal ( Figure 9) and vertical ( Figure 10) displacement rates with those obtained directly from GPS measurements acquired at each geodetic station. As performed for LOS comparisons, the accuracy was quantified through the root mean square error calculation. Table 3 summarizes the RMSE retrieved for each comparison.  These results were also validated by comparing the horizontal ( Figure 9) and vertical ( Figure 10) displacement rates with those obtained directly from GPS measurements acquired at each geodetic station. As performed for LOS comparisons, the accuracy was quantified through the root mean square error calculation. Table 3 summarizes the RMSE retrieved for each comparison.  Figure 9. Comparison between GPS horizontal displacement rates (blue triangles) and those retrieved through SBAS processing using Sentinel-1A data (red diamonds) and Envisat data (green diamonds).

Figure 10.
Comparison between GPS vertical displacement rates (blue triangles) and those retrieved through SBAS processing using Sentinel-1A data (red diamonds) and Envisat data (green diamonds).

Time Series Analysis
The processing of SAR data also provided the opportunity to generate the displacement time series relating to critical areas in order to follow the evolution of slope movements over time. In detail, the time series were generated on the area of maximum deformation identified by the totality of the elaborations, i.e., the area located on the northwestern side of Mt. Epomeo, south of the Lacco Ameno municipality (Figures 2 and  11).
In order to combine the displacements obtained by the three different datasets (Envisat, COSMO-SkyMed and Sentinel-1A), along their respective descending orbits, it was necessary to reproject the displacements according to the same incidence angle. For this aim, the incidence angle of the Sentinel-1 sensor (35.9°) was taken as a reference. The displacements measured along the descending LOS of the Envisat sensor were then divided by the cosine of the difference between the Sentinel-1 and the Envisat (23.6°) incidence angles. The same procedure was applied considering the COSMO-SkyMed descending displacements and the relative angle of incidence (27.4°). Moreover, the cumulative displacements over the entire period were obtained by estimating an average velocity in the absence of temporal coverage data. . Comparison between GPS horizontal displacement rates (blue triangles) and those retrieved through SBAS processing using Sentinel-1A data (red diamonds) and Envisat data (green diamonds).

Figure 10.
Comparison between GPS vertical displacement rates (blue triangles) and those retrieved through SBAS processing using Sentinel-1A data (red diamonds) and Envisat data (green diamonds).

Time Series Analysis
The processing of SAR data also provided the opportunity to generate the displacement time series relating to critical areas in order to follow the evolution of slope movements over time. In detail, the time series were generated on the area of maximum deformation identified by the totality of the elaborations, i.e., the area located on the northwestern side of Mt. Epomeo, south of the Lacco Ameno municipality (Figures 2 and  11).
In order to combine the displacements obtained by the three different datasets (Envisat, COSMO-SkyMed and Sentinel-1A), along their respective descending orbits, it was necessary to reproject the displacements according to the same incidence angle. For this aim, the incidence angle of the Sentinel-1 sensor (35.9°) was taken as a reference. The displacements measured along the descending LOS of the Envisat sensor were then divided by the cosine of the difference between the Sentinel-1 and the Envisat (23.6°) incidence angles. The same procedure was applied considering the COSMO-SkyMed descending displacements and the relative angle of incidence (27.4°). Moreover, the cumulative displacements over the entire period were obtained by estimating an average velocity in the absence of temporal coverage data. Figure 10. Comparison between GPS vertical displacement rates (blue triangles) and those retrieved through SBAS processing using Sentinel-1A data (red diamonds) and Envisat data (green diamonds).

Time Series Analysis
The processing of SAR data also provided the opportunity to generate the displacement time series relating to critical areas in order to follow the evolution of slope movements over time. In detail, the time series were generated on the area of maximum deformation identified by the totality of the elaborations, i.e., the area located on the northwestern side of Mt. Epomeo, south of the Lacco Ameno municipality (Figures 2 and 11).
In order to combine the displacements obtained by the three different datasets (Envisat, COSMO-SkyMed and Sentinel-1A), along their respective descending orbits, it was necessary to reproject the displacements according to the same incidence angle. For this aim, the incidence angle of the Sentinel-1 sensor (35.9 • ) was taken as a reference. The displacements measured along the descending LOS of the Envisat sensor were then divided by the cosine of the difference between the Sentinel-1 and the Envisat (23.6 • ) incidence angles. The same procedure was applied considering the COSMO-SkyMed descending displacements and the relative angle of incidence (27.4 • ). Moreover, the cumulative displacements over the entire period were obtained by estimating an average velocity in the absence of temporal coverage data. Figure 11. Detail of the maximum deformation zone and location of the area examined through the time series. Sentinel-1A vertical displacement map, epicentres [9] and faults [10] are also visible.
In Figure 12, the complete time series, including all the descending datasets available for a time interval of 17 years are shown, whereas in Figure 13a the vertical displacement time series generated for the Envisat and Sentinel-1 datasets are plotted. The plot in Figure  12 shows an almost constant negative displacement rate, measured along the descending LOS of all the datasets, as is also the case in Figure 13a. The only appreciable acceleration in both plots is identified in the time corresponding to the Mw 3.9 earthquake on 21 August 2017. This event is neatly recognized in the time series generated for vertical displacements and is shown in detail in Figure 13b.   [9] and faults [10] are also visible.
In Figure 12, the complete time series, including all the descending datasets available for a time interval of 17 years are shown, whereas in Figure 13a the vertical displacement time series generated for the Envisat and Sentinel-1 datasets are plotted. The plot in Figure 12 shows an almost constant negative displacement rate, measured along the descending LOS of all the datasets, as is also the case in Figure 13a. The only appreciable acceleration in both plots is identified in the time corresponding to the Mw 3.9 earthquake on 21 August 2017. This event is neatly recognized in the time series generated for vertical displacements and is shown in detail in Figure 13b.
Remote Sens. 2021, 13, 4253 10 of 15 Figure 11. Detail of the maximum deformation zone and location of the area examined through the time series. Sentinel-1A vertical displacement map, epicentres [9] and faults [10] are also visible.
In Figure 12, the complete time series, including all the descending datasets available for a time interval of 17 years are shown, whereas in Figure 13a the vertical displacement time series generated for the Envisat and Sentinel-1 datasets are plotted. The plot in Figure  12 shows an almost constant negative displacement rate, measured along the descending LOS of all the datasets, as is also the case in Figure 13a. The only appreciable acceleration in both plots is identified in the time corresponding to the Mw 3.9 earthquake on 21 August 2017. This event is neatly recognized in the time series generated for vertical displacements and is shown in detail in Figure 13b.

Discussion
All the displacement maps presented in this study and generated by processing different SAR data (Envisat, COSMO-SkyMed and Sentinel-1) allowed us to define the maximum deformation area within Ischia Island. In particular, LOS, horizontal and vertical ground displacement maps referring to the period November 2002-December 2019, identified the highest deformations localized in the northwestern slope of Mt. Epomeo. The central portion of this area, investigated through the time series analysis, was localized in the Lacco Ameno municipality, at a distance of about 1 km west from the epicentre of the Mw 3.9 earthquake that occurred on 21 August 2017, approximately at the western edge of the fault plane modelled in previous studies [9]. This zone and the Casamicciola Terme hilly area were the most damaged sectors of the island by this seismic event [25], whereas in the Fango locality many ruptures were observed immediately after the earthquake [8]. Moreover, considering the recent assessment of the epicentral environmental seismic intensity (ESI) [26] based on the distribution of the primary (surface ruptures and permanent displacement caused directly by the seismogenic source), secondary coseismic geological data (landslides, hydrological variations) and

Discussion
All the displacement maps presented in this study and generated by processing different SAR data (Envisat, COSMO-SkyMed and Sentinel-1) allowed us to define the maximum deformation area within Ischia Island. In particular, LOS, horizontal and vertical ground displacement maps referring to the period November 2002-December 2019, identified the highest deformations localized in the northwestern slope of Mt. Epomeo. The central portion of this area, investigated through the time series analysis, was localized in the Lacco Ameno municipality, at a distance of about 1 km west from the epicentre of the Mw 3.9 earthquake that occurred on 21 August 2017, approximately at the western edge of the fault plane modelled in previous studies [9]. This zone and the Casamicciola Terme hilly area were the most damaged sectors of the island by this seismic event [25], whereas in the Fango locality many ruptures were observed immediately after the earthquake [8]. Moreover, considering the recent assessment of the epicentral environmental seismic intensity (ESI) [26] based on the distribution of the primary (surface ruptures and permanent displacement caused directly by the seismogenic source), secondary coseismic geological data (landslides, hydrological variations) and collapse of drywall [5], the area where the SBAS analysis individuated the highest displacements corresponded to the zone with the maximum degree of intensity (VII ESI 2007 scale). Furthermore, the displacement maps obtained by the Sentinel-1 data processing, with unprecedented multitemporal coverage, highlighted significant movements also along the western and southern slopes of Mt. Epomeo. In fact, these sides are characterized by different geolithological compositions than that of the eastern flank of the Ischia Island, which is instead composed by volcanics younger than 10 ka [7] and where the lack of remote sensing deformation signals is recognized. Precisely, in the western side of Mt. Epomeo, corresponding to the observed slopes instability phenomena, low-strength volcanoclastic deposits occur along the Donna Rachele fumaroles area. In this zone, CO 2 and H 2 S gases rise through the NW-SE normal faults bordering the western side of Mt. Epomeo (Figure 2), these fractures connect highly fractured deep lavas (hosting the hydrothermal aquifer) with the surface [27]. The resulting circulation alters the shallower impermeable pyroclastic cover (Green Tuff and Citara Tuff) [27] weakening its mechanical properties [28][29][30]. Therefore, the displacements recorded by the Sentinel-1A sensor are presumably attributed to the volcanic system of the island, which is still active [10] and currently characterized by the emission of volcanic gas mainly in the western slope of Mt. Epomeo, corresponding to the fumarolic belt of Donna Rachele. Additionally, the observed deformations along the southern slope of Mt. Epomeo are most likely related to the complex relation between the tectonics and the morphology. In fact, by this side, the central resurgent block is marked by a double monocline fold facing toward south-southwest, which is interrupted by depressions, partially filled with debris avalanche, debris flow and volcanoclastic deposits, subsequently cut by deep canyons [31]. In particular, the volcanoclastic sequence of the southern slope of Mt. Epomeo comprises debris avalanche, debris flow and rockslide deposits interlayered with alluvial and shoreline sediments, extending from Mt. Epomeo's summit to the southern coast [31], according with the distribution of the displacements recorded by the Sentinel-1A sensor at the highly irregular and unstable southern flank of the resurgent block.
Combining the average velocities measured along the ascending and descending LOS of the Sentinel-1 and Envisat sensors it was possible to retrieve the vertical and horizontal displacement rates. In the northwestern slope of Mt. Epomeo, the maximum vertical displacement retrieved considering an almost continuous time interval of 17 years was 218 mm and the associated subsidence rate ranged from 10 to 20 mm/yr. These values are consistent with previous works, confirming the deflation trend of the area. In fact, these rates are in agreement with the SBAS results of C-band ERS-1/2 and Envisat data spanning the 1992-2010 time interval [32]. In addition, the high values of vertical downward deformation found in the northern and southern flanks of Mt. Epomeo are relative to two distinct areas characterized by an evident subsidence individuated by levelling surveys carried out between 1990 and 2003 [33] and high-precision levelling survey realized following the earthquake on November 2017 [34]. Such subsidence was explained through crack closure processes along two main ENE-WSW and E-W pre-existing faults, which represent the preferred degassing pathway of the hydrothermal system beneath Mt. Epomeo [1].
The totality of the results deriving from the SAR data processing were also compared with the ground measurements obtained by the global navigation satellite system (GNSS) stations located on the island. Taking into account the precision associated with each measurement, the comparisons between GPS and SAR data proved that our results were accurate. Indeed, the maximum derived RMSE for the comparisons between GPS and horizontal and vertical derived displacement rates was 2.7 mm/yr.
Finally, the analyses made on the displacement time series allowed us to evaluate the surface ground displacements over time in the maximum deformation zone detected through SBAS processing. The time series, representing the vertical evolution of the displacement, emphasize that the subsidence rate in the area of interest is constantly progressing with average values of 0.03 mm/day. Additionally, it is noted that the only significant acceleration begins immediately after the 21 August 2017 mainshock, speeding up the subsidence rate to 0.12 mm/day for 6 months. After about 180 days from the seismic event, in February 2018, the subsidence rate returned to pre-earthquake levels. This suggests that seismicity and related effects occurring on the island accelerate the subsidence rate specifically in the northern slope of Mt. Epomeo, according to recent findings [7,35], but also that the subsidence process of the resurgent block is in continuous progression, supporting the latest studies that attribute the subsidence phenomenon to the degassing of the magma body located at a depth of 2 km [36].

Conclusions
In this paper displacements on Ischia volcanic island ground surface were investigated over a continuous time interval of almost two decades, from November 2002 to December 2019. In particular, the SBAS multitemporal differential interferometry technique was applied to several SAR datasets with different temporal and spatial resolutions (X-band COSMO-SkyMed, C-band Envisat and Sentinel-1 data). Ground displacement maps, characterized by a wide land coverage, allowed us to identify the area on the island most affected by deformation, corresponding to the northwestern sector of Mt. Epomeo. In that zone, a detailed analysis of the ground displacements was performed and for the first time the effects over time caused by the 2017 seismic swarm were individuated. In particular, the ground deformation, following the mainshock, was accelerated by a subsidence rate of 0.12 mm/day that returned to pre-earthquake levels (0.03 mm/day) after 6 months from the event.
This study represents a first step toward an interdisciplinary approach [37] joining remote sensing elaborations, geolithological features and morphometric analyses aimed at the determination of predisposing features that help drive paroxysmal phenomena and identifying the prone areas where they could occur (e.g., landslides).
Finally, the peculiar tectonic, stratigraphic and hydrothermal conditions of the localized area along the northwestern, western and southern Mt. Epomeo's flanks, together with the determination of the mechanical properties of the outcropping materials, probably at the base of observed slope instability, could be the object of further on-field geophysical investigations, providing very useful information for the prevention and mitigation of such phenomena, recognized as the most frequent natural hazards on Ischia Island.

Funding:
The present work is supported by the FRASI-Integrated and multi-scale approach for the definition of seismic-induced landslide hazard in the Italian territory-research project, funded by the MATTM. Copernicus Sentinel-1 data are freely distributed by the European Space Agency and available at the Copernicus Open Access Hub. The Envisat data are provided free of charge by the European Space Agency. The COSMO-SkyMed data was obtained by the Italian Space Agency through the project id.213. Project carried out using CSK ® products, © ASI (Italian Space Agency), delivered under an ASI license to use.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.