Integration of Remote Sensing and Offshore Geophysical Data for Monitoring the Short-Term Morphological Evolution of an Active Volcanic Flank: A Case Study from Stromboli Island

: The Sciara del Fuoco (SdF) collapse scar at Stromboli is an active volcanic area affected by rapid morphological changes due to explosive/effusive eruptions and mass-wasting processes. The aim of this paper is to demonstrate the importance of an integrated analysis of multi-temporal remote sensing (photogrammetry, COSMO-SkyMed Synthetic Aperture Radar amplitude image) and marine geophysical data (multibeam and side scan sonar data) to characterize the main morphological, textural, and volumetric changes that occurred along the SdF slope in the 2020–2021 period. The analysis showed the marked erosive potential of the 19 May 2021 pyroclastic density current generated by a crater rim collapse, which mobilized a minimum volume of 44,000 m 3 in the upper Sciara del Fuoco slope and eroded 350,000–400,000 m 3 of material just considering the shallow-water setting. The analysis allowed us also to constrain the main factors controlling the emplacement of different lava ﬂows and overﬂows during the monitored period. Despite the morphological continuity between the subaerial and submarine slope, textural variations in the SdF primarily depend on different processes and characteristics of the subaerial slope, the coastal area, the nearshore, and “deeper” marine areas.


Introduction
Volcanoes are highly dynamic environments, whose evolution is affected by the alternance between constructive and destructive processes acting at different spatial and temporal scales [1][2][3].These processes continuously reshape volcanic flanks, making their monitoring and linking with eruptive and/or erosive phenomena through time fundamental for understanding both the morphological evolution of volcanic edifices and the genesis of specific morphological features, which can be used as "geomorphological markers" in other volcanic settings.Furthermore, morphological changes are often related to sudden and catastrophic phenomena (pyroclastic flows or landslides), which represent a main geohazard for the surrounding communities.This issue is particularly relevant in the case of insular volcanoes, where these phenomena can directly impact subaerial/submarine infrastructures or indirectly generate tsunami waves, as recently observed for the 2002 landslide at Stromboli [4,5] or the 2018 Anak Krakatau collapse [6].
Notwithstanding the direct link between the morphogenesis of subaerial and submarine volcanic flanks, integrated studies dealing with land-sea correlations are rare and have been mainly limited to large sector collapses involving the entire volcanic edifice [7][8][9].In the last few decades, some attempts have been carried out in this framework where shallow-water morphologies (e.g., insular shelf) were correlated with the facing island geology, providing useful insights on the island's vertical movements, the identification of paleo-centres or the mechanisms occurring when lava flows enter the sea [10][11][12][13].Even fewer studies have tried to integrate subaerial and submarine morphological changes during eruptive events affecting volcanic flanks [14][15][16][17][18].At Stromboli, the latter studies were mainly focused on the Sciara del Fuoco depression (Figure 1), representing the last of repeated lateral collapses that affected the NW flank in the last 13 ka [19,20].
temporal scales [1][2][3].These processes continuously reshape volcanic flanks, making their monitoring and linking with eruptive and/or erosive phenomena through time fundamental for understanding both the morphological evolution of volcanic edifices and the genesis of specific morphological features, which can be used as "geomorphological markers" in other volcanic settings.Furthermore, morphological changes are often related to sudden and catastrophic phenomena (pyroclastic flows or landslides), which represent a main geohazard for the surrounding communities.This issue is particularly relevant in the case of insular volcanoes, where these phenomena can directly impact subaerial/submarine infrastructures or indirectly generate tsunami waves, as recently observed for the 2002 landslide at Stromboli [4,5] or the 2018 Anak Krakatau collapse [6].
Notwithstanding the direct link between the morphogenesis of subaerial and submarine volcanic flanks, integrated studies dealing with land-sea correlations are rare and have been mainly limited to large sector collapses involving the entire volcanic edifice [7][8][9].In the last few decades, some attempts have been carried out in this framework where shallow-water morphologies (e.g., insular shelf) were correlated with the facing island geology, providing useful insights on the island's vertical movements, the identification of paleo-centres or the mechanisms occurring when lava flows enter the sea [10][11][12][13].Even fewer studies have tried to integrate subaerial and submarine morphological changes during eruptive events affecting volcanic flanks [14][15][16][17][18].At Stromboli, the latter studies were mainly focused on the Sciara del Fuoco depression (Figure 1), representing the last of repeated lateral collapses that affected the NW flank in the last 13 ka [19,20].Since its formation, the steep SdF slope has funneled towards the sea a large amount of volcaniclastic material produced by the persistent Strombolian activity [21] emitted by the summit crater terrace (Figure 1b).This persistent activity is occasionally interrupted by paroxysmal or effusive eruptions [22], during which the stress-strain conditions of the SdF slope can be modified, leading to the development of mass-wasting processes that range over different spatial scales [4,[23][24][25][26].During the 2002 eruptive crisis, two landslides (mobilized volumes in the order of 10 × 10 6 m 3 ) were triggered by dyke intrusion in the upper part of the edifice.These almost simultaneously affected both the submarine and subaerial part of the SdF and generated tsunami waves with run-up values up to 10 Since its formation, the steep SdF slope has funneled towards the sea a large amount of volcaniclastic material produced by the persistent Strombolian activity [21] emitted by the summit crater terrace (Figure 1b).This persistent activity is occasionally interrupted by paroxysmal or effusive eruptions [22], during which the stress-strain conditions of the SdF slope can be modified, leading to the development of mass-wasting processes that range over different spatial scales [4,[23][24][25][26].During the 2002 eruptive crisis, two landslides (mobilized volumes in the order of 10 × 10 6 m 3 ) were triggered by dyke intrusion in the upper part of the edifice.These almost simultaneously affected both the submarine and subaerial part of the SdF and generated tsunami waves with run-up values up to 10 m [4,26,27], highlighting the geohazard related to these processes.Since then, the morphological evolution of the SdF slope has been carefully monitored through repeated topo-bathymetric surveys, which are part of a more complex monitoring network encompassing ground-based devices as well as satellite monitoring systems [28].In this study, we show the results of the integrated land-sea morphological monitoring of the SdF slope in the period encompassing two different topo-bathymetric acquisitions in February 2020 and July 2021.During that period, Stromboli exhibited several phases of intense volcanic activity (Table 1), with strong explosions and frequent lava overflows [29].Some of these overflows (31 March 2020 and 19 May 2021) were associated with the collapse of the rim of the NE crater area (NEC hereafter), which produced avalanches of glowing rocks down the SdF.Of these collapses, the one with the largest volume (~0.5 × 10 6 m 3 ) occurred on 19 May 2021.A similar phenomenon was also observed at the beginning of the 12 January 2013 lava overflow [30] and on the 6 August 2014 event [31].The 19 May 2021 NEC rim collapse evolved into a small pyroclastic density current (PDC, hereafter).Once it reached the sea, the PDC generated a tsunami wave of approximately 20 cm in height [32], and its upper and more diluted part travelled over the sea surface for about 1 km [33].The PDC was followed by the emission of lava flows both in the following five days and in June 2021 (Table 1 and Figure 2).

Materials and Methods
The study was mainly based on the topo-bathymetric comparison of the subaerial and shallow-water submarine part of the SdF slope in the period between February 2020 and July 2021.All the available remote sensing and geophysical data were analyzed and This provided the possibility of comparing high-resolution digital elevation models (DEMs), obtained by photogrammetry using images collected by PLÉIADES satellites and Unmanned Aerial Vehicle (UAV, hereafter), with repeated multibeam (MBES, hereafter) bathymetries collected before and after the 2021 event.This data allowed us to observe the downslope evolution of these phenomena and to reliably estimate the eroded and accreted volumes along the entire volcanic slope (both subaerial and submarine).This information provided useful insights for hazard assessment because the volume and behaviour of these flows strongly control their tsunamigenic potential, as highlighted by the PDCs which entered the sea in July and August 2019, producing tsunami waves approximately 1.5 m in height [16,24].In this paper, we also show preliminary results of the first comparison between the COSMO-SkyMed Synthetic Aperture Radar (CSK-SAR) amplitude image of the SdF subaerial flank with MBES and Side Scan Sonar (SSS, hereafter) backscatter of the submarine flank with the aim of better constraining the main superficial textural changes associated with the emplacement of different volcanic deposits along the SdF slope.

