Characterising vent and crater shape changes at Stromboli: implications for risk areas

Active volcanoes are typically subject to frequent substantial topographic changes as well as variable eruption intensity, style and/or directionality. Gravitational instabilities and local accumulation of pyroclasts affect conditions at the active vents, through which gas-particle jets are released. In turn, the vent geometry strongly impacts the eruption characteristics. Here, we compare five high-resolution topographic data sets (<4 cm/pixel) of volcanic craters and vents from Stromboli volcano, Italy, that were acquired by unoccupied aerial vehicle (UAV) during five field campaigns between May 2019 and January 2020. This period includes two paroxysmal explosions (3 July and 28 August 2019) and exhibited significant changes on day-to-month timescales. Our results highlight changes to vent geometry and their strong control on the directionality of explosions. Recurrent UAV surveys enable the monitoring of temporal morphologic changes and aid the interpretation of observed changes in eruption style. Ultimately, this may contribute to repeatedly revised risk areas on permanently active volcanoes, especially those that are important tourist destinations.


Introduction
Vent evolution is a critical parameter for volcanic hazard assessment as shifts of vent geometry and position can be linked to shifts in eruptive mechanisms [Graettinger et al. 2015;Taddeucci et al. 2013a;Valade et al. 2016]. The geometrical evolution of craters can be correlated with processes of crater formation to enhance the understanding of active volcanic processes [Hanagan et al. 2020]. Direct observations detecting changes in the activity at persistently active volcanoes can provide insights into the shallow conduit system, which, in turn, also improves hazard assessment [Capponi et al. 2016;Salvatore et al. 2018;Simons et al. 2020]. Not only can the geometry be affected by the eruptive activity but vent geometry can also modulate the dynamics of volcanic explosions.
While vent geometry can be measured directly, our knowledge about conduit geometry has been constrained based on inactive fissures, eroded volcanoes, laboratory experiments or through indirect geophysical methods [e.g. Chouet et al. 2003;Keating et al. 2008;Parcheta et al. 2016;Zorn et al. 2020a]. Some morphological features are unlikely to be a direct proxy of the uppermost plumbing system before eruption/explosion.
For example, excavated craters/conduits following major Vulcanian/Plinian Corresponding author: markus.schmid@min.uni-muenchen.de explosions have probably widened due to the explosion intensity [Wilson et al. 1980] (giving a propensity for lithic components in pyroclastic deposits) and gravitational instabilities [Calvari et al. 2006], and drained lava lakes frequently exhibit funnel-like geometries [Patrick et al. 2019] that are probably the result of convection-driven thermal erosion. In contrast, spine growth (especially the cross sectional shape and extrusion direction) has been used to infer ascent velocities, magma properties and conduit geometry [e.g. Lacroix 1904;Vallance et al. 2008;Zorn et al. 2020a]. The effect of inclined conduits on eruption dynamics was first demonstrated experimentally by [James et al. 2004] who showed that gas slugs (constant starting volume and pressure differential) rising in inclined conduits will be less overpressured at burst compared to vertical conduits. Gas overpressure at the vent has been demonstrated experimentally to affect gas-only and gas-particle jet dynamics [Alatorre-Ibargüengoitia et al. 2011;Alatorre-Ibargüengoitia et al. 2010;Schmid et al. 2020]. Lagmay et al. [1999] and Major et al. [2013] linked asymmetric crater geometry to inclined eruption columns and, as a consequence, to a preferential distribution of proximal volcanic hazards. Taddeucci et al. [2013a] illustrated how the presence of a crater may change the dynamics of eruptive jets.
The link between vent geometry and jet dynamics has been investigated experimentally through rapid de-Vent and crater shape changes at Stromboli  Figure : Morphological variations of the crater terrace of Stromboli volcano, Italy, as seen from Pizzo between and . Variations are due to recurring subsidence/collapse events and accumulation of pyroclastic material and lava flows.
[A]- [E] show the entire crater terrace while [F] is a zoom to a spatter cone. [E] and [F] show the development of the spatter cone within a mere two days. Eruptive processes led to a geometry change of the spatter cone, visibly affecting the directionality of eruptive jets.

