Submarine and Subaerial Morphological Changes Associated with the 2014 Eruption at Stromboli Island

: Stromboli is an active insular volcano located in the Southern Tyrrhenian Sea and its recent volcanic activity is mostly conﬁned within the Sciara del Fuoco (SdF, hereafter), a 2-km wide subaerial–submarine collapse scar, which morphologically dominates the NW ﬂank of the ediﬁce. In August-November 2014, an effusive eruption occurred along the steep SdF slope, with multiple lava ﬂows reaching the sea. The integration of multisensor remote sensing data, including lidar, photogrammetric, bathymetric surveys coupled with SAR amplitude images collected before and after the 2014 eruption enabled to reconstruct the dynamics of the lava ﬂows through the main morphological changes of the whole SdF slope. Well-deﬁned and steep-sided ridges were created by lava ﬂows during the early stages of the eruption, when effusion rates were high, favoring the penetration into the sea of lava ﬂows as coherent bodies. Differently, fan-shaped features were emplaced during the declining stage of the eruption or in relation to lava overﬂows and associated gravel ﬂows, suggesting the prevalence of volcaniclastic breccias with respect to coherent lava ﬂows. The estimated volume of eruptive products emplaced on the SdF slope during the 2014 eruption, accounts for about 3.7 × 10 6 m 3 , 18% of which is in the submarine setting. This ﬁgure is different with respect to the previous 2007 eruption at Stromboli, when a large lava submarine delta formed. This discrepancy can be mainly related to the different elevation of the main vents feeding lava ﬂows during the 2007 eruption (around 400 m) and the 2014 eruption (around 650 m). Besides slope accretion, instability processes were detected both in the subaerial and submarine SdF slope. Submarine slope failure mobilized at least 6 × 10 5 m 3 of volcaniclastic material, representing the largest instability event detected since the 2007 lava delta emplacement.


Introduction
Stromboli is the north-easternmost island of the Aeolian Archipelago in the Southern Tyrrhenian Sea, representing the tip of a large, 3000 m-high volcanic edifice (Figure 1). Stromboli and Panarea volcanic edifices are part of a 45 km long volcanic area developed along a NE-SW regional extensional fault system (Figure 1b) above the thinned continental crust of the Calabrian Arc [1,2]. The evolution of Stromboli island occurred in the last 85 ka (apart from the older Strombolicchio edifice dated at about 200 ka) and can be divided in six main growth stages, with magmas ranging from calc-alkaline to potassic series, typical be divided in six main growth stages, with magmas ranging from calc-alkaline to potassic series, typical of an arc-volcanic setting, whose formation can be related to the subduction of the Ionian crust beneath the Calabrian Arc ( [2] and references therein).
The Stromboli edifice is characterized by a quasi-bilateral symmetry with respect to a tectonically controlled NE-SW axis (Figure 1c), which corresponds to the main orientation and distribution of dykes, eruptive fissures, secondary vents and active summit craters during most of the volcano evolution [1][2][3]. Dyke intrusion along this axis also promoted the development of recurrent lateral collapses along the NW and SE flanks of the volcano [2,[4][5]. Specifically, multiple sectors collapses affected the NW flank of the edifice in the last 13 ka, the last of them likely occurred during medieval times, leading to the development of the Sciara del Fuoco scar, SdF hereafter ( Figure 1c) [2].  [6]), with the location of the Aeolian Islands (black box); (b) zoom of the Aeolian archipelago in the Southern Tyrrhenian Sea, with the geographical location of Stromboli Island (in red) and the main SW-NE regional system (black dashed lines) controlling the Stromboli-Panarea alignment; (c) shaded relief map and contour lines (equidistance 500 m) of the Stromboli edifice (location in Figure 1b), with the indication of the 2002 tsunamigenic landslide scars (see also Figure 4a) and the main rift axis oriented along the SW-NE direction.
The SdF is a very steep subaerial/submarine partially filled depression that funnels eruptive material produced at the summit craters into the sea, where a large volcaniclastic fan, composed by debris avalanche deposits and overlying turbidite deposits, is recognizable to over 3000 m water depth (mwd, hereafter) [7]. This volcaniclastic system is mostly fed by the persistent activity, occurring at the summit craters (South-Western, Central and North-Eastern Craters, SWC, CC and NEC, respectively, in Figure 1a). This The SdF is a very steep subaerial/submarine partially filled depression that funnels eruptive material produced at the summit craters into the sea, where a large volcaniclastic fan, composed by debris avalanche deposits and overlying turbidite deposits, is recognizable to over 3000 m water depth (mwd, hereafter) [7]. This volcaniclastic system is mostly fed by the persistent activity, occurring at the summit craters (South-Western, Central and North-Eastern Craters, SWC, CC and NEC, respectively, in Figure 1a). This activity, mostly consisting of low-to mid-energy Strombolian explosions, shows intensity and Remote Sens. 2021, 13,2043 3 of 22 frequency fluctuations over time and is punctuated by: (a) high-energy explosive events, called paroxysms [8,9], often associated with pyroclastic flows as observed in 2019 [10], and (b) lava overflows from the summit crater area [11,12] and/or effusive flank eruptions, as recently occurred in 1985-1986, 2002-2003, 2007 and 2014 [13][14][15][16].
Dyke intrusions and effusive eruptions also induce stress changes on the volcanic flanks and can generate slope instability at different spatial scales, as observed during the initial phases of the last four flank eruptions [16][17][18]. Of these events, only the 2002-2003 eruption was associated with the occurrence of tsunamigenic submarine and subaerial landslides (red dashed line in Figure 1c; [19,20]). The tsunami waves related to the 2002 submarine and subaerial landslides severely impacted the Stromboli coastline with a maximum runup of 10 m [21], evidencing the high hazard associated with similar events that have repeatedly occurred also during the previous century [22,23]. Since then, SdF slopes have been carefully monitored through repeated topographic surveys for their submarine [24][25][26] and subaerial part [27]. These surveys are part of a more complex monitoring network including: ground-based thermal infrared data [3]; infrasonic monitoring [28], seismic network [29], borehole strainmeters [10], dilatometric sensors [30], ground tiltmeters [31], plume and soil gas geochemical monitoring [32,33], two ground-based interferometric synthetic aperture radar sites [34] and satellite thermal monitoring [35]. However, except for the study of the 2002-2003 eruption and related tsunamigenic landslides [36], no comprehensive subaerial and submarine studies of the eruptive events occurring at SdF were performed until now.
In this paper, we present an integrated analysis of multisensor remote sensing data for the characterization of the morphological changes associated with the 2014 effusive eruption occurred at Stromboli volcano, during which lava flows entered the sea in the NW shallow offshore of the island. In detail, the bathymetric comparison between multibeam surveys performed before and after the 2014 eruption have been integrated with lidar, photogrammetric and SAR amplitude images [16], aimed at fully understanding the morphological evolution of the SdF slope involved in this eruption. The integration of submarine and subaerial DEMs also allows one to estimate the volume of the eruptive products emplaced on the SdF slope during the 2014 eruption and compare it with the values previously estimated using satellite and ground based thermal observation [28,37], thus providing a crossvalidation of the two methodological approaches.
More generally, the results of this study are important to better understand the behavior and morphological evolution of lava flows penetrating into the sea, mainly in relation to the paucity of bathymetric monitoring of coastal areas during eruptive crisis [25,38], with observations mainly limited to the subaerial part of lava deltas [27,39,40] or scuba dives at water depth less than 50 m [41,42]. This is a very important issue considering the possible hazard associated with these events in active insular and coastal volcanoes.