Materials and Methods
The study was mainly based on the topo-bathymetric comparison of the subaerial and shallow-water submarine part of the SdF slope in the period between February 2020 and July 2021.All the available remote sensing and geophysical data were analyzed and integrated (Table 2 and Figure 3).

Topo-Bathymetric Monitoring of the SdF Slope and CKS-SAR Amplitude Image
The topographic change detection of the SdF subaerial slope was obtained by comparing DEMs generated both from PLÉIADES-1 tri-stereo satellite imagery [35] and a   The topographic change detection of the SdF subaerial slope was obtained by comparing DEMs generated both from PLÉIADES-1 tri-stereo satellite imagery [35] and a UAV survey.PLÉIADES-1 data were acquired on 7 April 2020, 4 May 2021, and 8 June 2021, while the UAV survey was performed on 7 July 2021 (Table 2 and Figure 3).To assess the accuracy of the heights and their horizontal position in the PLÉIADES-1 DEM, ground control points (GCPs) were collected (cartographic XY standard deviation: 0.15 m).
The photogrammetric UAV surveys were carried out with a Saturn-mini drone, patented by the Department of Earth Sciences of the University of Florence and specifically developed for extreme environmental uses [36].The drone was equipped with a GNSS receiver with a single frequency, multi-constellation, non-RTK capable device and carried a Canon IXUS 160 camera with 20 Megapixel resolution.The flight was performed on 7 July 2021 by acquiring 7 strips parallel to the slope at a height of about 80 m AGL (above ground level) with nadiral camera position.Previous flights carried out over the same area showed deformation of several meters in the centre of the SdF slope [37] due to intrinsic limitations of the photogrammetric technique and to the absence of GCPs in the central part of the slope.To overcome these issues and correct the photogrammetric surveys, 15 points from a LiDAR survey carried out in 2012 were used as virtual GCPs to co-register all the UAV surveys (see also Section 2.3).To further improve the results, point clouds obtained by the Structure-from-Motion technique were finally realigned to enhance the georeferencing accuracy using the Fine Registration tool, Iterative Closest Point (ICP), available within the CloudCompare software.After all the processing steps, the generated DEM had a resolution of 4.8 cm/pixel with a GCP georeferencing average error of 1.2 m.CSK X-band (3.1 cm of wavelength) SAR images were acquired five times from 18 December 2020 to 30 July 2021 in descending orbit (Table 3 and Figure 3) for both the analysis of SAR change detection (from December 2020 to June 2021) and the backscattering SAR signal (30 July 2021).Data from SAR sensors in volcanic environments have often been used to map areas affected by lithological and morphological changes, such as areas impacted by eruptive and post-eruptive (landslides or floods) phenomena [23,[38][39][40].
The SAR backscattering is determined by several factors, such as the local morphology and the surface microrelief related to the grain-size and the dielectric constant of the material at the surface [41,42].To define which roughness scale affects the backscatter properties, the Rayleigh criterion [43] was applied.In this way, the root mean squared height (hrms) variation on horizontal surfaces was evaluated following [41], and, for the CSK-SAR images used in this study, the hrms was 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 coarsegrained (i.e., blocks and bombs-dominated deposits, such as grain flows or 'a'ā lava surfaces).Local morphology produces irregularities with 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 morphological effects from the grain-size influence.Therefore, the term "roughness" is used to represent a combination of both factors [41].The CSK-SAR images were processed using the free SNAP software [44] both for the SAR change detection analysis and for the comparison with the SSS backscattering signal.The elaboration of the SAR images in SNAP consisted of: (i) the co-registration of all the images using the offset refinement based on the Shuttle Radar Topography Mission (SRTM) DEM at 1 arc second (made available by the United State Geological Survey) forming one unique stack and (ii) cropping around the target area.The cropped stack was then geocoded for correcting SAR geometric distortions using SRTM DEM to produce SAR orthorectified map-projected images.The backscattered intensity of each image was transformed into an amplitude image and then decibel-scaled, converting the data into a virtual band with the expression 10 × log10 (amplitude).In addition, speckle filtering was applied to improve the quality of the images, reducing the salt-pepper texture (speckle) [45].
RGB colour composites using multitemporal CSK-SAR images enabled the detection and interpretation of surface changes of the SdF slope.The use of the colour composition using the band and its ratio is a simple and fast approach to highlight strong changes [41].The colour composition was created taking into consideration a couple of images and inserting the oldest in the red band, the newer in the green band, and its ratio (old/new) in the blue band.

Bathymetric Monitoring of the SdF Slope and MBES/SSS Backscatter Data
Seafloor changes were computed through the difference between multibeam bathymetries collected in February 2020 and July 2021 (Table 2 and Figure 3; for details see Section 2.3).The first bathymetric survey was realized between 1 and 400 m water depth on 18 February 2020 by "Arena Sub SRL" onboard the small vessel "Valerio" using the multibeam system Reson Seabat 7125 working at a frequency of 200 kHz.The second bathymetric survey was realized between 5 and 420 m water depth on 27 July 2021 by "PRISMA SRL" onboard the 16 m-long "Euribia" vessel using the multibeam system Teledyne Reason SeaBat T50-P working at a frequency of 200 kHz.In both surveys, data were positioned through Trimble Applanix Pos Mv Wave Master II to guarantee a spatial positioning on the order of few decimeters at maximum.Each survey was preceded by the collection of sound velocity profile data through a Sontek CastAway probe.During the surveys, a mini SVS Valeport probe was used to update in real-time the sound velocity values close to the flat face of the multibeam transducer, because any error in this value would introduce an angular error both in the beam angle and the ray-tracing.For both surveys, tidal corrections were applied using Strombolicchio and Ginostra tide gauge data recorded by I.S.P.R.A, and data processing was performed using Teledyne PDS 2000 through the application of geometrical/statistical filters and the manual editing of fake soundings.Data were gridded with a cell size of 3 m for the entire dataset, and a grid with a cell size of 1 m was realized at depths <100 m.Gridded data were visualized as shaded relief maps, isobaths, perspective images, or slope gradient maps through the Global Mapper 15 software (Hallowell, ME, USA).
Backscatter data were collected both with MBES and SSS systems; the latter is a towed device, emitting an acoustic impulse from its side [46].In both systems, the backscattered energy is the result of the interaction between the emitted acoustic energy and the seafloor.Multi-directional scattering occurs at the interface, and some of the scattered energies are scattered back to the instrument.The backscatter strength depends on several factors, including the angle of incidence, the roughness of the seafloor, and the scattering behaviour of the material at the seafloor [47].Depending on the frequency of the acoustic signal, some acoustic energy is refracted into the sediments and scattered at deeper interfaces, resulting in the registration of volume backscatter.This effect increases with decreasing sonar frequency; thus, in our case, this effect can be considered negligible considering the very high frequency used in both systems, i.e., 200 and 600 kHz for MBES and SSS devices, respectively.
During the 2021 bathymetric survey, MBES backscatter data were collected and visualized through Global Mapper 15 in grey-colour scale.In addition, SSS data were collected in the first 200 m water depth along 4 lines using the towed EDGETECH 4125 system working at a frequency of 600-1600 kHz.The collected SSS device generally flew 10-20 m above the seafloor, with 125 m-wide channels on each side.SSS data were merged to obtain a single mosaic in geotiff format, then imported in Global Mapper 15 and co-registered to the bathymetry using homologous features (scarps, prominent blocks, and so on) identified in the two datasets.

DEM Co-Registration, Error Estimation, and Topographic Change Detection
Topographic change detection using multi-temporal DEMs was performed by differencing two DEMs of the same area derived from data acquired at different times to create elevation difference maps.This calculation is typically affected by errors associated with mismatches between the two DEMs, which leads to artefacts in elevation differences (∆h) [48].These errors can be quantified by measuring the DEM differences in areas where the two DEMs are expected to display no relevant change.In this work, the RMS displacement errors of the subaerial DEMs were calculated, taking as reference a 1-m-resolution DEM with high vertical and horizontal accuracy.The used reference DEM is a LiDARderived DEM elaborated from data acquired during an airborne survey that was carried out in May 2012 using a Leica ADS80 sensor, which has an instrumental vertical and horizontal accuracy of 0.10-0.20 and 0.25 m, respectively [23].The resulting RMS displacement errors (RMSE0 in Table 4) for the PLÉIADES-derived DEMs before the co-registration procedure were very high, being of 4.35 m, 4.81 m, and 4.01 m for the April 2020, the May 2021, and the June 2021 DEMs, respectively.The resulting RMS displacement error for the UAV-derived July 2021 DEM was 1.39 m.
These errors can be greatly reduced by accurately co-registering these DEMs to the 2012 reference DEM in the benchmark areas, i.e., areas without significant morphological changes.In this work, the co-registration was based on the minimization of the root mean square (RMS) difference between each DEM and the 2012 DEM by iteratively varying the three angles of rotation, the translation, and the magnification or reduction factor of the DEM by using a custom-made algorithm based on the MINUIT numerical minimization software library (e.g., [49]).The MINUIT co-registration procedure lowered the RMS displacement errors for the PLÉIADES DEMs from values above 4 m to values around 1 m, i.e., 0.61 m, 1.26 m, and 0.89 m for the April 2020, the May 2021, and the June 2021 DEMs, respectively (RMSE1 in Table 4).The RMS displacement error for the UAV-derived July 2021 DEM was reduced from 1.39 m to 0.9 m.In this work, the following differences between the co-registered DEMs were calculated for the subaerial part: (i) May 2021-April 2020 DEMs; (ii) June 2021-May 2021 DEMs; (iii) July 2021-June 2021 DEMs; (iv) July 2021-April 2020 DEMs, which covered the whole period of data acquisitions.The RMS difference in areas not affected by changes (σ ∆Z ) for these differences were, respectively, 1.07 m, 1.03 m, 1.52 m, and 1.93 m (Table 5).
The bathymetric data were processed using the same method, but the displacement error did not lower significantly after the co-registration procedure, so no corrections were applied.The bathymetric differences between July 2021 and February 2020 had a σ ∆Z of 0.41 m, which is lower than those of the co-registered subaerial DEMs (Table 5).The differences between two successive co-registered DEMs were used to detect and outline the extent of areas that were affected by topographic changes (Figures 4 and 5) and to calculate their volume.The volume (V) emplaced or lost between two acquisitions was calculated from the DEM difference according to the equation: V = Σ i ∆x 2 ∆z i , where ∆x is the cell size and ∆z i is the height variation within the grid cell i.These values were then summed 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., ErrV, high = A σ ∆Z , where A is the investigated area [48].

Topo-Bathymetric Changes
The availability of four DEMs generated from subsequent PLÉIADES-1 tri-stereo satellite imagery and UAV surveys (in April 2020 and May, June, and July 2021; see Table 2 and Figure 3) allowed the detection of the main changes that occurred along the SdF subaerial slope at the various time intervals.The first elevation difference map obtained between May 2021 (before the NE crater rim collapse) and April 2020 DEMs showed minor morphological changes throughout the SdF (Figure 4a).Changes mainly occurred in the central SdF coastal sector, with minor slope erosion (area SE1) or with two areas of accretion (SA1 and SA2), accounting for volumes of −30,000 ± 10,000 m 3 , +25,000 ± 8000 m 3 , +37,000 ± 13,000 m 3 , respectively (Table 6).In contrast, the elevation difference map obtained between June 2021 and May 2021 DEMs showed significant changes both in the outer rim of the NE crater and along the central part of the SdF slope (Figure 4b).In the first site, slope erosion (SE2) dominated with a maximum eroded thickness of approximately 30 m, accounting for a volume of −44,000 ± 3000 m 3 (Table 6).Along the central part of the SdF, the main morphological changes were associated with slope accretion (SA3), which extended from 630 m asl down to the coastline, accounting for +325,000 ± 83,000 m 3 .In detail, a narrow and thin (2-4 m) stripe of slope accretion was present in the upper part of the slope, thickening up to 12 m in the coastal sector and enlarging in a fanshaped geometry downslope of 200 m asl (Figure 4b).This slope accretion was also associated with the formation of a convex-seaward coastline, which prograded for approximately 35 m with respect to the previous (4 May) coastline.This convex-seaward coastline is made up of massive lava flows at the base, now partly eroded and overlain by brecciated lavas and coarse-grained volcaniclastic material, as visible from oblique photos (Figure 2c).Differently, the surrounding areas were dominated by small-scale fan-shaped features related to the emplacement of lava breccias and volcaniclastic coarse-grained material, as recognizable in the June 2021 DEM (Figures 5a and 6b).It is noteworthy that slope

Topo-Bathymetric Changes
The availability of four DEMs generated from subsequent PLÉIADES-1 tri-stereo satellite imagery and UAV surveys (in April 2020 and May, June, and July 2021; see Table 2 and Figure 3) allowed the detection of the main changes that occurred along the SdF subaerial slope at the various time intervals.The first elevation difference map obtained between May 2021 (before the NE crater rim collapse) and April 2020 DEMs showed minor morphological changes throughout the SdF (Figure 4a).Changes mainly occurred in the central SdF coastal sector, with minor slope erosion (area SE1) or with two areas of accretion (SA1 and SA2), accounting for volumes of −30,000 ± 10,000 m 3 , +25,000 ± 8000 m 3 , +37,000 ± 13,000 m 3 , respectively (Table 6).In contrast, the elevation difference map obtained between June 2021 and May 2021 DEMs showed significant changes both in the outer rim of the NE crater and along the central part of the SdF slope (Figure 4b).In the first site, slope erosion (SE2) dominated with a maximum eroded thickness of approximately 30 m, accounting for a volume of −44,000 ± 3000 m 3 (Table 6).Along the central part of the SdF, the main morphological changes were associated with slope accretion (SA3), which extended from 630 m asl down to the coastline, accounting for +325,000 ± 83,000 m 3 .In detail, a narrow and thin (2-4 m) stripe of slope accretion was present in the upper part of the slope, thickening up to 12 m in the coastal sector and enlarging in a fan-shaped geometry downslope of 200 m asl (Figure 4b).This slope accretion was also associated with the formation of a convex-seaward coastline, which prograded for approximately 35 m with respect to the previous (4 May) coastline.This convex-seaward coastline is made up of massive lava flows at the base, now partly eroded and overlain by brecciated lavas and coarse-grained volcaniclastic material, as visible from oblique photos (Figure 2c).Differently, the surrounding areas were dominated by smallscale fan-shaped features related to the emplacement of lava breccias and volcaniclastic coarse-grained material, as recognizable in the June 2021 DEM (Figures 5a and 6b).It is noteworthy that slope accretion shifted laterally to slope erosion (SE3), especially in the upper part of the slope (above 200 m asl), giving way to a 200 m-wide erosive stripe linked upslope with the main SE2 area (Figure 4b).The volume related to SE3 was difficult to estimate as the slope erosion was very thin (on average 1.16 m) slightly above the σ ∆Z (Table 6).
The elevation difference map obtained between the July 2021 and June 2021 DEMs (Figure 4c) showed very minor morphological changes with respect to the previous map and was mainly characterized by slope accretion both in the summit crater area (SA4) and along the central part of the SdF (SA5), especially close to the coastal sector.The volumes estimated were of +20,000 ± 6000 m 3 and +50,000 ± 39,000 m 3 for SA4 and SA5, respectively (Table 6).
Regarding the submarine part of the SdF, the elevation difference map obtained between the February 2020 and July 2021 bathymetries (Figure 5a) showed that the main morphological changes occurred in the first 200 m of the water depth (wd, herafter), dominated by seafloor erosion (ME) and alternated with smaller areas of seafloor accretion (MA).
In detail, two main erosive areas (ME1 and ME2) are present in the central part of the SdF: (i) ME1 affects a surface of approximately 48,500 m 2 between 5 and 175 m wd, with a maximum eroded material of 13 m, accounting for a volume of approximately 213,000 ± 20,000 m 3 , and (ii) ME2 affects a surface of approximately 28,700 m 2 between 10 and 160 m wd, with a maximum eroded material of 7 m, accounting for a volume of approximately 68,000 ± 12,000 m 3 (Table 6).Both erosive areas are related to the formation of downslope-elongated depressions on the 2021 bathymetry, which are partially floored by small-scale fan-shaped features (Figure 6b).The two erosive areas are divided by a small area of marked seafloor accretion (MA1), which encompasses a surface of approximately 13,000 m 2 between 3 and 180 m wd.MA1 shows a maximum accreted thickness of 10 m, accounting for a volume of approximately 35,000 ± 5000 m 3 .In the 2021 bathymetry, MA1 corresponded to the development of an 80-100 m wide morphological high, elongated downslope for approximately 300 m (Figure 6b).This ridge was characterized by relatively steep (ranging between 15 • and 30 • ) flanks, while its upper surface showed an uneven morphology, mostly made up of minor ridges (with local slope gradients up to 60 • ), metric blocks, and small fan-shaped lobes.The shallower part of the ridge was instead cut by a terraced feature (NSDT in Figure 6b), whose outer edge is located here around 4 m wd.Small erosive areas were also recognizable both in the SW (ME3 and ME4) and NE part (mostly ME5 and ME6) of the SdF submarine slope.They affect variable from about 7000 m 2 to 13,000 m 2 , with a maximum eroded thickness of 6 m.The computed mobilized volumes were approximately comprised of between approximately −10,000 and −20,000 m 3 (Table 6).Even smaller (1000-2000 m 2 ) and thinner (generally 1-3 m deep, at the limit of data resolution) erosive areas were present in the sector encompassed between ME1 and ME5-6 (minor erosive areas in Figure 5), mobilizing a total cumulative volume of −24,000 ± 12,000 (Table 6).All the erosive areas reflect the occurrence of shallow and smallscale landslide scars in the 2021 bathymetry, whose headwalls are mostly located around 10 m wd, affecting the edge of a terraced feature (NSDT in Figure 6b), recognizable along most of the submarine SdF slope.
At greater depths, the morphological changes are mostly associated with seafloor accretion, concentrated in three main areas (MA2, 3, and 4) located between 200 and 400 m wd in the central and NE part of the SdF; they developed downslope of the erosive areas ME1 and ME5-6.These sectors involve surfaces ranging from 42,000 to 58,000 m 2 , with a maximum accreted thickness of 5-8 m; estimated volumes for MA2, MA3, and MA4 account for +69,000 ± 17,000, +84,000 ± 24,000, and +115,000 ± 17,000 m 3 , respectively (Table 6).Two smaller areas (MA5-6) of seafloor accretion are also recognizable in the SW part of the SdF between 80 and 170 m wd, downslope of the erosive areas ME3 and ME4 (Figure 5a).MA5 and MA6 have surfaces of approximately 13,000-18,000 m 2 , an average thickness slightly over 1 m, accounting for a volume of +22,000 ± 8000 m 3 and +16,000 ± 5000 m 3 , respectively.Small erosive areas were also recognizable both in the SW (ME3 and ME4) and NE part (mostly ME5 and ME6) of the SdF submarine slope.They affect variable from about 7000 m 2 to 13,000 m 2 , with a maximum eroded thickness of 6 m.The computed mobilized volumes were approximately comprised of between approximately −10,000 and −20,000 m 3 (Table 6).Even smaller (1000-2000 m 2 ) and thinner (generally 1-3 m deep, at the limit of data resolution) erosive areas were present in the sector encompassed between ME1 and ME5-6 (minor erosive areas in Figure 5), mobilizing a total cumulative volume of −24,000 ± 12,000 (Table 6).All the erosive areas reflect the occurrence of shallow and small-scale landslide scars in the 2021 bathymetry, whose headwalls are mostly located around 10 m wd, affecting the edge of a terraced feature (NSDT in Figure 6b), recognizable along most of the submarine SdF slope.
At greater depths, the morphological changes are mostly associated with seafloor accretion, concentrated in three main areas (MA2, 3, and 4) located between 200 and 400 m wd in the central and NE part of the SdF; they developed downslope of the erosive areas ME1 and ME5-6.These sectors involve surfaces ranging from 42,000 to 58,000 m 2 , with a maximum accreted thickness of 5-8 m; estimated volumes for MA2, MA3, and MA4 account for +69,000 ± 17,000, +84,000 ± 24,000, and +115,000 ± 17,000 m 3 , respectively (Table 6).Two smaller areas (MA5-6) of seafloor accretion are also recognizable in the SW part of the SdF between 80 and 170 m wd, downslope of the erosive areas ME3 and ME4 (Figure 5a).MA5 and MA6 have surfaces of approximately 13,000-18,000 m 2 , an average thickness slightly over 1 m, accounting for a volume of +22,000 ± 8000 m 3 and +16,000 ± 5000 m 3 , respectively.

CSK-SAR Amplitude Image and MBES/SSS Backscatter Data
The investigation of SAR amplitude changes was performed considering the variation of amplitude between March 2021 and December 2020 (Figure 7a), May 2021 and March 2021 (Figure 7b), and June 2021 and May 2021 (Figure 7c) images.In the first period considered (image of March 2021 versus image of December 2020), it was possible to identify the overflow that occurred on 24-25 January 2021 in the central portion of the SdF, highlighted by a narrow stripe with a markedly rougher surface (Figure 7a).In the following period (image of May 2021 versus image of March 2021, Figure 7b), no further relevant changes were identified.In particular, the area affected by lava overflow in January 2021 was mostly unchanged, apart from a slight tendency to smoothing only in the upper portion.After that, between June 2021 and May 2021, the area affected by the January 2021 lava overflow showed no detectable changes (light grey colour in the amplitude image of Figure 7c) but is contoured by evidently rougher portions on both sides, due to the rolling of coarse-grained volcaniclastic material and partially superimposing on the lava overflow (see also Figure 2c).

CSK-SAR Amplitude Image and MBES/SSS Backscatter Data
The investigation of SAR amplitude changes was performed considering the variation of amplitude between March 2021 and December 2020 (Figure 7a), May 2021 and March 2021 (Figure 7b), and June 2021 and May 2021 (Figure 7c) images.In the first period considered (image of March 2021 versus image of December 2020), it was possible to identify the overflow that occurred on 24-25 January 2021 in the central portion of the SdF, highlighted by a narrow stripe with a markedly rougher surface (Figure 7a).In the following period (image of May 2021 versus image of March 2021, Figure 7b), no further relevant changes were identified.In particular, the area affected by lava overflow in January 2021 was mostly unchanged, apart from a slight tendency to smoothing only in the upper portion.After that, between June 2021 and May 2021, the area affected by the January 2021 lava overflow showed no detectable changes (light grey colour in the amplitude image of Figure 7c) but is contoured by evidently rougher portions on both sides, due to the rolling of coarse-grained volcaniclastic material and partially superimposing on the lava overflow (see also Figure 2c).4b,c).More interestingly, high CSK-SAR amplitude values were recorded all along the coastal strip of the SdF.This high-amplitude area matches a sandy/gravelly accumulation at the foot of the slope, with sparse sub-metric and metric blocks alternated with more coherent lava flows recognizable in the 7 July 2021 UAV optical image (Figure 8d) and oblique photos (Figure 2c).Lowamplitude values are generally found within topographic lows encompassed between the main lava flows and are commonly floored by relatively finer-grained volcaniclastic material, as recognizable by optical images collected both by the PLÉIADES-1 satellites and UAV survey.
Regarding the submarine part, MBES/SSS backscatter data show the presence of medium-high backscatter values all along the nearshore sector, within the first 10 m wd (Figure 8).This area is bounded downslope by a high backscatter stripe (Figure 8d), roughly oriented parallel to the coastline, matching the marked break-in slope associated with the edge of the submarine terrace (NSDT in Figures 6 and 8d).The flat terrace is characterized by overall smooth and medium-high backscatter values (Figure 8c,d), with the local accumulation of blocks ranging from sub-metric to metric in diameters or small lobate lava flows.4b,c).More interestingly, high CSK-SAR amplitude values were recorded all along the coastal strip of the SdF.This high-amplitude area matches a sandy/gravelly accumulation at the foot of the slope, with sparse sub-metric and metric blocks alternated with more coherent lava flows recognizable in the 7 July 2021 UAV optical image (Figure 8d) and oblique photos (Figure 2c).Lowamplitude values are generally found within topographic lows encompassed between the main lava flows and are commonly floored by relatively finer-grained volcaniclastic material, as recognizable by optical images collected both by the PLÉIADES-1 satellites and UAV survey.
Regarding the submarine part, MBES/SSS backscatter data show the presence of medium-high backscatter values all along the nearshore sector, within the first 10 m wd (Figure 8).This area is bounded downslope by a high backscatter stripe (Figure 8d), roughly oriented parallel to the coastline, matching the marked break-in slope associated with the edge of the submarine terrace (NSDT in Figures 6 and 8d).The flat terrace is characterized by overall smooth and medium-high backscatter values (Figure 8c,d), with the local accumulation of blocks ranging from sub-metric to metric in diameters or small lobate lava flows.
At greater depths, the quality of the SSS data increases with respect to the MBES backscatter, mainly because the former was acquired closer to the seafloor.On the other hand, the SSS data georeferencing was poorer and needed to be re-georeferenced with MBES bathymetry and backscatter through the matching of homologous morphologies.In the deeper areas, MBES and SSS data show the presence of areas with relatively medium backscatter above the main morphological highs, where downslope alignments of submetric or metric blocks are locally present (Figure 8c,d).The flanks of these morphological highs and the topographic lows/depressions are typically characterized by relatively higher backscatter values.Specifically, the two depressions related to ME1 and ME2 areas in the central part of the SdF (Figures 5a and 8c,d) are characterized by homogenous and very high-backscatter facies.At greater depths, the quality of the SSS data increases with respect to the MBES backscatter, mainly because the former was acquired closer to the seafloor.On the other hand, the SSS data georeferencing was poorer and needed to be re-georeferenced with MBES bathymetry and backscatter through the matching of homologous morphologies.In the deeper areas, MBES and SSS data show the presence of areas with relatively medium backscatter above the main morphological highs, where downslope alignments of sub-metric or metric blocks are locally present (Figure 8c,d).The flanks of these morphological highs and the topographic lows/depressions are typically characterized by relatively higher backscatter values.Specifically, the two depressions related to ME1 and ME2

