Mt. Etna Feeding System and Sliding Flank: A New 3D Image From Earthquakes Distribution in a Customisable GIS

High-resolution seismic imaging enables the reconstruction of ascending paths of magma and fluids, shallow molten accumulation and flank collapse areas, all crucial information for developing an efficient eruption forecasting strategy. Here, the Marching Cubes algorithm (MC - generally applied to medical visualization and three-dimensional (3D) modeling) is applied to 16 years of earthquake location data at Mt. Etna (Italy). The algorithm defines three-dimensional seismic clusters that take into account seismic location uncertainties and are embedded in a novel volcano-oriented Geographyc Information Systems (VolGIS) offering an interpretational environment comprising tomographic images and alternative geophysical models. The results show that a volume of very-low-seismicity is embedded in a high-velocity body, and acts as a zone of transition between transient magmatic events (west) and eastern deep seismicity related to the sliding eastern flank. The eastern cluster represents the 3D seismic signature of a deep (2–8 km below sea level) instability, affecting the portion of the eastern flank nearest to the feeding systems. This instability is likely caused by a combination of gravitational spreading and magmatic intrusions.


INTRODUCTION
Imaging a volcano with seismic waves is challenging but rewarding. A reliable image of the structural features of a volcano allows to reconstruct the regions where molten materials rise across the lithosphere (Jaxybulatov et al., 2014;Huang et al., 2015). These images can be linked with field observations, giving us a clear picture of ongoing volcano dynamics (Patanè et al., 2003). They can also discriminate zones of fluid and magma accumulation (Koulakov, 2013) as well as interfaces that release stress at depth , providing better constraints for seismic and volcanic hazard assessment. Regardless, the resolution of these maps is of the order of kilometres. Seismic locations, and particularly their "clustering," generally recorded across short periods of unrest, are often used to interpret tomographic anomalies, and from these the larger-scale dynamics leading to volcanism (e.g., Giampiccolo et al., 2020). It is thus necessary: 1) to develop more quantitative methods to define seismic clustering at a volcano; 2) to apply these methods to seismic locations recorded across years or decades, most likely to represent stationary or recurrent processes at a volcano.
Such an ideal dataset of seismic locations is available at Mt. Etna volcano. Mt. Etna is one of the most hazardous and monitored volcanoes in the world due to both its persistent eruptive activity and its proximity to highly-urbanized areas. The volcano is monitored by permanent and mobile networks and produces dense seismicity related to eruptive events and volcanic unrests (Alparone et al., 2015;Gruppo Analisi Dati Sismici, 2017). This seismicity has been used for decades to image shape, dimension and location of the volcano feeding systems with seismic tomography. The resolution of the resulting models, reaching a maximum depth of 20 km, has steadily increased (Sharp et al., 1980;Hirn et al., 1991;Cardaci et al., 1993;De Luca et al., 1997;Villaseñor et al., 1998). From the beginning of this century, researchers have interpreted the models jointly with remote sensing and field data, providing reliable interpretations of the shallow volcano dynamics (De Gori et al., 1999;Chiarabba et al., 2000;Laigle et al., 2000;Patanè et al., 2002;Patanè et al., 2003;Chiarabba et al., 2004;Patanè et al., 2006;Alparone et al., 2012;Díaz-Moreno et al., 2018;Giampiccolo et al., 2020). These studies have progressively recognized and improved the reconstruction of a high-velocity body (HVB) below both the central cone and the eastern sides of the volcano (De Gori et al., 1999;Alparone et al., 2012;Giampiccolo et al., 2020), generally interpreted as a solidified intrusive body Patanè et al., 2006;Alparone et al., 2012).
Unfortunately, seismic tomography models are affected by uneven resolution and instabilities caused by insufficient ray coverage (Rawlinson and Spakman, 2016). This makes it difficult, if not impossible, to model structures smaller than 1-2 km in an ever-changing volcanic environment (Koulakov, 2013). Due to the intrinsic limitations of the datasets used and the uncertainties related to ray-tracing algorithms in volcanic edifices, the highest resolution achieved by tomography on both the HVB and the nearby feeding systems is 1 km (De Gori et al., 1999;Alparone et al., 2012;Giampiccolo et al., 2020). The thickness of dikes is of the order of a few meters (Gudmundsson, 1983;Tibaldi, 2015), resolving them with travel time tomography is likely impossible. Also, there are strong uncertainties when detecting magma accumulation with tomography as demonstrated at Krafla volcano, where extensive geophysical imaging was available since the 90s (e.g., Julian et al., 1993). Based on these models, no large magma accumulation region was expected between depths of 3 and 7 km under the central caldera, at least until a deep geothermal well didn't drill into rhyolitic magma at these depths in 2009. Even novel tomographic imaging performed after drilling is unable to reconstruct magma at these depths at resolutions higher than 2 km 3 (Schuler et al., 2015).
Despite the dense seismicity and coverage offered by several seismic arrays, all recent tomographic works at Etna focus on specific time intervals related to volcanic eruptions (e.g., Patanè et al., 2002;Giampiccolo et al., 2020). In this way, the tomographic maps can be interpreted using geological and geophysical data produced by visible eruptive events. Regardless, such a dense seismicity should provide benchmark imaging models, necessary to understand the evolution of volcano dynamics across decades, including the relevant gravitational spreading observed across the wider eastern flank (Borgia et al., 1992;Urlaub et al., 2018). The shallow eastern flank is known to develop primarily aseismic creep (Rasà et al., 1996;Mattia et al., 2015;Bruno et al., 2017) and the wider flank is considered mostly aseismic. However, recent studies show that the deeper flank nearest to the feeding systems is associated with intense seismicity during unrest . At the edge of Mt. Etna's unstable sector intrusions encourage sliding during unrest (Alparone et al., 2020). Understanding if and how widely this seismicity has clustered across years and decades can thus help clarify the nature of sliding and its debated origin (Tibaldi and Groppelli, 2002;Urlaub et al., 2018). In particular, it can test numerical models proposing the existence of two sliding interfaces separating two domains: 1) a shallower flank, subject to gravitational instability and 2) a deeper flank, where compression is also caused by magmatic intrusions (Tibaldi and Groppelli, 2002;Apuani et al., 2013).
Forward methods that can properly account for data uncertainties can improve the resolution on, and interpretation of, seismic models. The "Marching Cubes" (MC -also known as "3D Contouring" o "Surface Reconstruction") is an algorithm for surface reconstruction that generates the 3D contour of a volume in space . It operates over a scalar field defined on a given volume. The algorithm approximates the data, which have to be related to physical quantities, to an isosurface, i.e., the surface of constant values in the scalar field. For example, the isosurface generated at the constant value greater than zero represents the enveloping surface around the volume where the scalar field is not null. The MC algorithm is used in medical imaging for 3D surface reconstructions of organs, tissues and anatomic parts and applied to datasets obtained by magnetic resonance tomography  as well as in pharmacology, chemistry, and meteorology (Newman and Yi, 2006). It is a fundamental tool in the framework of computational graphics and 3D modeling and has become an important technique for visual communication in the 3D animation and gaming industries. A review of the steady advances made by this approach can be found in Masala et al. (2013).
In seismology, the MC algorithm has been used for quick visualization of three-dimensional meshes and voxels. Typical examples are those proposed by Subramanian and Fussell (1990), who used it to make voxels for volume rendering in a ray tracing algorithm, and Ma and Rokne (2004), who developed a mesh propagation algorithm useful for the generation of seismic horizon surfaces. As of today, the MC algorithm has never been applied directly to seismic data and metadata, like earthquake locations, for the purpose of interpretation. With their uncertainties (hundreds of meters), the earthquake locations available from INGV across 16 years (Alparone et al., 2015;Gruppo Analisi Dati Sismici, 2017) can be used as the scalar field in the MC algorithm, at a resolution unavailable to tomography. The resulting isosurfaces will envelope volumes of seismic clustering and contour aseismic zones, whose discussion is always part of tomographic interpretations. The application of the MC algorithm to locations recorded across decades offers a unique method to define background (tectonic or gravitational) or repeated (likely magmatic) processes acting inside the volcano quantitatively.
To support interpretation of volcanic processes, any algorithm designed to interpret volcano dynamics must be implemented in an interactive visualization environment able to localize geological, geographical and geophysical data and models. Such an environment (VolGIS) has been developed to enable users to: 1) control visualization parameters in real-time; 2) locate the MC isosurfaces geographically, with precision sufficient to compare them with geophysical and geological data; 3) perform queries such as measure size and extension of the imaged structures (Guardo and De Siena, 2017). The integration of VolGIS and the MC algorithm is aimed at refining interpretations of volcano seismicity considering both imaged structures and dynamics modeled from geophysical signals. Once applied to seismic locations spanning years, it provides a tool to mitigate several limitations in our ability to interpret volcano dynamics, including: 1) the possibility of accounting for uncertainties in seismic locations when interpreting them with geophysical models; 2) a quantitative definition of "seismic clustering" and aseismicity in a geolocalized environment developed explicitly for the interpretation of seismic tomographic maps; 3) a seismic model of the volcano unrelated to specific volcanic unrest and representative of decadal-scale volcano dynamics.