Images [A]-[D] by Ulrich Kueppers, [E] and [F] courtesy of Angelo Cristaudo.
compression experiments where impulsive jets have been released from a vertical shock-tube setup. For instance, Cigala et al. [2017] explored the influence of four different radially symmetrical vent geometries on gas-particle jets and constrained the effect of vent geometry on residual overpressure at the vent, a parameter that contributes to jet expansion and particle dispersion. In similar experiments with the same setup, Schmid et al. [2020] focussed on the characteristics of gas-only jets released from vent geometries with reduced symmetry. In these experiments, six vent geometries were fabricated using two designs with the strongest impact (cylindrical, 15°diverging inner geometry of Cigala et al. [2017] combined with a slanted surface plane (5, 15, 30°). These experiments confirmed that the bilateral vent symmetry is a major controlling parameter for the expansion dynamics of gas-only jets, leading to jet inclination and asymmetric spreading angles despite a vertical conduit.
Understanding vent evolution and migration, as well as crater and/or vent asymmetry and their links to eruption dynamics and mechanisms, may improve hazard assessment. Tourist destination volcanoes (e.g. Stromboli, Italy; Villarica, Chile; Whakaari/White Island, New Zealand; Yasur, Vanuatu) are famous for the accessibility of observational points but infamous for unheralded strongly explosive events [e.g. Dempsey et al. 2020;Giordano and De Astis 2021;Viccaro et al. 2021]. Depending on eruption frequency, recurrent up-close surveys of the active crater area can reveal high(er) temporal resolution information on changes. Eruptive activity is known to vary on several scales, including -but not limited to -eruption frequency and height, erupted volume, grain size distribution, pyroclast temperature and jet directionality [e.g. Andronico et al. 2009;Harris and Ripepe 2007;Taddeucci et al. 2013b;Zanon et al. 2009]. Changes in eruption characteristics may be influenced by observable topographic changes of active vents (shape, size, depth, open/closed, rim height) [e.g. Capponi et al. 2016;Cole et al. 2015;Jessop et al. 2016;Solovitz et al. 2014;Suzuki et al. 2020].
Over the past 30 years, documentation of active craters has developed from increasingly complex sketches and photographs [e.g. Andronico et al. 2013;Calvari et al. 2014;Harris and Ripepe 2007] to aerial imaging and remote sensing data [James et al. 2020b;Turner et al. 2017;Zorn et al. 2020b]. Because of their versatility, UAVs have become a powerful tool for many geoscientific fields [Eltner et al. 2016;Niedzielski 2018]. Their ability to acquire high-resolution imagery, conduct measurements or collect samples in hazardous or inaccessible areas makes them an ideal tool for the volcanological community (see James et al. [2020b] for a comprehensive review of volcanological applications). Repeated observations of active vents by UAV Volcanica 4(1): 87 -105. doi: 1 .3 9 9/vol. 4. 1.871 5

Figure :
Strong morphological variations of N and N (vents and of the N vent area) of Stromboli volcano, Italy, due to recurring subsidence/collapse events and accumulation of pyroclastic material as seen by UAVs between 6 and . Arrows indicate North.
( Figures 1 and 2) and production of digital elevation models (DEMs) allows for precise, quantitative comparison of temporal changes. At persistently active vol-canoes such as Stromboli, UAVs provide an excellent opportunity to perform measurements several times a year and achieve a statistically robust morphological dataset over long periods, at high temporal and spatial scales. Turner et al. [2017] used UAVs to map active and inactive vents at Stromboli in May 2016. Here, we use UAVs to repeatedly acquire imagery for structure from motion (SfM) reconstruction in order to characterise and quantify changes of Stromboli's crater terrace, including the geometry of craters and vents at unprecedented spatiotemporal scales. We also use our photogrammetric models to quantify position, size and asymmetry of active craters and vents at Stromboli. In combination with observations of volcanic explosions we link the geometric variability of vents and craters to eruptive behaviour and the affected areas.

Activity and morphology at Stromboli volcano
Stromboli volcano, Italy, is perhaps best known for its continuous eruptive activity for the past 2000-2500 years [Rosi et al. 2000]. In the recent decades, Stromboli's activity has been characterised by mild, persistent explosive activity every few to tens of minutes ejecting ash, lapilli, and incandescent bombs up to heights of a few hundreds of meters above the vent [e.g. Andronico et al. 2013;Bertagnini et al. 1999;Patrick et al. 2007;Taddeucci et al. 2013b]. This 'normal' Strombolian activity is periodically interrupted by two types of more energetic explosions, known as 'major explosions' and 'paroxysms' as well as the eruption of lavas [Barberi 1993;Ripepe et al. 2008]. The largest events in the past century ejected metre-sized bombs and blocks as far as inhabited areas located ca. 2 km away [e.g. Calvari et al. 2011;Rittmann 1931;Rosi et al. 2000]. Such activity, albeit rare (most recently on 3 July 2019, 28 August 2019 and 19 July 2020), poses various hazards including pyroclastic density currents, ballistic impact, respiratory problems and vegetation fire [Brown et al. 2017]. The paroxysm on 3 July 2019 emitted an eruption column~5-8.4 km high and a pyroclastic flow that travelled down the Sciara del Fuoco [Giordano and De Astis 2021;Giudicepietro et al. 2020]. This eruption marked the start of a 2-month-long effusive phase. The 28 August 2019 paroxysm (~6.4 km high) was similar in style to the 3 July 2019 and was also accompanied by a pyroclastic flow and a lava flow [Giordano and De Astis 2021;Giudicepietro et al. 2020].
As a result of different eruptive styles, magnitude and frequency, the number and positions of vents as well as their geometry is affected [Calvari et al. 2014]. The active vents at Stromboli are located within the crater terrace, a break in slope at the top of the Sciara del Fuoco at about 800 m above sea level (asl), lying below the common lookout point Pizzo or Sopra la Fossa at around 918 m asl. The morphology of Stromboli's crater terrace has long been of great interest, as evi-denced by scientific descriptions, illustrations and photographs throughout the centuries. Washington [1917] reviewed 21 publications between 1768 and 1915 with a focus on the persistence of Stromboli's active vents. For the past several decades, the crater terrace has hosted three main vent areas: north-east (N), southwest (S) and the central (C), with S and C areas being often grouped together [Salvatore et al. 2018]. We use 'vent' as a term describing the opening in the ground from where gas and pyroclast jets are ejected, e.g. N1 for the vent number 1 in the north-east vent area, S2 for the vent number 2 in the south vent area . If there is a 'crater', i.e. a negative, subcircular volcanic landform around a vent, they are named after the associated vent (see Figure 3A and Figure 4A for an overview).
Recent publications [e.g. Gaudin et al. 2017;Harris and Ripepe 2007], as well as recurrent observations (at least once per year) of the crater terrace by several coauthors since 2005 have revealed significant variations in eruption style and frequency, and vent number (between 3 and 15, see Figure 1). The shallow plumbing system below the crater terrace has been investigated with tilt meters [Bonaccorso 1998], seismic networks [Chouet et al. 2003] and continuous GPS [Mattia 2004]. These studies suggest the presence of a NE-SW trending structural weak zone that coincides with the direction of dykes exposed in the edifice. Zanon et al. [2009] analysed explosions from Instituto Nazionale di Geofisica e Vulcanologia's (INGV) monitoring webcams and observed inclined jets which they linked to conduit geometry or conduit inclination. Calvari et al. [2014] proposed that morphology changes between 2002 and 2007, together with a massive collapse of the summit crater, have modified the shallow feeder conduit, leading to changes in the eruptive style (i.e. increasing the number of major explosive events and lava overflows).

Field campaigns
We contribute five aerial data sets of Stromboli's crater terrace over nine months (May 2019 and January 2020). Due to the fact that two paroxysms (on 3 July and 28 August 2019) occurred during this period [Giudicepietro et al. 2020], these data provide an opportunity to investigate both the 'normal' activity of Stromboli and also changes where deeper portions of the shallow plumbing systems were affected. During each campaign, several vents were active, with variable styles of 'normal' activity occurring (Figures 1, 2 and 8).
Field campaigns were conducted over 11-16 May 2019, 5-13 June 2019, 4-5 August 2019, 21-26 September 2019 and 25-27 January 2020. During this time, Stromboli volcano was erupting frequently, with several explosions occurring during the UAV mapping flights, resulting in gas or ash plumes as well as pyro-clastic ejecta up to 150 m above the vents. Therefore, all flights had to be conducted manually without predefined flight paths in order to be able to react quickly to prevent loss of, or damage to the UAVs. We mitigated systematic errors in our topographic models by following workflow and best-practice suggestions [e.g. Eltner et al. 2016;James et al. 2020a;James et al. 2019].
The flights were conducted at heights between 50 and 150 m above the main area of interest, with a double grid flight path and nadir to off-nadir camera angles. These low flight altitudes were chosen to accomplish a ground sampling distance (GSD) of a few centimetres. The August 2019 flights have a ground resolution of 5.8 cm/pix, the other flights (May, June, and September 2019 and January 2020) have an average resolution of 4.2 cm/pix (Table A1). During individual campaigns, up to 20 flights were performed over several days to ensure full coverage of the entire area of interest, under good light conditions and minimal obscuration by the degassing plume. Here, we present DEMs generated from images of individual flights with sunny to overcast sky and variable wind speed (Table A2). The camera was set to shutter priority with high shutter speeds (1/240-1/500 s). Due to the main focus of this study being the crater terrace and active vents, placing ground control points (GCPs) was not possible. Hence, the global navigation satellite system (GNSS) camera position information was used for georeferencing. In June we had the opportunity to fly on four out of five days (between 8 June and 12 June 2019) and focused on the short-term development of the N area. Four out of the five survey flights were performed by LMU staff (focused directly on the crater terrace and the geometry of vents and craters) and one flight by INGV Rome staff (August 2019). Additionally, we recorded UAV and ground-based video footage and imagery of the eruptive phenomenon to establish a direct link between vent geometry and the resulting eruption dynamics.

Hardware and software
We used two UAVs from DJI ® : Phantom 4Pro+ and Mavic 2Pro. The Phantom's camera has a 1 inch CMOS sensor and 8.8 mm focal length (equal to 24 mm as 35 mm equivalent) with a maximum resolution of 5472 × 3648 pixels and a mechanical shutter. The Mavic's camera uses the same size sensor but with a 10.2 mm (28 mm as 35 mm equivalent) focal length and an electronic shutter.
Different software packages were used for 3D reconstruction of the acquired aerial imagery and the analysis of the obtained models. The structure from motion (SfM) algorithm of Agisoft Metashape (Version 1.5.1 -1.6.1) was used to match image features to make a coarse 3D reconstruction of the surface. By comparing matching features across several images, the 3D position of the cameras can be calculated. Building on this, the multi-view stereo (MVS) algorithm of Agisoft Volcanica 4(1): 87 -105. doi: 1 .3 9 9/vol. 4. 1.871 5 Metashape uses the coarse cloud and the obtained camera parameters to perform the reconstruction of a dense cloud. The open source software CloudCompare (Version 2.10.2) and QGIS (Version 3.10) were used to perform cloud-to-cloud and cloud-to-mesh comparisons as well as DEM and orthomosaic comparisons, respectively.

Processing
Suitable images covering the area of interest were selected and imported into Metashape to check image quality (Agisoft Metashape Image Quality) and remove blurred images. As cut-off criteria, we used thresholds of 0.8 (May, June, August and September 2019) and 0.6 (January 2020). DJI ® UAVs are known to have accuracy problems in the flight height information stored in the image metadata that are beyond those of usual GPS inaccuracies. Therefore, a correction was applied to adjust for the in-flight vertical sensor offset . Areas with strong degassing and areas covered by the gas plume were masked to prevent artefacts during 3D reconstruction. The images were aligned to produce the sparse cloud, roughly representing the topography of the survey area. To improve image alignment, the camera model was optimised by including focal length (f ), affinity (b1), the centre of distortion (cx, cy) and both radial (k) and tangential distortions (p) of the lens within the bundle adjustment. The GNSS camera positions were included as control observations during the bundle adjustment with uncertainty estimates of 10 m in the three Cartesian directions.
Before further processing, the sparse cloud was cleaned by applying several filter criteria to remove points with weak geometry, large pixel matching errors and large pixel residual errors. The threshold for the reconstruction uncertainty was set to 15, points with a higher uncertainty were removed. The level for the projection accuracy was set to 3. The desired threshold for the reprojection error was 0.3 pixels. To reach this level, the threshold was set in a way that a maximum of 10 % of the total points was removed in every iteration until 0.3 was reached. Between iterations, the optimization (parameters: f , k1, k2, k3, cx, cy, p1, p2, b1) of the camera alignment was repeated, further decreasing the reprojection error to below 0.323 pix. The optimised sparse cloud was the basis for the creation of the dense cloud (high quality, mild depth filtering). The dense clouds were filtered by point confidence and points with a confidence level below 1 were removed. Where necessary, the dense point cloud was improved by manually deleting artefacts. From the dense cloud, all other products were calculated, e.g. meshed 3D models, tiled 3D models, DEMs and orthomosaics. The created DEMs have an average resolution of 8.4 cm/pix https://github.com/agisoft-llc/metashape-scripts/ blob/master/src/read_altitude_from_DJI_meta.py and https: //github.com/agisoft-llc/metashape-scripts/blob/master/src/ add_altitude_to_reference.py and a maximum resolution of 7.6 cm/pix (Table A2). To identify and locate vents, incandescence and fumaroles, orthomosaics were created by projecting the aerial images onto the 3D surface of the model or the DEM. As a result of this, the spatial resolution was increased to an average of 3.7 cm/pix, allowing a better recognition of ground features than from the DEMs alone.

Analysis
Visibility was best in June 2019. Accordingly, it was used as a reference model where prominent features were identified for the referencing of the other models (see Figure 3). Due to a lack of clear imagery, and because of significant topographic changes as a result of the two paroxysms, the September 2019 model was aligned to the already referenced August 2019 model. We manually picked individual points e.g. prominent rocks or pinnacles of rock within outcrops of bedrock or an exposed dyke ( Figure 3A). We assume that these features were stationary throughout the survey period but their appearance may have changed because of erosion, rockfall and pyroclastic deposition ( Figure 3B-E). As a result, the individual reference points may vary for each survey. In general, four to seven points were used as reference without additional check points. The number of suitable available reference points was limited because the large extent of the areas affected by the two paroxysms had not been anticipated during survey design.
CloudCompare was used to reference the models to each other and to create the rectified DEMs. The referencing yields root mean square (RMS) errors of 0.39 m (May-June), 0.76 m (June-August), 0.39 m (August-September) and 0.33 m (June-January). For the intersurvey comparison of this study, this sub-metre relative accuracy is sufficient and accurate global positioning was not required. Most of the analysis was carried out with QGIS, where it was possible to attain and compare surface elevation, slope angles, crater rim heights and crater and/or vent geometry. For all the craters and vents identified we measured area, aspect ratio, circumference, height of crater rims, difference between highest and lowest point of the crater rim (Δh rim ) and the orientation of the highest and lowest point around the crater. Based on the difference between the highest and lowest sector of the crater, we calculated the slant angle of the crater exit plane. QGIS was also used to address the evolution of crater and/or vent positions as well as the vent surface cover (open versus debris covered). We used the additional UAV video footage to identify active covered vents, and inactive vents during the survey period that were not visible on the DEMs or orthomosaics.  Figure .

Limitations
The aim of this study was to obtain the highest possible resolution of near-vent (within 100 m) topographic changes. Because of the proximity to the active vents, no ground control points (GCPs) could be placed within the area of interest. In an ongoing collaboration, our observations will be coupled with a study of topographic changes in a larger area above 700 m asl. Without reliable GCPs and check points, propagation of sys-Volcanica 4(1): 87 -105. doi: 1 .3 9 9/vol. 4. 1.871 5 tematic error and artefacts is difficult to quantify. Although our DEM comparisons are useful for providing a first order impression of uncertainties, we cannot assess systematic positioning or reconstruction errors.

Results
Our results show that the morphology of Stromboli's crater terrace and the geometry of craters and vents are transient on timescales of days to months.
3.1 Topographic changes of the entire crater terrace at the time scale of months The comparison between May and June 2019 and between September 2019 and January 2020 show the changes caused by 'normal' and elevated levels of activity (see Table A3 for levels of activity), while changes between June and August 2019 and between August and September 2019 were dominated by the two paroxysms and the effusive episode after the 3 July paroxysm.
During each campaign, we observed the conditions of crater topography (shape), fumaroles (strong or low degassing, 'hot' or 'cold') and active vents (open or closed) and any changes over day-to-month timescales. This allowed us to constrain quantitatively the impact of both constructive and destructive processes. Between May and June 2019 (~32 days) surface elevation changes ranged between´15 and +8 m in the survey area. In particular, negative height variations were restricted to the crater floors, while elevation gain occurred within small portions of the craters as well as in their surroundings. S1 was excavated by explosions, retrograde erosion and possibly subsidence. Gravitational instabilities affected the southern sector of S2 and led to elevation changes of up to´15 m (Figure 4). The paroxysm on 3 July 2019 removed at least 30 m (vertical) in the N and SC vent areas. The elevation changes between June and August (~53 days) were most apparent in the western portion of the SC area (between´10 and +16 m). On 23 September 2019, 26 days after the second paroxysm, up to 20 m of elevation was lost within S2 and the C vent area while elevation loss was´12 m in the N area. In January 2020,~124 days after the previous campaign (September 2019), only positive elevation changes were detected with the strongest increases west of S2 and around N1 and N2. A maximum of +32 m was gained in the SC area, and +20 m in the N area. Elevation loss related to the explosive excavation by the 3 July and 28 August 2019 paroxysmal events was partially masked by lava effusion and pyroclast deposition.

Geometric changes of vents at the time scale of days
The prevailing processes shaping the N crater area on a timescale of 4 days in June were the deposition of pyroclasts, retrograde erosion along scarps, growth of two circular features within N1 and possibly subsidence ( Figure 5). Pyroclast accumulation predominantly occurred within N1, while material loss from explosive excavation and/or subsidence dominated the crater floor of N2. Retrograde erosion was limited to a south-western section within N1 as well as the western and south-eastern portions of N2. The larger circular feature appeared to be a bank around the main vent within N1 outlined by erosion and/or the beginning of cone growth. On 1 July 2019 at 05.40 am local time, a hornito inside S2 crater showed energetic Strombolian explosions with abundant pyroclast and subordinate ash ejection. These explosions were emitting sub-vertical jets with a symmetrical dispersal of pyroclastic material. Two days later, on the morning of 3 July 2019 at 06.26 am local time, the hornito had lost several meters of its height and had developed a small notch in the rim (Figure 1E, F). Intermittent activity was at a similar level on both days [A. Cristaudo, pers. comm.] yet the emitted jets and pyroclasts were directd towards the notch (Figure 1E, F).

Vent and Crater positions and crater terrace morphology
The SC and the N areas were persistent sites of volcanic activity throughout the study period, although the activity at some vents ceased, then sometimes resumed at the same or at close-by positions. The spatiotemporal evolution of the vents and craters are depicted in Figure 6.

N vent area
From May to June 2019 the diameter of N1 was reduced and, as a consequence, its shape changed, while N2 remained unchanged except for the elevation changes of the crater floor. In the beginning of August 2019, the morphology of the crater terrace was completely altered. Instead of two pronounced craters (N1 and N2), a total of 12 new vents emerged, and six of them exhibited incandescence. These were arrayed along a curved line 13 m and 23 m north of the previous centres of N2 and N1 craters ( Figure 6A). By September, the activity in the N area was focused on N1 and N2 at their new locations 8 (N2) and 27 m (N1) from their pre-paroxysm location. Additionally, a new fumarole was active between N1 and N2 at a location where two vents were situated in August 2019. In January 2020, a larger vent (~3 × 4 m) formed a cone with a new crater and two smaller vents a) at the rim of  the crater and b) below the cone. A small (~0.5 m diameter) vent was active between N1 and N2 at a location where a fumarole was active in September 2019. The location of N2 did not change significantly from September 2019 to January 2020.