Reconstruction of the Main Eruptive and Erosive-Depositional Phenomena during the Monitoring Period
The topo-bathymetric changes observed by integrating remote sensing and geophysical analysis can be linked to specific eruptive or erosive-depositional processes that occurred in the SdF between February/April 2020 and July 2021.This link is quite clear for the subaerial slope, in relation to the high frequency of available DEMs (four in the investigated time span) and the direct observation of phenomena that occurred during the monitored period.
The limited morphological changes that occurred in the period April 2020-May 2021 (11 months) throughout the SdF are typical of the relatively slow morphological evolution of the SdF slope during "normal" eruptive stages with sporadic lava overflow events during phases of more intense/frequent activity (April 2022, January-February 2021), as already observed in previous monitoring activities [16].The large time gap that occurred between the two DEMs and their relatively low resolution did not allow clear detection of the lava overflow that occurred between 24 and 25 January 2021 along the central portion of the SdF slope.It is, however, recognizable by a roughness band through the comparison of CKS-SAR amplitude images collected during this period (Figure 7).Analysis of these images up to the 11 May 2021 (Figure 7a,b) also showed an almost unchanged texture for the surrounding areas, in agreement with the minimal morphological changes that occurred during this period.These findings highlight the advantage of integrating different remote sensing techniques and frequent surveys for a more accurate assessment of the morphological evolution of the Sciara del Fuoco, otherwise some eruptive or landslide phenomena would have gone undetected, such as the January 2021 lava overflow.
The large morphological changes recognized in the central part of the SdF between May 2021-June 2021 (1 month) can be related to the NEC rim collapse, which occurred on 19 May 2021, and the subsequent lava flows.The rim collapse was clearly evidenced by the maximum eroded thickness in the elevation difference map (SE2 in Figure 4b) as well as by the rougher area in the SAR change amplitude analysis (Figure 7c).This collapse event generated a peculiar type of PDC (i.e., a hybrid phenomenon between a landslide and a PDC, witnessed at volcanoes fed by mafic magmas, e.g., [33]), with an estimated volume of material of between 10,000 and 70,000 m 3 [32] that propagated down the SdF slope at a speed of approximately 50 km/h.Our volumetric estimation provided a minimum figure of approximately 44,000 m 3 for this event (SE2 in Figure 4b), taking into account that the post-collapse survey was realized approximately 20 days after the event, so that the crater rim collapse had already been partially filled by the subsequent volcanic activity.This filling continued over time, as evidenced by the SA4 accretion recorded during the next month (Figure 4c), where a total volume of approximately 20,000 m 3 was emplaced in the previous collapse depression.
According to [33], the ash-cloud deposit generated by the PDC was dominated by remobilized material ingested by the current during its movement down the SdF slope, thus indicating that significant erosion occurred.However, PDC erosion on the subaerial slope was barely recognizable on the June 2021-May 2021 elevation difference map (SE3 in Figure 4b).This can be again explained by the time gap between the collapse event and the DEM acquisition, which also recorded the infilling by lava flows and breccias in the subsequent five days (SA3 in Figure 4b) that propagated into the sea (MA1 in Figure 5a).Taking into account the spatial distribution of the slope erosion (SE3 in Figure 4b), we can tentatively reconstruct a 200 m-wide erosive strip, morphologically linked with the upslope crater rim collapse, but a reliable estimate of the eroded volume is very challenging and likely misleading, indicating the need for surveys immediately after the event.More interestingly, the spatial distribution of slope accretion, along with the textural analysis of optical and CKS-SAR images, provided insights for better constraining the emplacement of lava flows and volcaniclastic material in the SdF slope.The most interesting finding is the marked increase of slope accretion below 200 m asl, reaching its maximum at the coastline and forming an overall fan-shaped geometry (Figures 4b and 6b).Such a geometry could be attributed to a difference in slope gradients between the upper and lower slope of the SdF, but the comparison of along-slope sections performed in the pre-and post-19 May DEMs did not show significant variations in slope gradients, with an almost constant gradient of around 35 • for the entire slope.This observation could, however, be partially biased by the lack of a subaerial DEM just before the emplacement of the lava flows (possibly recording an original change in slope gradients due to the erosion exerted on the SdF slope by the NEC rim collapse-related PDC).It is (at least) suspicious that the post-collapse DEM recovered the same constant gradient profile observed in the pre-collapse DEM.An alternative explanation for the marked increase in thickness at the coastal area could be related to the abrupt decrease in velocity of the main lava flows once they reached the sea, thus enabling an increase in lava thickness due to inflation processes, similarly to that recognized in modern and ancient lava flows [50,51] as well as in analogue experiments [52].This process could have been particularly efficient in the sector where the largest coastline progradation was observed and where oblique photos showed that massive and brecciated lava flows were able to reach the coastline (Figure 2c).The marked progradation of the coastline in a very small sector could have also been favoured by the slope topography at that time, mainly due to the emplacement of the 24-25 January 2021 overflows that might have acted as a lateral ridge confining the May-June 2021 overflows.This inference is supported by the comparison of elevation difference maps (Figure 4b-d), CKS-SAR amplitude images (Figure 7b,c), and by the oblique photo taken on 27 July 2021 (Figure 2c).
In contrast, the surrounding areas are dominated by small-scale fan-shaped features and narrow elongated ridges, overall creating a composite volcaniclastic wedge (Figure 6b), a feature typically observed in the lower and middle part of the SdF slope during eruptions fed by high-elevation vents, as those that occurred in [2002][2003]2014, and 2019 events [15,16,53].This volcaniclastic wedge was fed by lava breccias related to lava flows or overflows, remaining mostly confined in the upper slope due to low effusion rates, and by a higher input of volcaniclastic material due to the brecciated NEC rim, increasing the frequency of debris/rock fall processes occurring along the SdF and making them able to reach the coastline (Figures 2c and 6b).The emplacement of this material is also responsible for the marked variation in SAR amplitude observed in Figure 7c as well as for high-amplitude values recorded in the central part of the SdF on the 30 July CKS-SAR image (Figure 8a).
The last elevation difference map, reconstructed for the period June 2021-July 2021 (1 month), shows the occurrence of minor morphological changes in the central part of the SdF.These can be related to minor lava overflows emitted during June as well as to less recurrent debris/rock fall events along the SdF slope due to the formation of a new crater rim.
The interpretation of the bathymetric changes is more complex because the two marine surveys were realized in February 2020 and July 2021, with a time gap of 1 year and 6 months, thus encompassing multiple eruptive and erosive-depositional phenomena.However, the integration of submarine data with the previous results from the subaerial slope and the 20 years' experience on submarine monitoring of the SdF allowed us to make reliable inferences on the main processes associated with the observed seafloor variations.Indeed, previous monitoring activities showed that the main morphological changes in the submarine SdF were observed during an eruptive crisis or increased periods of Strombolian activity.The "ordinary" slope evolution was, however, limited to minor seafloor accretion related to the partial dismantling of the volcaniclastic wedge after the end of stronger eruptive periods [16] and/or small landslide/erosive processes affecting the shallowest part of the SdF.In the 2020-2021 difference map (Figure 5), the latter processes were testified by the small erosive areas recognized along the SW (ME3-4) and NE (ME5-6 and minor erosive areas) sectors of the SdF, where the slope failures were likely triggered by major storms that struck the NW flank of Stromboli during the winter period.The shallow part of the SdF submarine slope is mainly made up by a chaotic arrangement of coarse-grained volcaniclastic material (volcaniclastic sand and gravel, lava breccias, and blocks, [54]), and stability analyses showed that this slope is prone to small-scale shallow landslide processes [25,55].Particularly, the nearshore submarine depositional terrace edge (NSDT in Figures 6 and 8d), located around 10 m water depth, is an area prone to the development of small-scale slope instabilities, similarly to what was observed along the shallow submarine part of other volcanic islands [56,57].This nearshore terrace is a common feature developed along volcanic islands, representing the morphological expression of a depositional body with prograding geometry, whose edge depth can be related to the local storm-wave base level [2,[58][59][60].The terrace edge represents an area where slope gradients markedly increase and a change from sediments repeatedly reworked by wave action to a slope dominated by avalanching process, thus favouring the development of small-scale slope instabilities.It is noteworthy that the runout of the mobilized material is generally low (see MA2-6 in Figure 5a) considering the steep submarine slope, but it is likely related to the fact that this gravity instability mobilized the shallow part of the SdF, made up by coarse-grained volcaniclastic material, with movement ceasing when the slope recovers its angle of repose, similar to what was observed in the previous monitoring activities at Stromboli and at the submarine slope of Hawaiian 'a'ā lava-fed deltas [18].
The two larger and deeper erosive areas (ME1-2) observed in the central part of the SdF exactly match the topographic changes of the subaerial slope associated with the 19 May 2021 NEC rim collapse (Figure 5b), highlighting a cause-effect relationship between them and supporting the erosive potential of these processes, as indirectly suggested by [33].This relationship is also in agreement with the lack of a submarine depositional terrace in the ME1-2 depressions down to 6 m of water depth (limit of the bathymetric survey), in contrast to the surrounding areas, which were not affected by localized erosion (Figures 6b and 8c,d).Marked erosion was also generated in the central part of the SdF due to the entry into the sea of the 3 July 2019 PDC [16].In the latter case, the survey was realized just a few days after the event, clearly recording the original erosive signature of the current once it entered the sea.In the present study, the marine survey was undertaken two months after the event, thus also recording the morphological changes associated with the emplacement of subsequent lava flows and breccias, as well as of recurrent small-scale debris/rock fall events.This is clearly evidenced by the seafloor accretion MA1 interposed between the two erosive areas ME1-2, the former exactly matching the main coastline progradation due to the May 2021 lava flows.It is noteworthy that the MA1 small ridge (Figure 6b) is similar to other morphological highs formed during the 2007 and 2014 lava entrance to the sea and lava delta formation [15,61], made up by variable percentages of coherent lavas and chaotic breccias.In the 2021 case, the small size of the MA1 ridge coupled with its relatively smooth morphology and slope gradients of up to 30 • (except for very narrow ridges with gradients up to 60 • ) suggests the prevalence of lava breccias and volcaniclastic materials rather than coherent lava flows.It is also noteworthy that a submarine terrace developed on top of the MA1 ridge at water depths (around 4 m) shallower than the surrounding terrace edge, and coastal platforms developed at the SdF shoulders.This lower depth is likely in agreement with the lack of significant storm events between the ridge formation and the survey.The shallower slope erosion in ME2 with respect to ME1 could be explained by the accretion recorded in the subaerial sector between the June and July 2021 DEMs (Figure 4c).Particularly, the cumulate elevation difference map obtained between the first and last available DEMs (February 2020-July 2021) clearly highlights a preferential accretion in the westernmost part of the central SdF (Figures 4d and 5b).
This interpretation suggests a possible significant underestimation of the total volume mobilized in the central part of the SdF during the 2021 collapse event, which could be estimated as between 350,000 and 400,000 m 3 of mobilized material against the 280,000 computed from the elevation difference map.This value was obtained by reconstructing a post-collapse surface without the subsequent slope accretion and enlarging the erosive area coastward.On the other hand, the integration of subaerial and submarine data enables us to provide a minimum volume of volcanic material emplaced during the May overflows, accounting for approximately 360,000 m 3 .This value does not consider both the volcanic material that filled the two erosive depressions ME1 and ME2 and the material emplaced between the coast and the limit of the bathymetric survey.However, the volume emplaced in the submarine setting is much lower than the one emplaced in the subaerial slope, as typically occurs in cases of effusive eruptions fed by a high-elevation vent and with rapidly declining eruptive rates [15 and reference therein].