DATA AND METHODS
In the present application, the MC algorithm uses as data the earthquakes nucleated between the years 2000 and 2016 at Mt. Etna, in the framework of the VolGIS. The first ten years of data are described and analyzed by Alparone et al. (2015). We added six years of location data (from 2010 to 2016), which were provided by the Gruppo Analisi Dati Sismici (2017). These earthquakes are located in an area that spans between 14.707 and 15.295°E and from 37.509 to 37.900°N, in a depth range from −3 to 10 km b.s.l. We selected the "Date," "Latitude," "Longitude," "Depth," "Magnitude," "RMS," "ERH" and "ERZ" from the event location file. The last two fields, ERH and ERZ, indicate the horizontal and vertical uncertainties in localization (in kilometres), respectively.

Geographical Information System Setup
After setting the data-frame coordinate system, both a DEM of Etna (Bisson et al., 2016) and all the seismic events were loaded to the GIS-Workspace. Without any selection dependent on location uncertainties, both the top view (at different scales, i.e., 1:75,000, 1:125,000) ( Figures 1A,B) and W-E cross-section view ( Figure 1C) of the earthquake spatial distribution highlight a volume of low-to-zero seismicity for the selected period. The methods proposed allow to assess if this anomaly is an effect of graphic visualization or a feature related to the volcano dynamics. VolGIS incorporates the Marching Cubes algorithm (MC) and analyses the earthquake distribution in order to create a quantitative 3D model of the high-seismicity volumes. Geolocalized tomographic maps and models of deformation are then included in the GIS for interpretation.