SC vent area
The crater diameter of S2 increased from May to June 2019 while the C vent remained unaltered. In early August 2019 the location of S1 shifted 13 m to the north and a new elongated vent (S3) became active on the western side of the crater terrace, around 46 m west of the position of S2 in June ( Figure 6). By September 2019, S1 had shifted another 11 m to the west and S2 formed an elongated crater system together with two C craters. C1 crater was at the same location as it was in May and June, while a new fourth crater was visible 28 m north of the centre of the S2+C crater system. S3 changed from an~30-m-long fissure to a circular crater with an~18 m diameter. In January 2020, three vents were visible in the SC vent area (S1, S2, and C1). S1 was at the same location as in September, S2 had a similar crater size as in May but was at a new location. C1 was isolated from the S2 crater, close to C1's location in May, while C2 was not visible in January 2020 anymore. The circumference of individual vents and craters was between and 46 and 227 m, with aspect ratios between 0.23 and 0.96. The height of the crater rim around the crater varied for each vent, in some cases considerably. N1 crater in June showed the biggest difference in Δh rim . The NW side was 21 m lower than its highest point in the NE. With around 4 m in September, S3 exhibited the lowest difference in crater rim height. In combination with the crater diameter, for each crater we calculated the slant angle of the theoretical surface plane of the crater rim. These range from 10°(N2 January 2020) to 39°(S1 September 2019) and had an average slant angle of~16°calculated from n "22 measurements.
We tracked the changes of N1, N2, S1, and S2 along four transects through the crater terrace (see transects in Figure 4A). The changes of N1 are illustrated in Figure 7A, where the edifice around the crater was growing over time, accompanied by crater enlargement between May and June 2019 (circumferences: 163 and 168 m). After the paroxysm on 3 July 2019, the location of N1 shifted towards the north and the circumference of the newly built crater was considerably smaller (131 m in September 2019 and 60 m in January 2020). Even though vent location and crater shape was modified by the events in July and August 2019, the southern crater rim of N1 was always higher than the northern crater rim. The difference between the highest and lowest point of the crater rim (Δh rim ) was 21 m in June 2019, changing to Δh rim of 6 m in January 2020 (Table A1). N2 ( Figure 7B) showed a similar evolution (circumference: 150 m in May and June 2019, 67 m in Septem- The S1 transects of May and June 2019 show change from cone to crater that was caused by a strong explo- Table : Geometrical parameters of N , N , S , S , and S craters from May to January detailing area (A), circumference (C), long axis (a), short axis (c), the vertical difference between highest and lowest point of the crater rim (Δh rim ), the azimuth of the lowest point of the crater rim (AZM min ) and the dip angle of the theoretical surface plane between highest and lowest point of the crater rim (θ surface ).

Date
Name A  sion on 15 May 2019. The transects through S1 also cut through the C crater ( Figure 7C) and shows that the location of the C crater was stationary. Figure 7D shows transects through S2 (in May, June, September 2019 and January 2020). Strong elevation changes can be seen, including the removal of the crater rim (visible in May and June 2019 transects) and subsequent growth by pyroclastic accumulation (August 2019 -January 2020) at a new location.
We compared our models (2019 and 2020) and photographs (2013-2019) with 1) sketches from 1994-2004 , 2) photographs from 2007-2012 [Calvari et al. 2014] and 3) a model from 2016 [Turner et al. 2017]. These observations represent only snapshots in time and therefore interpreting a trend for small scale features development was not possible. However, interpolation between repeat surveys in a larger context might be useful. For example, from 1994 to 1997, the number of active vents in the N portion of the crater terrace decreased from 7 vents (5 hornitos, 2 with craters) to 2 craters with 1 vent each (N1 and N2). Also, from May 2000 to May 2001 the activity within C shifted as a hornito built up to the north of its pre-2000 location, where activity ceased in 2002. During the same timespan, up to 4 vents were active within a locally stable crater in S.
In 2004, probably as a result of the paroxysmal episode in 2003, the locations of N and S shifted eastward.
In 2007 the entire crater terrace was a deep depression affecting all vent areas. Between 2008 and early 2011, the crater in S was increasing in size and a hornito in C transformed into a crater. This was occurring alongside the refilling of the crater terrace by pyroclastic deposits. With continued deposition, N migrated towards the scarp of the Sciara del Fuoco.
In September 2011, only~2 vents in N and one vent in C were visible. This changed by September 2012 when there was a large crater in S formed by explosive excavation and an increased number of active vents in C. From 2013 to 2017, the changes were dominated by the enlargement and deepening of the crater in S and both craters in N probably due to a combination of explosive excavation and retrograde erosion of the crater rims. In May 2016 no vent was active in C, but activity resumed in 2017 as part of a larger S crater complex.