Textural Variations along the SdF Slope by Integreating CKS-SAR Image and Marine Backscatter Data
The qualitative comparison between CKS-SAR amplitude images and marine backscatter data carried out for the first time in this study indicates that, despite the morphological continuity between the subaerial and submarine slopes observed in the previous section, textural variations along the SdF slope primarily depend on the different processes and characteristics of the subaerial slope, the coastal area, the nearshore, and "deeper" marine areas.
The textural variations of the subaerial slope mainly depend on the type, frequency, and intensity of volcanic activity.Low CKS-SAR amplitude values are commonly observed during periods of low volcanic activity at the summit craters, when relatively fine-grained volcaniclastic material tend to floor the main topographic low, as was locally recognizable in the 30 July 2021 CKS-SAR image (Figure 8a,b).The only exception is the central depressed area of the SdF, which is aligned with the brecciated NE crater rim, where most of the coarse-grained material emitted by this crater was directly funneled seaward through the steep SdF slope.Similarly, the superficial texture of the lava flows, which is, in turn, dependent on the vent location, effusive rates, local steep gradients, and successive erosivedepositional processes, affects the final amplitude response.In general, at Stromboli, lava flows are mainly characterized by rugged superficial textures; thus, they are generally characterized by high-amplitude values, as observed in the SW and NE part of the SdF on the 30 July CKS-SAR image (Figure 8a).
The coastal area is surprisingly characterized by an almost continuous and homogenous stripe of high-amplitude values.This behaviour could be caused by several factors and their combinations.Indeed, the SAR amplitude value can be influenced by [62]: (i) the height of the sea wave, with the X-band being more sensitive with respect to longer wavelengths; (ii) the material size, with the coarser material showing a higher backscattering potential mainly in shoreline detection; (iii) the polarization, with the HH-polarized images (as the used CSK-SAR) usually generating larger contrast along the shoreline; and (iv) incidence angle with respect to the target.In addition, the dielectric constant enhances the radar backscattering signal in non-submerged areas [63].In fact, in a soil-water mixture, such as the shoreline, a dielectric constant is a function of the soil grain dimension and the volume of the moisture content [64].The highest amplitude values along the shoreline can be observed in correspondence with the areas where massive/brecciated lava flows or the accumulation of larger clasts and blocks are present.
The nearshore area is quite variable in terms of backscatter intensity (i.e., somewhat corresponding to the SAR amplitude), being a dynamic area resulting from the continuous interaction between subaerial and marine processes.This sector is mainly affected by wave-action processes that largely rework the volcanic material present, giving rise to a submarine depositional terrace all along the coastal area (Section 4.1).This terrace is generally characterized by smooth and medium-high backscatter areas, related to sandy and gravelly sediments, as ground-truthed by ROV footage (unpublished data), alternating with areas characterized by higher backscatter values linked with the local accumulation of metric and sub-metric blocks or remnants of lava flows.The blocks are mainly due to brecciation or successive erosion of the most recent or largest lava flows (Figure 8c,d) or, more sporadically, to larger avalanching processes occurring on the subaerial slope.In some cases, the largest events can temporarily dismantle part of the terrace, as was also observed during the 19 May 2021 NEC rim collapse, which formed two erosive depressions characterized by homogenous and high backscatter values (ME1 and ME2 in Figure 8c,d) representing a by-pass area for successive avalanching processes.
Textural variations in the deeper area are more variable and strongly dependent on the topography.They mainly reflect the interaction at different temporal and spatial scales between local dynamics associated with small-scale subaerial/submarine avalanching processes and larger events (i.e., effusive eruptions, such as those observed in 2007 and 2014, PDC generated from strong paroxysms in 2019, or crater rim collapses in 2021).These can significantly modify the submarine morphology in the first few hundreds of meters.This area is characterized by medium-high backscatter sectors, with the local accumulation of blocks alternating with smooth and higher backscatter sectors.The former sectors correspond to the morphological ridges formed at the entry point of lava flows into the sea, i.e., remnants of lava deltas, which are often dismantled in metric or sub-metric blocks by lava brecciation and are progressively covered by finer-grained material with respect to the adjacent topographic lows.These latter areas funnel downslope the coarse-grained volcaniclastic material remobilized from coastal and near-shore areas during storm-waves or that are directly transported into the sea by avalanching processes along the subaerial slope, especially during increased periods of Strombolian activity.