Marching Cubes
The MC samples a point cloud set by taking eight vertices simultaneously from a provided spatial grid, marking a cube. For each cube, the algorithm computes 2 8 256 possible polygons configurations values that, due to symmetries, are reduced to 15 possible polygons configurations (Supplementary Figure 1A). The MC then queries a list of pre-calculated geometries to pick the combination of polygons that best represents the isosurface passing through the cube, given a predefined grid spacing and density isovalue set by VolGIS. Once the user inputs these parameters, the software applies the algorithm and generates the isosurface according to the spherical spreading law: Our point cloud-set is the dataset of earthquake locations. After defining a grid, the contribution of each earthquake location (j) in the cloud to the vertices (i) of the grid is thus computed using the inverse square of the distance dependency (d(i, j) 2 ) between each location and the vertex. The result is a scalar field defined over the grid vertices by adding up all the inverse distances j. From Eq. 1 we obtain a value W(i) for each vertex of the cube (Supplementary Figure 1B); the value decreases when the number of locations near to a vertex increases.
For a given dataset, the isovalue establishes what percentage of earthquake locations near to each vertex is required for the vertex to be inside the surface. The i-th vertex is included if W(i) has a value less than or equal to the isovalue, otherwise it falls outside of the surface (Supplementary Figure 1C). Another way to achieve the same effect is to uniformly increase or decrease all the vertices values (Lopes and Brodlie, 2003).