Morphological changes of the crater terrace
We have used our data to evaluate the persistence of craters and vents at Stromboli's crater terrace. We found that the prevailing classifications of vents into two-to-three main centres of activity was applicable throughout the timespan of our observations. The position of craters and vents seem to be structurally controlled and we suggest that the S2 and C vents are aligned along a NE-SW trending feature that appears to be parallel to an old dyke west of the crater terrace (see Figure 3A). The paroxysm of 3 July 2019 excavated around 30 m of material, with explosions seemingly being produced by S2 and N2. This value is a conservative minimum, because lava effusion and explosive activity built up the crater terrace during 32 days between the paroxysm and data collection on 4 August 2019. The two paroxysms in 2019 have destroyed the uppermost (few tens of meters) plumbing system and replaced it by a zone of variably sized and cohesive fall-back material as already described by Calvari et al. [2014] for the period between 2007 and 2012. This heterogeneous zone provided additional pathways for magma to reach the surface and enabled the occurrence of the multitude of active vents in the N vent area visible in August 2019. However, it appears that the deeper plumbing system beneath the uppermost zone has not changed substantially, because the activity was soon re-focused into fewer active vents inside N1 and N2 craters. We suggest that the northward migration of N1 is surface morphology-controlled rather than being due to subsurface structure. Since a large amount of material was removed in the area that now hosts N1 and the estab-lished conduit was partially destroyed, the new uppermost conduit of N1 migrated approximately 35 m north of its original location during the refilling of the crater terrace. We suggest that changes to the southern wall of N1 and to the western wall of N2 represent retrograde erosion of parts of the steep, and possibly overhanging crater walls.
During the last 26 years the crater terrace has been subjected to opposing constructive and destructive processes. It seems that changing vent locations and openings of new vents promoted hornito growth. These hornitos were destroyed by explosions and/or collapse and the resulting craters were subsequently enlarged by explosive activity and erosion until a larger event eventually overprints the morphology.