Conclusions
This study clearly demonstrated the importance of a multi-methods and multi-temporal land-sea approach through frequent surveys of highly dynamic environments, such as the SdF collapse scar at Stromboli, to reconstruct the submarine/subaerial morphological slope evolution, as well as to estimate the erupted/remobilized volumes with a significant level of detail.This approach has important implications both for better understanding the main eruptive and erosive/landslide processes that continuously reshape active volcanic flanks and in providing insights for geohazard assessment.In particular, considering that these phenomena can generate local but significant tsunami waves and that data on tsunamigenic flows are scarce and difficult to collect, the present work gains further relevance.In fact, in the case of Stromboli, the presence of different buoys in front of the Sciara del Fuoco allows the detection of possible tsunamigenic waves.This combined with the accurate volume assessment carried out here for events of different sizes constitutes a rather unique dataset, which is of the utmost importance in the validation and calibration of tsunami models.This approach also allowed for the assessment of the impacts in terms of seafloor erosion and tsunamigenic potential caused by the hybrid PDC associated with the NEC rim collapse in 2021 and the ones due to the 2019 PDC generated during a strong paroxysmal event.The 2021 event eroded a volume (350,000-400,000 m 3 ) significantly lower than the 2019 PDC (1,300,000 m 3 [16]).In accordance, the larger 2019 PDC generated tsunami waves up to 1.5 m, and they were less than 1 km for the Sciara del Fuoco event in 2019 [16], whereas negligible tsunami waves (maximum height of 0.20 cm, [32]) were detected for the 2021 event.Furthermore, this study has highlighted the need for frequent (i.e., from 1 to 3 surveys per year, based on the state of activity of the volcano), timely (i.e., just after the main eruptive or mass-wasting events), and simultaneous collection of topo-bathymetric surveys as well as their integration with other remote sensing data (CKS-SAR amplitude images or side scan sonar backscatter data).These should be combined with the continuous monitoring of ground deformation along subaerial volcanic flanks through seismic networks, time-lapse SAR imageries, Interferometric SAR data, and other geodetic methods [65].On submarine flanks, the bathymetric surveys can be combined with different monitoring systems for the detection of eruptive activity and/or mass-wasting processes, such as ocean bottom seismometers, hydroacoustic arrays, and smart cables [66,67].
In conclusion, such an accurate and holistic characterization can only be possible if there is the collaboration of an interdisciplinary group of researchers combined with the use of multiple traditional and innovative techniques.In the present case, this holistic approach was fully adopted and implemented, leading to a more comprehensive understanding of the main morphological and textural variations that occur on the entire SdF slope and ultimately to highlight further the complexity of this hazardous and fascinating volcano.