The 2014 Eruption: Chronology and Subaerial Morphological Analysis
The 2014 effusive eruption began on the morning of 7 August 2014 and ended on 13 November 2014. The eruption was preceded by two months of more intense activity and anomalous values in geophysical and geochemical parameters [28,32,43], with stronger and more frequent Strombolian explosions, overflows from the crater terrace and small landslides along the subaerial SdF triggered by crater-rim collapse or by remobilization of volcaniclastic material by overflows (Figure 2a, [44]). The most frequent explosive activity and the lava overflows were generated mainly by the NEC, with the accumulation and remobilization of volcaniclastic material along the central part of the SdF [16,44]. On 6 August 2014, starting in the early afternoon, a series of overflows from the NEC area reached the sea. These events were associated with crater-rim collapse landslides that occurred in the same sector (Figure 2a; [16]). The actual onset of the eruption (between 3 and 5 a.m. UTC on 7 August 2014) was marked by the propagation of magma from the crater terrace towards the NEC, with the development of an eruptive fracture that fed a The propagation of the fracture triggered the collapse of the volcaniclastic talus below the NEC [16]. The first days of the 2014 eruption were characterized by a high effusive rate, which decreased dramatically in the following days, and then remained low for the duration of the eruption [28,45]. During this initial phase, a series of lava flows reached the sea in different points, always located in the northernmost area of the SdF (Figure 2b-d; [46]). Subsequently, with the decrease of the effusive rate, the lava flows did not propagate far from the vent, with lava fronts that stood at around 400 m a.s.l.. This produced the typical morphology of Stromboli lavas with vent positioned at high altitudes [14], with the development of a proximal shield, a medial zone fed by small flows, and characterized by frequent lava crumbling down slope and producing a debris field, with debris moving downslope, to the coastline (Figure 2c,d; [46]).
The volume of lava emitted by the 2014 eruption has been estimated to be 7.4 × 10 6 m 3 by [37] and 5.5 × 10 6 m 3 by [28], with 2.697 ± 0.190 × 10 6 m 3 of lava emplaced on the subaerial slope (black dashed line in Figure 1b; [47]). After the 2014 eruption, at least for the period 2015-2016, the eruptive activity remained at very low levels [29], mainly characterized by sporadic, low intensity Strombolian explosions. In this period the volcaniclastic sedimentation from the crater area towards the SdF was significantly reduced, and erosion of the 2014 lava field mainly occurred, with small landslides mobilizing volumes in the order of 10 3 -10 4 m 3 [34].

Data and Methods
Data used for this research mainly rely on the analysis of morphological changes detected by using high-resolution topo-bathymetric surveys collected before and after the 2014 eruption along the SdF collapse scar and time-lapse SAR (synthetic aperture radar) amplitude images collected between May and August 2014 ( Figure 3 and Table 1). The propagation of the fracture triggered the collapse of the volcaniclastic talus below the NEC [16]. The first days of the 2014 eruption were characterized by a high effusive rate, which decreased dramatically in the following days, and then remained low for the duration of the eruption [28,45]. During this initial phase, a series of lava flows reached the sea in different points, always located in the northernmost area of the SdF (Figure 2b,c,d; [46]). Subsequently, with the decrease of the effusive rate, the lava flows did not propagate far from the vent, with lava fronts that stood at around 400 m asl. This produced the typical morphology of Stromboli lavas with vent positioned at high altitudes [14], with the development of a proximal shield, a medial zone fed by small flows, and characterized by frequent lava crumbling down slope and producing a debris field, with debris moving downslope, to the coastline (Figure 2c,d; [46]).
The volume of lava emitted by the 2014 eruption has been estimated to be 7.4 × 10 6 m 3 by [37] and 5.5 × 10 6 m 3 by [28], with 2.697 ± 0.190 × 10 6 m 3 of lava emplaced on the subaerial slope (black dashed line in Figure 1b; [47]). After the 2014 eruption, at least for the period 2015-2016, the eruptive activity remained at very low levels [29], mainly characterized by sporadic, low intensity Strombolian explosions. In this period the volcaniclastic sedimentation from the crater area towards the SdF was significantly reduced, and erosion of the 2014 lava field mainly occurred, with small landslides mobilizing volumes in the order of 10 3 -10 4 m 3 [34].

Data and Methods
Data used for this research mainly rely on the analysis of morphological changes detected by using high-resolution topo-bathymetric surveys collected before and after the 2014 eruption along the SdF collapse scar and time-lapse SAR (synthetic aperture radar) amplitude images collected between May and August 2014 ( Figure 3 and Table 1).