RESULTS
The MC algorithm solves a forward problem. It is thus necessary to test different datasets in order to assess the robustness of the imaged structures. The MC was applied to three different datasets ( Figure 2; Supplementary Figure 2, 3), created according to three different horizontal and vertical uncertainties (ERH and ERZ). These applications test the robustness of the reconstructed structures, providing a control on visualization based on modeled uncertainties. We also test the stability of structure locations and shapes applying a bootstrap test.

Dataset 1: 6,907 Earthquakes
The first dataset comprises 6,907 earthquakes with a root-meansquared (RMS) average value of 0.135. We took into consideration that the average ERH and ERZ are equal to 548 and 712 m, respectively. By dividing the entire area (51,779,297 m for each side) by the ERZ (712 m), we obtained the maximum amount of cells that can be included in the analysis (72.7). The software rounded down (72) to respect the ERH-ERZ constrain. Then the area was separated into 72 cells, obtaining a value equal to 719,15 m per each cell-side.
The MC algorithm produced three isosurfaces, or clusters, that envelope the high seismicity volumes and constrain a volume of low (almost absent) seismicity (Supplementary Figure 2A): • The first two clusters, located beneath the summit craters (C1) and south to it (C2), have a vertical extent of 2.79 km (from 0.24 km above sea level to 2.56 km below it) and 3.6 km (from 0.56 to 4.16 km b.s.l.), respectively. • The third cluster (C3) is located under the eastern sector of the volcano. It has an approximated thickness of 5 km, dipping eastward from its shallow point at 1.4-7.39 km b.s.l. (Supplementary Figure 3B). • The algorithm allows the imaging of a fourth cluster, located between the second and third one. However, its dimension is smaller than the earthquake vertical uncertainty, so it is not interpreted in this first analysis.

Dataset 2: 2,824 Earthquakes
Here, we set the grid spacing to 505.45 m, since the used events have both ERH and ERZ equal or less than 500 m. In this case the total earthquakes and the average RMS are equal to 2,824 and 0.149, respectively. The results show the same three clusters obtained using the first dataset, together with three smaller ones of measurable dimensions (Figure 2A). The clusters  obtained with this dataset have a smaller extension compared to those obtained using Dataset 1, but the variations are within the dataset uncertainty. Both shape and position of the clusters remain the same. This result is due to the lower uncertainties in source location.
• The cluster C1 has a vertical extent of 3.28 km (from 550 m a.s.l. to 2.73 km b.s.l.). The second cluster (C2) extends for 3.05 km, spanning from 0.12 to 3.17 km b.s.l. • Cluster C3 keeps both the same dimensions and position of the previous analysis. • Located between C2 and C3, cluster C4 is now measurable, with a focal point at about 1.19 km b.s.l. and a vertical extent of 1.67 km. • Two additional clusters are located above (C5) and below (C6) the second one (C2). They have both an average vertical extent of about 1 km and they range from 1.89 to 0.71 km a.s.l. and from 3.52 to 4.52 km b.s.l., respectively ( Figure 2B).

Dataset 3: 1,937 Earthquakes
In the third and last analysis we selected earthquakes with both ERH and ERZ equal to or less than 400 m, for a total of 1,937 events. The RMS average value is again 0.149 and the grid spacing is set to 402.73 m. The clusters imaged at this resolution are five in total (Supplementary Figure 3A) -with C6 (defined in the previous analysis - Figure 2) having dimensions smaller than the uncertainties. In this case, one of the clusters (C1) has a vertical extent higher than in the previous cases (3.72 km compared to 3.28 km for Dataset 2). This variation is higher than the uncertainty and changes the cluster shape. The other four clusters keep their average thickness and position, but their shape changes due to the limited number of earthquakes taken into account (Supplementary Figure 3B).

Robustness of the Isosurfaces Relative to the Amount of Data
Dataset 2 (ERH and ERZ ≤ 500 m) shows the highest numbers of clusters and will be used for both the bootstrap test, which is independent of the underlying spatial distribution of the data, and further discussions. The isosurfaces are defined robust in the sense that they are visually continuous as the data change in value (Lopes and Brodlie, 2003). We subtracted 10% (2,542 earthquakes), 20% (2,259) and 40% (1,694) from Dataset 2 randomly and the procedure was then repeated 10 times for each percentage value. The MC algorithm reconstructs isosurfaces almost identical to those generated using the original dataset (Supplementary Figure 4,5), confirming that Dataset 2 provides a robust spatial distribution of isosurfaces with reductions up to 40%.