Figure 1 .
Figure 1.(a) Shaded relief and contour map (equidistance 500 m) of Stromboli volcano, with location of the Sciara del Fuoco using the April 2020 and February 2020 digital elevation models for the subaerial and submarine parts, respectively; the inset shows the location of Stromboli volcano (red point) in the Tyrrhenian Sea.(b) This is a 3-D image of the Sciara del Fuoco collapse scar, which can be morphologically divided in the north-eastern (NE), central (C), and south-western (SW) sectors and in the subaerial (coloured) and submarine (grey) parts.

Figure 1 .
Figure 1.(a) Shaded relief and contour map (equidistance 500 m) of Stromboli volcano, with location of the Sciara del Fuoco using the April 2020 and February 2020 digital elevation models for the subaerial and submarine parts, respectively; the inset shows the location of Stromboli volcano (red point) in the Tyrrhenian Sea.(b) This is a 3-D image of the Sciara del Fuoco collapse scar, which can be morphologically divided in the north-eastern (NE), central (C), and south-western (SW) sectors and in the subaerial (coloured) and submarine (grey) parts.

23 Figure 2 .
Figure 2. (a) 19 May 2021 oblique photo of the SdF showing the lava flows that reached the sea after the start of the overflows that caused the collapse of the NE crater-rim; (b) 7 July 2021 photo of the SdF showing the accumulation of lava and volcaniclastic deposits related to the eruptive activity of May and June 2021; (c) 27 July 2021 oblique photo of the central part of the SdF showing the remnants of the 24-25 January 2021 and the 19 May 2021 overflows, partially draped by the successive May-June 2021 volcaniclastic material.