Lidar, Photogrammetric and Bathymetric Surveys
The difference map, obtained by comparing repeated digital elevation models (DEMs, hereafter), enable us to quantify the morphological change occurred on the SdF slope during the time interval between the surveys, with positive and negative values associated with slope accretion and erosion, respectively. The volumes associated with slope changes were obtained by integrating the difference in depth over the area of interest through the software Global Mapper. The reliability of computed volumes depends on the resolution and accuracy of the DEMs, as described below both for the submarine and subaerial slope.
For the submarine slope, seafloor changes were computed through the difference between multibeam bathymetries collected in 2013 and 2016 onboard the R.V. Urania and Minerva 1 (CNR). On 15 February 2013, bathymetric data were collected between 25 and 700 m water depths using the multibeam echosounder Kongsberg Simrad EM710 working at a frequency of 70-100 kHz. On 1 January 2016, bathymetric data were collected in the depth range 10-700 m using Teledyne Reson 7125 and 7160 multibeam echosounders working at a frequency of 400 kHz and 45 Hz in shallow-(<100 mwd) and deep-water (>100 mwd), respectively. In these surveys, data were DGPS-positioned and processed with hydrographic software (Caris Hips and Sips Professional), using daily sound speed profiles and patch test of transducers in the survey zone. Furthermore, a hull-mounted sound speed sensor was used to update in real-time the sound velocity values close to the flat face of the multibeam transducer. Tidal corrections were performed using data from the nearby tide gauge stations (www.mareografico.it). Random spikes and organized noise (multiple) produced by local very steep slopes were removed by applying geometrical and statistical filters. Cleaned multibeam data were gridded using a weighted averaging algorithm to produce the DEMs at variable resolution, from 1 m in shallow-water (<100 mwd) to 5 m at greater depths. DEMs used for bathymetric comparison were gridded with cell-size of 3 m and limited at depths <500 m, because at greater depths acoustic noise largely increased, making unreliable the computed changes. A vertical accuracy range of ±0.5 m was estimated for the difference map by comparing the difference in depth of stable benchmarks between pairs of successive bathymetric sets. Small depth changes below this error range in the difference map were not considered in volume computations.
For the subaerial slope, topographic changes were estimated by comparing two DEM reconstructed starting from a LiDAR airborne survey (4-8 May 2012) and PLÉIADES-1 tri-stereo satellite imagery (26 May 2017). The LiDAR-DEM was obtained through the processing of the 3D point cloud that was acquired by using the Leica ADS80 sensor, which has instrumental vertical and horizontal accuracy of 0.10-0.20 and 0.25 m, respectively [46,47]. The acquired point cloud has a mean point density of 8 pt/m 2 . The 2017 DEM derived from the tri-stereo optical imagery acquired by the PLÉIADES-1 satellites, provided by Airbus with a nominal xy resolution of 1 m × 1 m. The PLÉIADES-1A (PHR1A) and PLÉIADES-1B (PHR1B) satellites can sense different synchronous images of the same area, with an angle variable between 6 and 28 • . The tri-stereo approach consists of the acquisition of three nearly simultaneously images (one backward looking, one forward looking and a third near-nadir image) allowing the DEM reconstruction through the photogrammetric processing of the three stereo images [48]. Although the acquisitions used in this study are 100% cloud-free, the presence of the gas plume from the crater terrace prevented the reconstruction of the topography in this area. 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). A block adjustment including all the satellite scenes was performed. The block adjustment was validated when the following accuracy was achieved: (i) an image's pixel xy bias smaller than 0.3 pixels, (ii) an image's pixel xy standard deviation smaller than 0.3 pixels and (iii) an image's pixel xy maximum residuals smaller than 2 pixels (for details see [16,47]). The z standard deviation was about 1.5 m (47).

Change Detection with SAR Amplitude Images
Data from SAR sensors have often been used to maps areas that have been affected by lithological and morphological changes, i.e., to identify areas that have been impacted by eruptive and post-eruptive (landslides or floods) phenomena [49][50][51]. The SAR backscattering is determined by several factors, as the local morphology and the surface microrelief related to the grain-size, and the dielectric constant of the material at the surface [52,53]. To define which roughness scale affects the backscatter properties, the Rayleigh criterion was applied. In this way, the root mean squared height (hrms) variation on horizontal surfaces has been evaluated following [52] and, for the COSMO-SkyMed SAR (CSK-SAR) images used in this study, the hrms is approximately 4 mm. In an area with volcaniclastic sedimentation, this indicates that it is possible to identify the variation of the sedimented material between fine-grained (i.e., ash-dominated deposits) and coarse-grained (i.e., blocks and bombs dominated deposits, such as grain flows or 'a'ā lava surfaces). Local morphology produces irregularities having wavelengths at least twice as large as the satellite resolution cell due to changes in the local incidence angles, and it is not directly possible to separate the local morphology effects from the grain-size influence. Therefore, the term "roughness" is used to represent a combination of both factors [52].
To detect and interpret changes in land cover in connection with the SdF slope, RGB color composites are used [44,52]. The composite is created using these combinations: The ratio is used because it depends on the relative average radar backscattering changes between two images, and therefore it does not depend on the pixel's intensity level [54].
A dataset comprising 7 COSMO-SkyMed SAR (CSK-SAR) images, acquired in descending orbit between 8 May 2014 and 12 August 2014, was used for this study ( Figure 3 and Table 1). The products were coregistered, using the offset refinement based on the shuttle radar topography mission (SRTM) digital elevation model (DEM) at 1 arc second forming one unique stack, which is cropped around the target area. The cropped stack is geocoded by correcting SAR geometric distortions using SRTM DEM, producing SAR orthorectified map-projected images. Backscattered intensity of each image is transformed in amplitude image and then decibel scaled, converting the data into a virtual band with the expression 10*log10 (amplitude). Finally, the quality of the images was enhanced using a multitemporal speckle filter that reduces the "salt and pepper-like" texturing (speckle) of the CSK-SAR data.

