The improvement of the ground deformation monitoring system from 2001 on (ref. 21) has allowed to perform detailed modelling of acting pressure sources. Many authors (e.g. refs 22, 23, 24, 25) clearly delineated a pathway of the magma rising towards the surface, which borders the western side of the high seismic velocity body underneath the summit craters, visible in several tomographic studies (e.g. ref. 26). Along this pathway, the magma can accumulate and differentiate at various levels. In particular, ref. 24 proposed the presence of inflation sources located in a volume at a depth between about 5 and 8 km b.s.l (see inset in Fig. 2).
Regarding the ground deformation pattern affecting the volcanic edifice, we focus on distance measurements between benchmarks located on the western flank of the volcano, in order to identify the inflation/deflation phases. We prefer the baselines on the western flank to the baselines deployed on the eastern one as the former are much less affected by major faulting and eruptive fracturing (see e.g. ref. 23). In fact, the behaviour of the western flank is typically considered as a clear indicator of the volcano state related to the processes ongoing in the plumbing system (e.g. refs 27, 28). By default, we focused on the baseline EMGL-EMAL; in addition, also the baseline EMEG-ESLN was considered (see inset in Fig. 2). During times when both baselines could be measured, we noticed that their trend over time was almost identical, given their similar orientation and location. We therefore used EMEG-ESLN when EMLG-EMAL was unavailable.
The ground deformation curve (Fig. 2) shows a general trend of inflation over the 20 years. On occasion of major eruptions this trend is interrupted by rather short-lived periods of deflation; for instance, we note about 1 year of deflation in 2005 and also from 2013 to 2015 (ref. 24). Rapid inflation occurred from 2002 until 2004, on occasion of the strong 2002/03 flank eruption, until the 2004 eruption (ref. 29). Inflation accelerated from 2015 until 2021. This phase accounts for almost 50% of the inflation encountered in the years 2002–2022. In 2021, we witnessed a major deflation, linked to the paroxysmal activity with about 60 lava fountains (ref. 25).
Seismicity: Overall Picture
Figure 3 shows the distribution of earthquake foci recorded in the area of Mt. Etna since 2002. The events, about 19,000 with magnitude mostly between 1 and 2.5 (Mmax = 4.8), were located using the “tomoDDPS” software (ref. 30). Compared to standard location techniques, tomoDDPS improves the accuracy of hypocentre location through a combination of absolute and differential arrival-time readings between couples of closed-spaced earthquakes. Besides, the code allows computing the seismic ray-tracing in a 3D velocity model; here we used a model derived by the integration of passive and active seismic data (ref. 31).
The seismicity on Etna is rather heterogeneously distributed. As a general trend, hypocentre depth increases from E, ESE towards NW. Shallow foci, with a depth of less than 5 km, are frequent in the central part of the volcano, as well as on its eastern and southern flanks. A similar distribution is found for the earthquakes with depths of about 5–10 km. Foci with depth between 10–15 km are frequent in the SW part of the area, but also to the N and NNE of the summit craters, while deeper events (depth in the range of 20–35 km), are located in the NW part of the volcanic edifice. A fairly clear boundary between domains with shallow and deeper foci can be outlined; it crosses the volcano from SW to NE, intersecting the summit area.
The distribution of the hypocentres reflects the interaction of local effects due to the dynamics of the volcano, and the regional tectonic framework of the area. In this picture, the large amount of seismicity in the central area point to the mid-shallow plumbing system, while the progressive increase of focal depth towards NW can be understood as an expression of the crustal convergence that takes place beneath the thin volcanic cover; i.e., under the volcano we encounter the NW subduction of the Hyblaean foreland crust below the fold and the thrust system of the SE verging Sicilian chain (see refs 14, 32, 33). Likewise, the large number of hypocentres of the eastern sector, especially in-depth levels below 5 km, can be related to the activation of the numerous faults characterising that flank, supposedly linked to a regional transtensional fault system extending offshore (refs 15, 34).
Seismogenic Volumes And Strain Rate
We quantified the seismic activity over time in calculating the cumulative curves of seismic strain within some key volumes of the volcanic area. The seismic strain curves were obtained using Richter’s relation for the seismic energy (E in erg), i.e., log(E) = 9.9 + 1.9ML − 0.024ML2 (ref. 35).
For the definition of the key volumes, we started from the depth distribution of the foci. In the distribution of focal depth (inset in Fig. 3), we can see that about 5,000 events are rather shallow, less than or equal to 3 km b.s.l. Further changes of the curvature in the distribution can be noted around 7–8 km and 17–18 km. We therefore decided to distinguish four main depth ranges, i.e. < 3 km, 3–7 km, 7–17 km and > 17 km. As we shall see later, this choice is appropriate also for the events for which fault plane solutions were calculated. To identify the horizontal boundaries of the key volumes we mapped the seismic energy released annually in the volcanic area. Specifically, we defined a three-dimensional grid, with a horizontal step of about 1 km and vertical spacing corresponding to the depth limits defined above. We calculated the sum of the energy released by the earthquakes located within the grid cells. Focusing on volumes with a high release of seismic energy, we finally identified 11 volumes. For each of them we calculated the seismic strain release curves and compared them to the ground deformation curve mentioned earlier. (Fig. 4).
In the shallowest layer we distinguish four key volumes; they show major release of seismic strain during the flank eruptions characterized by intense seismic swarms, i.e. the eruptions of 2002/03 and 2018. This is especially true for the V2 curve which refers to events localized in the eastern sector, while the curve V3, encompassing the summit craters area, shows a strong increase even during the 2008 eruption. V4, which refers to the NE Rift and the Pernicana fault sector, i.e. the northern edge of the volcanic block subject to eastward sliding, reveals a certain degree of activity even in non-eruptive periods.
In the layer between 3 and 7 km we have distinguished 3 key volumes. Volume V5 encompasses foci in the SW part of the area, V6 lies on the eastern flank of the volcano, while V7 includes the summit craters and adjacent areas to the west. The similarity of seismic strain release and ground deformation curves is strongest for the volume V6, i.e. the events on the eastern flank. Moderate increments of the V7 are detected before the eruptions of 2002/03 and 2018. On the other hand, seismic strain release of V5 is strongly affected only by the December 2018 eruption.
In the layer between 7 and 17 km we consider the volumes V8 with foci located SW of the summit, and V9, which is situated E and SE of the central craters. Whereas the seismic strain release in V9 is rather smooth, several steps can be observed in the curve for V8 which, however, do not show a particular link to the volcanic activity.
In the basement with depth larger than 17 km we identified two volumes, V10 situated SW and V11 NW of the summit area. V11 is among the seismically most active volumes. Its curve shows strong steps just before some eruptive activities (i.e. 2006, 2011 and 2013), however the same happens in inter-eruptive periods; therefore, a clear relation is lacking to the ground deformation curve or volcanic activity. V10 shows a lower degree of seismic activity, again poorly linked to ground deformation or volcanic activity.
We defined a measure of similarity between the normalized curves of ground deformation and strain release for the 11 key volumes. For this purpose, we calculated the differences between the normalized curves and integrated the differential curves, which we can call “accumulated divergence”. It should be small when the two curves are similar – such as for the pair deformation/seismic strain V6 – or higher for strongly differing curves, such as can be seen for the pair deformation/seismic strain V5. Taking the absolute maxima of the accumulated divergences for all 11 key volumes, we obtain the values reported in Table 2.
Table 2
Similarity between the ground deformation and strain release curves for 11 key volumes.
V1 | V2 | V3 | V4 | V5 | V6 | V7 | V8 | V9 | V10 | V11 |
0.286 | 0.146 | 0.060 | 0.160 | 0.32 | 0.013 | 0.208 | 0.108 | 0.136 | 0.253 | 0.11 |
From this criterion the ground deformation is most similar to seismic strain release of V6, i.e. the volume to the east, situated at depth between 3 and 7 km. The second most similar one is the volume 3, which encompasses the shallow level of the summit area.
Focal Mechanism Classification And Clustering
The radiation pattern of seismic energy released during an earthquake provides useful information concerning the forces acting in the seismic source and the orientation of the rupture. Using Anderson's classical scheme for pure shear fracture (ref. 36), we can distinguish three principal types of faulting: “Normal”, “Reverse”, and “Horizontal Strike-Slip”. A more comprehensive classification scheme is the one proposed by ref. 37, which includes mechanisms not attributable uniquely to any of the referred classes. These classifications do not include any information about the strike of the P- and T-axes; this can be limiting, for example in studies such as this one where we investigate the orientation of the stress field controlling the tectonic processes of an area.
Following a strategy proposed by ref. 38 (see also ref. 39), we used unsupervised pattern recognition. It is based on an a-priori defined metric of similarity between patterns rather than exploiting a-priori defined target information furnished by an expert (see Supplementary Information for more details). We applied these unsupervised learning techniques to moment tensors derived from fault plane solutions of Mt. Etna earthquakes. We focused on 577 focal solutions of events which were recorded in the time span 2002–2021 and had magnitudes ranging from 2.1 to 4.8. The mechanisms were obtained by calculating the ray-tracing in a 3D velocity model (ref. 31), using information regarding locations, arrival times and polarities taken from some datasets and published studies (refs 16, 17, 33, 40, 41, 42, 43).
In order to divide the focal mechanisms into homogeneous groups, we used an updated version of the software KKAnalysis (refs 44, 45), which offers different types of unsupervised classifiers. As a first step, the dataset was processed by applying the so-called Kohonen Maps, also known as “Self-Organizing Maps” (SOM, see Supplementary Information). SOM enable reducing the dimensionality of the feature vectors, allowing an efficient representation of pattern characteristics in space and time. Similar to ref. 39, we adopted moment tensor components as features. They can be handled with centroid-based clustering, allowing the definition of prototypes, i.e., objects which represent a number of patterns with reasonable fidelity. The algorithm grouped all the 577 focal solutions in 114 nodes, each of which can be understood as micro-clusters, whose centroid vectors form prototypes. We can substitute the original feature vectors by the prototypes which means a smoothing of the data; however, we are still dealing with a large number of groups (i.e. 114). Thus, in a second step, we performed the K-means clustering by exploiting the results obtained from the SOM. Specifically, the original moment tensors of the dataset were replaced by the ones of the prototypes, being representative of the original items. This significantly reduced the scatter of the initial data; indeed, instead of having 577 different moment tensors the variability reduces to only 114 versions of feature vectors. For K-means we can use simple rules for establishing the adequate number of clusters, namely the “Davis-Bouldin” index (see Supplementary Information). Using this index 8 clusters are preferred and we deal with a still smaller number of representative feature vectors – i.e., eight instead of the 114 SOM-nodes or the 577 original ones. The overall characteristics of each group were summarized by identifying the prototypes, i.e. patterns that are similar to all members of a cluster. This was done using the so-called Kagan angles (ref. 46). The metric is based on the angle by which a system P, T, B must be rotated so that the axes fall in the same position as P’, T’, B’ (see Supplementary Information). The best representation is found by taking the cluster members one by one and calculating the sum of Kagan angles with respect to all other members. The representative sample is the one which has the lowest sum of Kagan angles (see ref. 38 for further details).
In Fig. 5, we present the distribution of the focal mechanisms in the area of Mt. Etna. The colouring of the beach balls represents their membership to the clusters; the prototype beach balls of the clusters are shown in the upper right corner of each panel. Clusters 1, 2, 3, 4 and 6 correspond essentially to horizontal strike slip mechanisms. They are distinguished by differences in the direction of the P-axes, which is N for members of cluster 1, NNW for cluster 4, E and NE for clusters 2, 3 and 6. Cluster 5 consists of clear normal faulting events, even though the strike directions of the P- and T-axes show a considerable scatter. Clusters 7 and 8 are mainly of reverse type; P axes are directed NW-SE for C7, while they spread from ENE-WSW to ESE-WNW for C8.
The distribution of the clusters as a function of time and depth is shown in Fig. 6. Of note is the concentration of certain groups in distinct depth levels. In particular, the clusters C1, C5 and C8 are primarily representative of shallow events (median depth less than 4 km). They occur in the summit area or along the main fault systems of the volcano, namely the southern and north-eastern rift zones (close to the central craters) and the Pernicana fault (NE of the summit). Clusters C2 and C6 are located mainly in the central-eastern sector, in the depth range 3–7 km. The elements of the cluster C3 are located NE or SW of the summit area, at depths from very shallow to about 15 km. The foci belonging to clusters C4 and C7 are found in the NW sector, at depths over 17 km; furthermore, they tend to cluster in time, forming swarms, such as in June 2006 (before the 2006 eruption), in December 2009 or in May 2011 (before the lava fountain activity). The foci of the clusters 5 and 8 often coincide with periods where strong volcanic activity is observed. During the eruptions 2002/03 and December 2018, the two event types occurred at very similar times. This is a surprise at the first moment: cluster 5 are mainly normal faulting mechanism, whereas the mechanism of cluster 8 corresponds to horizontal strike slip/reverse faulting. This is probably an effect of the complex stress redistribution during the process of dike intrusion and flank movement.
For the interpretation of the seismo-tectonic stress field we plotted the orientation of the P- and T-axes of each focal mechanism, distinguishing the four main depth ranges mentioned above (see Fig. 7). For the shallow group we find strongly varying directions of P-axes, in particular in the summit area. For many of these events, however, the fault plane solutions correspond to normal faulting mechanisms with almost vertical P-axes. The horizontal orientation of the T-axes reveals a more stable picture for these events. The pattern of the T-axes is consistent with a SE and ESE movement of the SE flank of the volcano. Such movements are also revealed by INSAR imagery (ref. 47).
The events populating the layer in the depth range 3–7 km are found mainly in the eastern and southern part of the area. Again, we note a strong variation of the P-axes, with a prevalent strike direction in E and ESE in the eastern part, S in the southern part and again E in the western part of the area. On the whole the directions of the P-axes are oriented in a more or less radial geometry, pointing to a centre close to the summit area of the volcano. The T-axes make up concentric rings, again around a centre close to the summit.
In the layer between 7–17 km depth the events are concentrated in the eastern and south-western quadrants. The orientation of both P- and T-axes confirm the picture outlined above, i.e. P-axes pointing to a centre close to the summit, and T-axes delineating concentric rings around this centre.
Events deeper than 17 km are found mainly in the north-western sector of the area, with P-axes striking NW. Some events located in the western side have P-waves striking E and NE. The strike direction of the T axes is essentially perpendicular to the orientation of the P-axes. On the whole, the picture created by the orientation of the axes matches the ones found at 3–7 km and 7–17 km. At the same time, the orientations found for the events located in the NW sectors coincide with the regional stress field encountered for the central and northern part of Sicily (refs 43, 48).