Figure 2 .
Figure 2. (a) 19 May 2021 oblique photo of the SdF showing the lava flows that reached the sea after the start of the overflows that caused the collapse of the NE crater-rim; (b) 7 July 2021 photo of the SdF showing the accumulation of lava and volcaniclastic deposits related to the eruptive activity of May and June 2021; (c) 27 July 2021 oblique photo of the central part of the SdF showing the remnants of the 24-25 January 2021 and the 19 May 2021 overflows, partially draped by the successive May-June 2021 volcaniclastic material.

Figure 3 .
Figure 3. Chronology of data acquisition between February 2020 and July 2021 and main eruptive events in the same period (see also Table2).

Figure 3 .
Figure 3. Chronology of data acquisition between February 2020 and July 2021 and main eruptive events in the same period (see also Table2).

Figure 4 .
Figure 4. Sciara del Fuoco topographic changes recognized between April 2020 and July 2021 through DEM comparisons; SE and SA stand for subaerial erosion and accretion, respectively.(a) May 2021 vs. April 2020 elevation difference map draped over the May 2021 DEM; (b) June 2021 vs. May 2021 elevation difference map draped over the June 2021 DEM.The main morphological changes were due to the 19 May 2021 NEC-rim collapse PDC and subsequent lava flows; (c) July 2021 vs. June 2021 elevation difference map draped over the July 2021 DEM; (d) July 2021 vs. April 2020 comparison shows the morphological changes of the subaerial SdF during the complete investigation period.The frames (a-c) were cleaned for the residual mismatching between DEMs after the co-registration.Frame (d) was not cleaned to show the distribution and the magnitude of the residual mismatching between the first and the last DEMs.