DISCUSSIONS
The MC algorithm applied to the seismicity of Mt. Etna reveals three-dimensional high-resolution seismically-dense structures, which visually and spatially constrain a volume with low seismicity (VLS). A comparison of the retrieved structures with studies related to the 2000-2016 activity of the volcano (Bonaccorso et al., 2002;Bonforte et al., 2008;Currenti et al., 2008;Alparone et al., 2012;Bonforte et al., 2013;Carbone et al., 2014;Bruno et al., 2017): 1) shows that the VLS is enclosed in the high-velocity body (HVB) depicted by previous tomographic models; 2) allows to characterize the VLS as a portion of the feeding system; and 3) images the persistent seismicity associated to the sliding of the combination of gravitational and magmatic forces affecting the deeper eastern flank, closest to the feeding system.

Spatial Relation with Seismic and Deformation Models
The comparison between both the tomographic maps and geodetic models with the seismic clusters allows to define a preliminary seismic zonation map of the volcano without separating the seismicity in different time periods. It is hence possible to characterize the clusters as produced by either gravitational or magmatic processes, such as dikes' intrusion.
The HVB thus consistently comprises all clusters and the VLS except for C3 (red contour). The VLS is thus a smaller-scale feature inside the high-velocity pattern that separates the main western clusters C1 and C2 from C3, under the eastern part of the volcano. The obvious questions are if the clusters are caused by different dynamics (magmatic, tectonic and gravitational) and if the VLS is part of the feeding systems of the volcano.
We also compare clusters C1 and C2 with the location of the dike and crack models obtained from geodetic and GPS analyses performed between 2000 and 2016 (Figure 3).
• The locations and main orientations of clusters C1 and C2 match the crack projection modeled by Bonaccorso et al. (2002) relative to the lateral eruption of July 2001 ( Figure 3A). The source of ground deformation modeled eruptions, respectively . Both models are located less than 1 km west of the VLS -the first below cluster C1, the second on top of it ( Figure 3B). • The model proposed by Bruno et al. (2017), relative to the May 2015-September 2016 period, is located in the first cluster, in a range that spans between 2 and 6 km b.s.l., about 2 km west of the VLS ( Figure 3C). • The deformation sources modeled using continuous GPS data (Cannata et al., 2018) are related to several eruptive events. They highlight a more complex ground deformation patterns then those obtained by measurements performed Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 589925 during a specific eruptive event (e.g., compare the deformation sources from  ( Figure 3B) and (Palano et al., 2017) (Figure 3D), for the period July 2004-July 2005). However all these ground deformations are primarily located inside or around cluster C1 ( Figure 3D). • Clusters C1 and C2 are adjacent to the fracture-weaknesszone (FWZ) modeled by Carbone et al. (2014). This zone, located just 1.6 km SW from the VLS, is interpreted as a deeper pressure source filled with highly-pressurized gases ( Figure 3E). • The surface dike projection modeled by Pezzo et al. (2020) comprises most of C1, C2 and most of the VLS ( Figure 3F). • The spatial relation between these models and the clusters west of the VLS (C1 and C2) shows that the corresponding seismicity is produced by transient magmatic activity, marking the well known "South-rift" Pezzo et al., 2020). The weak zone modeled by Pezzo et al. (2020), considered as the superficial evidence of a deeper dike intrusion, is the only deformation anomaly crossing the VLS consistently.

