Catching Geomorphological Response to Volcanic Activity on Steep Slope Volcanoes Using Multi ‐ Platform Remote Sensing

: The geomorphological evolution of the volcanic Island of Stromboli (Italy) between July 2010 and June 2019 has been reconstructed by using multi ‐ temporal, multi ‐ platform remote sensing data. Digital elevation models (DEMs) from PLÉIADES ‐ 1 tri ‐ stereo images and from Light Detection and Ranging (LiDAR) acquisitions allowed for topographic changes estimation. Data were comprised of high ‐ spatial ‐ resolution (QUICKBIRD) and moderate spatial resolution (SENTINEL ‐ 2) satellite images that allowed for the mapping of areas that were affected by major lithological and morphological changes. PLÉIADES tri ‐ stereo and LiDAR DEMs have been quantitatively and qualitatively compared and, although there are artefacts in the smaller structures (e.g., ridges and valleys), there is still a clear consistency between the two DEMs for the larger structures (as the main valleys and ridges). The period between July 2010 and May 2012 showed only minor changes consisting of volcanoclastic sedimentation and some overflows outside the crater. Otherwise, between May 2012 and May 2017, large topographic changes occurred that were related to the emplacement of the 2014 lava flow in the NE part of the Sciara del Fuoco and to the accumulation of a volcaniclastic wedge in the central part of the Sciara del Fuoco. Between 2017 and 2019, minor changes were again detected due to small accumulation next to the crater terrace and the erosion in lower Sciara del Fuoco.

In this work, geomorphological variations on high-gradient volcano flanks in response to shifts in volcanic activity have been identified by means of remote sensing techniques. The Sciara del Fuoco (SdF) depression on the island of Stromboli (Italy) was studied, as it is the optimal test-site for monitoring the effect of volcanic eruption on steep-slope volcano flank because i) it is affected by persistent volcanic activity, ii) it is prone to mass-wasting phenomena, and iii) it is one of the best studied and, among all, most monitored volcano sites on Earth, providing exceptional validation data and ground-truth constraints. A large set of multi-temporal data of the SdF was acquired, including Light Detection and Ranging (LiDAR) data, tri-stereo PLÉIADES-1 imagery, high-spatial-resolution (HSR) optical imagery (QUICKBIRD and PLÉIADES-1), and moderate spatial resolution (MSR) SENTINEL-2 multi-spectral instruments (MSI) imagery. Multi-temporal data permitted the mapping of areas affected by major lithological and morphological changes, as well as the volumes of deposited/eroded material. Moreover, the results of this study were strengthened by the integration of previous geomorphological reconstructions [23,28,41,43,44]. The results led to the identification of topographical variations and geomorphological processes that occurred in response to the variation in eruptive intensity.

Study Area
The 916 m-high Stromboli Island is the emerged portion of a ~3000 m-high stratovolcano located in the Tyrrhenian Sea off the southern coast of Italy ( Figure 1). The volcano experienced several large mass-wasting phenomena, which formed two scars, one on its NW flank (Sciara del Fuoco; SdF) and the other one on its SE flank, with a consequence bilateral flank instability [45][46][47]. The SdF depression is filled with volcaniclastic deposits and lavas [48][49][50] that are emitted from a summit crater terrace located at ≈ 750 m a.s.l., and from ephemeral vents within the SdF [51]. The distinctive persistent Strombolian activity is characterized by intermittent explosions from three vent areas (NE, SW and central) that are located in a summit crater terrace [52]. This activity, showing intensity and frequency fluctuations over time, is often punctuated by lava overflows from the crater terrace [44], and/or by flank eruptions, with the outpouring of lava flows from ephemeral vents [32]. Recent eruptive activity and slope instability events are summarized in Table 1. The last flank eruption, characterized by lava emission from an ephemeral vent located at 650 m a.s.l., started on 7 August 2014, was preceded by two months of increased Strombolian activity and several lava overflows from the craters, and lasted until 13 November 2014 [41,[53][54][55][56][57][58][59][60]. Both in December 2017 and 2018, the eruptive activity increased (spattering and overflows), mainly in the North-Eastern Crater (NEC) area, which led to the growth of the NEC cone outside of the crater terrace [1,2]. This activity was related to frequent, small scale, instability phenomena, such as rockfalls and gravel slides that evolved into gravel flows [2].  Slope instability phenomena at Stromboli are classified, on the base of their size and movement, into three types [2,62]: 1) "deep-seated gravitational slope deformations" evolving to "rock or debris avalanches" from the SdF (volumes > 10 6 m 3 ); 2) "rock (rotational or planar) slides" evolving to "rock avalanches" from the SdF (volumes ≈ 10 6 m 3 ); and 3) "rock falls" or "gravel/debris slides" evolving to "gravel/debris flows" (volumes ≈ 10 5 m 3 ). The greatest hazard associated with large-scale landslides (type-1) at Stromboli is their ability to generate tsunamis, whose effects can propagate far away from the source areas [63,64].