Link between vent geometry and eruption dynamics
We observed changes to individual vents on timescales of days. The changes of vent geometry were accompa- nied by changes in style and direction of explosions. For example, until partial destruction of S1 in May 2019, the explosions were dominated by long-lasting (up to 50 s), sub-vertical gas-rich jets with incandescent pyroclasts. After this event, explosions were heavily laden with brownish ash, probably related to debris coverage of the vent due to backfall from explosions from other vents and erosion of the walls as described before by Capponi et al. [2016] and Simons et al. [2020]. The explosions and ballistic trajectories were directed towards the newly formed low side of S1 which may be relevant for risk assessment/hazard management (Figure 8). The partial destruction of the cone of S1 in May 2019 exposed the uppermost 12 m of a formerly cylindrical sub-vertical conduit. This confirms the impact of asymmetrical exit geometry on jet characteristics released from a vertical conduit.
The modification of S2 in early July 2019, shown in Figure 1E, F, is another example where short-term changes of the vent geometry led to directed explosions. This change occurred over few days and supports the assumption that the new direction of the explosions was a result of a modified exit geometry of the vent ( Figure 1E, F). Previous authors have linked inclined explosions at Stromboli to conduit inclination [Zanon et al. 2009]. We suggest that, for the two cases presented here (S1, May 2019; S2, July 2019), conduit incli-nation played no role in the change of direction of the explosions and that the inclined explosions were solely based on the asymmetric crater and/or vent geometry. Further evidence of this is that the timescales over which the explosion behaviour changed (days) were too short to be a consequence of tilting (due to inflation or gravitational creep) of the crater terrace and the underlying plumbing system. The experiments of Schmid et al. [2020] confirmed that bilateral vent symmetry, i.e. one side higher than the other (Δh rim ), produced inclined jets and asymmetric spreading angles despite a vertical conduit. These gas-only experiemts provided a link to the underlaying pure fluid dynamical proceses related to complex vent geometries. The pyroclastrich events with variable-sized ejecta at Stromboli add complexities related to the coupling between particles and the gas phase, but the underlying principles remain valid.
Vents covered by debris may also impact the directionality of jets released from explosive eruptions. These processes have been documented by experimental studies of explosions of known energy, explosion depth and covering lithology by Graettinger et al. [2015] and Taddeucci et al. [2013a]. At Stromboli, both open and covered vents have been observed.