Morphological Evolution of the Submarine SdF between 2013 and 2016
The DEMs used to constrain seafloor changes associated with the August-November 2014 eruption were collected in February 2013 and January 2016, respectively ( Figure 3). Despite the time gap before and after the eruption, our 20 years' experience on submarine monitoring of the SdF has evidenced that the main morphological changes occur during eruptive crisis or increased period of Strombolian activity, while the "ordinary" evolution of the slope is limited to small readjustments of the shallowest part of the SdF [20,[24][25][26].
The pre-2014 eruption submarine SdF slope can be morphologically divided in the SW and NE sectors, whose boundary roughly corresponds to the SW limit of the 2002 landslide scars (dashed red lines in Figure 4a). Considering also previous surveys, the SW sector in the time frame 2002-2017 was affected by minor morphological changes, whereas significant changes occurred in the NE sector both due to the morphological readjustment of the 2002 landslide scar, and to the emplacement of the 2007 lava delta (dotted blue lines in Figure 4a). These events left a depressed area in the southern part of the NE sector of the SdF slope, hereafter defined as the central part of the SdF.
The difference between pre-and post-2014 bathymetries shows that the main morphological changes on the NE sector of the SdF occur within the first 250 mwd, even if they locally extend down to 500 mwd (seaward limit of the difference map, Figure 4b). Seafloor accretion largely overwhelms erosion, accounting for a total estimated volume of +1.75 × 10 6 m 3 and −3.5 × 10 5 m 3 , respectively ( enhanced using a multitemporal speckle filter that reduces the "salt and pepper-like" texturing (speckle) of the CSK-SAR data.

Morphological Evolution of the Submarine SdF Between 2013 and 2016
The DEMs used to constrain seafloor changes associated with the August-November 2014 eruption were collected in February 2013 and January 2016, respectively ( Figure 3). Despite the time gap before and after the eruption, our 20 years' experience on submarine monitoring of the SdF has evidenced that the main morphological changes occur during eruptive crisis or increased period of Strombolian activity, while the "ordinary" evolution of the slope is limited to small readjustments of the shallowest part of the SdF [20,[24][25][26].
The pre-2014 eruption submarine SdF slope can be morphologically divided in the SW and NE sectors, whose boundary roughly corresponds to the SW limit of the 2002 landslide scars (dashed red lines in Figure 4a). Considering also previous surveys, the SW sector in the time frame 2002-2017 was affected by minor morphological changes, whereas significant changes occurred in the NE sector both due to the morphological readjustment of the 2002 landslide scar, and to the emplacement of the 2007 lava delta (dotted blue lines in Figure 4a). These events left a depressed area in the southern part of the NE sector of the SdF slope, hereafter defined as the central part of the SdF.
The difference between pre-and post-2014 bathymetries shows that the main morphological changes on the NE sector of the SdF occur within the first 250 mwd, even if they locally extend down to 500 mwd (seaward limit of the difference map, Figure 4b). Seafloor accretion largely overwhelms erosion, accounting for a total estimated volume of +1.75 × 10 6 m 3 and -3.5 × 10 5 m 3 , respectively ( Table 2).  The pre-/post-eruption comparison for the subaerial part of the SdF shows only accretion that largely overwhelms erosion (limited to a few sectors here indicated by arrows, see also (adapted from [47]). The dashed red line delimits the area with main morphological changes associated with the 2014 lava flows entering the sea. Table 2. Volumetric estimation related to (a) submarine morphological changes (depth range 10-500 mwd) occurred between 2013 and 2106, (b) accretion and erosion of the entire SdF slope associated with the 2014 effusive eruption, (c) volcaniclastic material emplaced in the central part of the SdF, mainly before the 2014 eruption (see text for detail). The red and light-blue colors are referred to the subaerial and submarine flanks, respectively; the black color is used for the total volume. * Volume inferred in the gap area between the coastward limit of the bathymetric survey (around 10 mwd) and the coastline. Note that the error on z regard to the submarine and subaerial DEMs is on the order of 0.5 m and 1.5 m, respectively. The symbol ≈ is used for approximately.  (Figures 6 and 7) and a previously existing basement high (BH) are also indicated, along with the trace of the bathymetric profiles (P1 to P5, blue dotted lines) shown in Figure 8.
A smaller erosive area (E2 in Figure 5), with thickness lower than 5 m, was recognizable just to the north of E1. These erosive areas correspond to the formation of channelized features on the 2016 DEM, bounded by steep-sided ridges (R1 and R2 in Figures 5  and 6b).  (Figures 6 and 7) and a previously existing basement high (BH) are also indicated, along with the trace of the bathymetric profiles (P1 to P5, blue dotted lines) shown in Figure 8.
A smaller erosive area (E2 in Figure 5), with thickness lower than 5 m, was recognizable just to the north of E1. These erosive areas correspond to the formation of channelized features on the 2016 DEM, bounded by steep-sided ridges (R1 and R2 in Figures 5 and 6b). Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 23 The E1 channelized area ( Figure 5 and Ch in 6b) roughly matches the location of a large morphological high (MH in Figure 6a) on the 2013 DEM; the smaller E2 erosive area corresponds, instead, to an area previously characterized by small-scale coalescing fan-shaped features (Figure 6a).
Seafloor accretion, as resulting from the comparison between pre-and post-2014 eruption, can be divided in four main areas (A1-A4 in Figure 5): A1-A3 were mainly located in the first 200 mwd, while A4 was located at depths greater than 200 mwd.  Figure 5). In (b), the subaerial difference map between pre-and post-2014 eruption is also draped over the coastal topography (for the color scale refer to The E1 channelized area ( Figure 5 and Ch in Figure 6b) roughly matches the location of a large morphological high (MH in Figure 6a) on the 2013 DEM; the smaller E2 erosive area corresponds, instead, to an area previously characterized by small-scale coalescing fan-shaped features (Figure 6a).
Seafloor accretion, as resulting from the comparison between pre-and post-2014 eruption, can be divided in four main areas (A1-A4 in Figure 5): A1-A3 were mainly located in the first 200 mwd, while A4 was located at depths greater than 200 mwd.
The A4 area covers a surface of 1.8 × 10 5 m 2 down to 500 mwd (seaward limit of the difference map) for a total estimated volume of +6 × 10 5 m 3 ; it morphologically corresponds to the development of a large fan-shaped feature, as visible on the 2016 DEM (Figure 6b), with its apex located within the E1 channelized area (see P3 in Figure 8).
The A1 area covers a surface of 1.25 × 10 5 m 2 , with a thickness of few tens of meters ( Figure 5 and P1 in Figure 8) for a total estimated volume of +8.35 × 10 5 m 3 ; it morphologically corresponds to the development of small coalescing fan-shaped features alternated with downslope elongated ridges (R3, R4 and R5 in Figures 5 and 7b), characterized by a rough morphology and steep lateral gradients (over 60 • , P1 in Figure 8)  The A4 area covers a surface of 1.8 × 10 5 m 2 down to 500 mwd (seaward limit of the difference map) for a total estimated volume of +6 × 10 5 m 3 ; it morphologically corresponds to the development of a large fan-shaped feature, as visible on the 2016 DEM (Figure 6b), with its apex located within the E1 channelized area (see P3 in Figure 8).
The A1 area covers a surface of 1.25 × 10 5 m 2 , with a thickness of few tens of meters ( Figure 5 and P1 in Figure 8) for a total estimated volume of +8.35 × 10 5 m 3 ; it morphologically corresponds to the development of small coalescing fan-shaped features alternated with downslope elongated ridges (R3, R4 and R5 in Figures 5 and 7b), characterized by a rough morphology and steep lateral gradients (over 60°, P1 in Figure 8 Figure  8).
The A2 accretion area also corresponded to a well-defined ridge (R2 in Figure 6b) on the 2016 DEM, having steep flanks (slope gradients up to 60°, P2 in Figure 8); this ridge extended over an area of 2.6 × 10 4 m 2 down to 200 mwd, with a maximum thickness of 20 m (Figures 5 and 6b) for an estimated volume of +1.7 × 10 5 m 3 . R2, together with smaller ridge (R1 in Figures 5 and 6b), bound the E1 erosive area. R2 was not present on the 2013 DEM, roughly matching the location of a previous channelized feature (Ch in Figure 6a), while R1 was already present before the eruption as part of a morphological high (MH in Figure 6a) and did not correspond to an accreted area in the difference map of Figure 5.
Finally, A3 was located in the NE-most part of the study area, covering a surface of 30.000 m 2 down to 200 mwd, with an average seafloor accretion of approximately 10 m ( Figure 5 and P2 in Figure 8) for an estimated volume of about +1.5 × 10 5 m 3 . Seafloor accretion here morphologically corresponds to the development of irregular and smooth  Figure 8).
The A2 accretion area also corresponded to a well-defined ridge (R2 in Figure 6b) on the 2016 DEM, having steep flanks (slope gradients up to 60 • , P2 in Figure 8); this ridge extended over an area of 2.6 × 10 4 m 2 down to 200 mwd, with a maximum thickness of 20 m (Figures 5 and 6b) for an estimated volume of +1.7 × 10 5 m 3 . R2, together with smaller ridge (R1 in Figures 5 and 6b), bound the E1 erosive area. R2 was not present on the 2013 DEM, roughly matching the location of a previous channelized feature (Ch in Figure 6a), while R1 was already present before the eruption as part of a morphological high (MH in Figure 6a) and did not correspond to an accreted area in the difference map of Figure 5.
Finally, A3 was located in the NE-most part of the study area, covering a surface of 30.000 m 2 down to 200 mwd, with an average seafloor accretion of approximately 10 m ( Figure 5 and P2 in Figure 8) for an estimated volume of about +1.5 × 10 5 m 3 . Seafloor  Figure 8). These fan-shaped features were emplaced above a previous channelized feature recognizable on the 2013 DEM (Ch in Figure 6a), carving the NW flank of a basement high (BH in Figures 5 and 6) recognizable along the northern shoulder of the SdF collapse scar.
In order to compare submarine and subaerial morphological changes associated with the 2014 eruption, a topographic profile (P5 in Figure 8) crossing the entire NE sector of the SdF (location in Figure 5)  fan-shaped features on the 2016 DEM, developed above and below of a marked slope break (Figure 6b and P4 in Figure 8). These fan-shaped features were emplaced above a previous channelized feature recognizable on the 2013 DEM (Ch in Figure 6a), carving the NW flank of a basement high (BH in Figures 5 and 6) recognizable along the northern shoulder of the SdF collapse scar. In order to compare submarine and subaerial morphological changes associated with the 2014 eruption, a topographic profile (P5 in Figure 8) crossing the entire NE sector of the SdF (location in Figure 5) is reported, showing the slope accretion occurred between 2012 and 2017 surveys.  Figures 5, 6 and 7) showing seafloor accretion and erosion (red and light blue polygons, respectively) between pre-and post-2014 eruption topography. Acronyms as in the previous figures. Below the profiles, the graph of slope gradients is referred to the post-eruption topography. The location of the apex fan recognizable on the hillshade bathymetry of Figure 6b is reported on profile P3, indicating that the scar is partially filled by landslide deposits.