Figure 4 .
Figure 4. Sciara del Fuoco topographic changes recognized between April 2020 and July 2021 through DEM comparisons; SE and SA stand for subaerial erosion and accretion, respectively.(a) May 2021 vs. April 2020 elevation difference map draped over the May 2021 DEM; (b) June 2021 vs. May 2021 elevation difference map draped over the June 2021 DEM.The main morphological changes were due to the 19 May 2021 NEC-rim collapse PDC and subsequent lava flows; (c) July 2021 vs. June 2021 elevation difference map draped over the July 2021 DEM; (d) July 2021 vs. April 2020 comparison shows the morphological changes of the subaerial SdF during the complete investigation period.The frames (a-c) were cleaned for the residual mismatching between DEMs after the co-registration.Frame (d) was not cleaned to show the distribution and the magnitude of the residual mismatching between the first and the last DEMs.

Figure 5 .
Figure 5. (a) Elevation difference map (cleaned for the residual mismatching between DEMs) between the July 2021 and February 2020 bathymetries draped over the 2021 DEM; ME and MA for marine erosion and accretion, respectively.(b) Elevation difference map for the subaerial (July 2021-April 2020) and submarine (July 2021-February 2020) SdF draped in semi-transparency over the July 2021 submarine and subaerial DEMS; note that the elevation difference maps were not cleaned to show the distribution and the magnitude of the residual mismatching between the first and the last DEMs.

Figure 5 .
Figure 5. (a) Elevation difference map (cleaned for the residual mismatching between DEMs) between the July 2021 and February 2020 bathymetries draped over the 2021 DEM; ME and MA for marine erosion and accretion, respectively.(b) Elevation difference map for the subaerial (July 2021-April 2020) and submarine (July 2021-February 2020) SdF draped in semi-transparency over the July 2021 submarine and subaerial DEMS; note that the elevation difference maps were not cleaned to show the distribution and the magnitude of the residual mismatching between the first and the last DEMs.

Figure 6 .
Figure 6.(a) The 2020 shaded relief map and contours (equidistance 50 m) of the central part of the SdF, where a main depression is recognizable, partially filled by small ridges and fan-shaped features; NSDT: nearshore submarine depositional terrace and related average depth of the outer edge.(b) The 2021 shaded relief map and contours (equidistance as in (a)); black and red contour lines refer to the 2020 and 2021 DEMs, respectively.MA1 and ME1-2 as in Figure 5.

Figure 6 .
Figure 6.(a) The 2020 shaded relief map and contours (equidistance 50 m) of the central part of the SdF, where a main depression is recognizable, partially filled by small ridges and fan-shaped features; NSDT: nearshore submarine depositional terrace and related average depth of the outer edge.(b) The 2021 shaded relief map and contours (equidistance as in (a)); black and red contour lines refer to the 2020 and 2021 DEMs, respectively.MA1 and ME1-2 as in Figure 5.

Figure 7 .
Figure 7. RGB colour composite of amplitude CSK images: (a) 24 March 2021-18 December 2020; (b) 25 May 2021-24 March 2021; (c) 28 June 2021-25 May 2021.The analysis of the 30 July 2021 CSK-SAR image (Figure 8a) showed the alternation between areas with low and high amplitudes throughout the SdF slope.The high-amplitude areas (lighter tones) correspond to the following locations: (i) the summit crater terrace, (ii) the 2007-2014 lava flow field in the NE part of the SdF, (iii) the various lava flows emplaced along the SW part of the SdF during the last decades, and (iv) the central part of the SdF, where they match the slope accretion area recorded between June 2021-May 2021 and July 2021-June 2021 elevation difference maps (Figure4b,c).More interestingly, high CSK-SAR amplitude values were recorded all along the coastal strip of the SdF.This high-amplitude area matches a sandy/gravelly accumulation at the foot of the slope, with sparse sub-metric and metric blocks alternated with more coherent lava flows recognizable in the 7 July 2021 UAV optical image (Figure8d) and oblique photos (Figure2c).Lowamplitude values are generally found within topographic lows encompassed between the main lava flows and are commonly floored by relatively finer-grained volcaniclastic material, as recognizable by optical images collected both by the PLÉIADES-1 satellites and UAV survey.Regarding the submarine part, MBES/SSS backscatter data show the presence of medium-high backscatter values all along the nearshore sector, within the first 10 m wd (Figure8).This area is bounded downslope by a high backscatter stripe (Figure8d), roughly oriented parallel to the coastline, matching the marked break-in slope associated with the edge of the submarine terrace (NSDT in Figures6 and 8d).The flat terrace is characterized by overall smooth and medium-high backscatter values (Figure8c,d), with the local accumulation of blocks ranging from sub-metric to metric in diameters or small lobate lava flows.At greater depths, the quality of the SSS data increases with respect to the MBES backscatter, mainly because the former was acquired closer to the seafloor.On the other hand, the SSS data georeferencing was poorer and needed to be re-georeferenced with MBES bathymetry and backscatter through the matching of homologous morphologies.In

23 Figure 8 .
Figure 8.(a) 30 July 2021 CSK-SAR amplitude image of the subaerial SdF slope combined with the 27 July 2021 multibeam (MBES in (a)) and side scan sonar (SSS in (b)) backscatter data of the shallow submarine portion of the SdF; black lines are contours (equidistance 100 m), and blue dashed lines are the SSS track lines.(c) Zoom of MBES backscatter and CKS-SAR image; (d) Zoom of the SSS backscatter with the 7 July 2021 AUV optical image; ME1, ME2, and NSDT are as in in the previous figures.Note that SSS data ensonified seafloor areas more proximal to the coast with respect to MBES data due to the different acquisition geometry of the two systems.

Figure 8 .
Figure 8.(a) 30 July 2021 CSK-SAR amplitude image of the subaerial SdF slope combined with the 27 July 2021 multibeam (MBES in (a)) and side scan sonar (SSS in (b)) backscatter data of the shallow submarine portion of the SdF; black lines are contours (equidistance 100 m), and blue dashed lines are the SSS track lines.(c) Zoom of MBES backscatter and CKS-SAR image; (d) Zoom of the SSS backscatter with the 7 July 2021 AUV optical image; ME1, ME2, and NSDT are as in in the previous figures.Note that SSS data ensonified seafloor areas more proximal to the coast with respect to MBES data due to the different acquisition geometry of the two systems.

Table 2 .
Timetable of the topo-bathymetric surveys used for monitoring the morphological evolution of the SdF during the February 2020-July 2021 period; mwd: metre water depth.

Table 3 .
Information about the COSMO-SkyMed SAR images and their use in the morphological evolution analysis, with the date of related images.

Table 4 .
Results of the PLÉIADES and UAV survey DEM adjustments, with respect to the 2012 LIDAR DEM.RMSE0 and RMSE1 are root mean square errors computed between the 2020-2021 DEMS and the 2012 DEM before and after the co-registering procedure, respectively.

Table 5 .
Cell-size and σ ∆Z of the different elevation difference maps used for this study.

Table 6 .
Area, volume, average thickness, and associated σ ∆Z for the main topo-bathymetric changes recorded on the Sciara del Fuoco slope during the overall monitoring period, including the subaerial and the submarine analysis (July 2021-February 2020).