Towards a morphological monitoring of persistently active volcanoes
Spatiotemporal topographic data provides unique information that should be incorporated into multiparametric volcano monitoring, both to increase our understanding of eruptive processes and to better assess and mitigate specific hazards. Tourist-destination volcanoes are of prime importance for the local economy [Erfurt-Cooper et al. 2015]. Several such volcanoes have decade-long observations and there is some understanding of the frequency and magnitude of explosive eruptions that put the local population and tourists at risk [Bertolaso et al. 2009;Erfurt-Cooper et al. 2015]. If the frequency of such events-to-be-avoided is low, permanent access bans may be difficult to enforce. The probability of ballistic impacts or inundation by pyroclastic density currents are two ways to define exclusion zones on active volcanoes [e.g. Alatorre-Ibargüengoitia et al. 2012;Alatorre-Ibargüengoitia et al. 2016;Lavigne et al. 2017;Romero et al. 2017;Toyos et al. 2007], frequently defining large portions of the volcano's flanks as off limits. On volcanoes where 1) continuous activity is believed to lead to limited risk of above-average magnitude explosions, 2) this number is small enough to be considered tolerable in the year-long average by the local civil protection and 3) with high agricultural and/or tourism pressure, topographic analysis of active vents at high spatial and temporal resolution may contribute to define risk areas or exclusion zones of substantially smaller extent and at potentially more frequent revision intervals. It is up to the local authorities if the resources for monitoring, observation and interpretation outweigh the economic interests of parts of the population to allow for such a symbiosis of human activity in potentially varying parts of a volcano. It goes without saying that access to agricultural land or tourist viewpoints is subordinate to public safety.