CSK-SAR Data Collected between May and August 2014
CSK-SAR images, collected between 8 May 2014 and 24 May 2014, show a decrease in the amplitude ratio (i.e., loss in SAR backscattering between the two consecutive images) in the central portion of the SdF (Figure 9a), testifying the erosion of volcaniclastic  Figures 5-7) showing seafloor accretion and erosion (red and light blue polygons, respectively) between pre-and post-2014 eruption topography. Acronyms as in the previous figures. Below the profiles, the graph of slope gradients is referred to the post-eruption topography. The location of the apex fan recognizable on the hillshade bathymetry of Figure 6b is reported on profile P3, indicating that the scar is partially filled by landslide deposits.

CSK-SAR Data Collected between May and August 2014
CSK-SAR images, collected between 8 May 2014 and 24 May 2014, show a decrease in the amplitude ratio (i.e., loss in SAR backscattering between the two consecutive images) in the central portion of the SdF (Figure 9a), testifying the erosion of volcaniclastic deposits. Afterwards, the progressive increase in the explosion rate (16), and the emplacement of lava overflows and rock-falls evolving in gravel flows, produced the volcaniclastic sedimentation in the central part of the SdF (Figure 9b,c). deposits. Afterwards, the progressive increase in the explosion rate (16), and the emplacement of lava overflows and rock-falls evolving in gravel flows, produced the volcaniclastic sedimentation in the central part of the SdF (Figure 9b,c). This sedimentation (and related superficial roughness) clearly increased in the images collected before the onset of the eruption (11 and 27 July 2014; Figure 9d,e). The This sedimentation (and related superficial roughness) clearly increased in the images collected before the onset of the eruption (11 and 27 July 2014; Figure 9d,e). The events of 6 August 2014, with a series of overflows associated with frequent landslides, significantly increased the roughness in the central part of the SdF (Figure 9f). In the same image, the variation in amplitude on the coastline is due to the entrance into the sea of the lava flows. The lava effusion that began on 7 August 2014 instead produced no remarkable changes in SAR backscattered amplitude, probably due to the similar texture between the pre-effusive and the early-lava flows surfaces. Instead, it is possible to note how the area of proximal volcaniclastic accumulation (i.e., NEC talus) is characterized by accumulations during the whole period considered, except for the last image where an erosive process is highlighted (reduction in the backscattered amplitude in Figure 9f), consistent with the landslide occurred in the early stages of the eruption.