The Volume with Low Seismicity as Portion of the Feeding System
The VLS is in an extensional regime between the eastern and western flanks of the volcano Alparone et al., 2011;Carbone et al., 2014;Murray et al., 2018). A comparison between the VLS and different Coulomb stress distribution models (Currenti et al., 2008) (Supplementary Figure 7) shows that the VLS falls inside a wide volume, extending to a depth of 6-7 km b.s.l., characterized by positive Coulomb stress changes (values equal and higher than 0.8 MPa). Here failure on a fault is favored (Cocco and Rice, 2002). If the VLS falls in a fragile HVB, why doesn't it show seismicity?
In active volcanoes earthquake-free zones are typically interpreted as the result of magma storage (Decker, 1984;Scandone and Malone, 1985;Parfitt and Wilson, 2008). At Mt. Etna such a seismic gap was previously recognized and considered as a magma supply in the shallow HVB (Chiarabba et al., 2000). If the origin of the HVB are cooled magmatic rocks they can thus pass from a brittle fracture behavior to a plastic deformation one after dikes' intrusion (Parisio et al., 2019). Such a dike intrusion has been modeled for the 2018 eruption (Aloisi et al., 2020;Pezzo et al., 2020) inside the VLS. Rocks near these dikes intrusions can deform in a ductile manner (aseismically) while the cooler rocks far from the dike will fail in brittle mode (Parisio et al., 2019). This justifies the drastic seismicity decrease FIGURE 4 | W-E cross-section highlighting (in gray) the area with lowseismicity above and below the VLS (in dashed blue-red).
FIGURE 5 | A W-E cross-section view of Cluster C3 (red outline). The previously modeled detachment planes are represented in yellow (Bonforte and Puglisi, 2006) and light blue .

Cluster East of theVolume with Low Seismicity: Seismicity Underlying Flank Instability
Cluster C3, located east of the VLS is oriented W-E and has a thickness of 5 km and a length of about 7.4 km ( Figure 5). The top surface of C3 matches the element known in literature as "décollement" or detaching plane Tibaldi and Groppelli, 2002;Bonforte and Puglisi, 2003;Lundgren et al., 2003, Lundgren et al., 2004Walter et al., 2005;Bonforte and Puglisi, 2006;Puglisi et al., 2008;Aloisi et al., 2011;Palano, 2016). This plane is located at about 3 km depth, dipping from 1.4 km b.s.l. to almost 4 km b.s.l., from top to bottom, respectively. It is in overall agreement with previous models Bonforte and Puglisi, 2003;Bonforte and Puglisi, 2006;Puglisi et al., 2008;Azzaro et al., 2013) (Figure 5) despite being obtained with seismic locations distribution instead of deformation data. It is thus possible to consider the top of the 3D surface of this cluster as the boundary zone between rocks with different mechanical properties (Chiarabba et al., 2000;Laigle et al., 2000;Alparone et al., 2011;Palano, 2016), on which the superficial mass of the volcano flank slides due to gravitational instability Tibaldi and Groppelli, 2002;Bonforte and Puglisi, 2003;Lundgren et al., 2003, Lundgren et al., 2004Walter et al., 2005;Apuani et al., 2013) (Figure 6). On the other hand the bottom surface of cluster C3 dips from 4 to 7.39 km b.s.l. This surface fits with a magma intrusion-induced slip surface hypothesized by previous studies (Borgia et al., 1992, Borgia et al., 2000Tibaldi and Groppelli, 2002;Lundgren et al., 2004;Neri et al., 2004;Walter et al., 2005;Apuani et al., 2013). There is a general agreement that the wider eastern flank, especially the portion in the coastal area, is sliding primarily due to gravitational forces (Urlaub et al., 2018). However it is unlikely that a 5 km-thick seismic structure (C3) is produced by gravitational collapse only. If our interpretation of the VLS is correct, a portion of the western seismicity inside C3, with volume equivalent to C1 and C2, could be equally caused by dikes' intrusions. It is thus possible to explain the seismicity across the eastern flank using a simple conceptual model of two flanks: a shallower and a deeper one (Tibaldi and Groppelli, 2002) ( Figure 6). The shallower one, above the cluster C3, is mostly affected by gravitational collapse and is sliding along the C3 top surface. The deeper flank is constrained by cluster C3 and slides along its bottom surface. The corresponding seismicity is due to a combination of gravitational and magmatic forces, most likely from an intermediate storage zone like that modeled by Aloisi et al. (2020) between 4 and 7 km under the central cone. The VLS represents the shallowest evidence of magmatic intrusions (Neri et al., 2018;Aloisi et al., 2020) contributing to the sliding of the deeper eastern flank.
In our interpretation C3 thus represents the first 3D seismic image of the deeper portion of the eastern flank, in contact with the feeding systems and sliding due to both magmatic and gravitational forces. The seismicity inside the cluster has already been associated with magmatic or hydrothermal Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 589925 activity during the 2018 unrest : we prove that this seismicity is a consistent feature before this event, likely related to gravitational sliding (at the top of the cluster) and repeated magmatic intrusions. The existence of a ductile VLS where dikes can intrude thus better agrees with a reciprocal feedback between: 1) gravitational collapse encouraging the uprizing of magma and 2) dikes' intrusions favoring sliding in the deeper part of this portion of the eastern flank (Tibaldi and Groppelli, 2002;Bonforte and Puglisi, 2003;Aloisi et al., 2011;Apuani et al., 2013;Bruno et al., 2017). Previous authors have underlined that "eruptions do not trigger catastrophic flank collapses" analyzing displacement of aseismic shallow faults related to the flank kinematics (Urlaub et al., 2018). Our results show that deep magmatic activity produces persistent (at least decadal-long) seismic signatures that strongly differ from those at surface. Deeper magmatic activity or flank eruption (Alparone et al., 2020) are thus a likely trigger of shallower sliding and, consequently, catastrophic collapse of the shallower portion of the flank. For a complete understanding of how catastrophic this pattern would be, it is necessary to consider the relation between local tectonic and eruptive events (Bonforte et al., 2019;De Novellis et al., 2019;Aloisi et al., 2020;Pezzo et al., 2020), i.e., conducting a time dependent analysis.