Closing remarks
Features of volcanic vents and craters exert a prime control on explosive volcanic activity. The data presented here show the development of the craters and vents at Stromboli volcano in unprecedented detail. The high temporal resolution allowed for a distinction between the effects of 'normal' activity and the elevated activity during the two paroxysms. In addition to qualitative description, it was possible to quantify geometry and morphology changes due to the high resolution of DEMs calculated from images taken during UAV surveys. Such surveys can be accomplished at short repeat intervals and at tolerable risk exposure for the pilot and their observer, which is important since the geometries can change on short timescales. High temporal and spatial resolution may allow quantification of as-yet unconstrained eruption parameters e.g. erupted mass, to date only indirectly-and crudely-known via measur-ing the degassing behaviour. At present, the limiting factor is the comparatively high error from the alignment of the models that could be improved by adding high quality GCPs or UAVs with real-time kinematic or post-processed kinematik (RTK/PPK) capabilities. Long-term observations showed that crater and vent geometry can be stable or transient over periods of weeks to months ('normal' activity) but strongly altered during the two paroxysms in July and August 2019. Nevertheless, deeper portions of the shallow plumbing systems were seemingly unaffected as successive activity soon focused back to the three centres of activity as before the paroxysms. Moreover, 'normal' eruptive activity with predominantly near-vent deposition of erupted material commonly rebuilt volcanic landforms resembling the pre-paroxysm configurations within weeks to few months.
Additionally, we observed the paramount impact of crater and vent geometry on pyroclast ejection characteristics, a fact that has strong implications for areas potentially affected by bomb impact and fall of pyroclasts. UAV photogrammetry can acquire unbiased data sets that have a much higher information content and comparison potential than photographs or sketches. DEMs enable high-quality measurements of predominant geometric features. Laboratory experiments showed that crater and vent geometry influence the directionality of volcanic jets. This parallels observations of changed eruptive behaviour and directionality following fairly sudden (hours to few days) geometric changes. Repeated UAV surveys can be used to evaluate risk areas, with low risk for the operators and with standard computational capabilities. A quantitative comparison of vent and crater geometry as well as crater terrace morphology was not possible due to the different media throughout the years. As Stromboli is frequently visited by scientists and UAVs are available in many working groups, large collective timeseries at high temporal and spatial resolution are achievable. Therefore, we made our DEMs, orthomosaics and processing reports publicly available (GFZ Data Repository: https: //doi.org/1 .588 /fidgeo.2 21. 15).