Seafloor Accretion Associated with the 2014 Eruption
The seafloor accretion observed in areas A2 and A3 ( Figure 5) fits well the entry points of the main lava flows reaching the sea on the NE-most sector of the SdF during the 2014 eruption (outlined by the dashed black lines in Figure 1b), emplaced above the remnants of the 2007 lava delta (Figure 6a; [16]). Once entered the sea, these lava flows built 40-60 m wide and some hundreds m long lobes (for instance R2 in Figure 6b), similarly to what was observed during the emplacement of the 2007 lava flows at Stromboli [25]. However, some morphological differences can be noted between the lobes emplaced in A2 and A3: while the former one corresponds to a well-defined and steep-sided ridge (R2 in Figure 6b), a more irregular, fan-shaped feature is present in the area A3 (see P2 in Figure 8 for comparison). This difference in morphology could be ascribed to the different amount of coherent lava flows with respect to chaotic breccias, using as reference the scuba observations made for the 2007 lava delta, where this kind of products were recognized [25]. Similarly, ancient analogous of 'a'ā lava-fed deltas outcropping in Antarctica were found to be made up by the coalescence of "hyaloclastic lobes", with variable percentages of coherent lavas and chaotic breccias [57].
The higher slope gradients (locally over 60 • , P2 in Figure 8) observed along the flanks of R2 can be related to the high percentage of coherent lava flows forming this ridge. Lava flows feeding R2 were, in fact, active along the slope during the first days of the eruption, when the effusion rates were high (up to 20 m 3 /s according to the model of [28]). The high effusion rates can also justify the thicker accretion recognizable on the difference map in the coastal sector of the subaerial slope (red and violet colors in Figure 5 just above A2). Such conditions would have favored the penetration of lava flows into the sea as coherent bodies, as suggested by several laboratory and field studies [55,56]. On the meantime, the lower thickness observed for the 2014 lava flows in the coastal sector facing A3 (yellow colors in Figure 5) would justify the development of the irregular fan-shaped feature in this area, due to the prevalence of volcaniclastic breccias on the slope. The volcaniclastic contribution appeared predominant also moving downslope, where the emplacement of a smoothed fan is observed in A3, down to almost 200 mwd (Figure 6b). This composite fan likely formed during the declining stage of the 2014 eruption, when effusion rates reached a steady value of 0.4 m 3 /s, and the entrance of well-fed, coherent lava flows into the sea was no longer observed. Moreover, the steep slope break at around 100 mwd could have favored the brecciation of the less sustained submarine lava flows that reached this area, as testified by the thinned seafloor accretion over the slope break (see P4 in Figure 8).
The overall comparison between the pre-and post-bathymetries thus suggests a key role played by the paleomorphology of the SdF slope in controlling the location of the different lava flows and volcaniclastic lobes, which were emplaced within previous channelized features (Figure 6; P2 in Figure 8). Beyond the SW limit of the main 2014 lava flows (as reported by [16], indicated with the dashed line in Figure 4b), a series of steep-sided ridges (R3, R4 and R5 in Figures 5 and 7b) form the accretion area A1. From the comparison with onshore data (Figure 10) they can be interpreted as the submarine extension of lava flows that were mainly emitted during the first days of the 2014 eruption, when effusion rates were very high, and multiple flows reached this coastal sector (Figure 2b,d).
the comparison with onshore data (Figure 10) they can be interpreted as the submarine extension of lava flows that were mainly emitted during the first days of the 2014 eruption, when effusion rates were very high, and multiple flows reached this coastal sector (Figure 2b,d). Besides these ridges, the seafloor accretion in A1 is also made up of coalescing and overlapping small fan-shaped features (Figure 7b) that can be interpreted as the result of erosive-depositional processes associated with gravel flows entering the sea, fed by small-scale instabilities on the steep (>30°) subaerial and submarine slope of the SdF. These small fan-shaped features are a common, but ephemeral feature, often observed in many bathymetric surveys performed since the 2002. They are associated with the morphological adjustment of the SdF slope that commonly occurs due to severe winter storms affecting the NW flank of the island. A good match is observed between subaerial and submarine accretion in the central sector of the SdF (Figures 5 and 10a), due to the progressive infilling of this depressed area (P1 and P5 in Figure 8). This infilling can be mainly related to the abundant volcaniclastic material derived from (a) overflow-induced landslides, (b) the NEC-talus collapse (Figure 2a)  Besides these ridges, the seafloor accretion in A1 is also made up of coalescing and overlapping small fan-shaped features (Figure 7b) that can be interpreted as the result of erosive-depositional processes associated with gravel flows entering the sea, fed by small-scale instabilities on the steep (>30 • ) subaerial and submarine slope of the SdF. These small fan-shaped features are a common, but ephemeral feature, often observed in many bathymetric surveys performed since the 2002. They are associated with the morphological adjustment of the SdF slope that commonly occurs due to severe winter storms affecting the NW flank of the island. A good match is observed between subaerial and submarine accretion in the central sector of the SdF (Figures 5 and 10a), due to the progressive infilling of this depressed area (P1 and P5 in Figure 8). This infilling can be mainly related to the abundant volcaniclastic material derived from (a) overflow-induced landslides, (b) the NEC-talus collapse (Figure 2a) and (c) erosion of the upper part of the SdF (arrow in Figure 4b), occurred before, during and after the 2014 eruption.
CKS-SAR data show a relevant accumulation of volcaniclastic material within the depressed area encompassed between the 2007 lava delta and the SW limit of the 2002 landslide scar some months before the starting on the 2014 eruption ( Figure 9). Differently, relatively mild slope dynamics, mainly characterized by weak erosive-depositional processes, occurred in the SdF in the period late 2014-early 2017 (Figure 3, [34]), when the post-2014 submarine/subaerial surveys were realized. A possible limit between the volcanic material emplaced in the submarine slope before and during the 2014 effusive eruption has been inferred by through the analysis of the bathymetric cross-sections by reconstructing a hypothetical surface that approximately matches the base of the ridges before mentioned (R2 to R5, dashed black line indicated as "inferred 2014 base" in P1 and P2 of Figure 9).