CONCLUSIONS
The earthquakes nucleated between 2000 and 2016 at Mt. Etna have been used as data in the Marching Cube algorithm, a 3D modeling and visualization analysis. The MC has been framed in an experimental GIS (VolGIS), which allows the users to obtain a three-dimensional seismic image in a few seconds and with a high quality graphics output. The system offers unique support in the imaging of relatively small-sized volumes within tomographic and deformation models resolved over kilometres. The framework allows to establish benchmark seismic zonation-mapping at volcanoes including uncertainties from seismic locations. It is specifically designed for volcanoes that have been mapped with seismic tomography and where extensive seismic networks and highrate seismicity are available. These are today standards at many volcanoes worldwide, like Kilauea, La Réunion and those of the Canary Island, all characterized by high-rate seismicity and flank instability.
The results at Etna constrain a volume of low (almost-absent) seismicity, persistent across a 16 years period. We interpret this volume as a portion of the feeding system of the volcano embedded in a high velocity body. In this interpretation, the VLS represents the zone of transition between seismicity caused only by transient dike-related activity (west of the VLS) and seismicity caused by mixed magmatic and gravitational components. The model provides the first 3D seismic evidence of the different dynamics affecting the deeper and shallower portion of flank closest to the feeding systems, with important implications on the estimation of the volumes affected by sliding and the forces causing it. Only time-dependent analyses will provide exact spatial correlations of the VLS with source deformation models, insight into the temporal instability of the eastern flank and whether the VLS is a stationary feature or can be reconstructed only after a major eruption of the volcano.

DATA AVAILABILITY STATEMENT
The used software (VolGIS) provides graphic outputs that coincide with those shown in the article and in the Supplementary Material. The earthquake location data can be requested to the INGV-Osservatorio Etneo (Gruppo Analisi Dati Sismici, 2017).

AUTHOR CONTRIBUTIONS
RG designed the software and directed its development, built the 3D models, generated the maps and wrote the first draft text. LS guided the overarching research direction and advised on the geophysics and volcanology. CD coordinated the technical work on the software. All authors contributed to revising the text.

ACKNOWLEDGMENTS
This work was made possible by the outstanding efforts of the staff at the INGV-Osservatorio Etneo. We are particularly grateful to A. Ursino, S. Alparone, and G. Currenti, who provided us the seismic dataset and helped us to correctly interpret the variations of Coulomb stress at Etna. We thank Jurgen Neuberg for suggesting the comparison with the Coulomb stress and Guido Ventura for the insightful discussions on the volcanological significance of our results. We also thank the computer scientists of the LVCC for the support in the software development and A. Colubri for suggesting the application of the Marching Cube algorithm.