Geomorphological Mapping
In the last few decades, different remote sensing techniques have been used to catch environmental responses to volcanic activity (Table 2). Here, the geomorphological mapping of the SdF was based on the integration of field surveys, multi-temporal, satellite-borne HSR optical imagery, morphometric analysis and topographic change detection (TCD). The topographic data used were four different DEMs that were derived from 2010 and 2012 LiDAR data and 2017 and 2019 PLÉIADES-1 images. The HSR-orthorectified imagery dataset (Table 3 and     Moreover, multi-temporal SENTINEL-2 MSI images were used to constrain the main geomorphological events in a higher temporal resolution (Table 4 and Figure 3). The spatial resolution of SENTINEL-2 MSI bands is 10-60 m and has a coverage between −56 and +84 degrees latitude with a 290 km swath width. The minimum revisit time at the equator is 5 days, which increases a when maximum cloud cover is set (max. = 5% for this study). Several band combination and ratios have been used to enhance contrasts between features, as well as to reduce the variations in topographic illumination.

Topographic Source Data and Methods
The detection of topographic changes measured by differentiating pre-, sin-, and post-eruptive DEMs is nowadays considered the most suitable method to accurately quantify the volume of material that is emplaced or removed during volcanic eruptions (Table 1). In this study, the topographic change detection of Stromboli Island from 2010 to 2019 was performed by using DEMs that were derived from data that were acquired by airborne LiDAR data and PLÉIADES-1 tri-stereo satellite imagery.
Airborne LiDAR is an active remote sensing system that integrates data from a laser ranging and scanning unit and a position and orientation system (POS), that is generally mounted on a small aircraft, and that is used to acquire a dense cloud of points with known coordinates of the imaged surface. The point cloud is then generally used to generate an HSR (1 m) DEM. In this work, the first LiDAR data of the Stromboli volcano were acquired on 22 July 2010 by using an Optech Airborne Laser Terrain Mapper (ALTM) GEMINI laser altimeter (http://optech.on.ca). The nominal instrumental planimetric and vertical accuracies were in the ranges of 0.8-1.35 and 0.25-0.35 m, respectively. The point cloud had a mean point density of 1 pts/m 2 , which allowed us to generate a 1-m resolution DEM. The second LiDAR-DEM was obtained by elaborating the 3D point cloud that was acquired during an airborne survey that was carried out on May 2012. The data were acquired by using the Leica ADS80 sensor, which has instrumental vertical and horizontal accuracy of 0.10-0.20 and 0.25 m, respectively [41,65]. The acquired point cloud had a mean point density of 8 pt/m 2 .
The 2017 and 2019 DEMs were derived from the tri-stereo optical imagery that was acquired by the PLÉIADES-1 satellite. The PLÉIADES-1 constellation is composed by two satellites, PLÉIADES-1A (PHR1A) and PLÉIADES-1B (PHR1B). These satellites can sense three or more synchronous images of the same area, with an angle variable between ~6° and ~28°. The stereoscopic triplet is composed of three nearly simultaneously acquired images, one backward looking, one forward looking, and a third near-nadir image [35,41,66]. The tri-stereo images of this work were 100% cloud free, but the presence of smoke emitted that was by the vents prevented the correct acquisition of the crater terrace. The total areal acquired was 58 km 2 . This acquisition mode permitted the obtainment of a DEM through the photogrammetric processing of the three stereo images (tri-stereo mode). The dataset also comprised panchromatic and multispectral very high resolution (VHR) optical imagery with spatial resolutions of 0.5 and 2 m, respectively. To assess the accuracy of the heights and their horizontal position in the PLÉIADES-1 DEM, ground control points (GCPs) were collected on the map database (cartographic XY standard deviation: 0.15 m). Tie points were automatically collected in the images (Table 5 and Supplementary Material). A block adjustment including all the satellite scenes was performed. The block adjustment was validated when the following accuracy was achieved: (i) a pixel xy bias smaller than 0.3 pixels; (ii) a pixel xy standard deviation smaller than 0.3 pixels; and (iii) a pixel xy maximum smaller than 2 pixels.

DEMs Co-Registration and Topographic Change Detection
Topographic change detection (TCD) using multi-temporal DEMs was performed by differencing two DEMs of the same area that were derived from data that were taken at different times. This calculation is often affected by errors depending on the presence of a mismatching between two DEMs of the same area, which leads to artefact h [30]. This error can be detected and reduced by the measurement and minimization of the DEM differences in areas where the two DEMs are supposed to be equal, i.e., those areas that are not reasonably affected by relevant natural changes.
In this work, DEM to DEM coregistration was based on the minimization of the root mean square (RMS) error between one DEM to the other by iteratively varying the three angles of rotation, the translation, and the magnification or reduction factor of one DEM [39,67,68] by using a custom-made algorithm based on the MINUIT minimization library [69]. MINUIT is a tool that can be used to find the minimum value of multiparameter functions and that can be freely downloaded (http://www.cern.ch/minuit). Coregistration was performed by following the key steps listed below: 1. Preliminary DEMs comparison: Two DEMs were compared, and their differences were visualized in map in order to detect the areas that were affected by topographic changes due to volcanic activity and/or natural phenomena as well as the distribution of errors outside the areas affected by real changes. 2. Detection of target areas and initial RMS displacement calculation: A given number of areas without relevant natural changes around the region of interest, in this case the Sciara del Fuoco scar, were selected for DEMs coregistration and the initial RMS displacement between the two DEMs calculated for these areas. 3. Calculation of the minimization parameters: The minimization algorithm based on MINUIT was launched and the minimization parameters, and the final RMS displacement was calculated for the target areas. 4. DEMs coregistration: Minimization parameters were used for coregistering one DEM (slave) on the DEM previously chosen as a reference (master). 5. Selection of independent areas and calculation of RMS displacement error : Since the areas used for calculating were targeted by the minimization process, could not be used for the error calculations of thicknesses and volumes. Independent areas, not affected by natural changes and that were as close as possible to the region of interest, were used to check the residual mismatching between the two DEMs. Note that in this case, the resulting RMS displacement error was the RMS residuals (in elevation) between one DEM and another arbitrary one chosen as a reference, rather than a true absolute error [68]. The differences between two successive properly coregistered DEMs were used to detect and outline the extent of the areas that were affected by topographic changes (Figure 4) and to calculate the volume and thickness variation inside them (Table 6). Moreover, the perimeters of the observed changes were double checked by using the orthorectified (airborne and PLÉIADES-1) images. The areas where the volumes gained and lost by the SdF were calculated, and the selected profiles are shown in Figures 4 and 5.    Figure 4d. Labels describe the observed phenomena, and the sectors shown in Figure 4 are reported in brackets. Erosion is illustrated in cold tones, and deposition is illustrated in warm tones. The most recent surface, i.e., the 2019 DEM, is marked by a squared pattern.
The volume (V) emplaced or lost between two acquisitions was calculated from the DEM difference according to the equation: ∑ Δ Δ [30], where x is the grid step and zi is the height variation within grid cell i. These values were then summed up for all the cells in the selected areas in which the volume changes were calculated. An upper bound on the error for the volume estimate was given by assigning to each pixel the maximum possible error, i.e., , , where A is the investigated area. Thus, the errors provided here are reasonably overestimated [30].

DEMs Comparison
A preliminary comparison between 2010 and 2012 LiDAR DEMs had = 0.46 m, which was considered low enough to not require any further correction between the two LiDAR DEMs. In this case, even if no correction was performed, was lower than , as it had a value of 0.35 m ( Table 6) Table 6). The 2012-2017 LiDAR vs PLÉIADES coregistration was already performed by using a similar method in [41], which provided some slightly different values than this work, i.e., = 3.37 m, = 1.27 m. This difference was due to the different areas that were selected for calculating RMS displacement that included almost the entire island with the exclusion of the SdF. In [41], was used as . The 2010-2019 coregistration provided the following values: = 3.54 m, = 0.61 m, = 0.75 m ( Table  6). The 2019-2017 comparison was calculated by comparing the PLÉIADES DEMs before and after their coregistration with LiDAR DEMs. The following RMS displacement areas were calculated: = 1.06 m, = 0.46 m, = 0.53 m. Though both PLÉIADES tri-stereo and LiDAR DEMs had 1 m of spatial resolution, the difference in accuracy can be easily perceived visually. Figure 6a shows a comparison between the hillshading derived from both DEMs of an SdF sector. Due its slope and the type of deposits, this area is not easy to be correctly sensed [70]. The LiDAR modelled the surface without any macroscopic artefacts, allowing for the detection of the complex system of small valleys and ridges inside it. On the other hand, the PLÉIADES's DEM showed a diffuse noise, which did not allow for us to distinguish the internal structure of the SdF. This noise was better highlighted by the sky view factor 2 (Figure 6c,d) and the openness down map (Figure 6e,f). Sky view factor 2 (SVF2) is a variant of sky view factor (SVF). The SVF was described by a solid angle (Ω) open to the sky and is expressed in terms of the total sky view possible from any given pixel, i.e., SVF = Ω/2π [71]. It is a fraction of the sky visible from each pixel and ranges from one to zero. SVF2 also accounts for negative values, allowing for the presence of angles under the horizon, and it can reach a maximum possible value of 2 [71]. SVF2 enhances the perception of the relative height of surface elements, thus helping to detect crests and ridges. The LiDAR DEM clearly showed these structures inside the SdF, thus permitting us to follow their downslope directions (Figure 6c). In the PLÉIADES's DEM, several artefacts cross-cut the SdF, thus masking the smaller structures (Figure 6d). Openness down is a measure of belowground openness [71][72][73]. It has high values inside valleys, gullies, and craters, effectively resulting in the detection of cracks and fractures [68]. An openness down map derived from LiDAR showed welldefined valley directed downslope (Figure 6e), while the PLÉIADES one highlighted the artefact "cobble" texture inside the SdF (Figure 6f). These data agree with a previous analysis that was carried out in [74,75], whose authors found that the elevation differences between LiDAR and PLÉIADES data in urban areas and farmlands are largely on the order ± 1 m.
At the same time, by observing the profiles of Figure 5, it is possible to appreciate the coherence among the PLÉIADES and LiDAR DEMs, where no natural morphological changes were expected (i.e. at the northwestern sectors of each profile). A good agreement between the PLÉIADES and LIDAR data was described by [74,75]. However, the latter also reported the presence of a noise in the PLÉIADES tri-stereo DEMs, suggesting its use only for features greater than 2 m.

Morphological Changes
The integration of the TCD and the analysis based of the optical images were the basis for the definition of the short-term (between July 2010 and June 2019) geomorphological evolution of the SdF steep-slope, which is frequently affected by an accumulation of volcanic material, erosion, and masswasting phenomena. Major variations are illustrated in Figures 4 and 5, and volume changes are shown in Table 6.
The period between July 2010 and May 2012 showed changes (Figures 4a and 5) only related to phases of more intense eruptive activity (stronger/frequent explosions and overflows). The lava overflowed from the South-Western Crater (SWC) area on 11-12 December 2011 [44]. It diverged in two branches. The southernmost was 400 m long and reached 560 m a.s.l., and the northernmost was longer at 730 m and reached a minimum elevation of about 207 m a.s.l. (Figure 4a and 5, region 1). From the same crater, a small overflow also occurred on 18 August 2011, which emplaced over the 11-12 December overflow [44]. Both overflow emplaced a volume of 0.118 ± 0.022 × 10 6 m 3 of lava outside the crater terrace, with a mean thickness of 1.850 m. Overflow also occurred on 1-2 August 2011 from the NEC (Figure 4a, region 3) [54]. The lava run over a pre-existing debris talus accumulated below the NEC (Figure 4a, region 2 and Figure 5, P1). The lava front reached a minimum elevation of 530 m after having flowed a distance of 470 m. The deposits forming the debris talus later covered the feeding channel of the 1-2 August overflow, which was no longer visible at the time of the 2012 LiDAR survey. Though it was not possible to distinguish the debris talus from the lava where they overlapped, their volumes, which must be considered as an approximation of the real values, were provided,. Region 2 of Figure 4a, which was most related to the debris talus, had a volume of 0.289 ± 0.019 × 10 6 m 3 for a mean thickness of 5.465 m. The terminal part of the lava flow, which included the distal sector of the feeding channel and the dispersed flow (Figure 4a, region 3), had a volume of 0.055 ± 0.005 × 10 6 m 3 with a mean thickness of 3.870 m. In profile 1 of Figure 5, it is also possible to observe a small area that was strongly affected by erosion, which caused the loss of a rock wall of about 14 m. This was due to the disaggregation of the March 2002 lava flow to just north of the 2011 overflow, which was affected by headward erosion.
During the period between May 2012 and May 2017, large topographic changes occurred. The largest variations were related to the emplacement of the 2014 lava flow field in the NE part of the SdF (Figure 4b, region 4). The 2014 lava flow volume that was calculated here was 2.697 ± 0.190 × 10 6 m 3 . This value differs from the volume that was calculated in [41], which was 3.07 ± 0.37 × 10 6 m 3 , because the NEC talus (sector 5) was not included here. In any case, part of this talus was included in the volume calculation because it was not possible to distinguish between the talus and the proximal sector of the lava (Figure 5, P2). The lava average thickness was 11.964 m. In the median portion of the SdF, the lava filled pre-existing depressions (  (Figure 5, P3). During the same time period, the SdF has been affected by massive erosive phenomena. The older lava has been affected by headward erosion, which is visible in profile 1 of Figure 5, carried on to lose material. Besides this, two large areas have been affected by erosion: one constrained between the 2014 lava and the volcaniclastic wedge ( Figure  4b, region 8) and the other south of the volcaniclastic wedge (Figure 4b, sector 9). Erosion in region 8 caused the loss of 0.153 ± 0.044 × 10 6 m 3 of volcaniclastic material for a mean thickness loss of just below 3 m. Locally, the erosion had high impact, arriving at 10 m of thickness lost ( Figure 5, P2). The erosion rate in this sector was about 2.5 × 10 3 m 3 /month. Region 9 was wider than region 8; it lost 0.425 ± 0.104 × 10 6 m 3 of material for a mean thickness loss of just below 3.5 m, arriving locally even just over 10 m. The erosion rate in this sector was about 7.1 × 10 3 m 3 /month. The material eroded in region 9 partially accumulated downslope for a volume of 0.031 ± 0.010 × 10 6 m 3 (Figure 4b, sector 7).
Between 2017 and 2019, Strombolian activity kept increasing the volume of the NEC talus of 0.061 ± 0.008 × 10 6 m 3 for a mean increase in thickness of 3.9 m (Figure 4c, region 10). Besides this, two eroded areas in the lower SdF were detected (Figure 4c, sector 11 and 12) that corresponded to sectors 6 and 7 of Figure 4b

Discussion
The 2010-2019 volcanic activity fluctuated in the intensity and eruptive style from mildexplosive (Strombolian) to effusive [2,41,44,53,[55][56][57][58][59]. This led to an oscillating production of volcaniclastic material and the emplacement of a lava flow filed and several lava overflows [23,41]. During periods characterized by a low intensity/frequency of Strombolian activity, the production of materials ejected from the crater terrace towards the SdF was generally low, and erosion was the prevailing process, mainly affecting the central sector of the SdF. After the 2014 flank eruption, the gravitational re-adjustment of the lava breccia sometimes produced rockfalls [23,41].
The data presented here allow for the reconstruction of the 2010-2019 geomorphological evolution of the SdF slope ( Figure 7; Table 7). The main source area of volcanic material is the crater area, where the deposition of scoria bombs-to-ash deposits raised the crater terrace floor, constructing scoria cones within it, and accumulating spatter agglutinates, and forming hornitos. The (volcaniclastic/lava) material temporarily stored in the proximal and medial areas was usually remobilized due to gravity (in high-angle zones), such as erosion or landslides. In medial areas, characterized by the presence of volcaniclastic deposits (i.e., the central part of the SdF), the flowing of lavas or grain flows could have produced erosion and remobilization of material. The distal area was also a temporary storage area where the phenomena of accumulation from higher elevation counteracted the sea erosion. With the addition of the results of previous works on the submarine part [47,48,60,[75][76][77], it was possible to identify other main temporary storage areas (the near-coast volcaniclastic apron identified by [60,76,77]) up to the most distal zones, which have been reached sporadically by the most energetic phenomena [47,48].   The combination of this work with the results of [23] provides an evolutionary framework, at least for the period after the 2007 eruption. This work shows that the emplacement of the lava field of 2014 was preceded and accompanied by the accumulation of a volcanoclastic wedge from the midslope to the coastline (and also below) [60,78,79]. In this case, the positive anomalies detected by the TCD must have been produced by the lava flow field and also by the volcanoclastic materials that were deposited during periods of high explosive activity. These results are partly in contradiction with the statement made in [32]. The authors of that article suggested that the accumulation of material on the SdF depends only on the effusive phases, and, therefore, the positive anomalies in the TCDs should be ascribed only to the emplacement of lava flows. The quantification of this volcanoclastic wedge is relevant, as is the submarine volcanoclastic apron, because it was composed of the same material that was involved in the 30 December, 2002 tsunamigenic landslide [77] in addition to being located in the same area.

Conclusions
In this work, the geomorphological responses to the volcanic activity at the Stromboli steepslope volcano were identified by means of DEM comparisons. The results were strengthened by the analysis of medium and high-resolution satellite optical images. Multi-temporal lithology and landform variations depended on the shift in spatial and temporal distribution of processes. Stromboli, as many other frequently active volcanos, is characterized by periods of low volcaniclastic productivity (no or low-level volcanic activity) when erosion prevails, and this erosion if mainly concentrated on the loose, older volcaniclastic deposits. On the other side, the increase in the explosive activity produces more volcaniclastic material with a consequent increase in sedimentation along the slope. The main source area of volcanic material is the crater area, where volcanic activity is used to accumulating spatter agglutinates, forming hornitos. Lava accumulation can occur in different positions, depending on the location of the eruptive vent and the effusive rate.
This work shows that the differences between DEMs, if not acquired at high temporal frequency, should be accompanied by other observations (such as those based on optical images) in order not to incur the error of attributing topographical variations to a single event/a single lithology. In the case of Stromboli, the topographic changes between May 2012 and May 2017 are not due only to the emplacement of the lava field 2014 but also due to the accumulation of the volcanoclastic wedge. The presence of material with different technical characteristics is fundamental in a volcano flank that is prone to landslide, since it allows for the identification of more unstable volumes.
The comparison between the accuracy levels of PLÉIADES and LiDAR DEMs that was performed here should be considered preliminary and not exhaustive. Further analysis must involve the source data, the used modelling methods, and data that were acquired with various methods at different times. Having stated this, a first indication on the differences between two largely used methods for modelling topography in volcanic areas is provided here. This analysis can be used by the volcanological community and the civil protection authorities in case of a cost-benefit analysis for planning the best method for updating the topography and quantifying the morphological changes of an active, steep-slope volcano.
This study demonstrates how remote sensing data is useful for monitoring eruptive and geomorphological phenomena in active volcanic areas. This approach can be applied to other volcanic areas that are characterized by similar phenomena, but at the same time it can be generalized and applied in other geological contexts that are characterized by steep slopes.