Seafloor Erosion Associated with the 2014 Eruption
During the phase preceding the onset of the eruption, there were frequent gravel flows both from the edge of the crater terrace and along the SdF, mainly triggered by overflows. The volume of these gravel flows never reached 10 6 m 3 , remaining in the order of magnitude 10 4 -10 5 m 3 of remobilized material [16]. In the submarine slope, two failure events (E1 and E2 in Figures 5 and 10a) occurred where the main 2014 lava flows entered the sea, indicating a strong link between slope failures and 2014 eruptive dynamics. It is noteworthy that E1 is the larger submarine instability event recorded through the bathymetric monitoring along the SdF since the 2007 lava delta emplacement.
The computed volume of this event is around −3.5 × 10 5 m 3 ( Table 2), but this figure is very likely underestimated because: (a) bathymetric data are limited to water depth of 10 m, and the difference map indicates that the scar likely cuts back coastward and (b) the scar is partially filled both by landslides deposits in its distal part (see Figure 6b and P3 in Figure 8) and by small-scale fan-shaped features related to successive erosive-depositional process in its proximal part (Figure 6b). A more realistic estimate of the mobilized volume can be obtained by taking into account the A4 seafloor accretion area (Figures 4b, 5 and 10a), interpreted as the landslide deposits recognizable at the base of the scar and accounting for approximately +6 × 10 5 m 3 ( Table 2). The low runout of the mobilized material is likely due to a rapid dissipation of pore pressures during the slide, similarly to what was reported for small submarine slope failures affecting the 2018 'a'ā lava delta off Hawaii [38]. Again, this figure is likely underestimated because (a) the mapping of the landslide deposits is limited down to 500 mwd ( Figure 5), and (b) the distal part of the moving landslide could be evolved downslope as sedimentary gravity flow, due to the steep gradients characterizing the SdF slope.
Regardless, the estimated volume of the E1 slope failure is two order magnitude larger than any instability processes occurred since the emplacement of the 2007 lava delta, and it also affected the more stable part of the delta according to the previous bathymetric surveys [26]. Such findings support our inference that triggering mechanisms for this larger slope failure are likely different from those common involved in the "ordinary" readjustment of the submarine slope, such as the loading of the 2014 lava flows on the submarine slope or the seismic activity associated with the 2014 eruptive event. However, the combined action of storm-waves cannot be excluded among the triggering mechanisms due to the shallow depths.

Volumetric Considerations on the 2014 Eruption and SdF Infilling
The estimation of the submarine volume emplaced during the 2014 eruption is difficult, because the subaerial and submarine surveys were not performed just after and before the eruption, as for instance realized during the 2007 eruption [15,25]. Data were collected one year before and one after the eruption for the submarine part and a slightly longer time frame (2012-2017) for the subaerial part. Despite such limits, some volumetric constraints on the overall material emplaced in the SdF during the 2014 eruption can be derived.
According to [16] the main lava flows emitted during the 2014 were emplaced in the NE-most sector of the SdF, encompassed by the dashed black line in Figures 1 and 9a, with a volume of approximately +2.7 × 10 6 m 3 ( Table 2). Seafloor accretion in A2 and A3 (Figures 4b and 10) can be confidently associated with these lava flows, accounting for approximately +3.3 × 10 5 m 3 (Table 2). However, taking into account a southwestern, further extension of lava flows in the central part of the SdF on the base of submarine evidence (from R3 to R5 in Figures 5, 7, 8 and 10) we proposed in this work that further +1.5 × 10 5 m 3 can be added as products of the 7 August-13 November 2014 eruption (Table 2). Moreover, this figure is surely underestimated because the bathymetric com-parison pre-/post-eruption is limited at depths greater than 10 m. In order to compute the volume not considered, an average thickness of 5 and 10 m has been assumed for the SW and NE sectors, respectively, based on the trend derived from the difference map and integrated over the coast. The estimated volume is +2 × 10 5 m 3 of material that added to the previous two figures accounts for a total estimated volume of ≈6.8 × 10 5 m 3 in the submarine part of the SdF (Table 2).
A similar approach has been also applied for the central part of the subaerial SdF facing A1, where the volume possibly related to the 2014 eruption was distinguished with respect to the underlying volcaniclastic material by reconstructing a surface interpolated through the base of the subaerial ridges facing the submarine ones (P5 in Figure 8). This volume was estimated to about +3.5 × 10 5 m 3 , which summed to the previous +2.7 × 10 6 m 3 accounts for a total volume of +3.05 × 10 6 m 3 emplaced on the subaerial slope ( Table 2). It is noteworthy that the volume emplaced in the surveyed submarine area was markedly lower (about 5 times) than the volume emplaced on the subaerial SdF ( Figure 11 and Table 2). This is different from what observed during the previous 2007 eruption, where the volume of the erupted material estimated for the submarine part (where a large lava delta in the order of 7 × 10 6 m 3 of volume was formed, see [25]) was approximately three times larger than the subaerial one ( Figure 11 and Table 2).
According to [16] the main lava flows emitted during the 2014 were emplaced in the NE-most sector of the SdF, encompassed by the dashed black line in Figures 1 and 9a, with a volume of approximately +2.7 × 10 6 m 3 ( Table 2). Seafloor accretion in A2 and A3 (Figures 4b and 10) can be confidently associated with these lava flows, accounting for approximately +3.3 × 10 5 m 3 (Table 2). However, taking into account a southwestern, further extension of lava flows in the central part of the SdF on the base of submarine evidence (from R3 to R5 in Figures 5, 7, 8 and 10) we proposed in this work that further +1.5 × 10 5 m 3 can be added as products of the 7 August-13 November 2014 eruption ( Table 2). Moreover, this figure is surely underestimated because the bathymetric comparison pre-/post-eruption is limited at depths greater than 10 m. In order to compute the volume not considered, an average thickness of 5 and 10 m has been assumed for the SW and NE sectors, respectively, based on the trend derived from the difference map and integrated over the coast. The estimated volume is +2 × 10 5 m 3 of material that added to the previous two figures accounts for a total estimated volume of ≈6.8 × 10 5 m 3 in the submarine part of the SdF (Table 2).
A similar approach has been also applied for the central part of the subaerial SdF facing A1, where the volume possibly related to the 2014 eruption was distinguished with respect to the underlying volcaniclastic material by reconstructing a surface interpolated through the base of the subaerial ridges facing the submarine ones (P5 in Figure  8). This volume was estimated to about +3.5 × 10 5 m 3 , which summed to the previous +2.7 × 10 6 m 3 accounts for a total volume of +3.05 × 10 6 m 3 emplaced on the subaerial slope ( Table 2). It is noteworthy that the volume emplaced in the surveyed submarine area was markedly lower (about 5 times) than the volume emplaced on the subaerial SdF ( Figure  11 and Table 2). This is different from what observed during the previous 2007 eruption, where the volume of the erupted material estimated for the submarine part (where a large lava delta in the order of 7 × 10 6 m 3 of volume was formed, see [25]) was approximately three times larger than the subaerial one ( Figure 11 and Table 2).  Table 2) and percentages between volcanic material accreted on the subaerial (light brown) and submarine (light blue) slope of the SdF in correspondence of different eruptive phases.  Table 2) and percentages between volcanic material accreted on the subaerial (light brown) and submarine (light blue) slope of the SdF in correspondence of different eruptive phases. This discrepancy can be mainly associated with the different elevation of the main vents feeding lava flows during the 2007 eruption (around 400 m a.s.l.) and the 2014 eruption (around 650 m a.s.l). Considering that, during the 2014 eruption, effusion rates drastically fell from 20 to 0.5 m 3 /s just two days after the eruption [28], lava flows erupted by a high-elevation vent during the declining stage of the eruption could have remained mostly confined to the median part of the SdF slope, without reaching the sea. Another interesting finding is that the total volume estimated for the 2014 eruption through the above-described computations (approximately +3.7 × 10 6 m 3 , Table 2) is markedly lower than the total bulk volume of +7.4 × 10 6 m 3 estimated using data derived by the new satellite Technology Experiment Carrier-1 (TET-1) [37], whereas it is more comparable to the total bulk volume of +5.5 × 10 6 m 3 obtained by thermal images from satellites using the moderate resolution imaging spectroradiometer (MODIS) sensor [28]. This large difference can be mostly explained by the different methods used for constraining the total volume, as in the case of MODIS approach, where the error is usually estimated at 50% [58].
Nevertheless, our approach relies on topographic differences realized not contextually at the start and end of eruption, so hindering a clear discrimination between material the emplaced during the eruption with respect to the volcaniclastic material emplaced (mostly) before and after the eruption, especially in the central part of the SdF. This area, in fact, is a morphological low with respect to the surrounding NE and SW sectors, representing a main pathway for the funneling to the sea of volcaniclastic material produced by persistent strombolian activity and erosive processes acting in the upper SdF slope (Figures 4b and 10), and this makes difficult to discriminate between such deposits and the previous 2014 overflows. As a whole, our results show a good matching between the subaerial and submarine slope accretion in the central part of the SdF (Figures 4b, 5, 9 and 10a). As regards the subaerial slope, this infilling is mostly limited between 500 m and the coastline, resulting in a volume of approximately +1.35 × 10 6 m 3 ( Table 2). For the submarine slope, the infilling is mostly limited down to −250 m, with an average accretion of 10 m for an estimated volume of approximately +7 × 10 5 m 3 , 4.9 × 10 5 m 3 of which were directly measured through the difference map (limited to depths greater than 10-15 m). The cumulative value between submarine and subaerial slope accounts for approximately 1.84 × 10 6 m 3 accumulated during few years, representing not only eruptive dynamics but also gravity instability processes reworking the SdF slope. It is noteworthy that this infilling is also responsible for the development of a steep submarine volcaniclastic apron in the proximal sector of the submarine SdF that was one of the main factors controlling the geometry of the 2002 tsunamigenic landslide [20].

Conclusions
The integration of repeated bathymetric, lidar and photogrammetric surveys along with SAR amplitude images was a fundament tool to understand the morphological evolution of the SdF slope induced by the 2014 eruption and its pre-eruptive stage, which was characterized by two months of intense Strombolian activity. Difference maps showed that slope accretion largely overwhelmed erosion along most of the SdF slope both before and during the 2014 eruption. The total volume of erupted material during the 2014 eruption was approximately estimated in +3.7 × 10 6 m 3 , ≈80% of which emplaced in the subaerial slope. This figure is different from what depicted for the previous 2007 eruption, where a large submarine lava delta formed, representing ≈75% of the total volume. This difference can be explained by the location of the main vents feeding the 2007 (around 400 m a.s.l.) and 2014 eruption (around 650 m). Before the 2014 eruption, the total volume estimated is about +1.84 × 10 6 m 3 , the 60% of which emplaced in the subaerial slope between 500 m and the coastline. The remaining 40% was emplaced in the submarine slope down to 250 mwd, promoting the formation of a steep volcaniclastic apron, a feature that increase the tsunamigenic potential of the SdF, as it was considered to have controlled the geometry of the submarine tsunamigenic landslide occurred in 2002.
The 2014 eruption also promoted significant slope instabilities along the SdF slope. In the subaerial slope, landslides mobilized volumes in the order of −10 4 /10 5 m 3 , mostly affecting the upper part of the SdF during the early stages of the eruption. In the submarine slope, a main landslide scar was detected in the NE part of the SdF within the first 200 mwd, just off the entry points of the main 2014 lava flows. Landslide deposits are recognizable at the base of the scar, extending down to 500 mwd, corresponding to the limit of the bathymetric survey. This landslide mobilized a volume of at least −6 × 10 5 m 3 , representing the largest slope failure affecting the submarine part of the 2007 lava delta since its emplacement.
DEM analysis also showed a good correlation between the subaerial and submarine slope and a variability in submarine morphologies associated with the entrance into the sea of lava flows (i.e., steep-sided ridges and fan-shaped features), reflecting a different ratio between coherent lava flows and volcaniclastic breccias that, in turn, was likely controlled by changes in effusion rates and paleo-topography. This observation is very important for the interpretation of seafloor volcanic morphologies around insular and coastal volcanoes, mainly in relation to the paucity of marine studies monitoring the behavior of lava flows penetrating into the sea.
The presented results highlight the importance of integrated submarine and subaerial studies to monitor active volcanoes, providing a comprehensive view of the main processes (constructive vs destructive) associated with eruptive dynamics, with significant implications also for the related geohazard assessment.