Frequency-magnitude distribution of earthquakes at Etna volcano unravels critical stress changes along magma pathways

The high seismic productivity of volcanic areas provides the chance to investigate the local stress conditions with great resolution, by analysing the slope of the frequency-magnitude distribution of earthquakes, namely the b-value. Here we investigated the seismicity of Mt. Etna between 2005 and 2019, focusing on one of the largest known episodes of unrest in December 2018, when most of the intruding magma aborted, rather oddly, its ascent inside the volcano. We found a possible stress concentration zone along magma pathways, which may have inhibited the occurrence of a larger eruption. If the origin of such hypothetical loaded region is related to tectonic forces, one must consider the possibility that geodynamic processes can locally result in such rapid crustal strain as to perturb the release of magma. Strong b-value time-variations occurred a few days before the unrest event, suggesting new possibilities for investigating the volcano state and impending eruptions. Stress concentration zones along magma pathways may locally vary the slope of the frequency-magnitude distribution of earthquakes and temporarily inhibit the occurrence of large eruptions, according to an analysis of the seismicity of Mount Etna between 2005 and 2019.

M easuring changes in crustal stress associated with the ascent of magmatic bodies may be crucial in understanding a volcano's state and could make the assessment of incoming volcanic eruptions more effective. On 24 December 2018, Mt. Etna (eastern Sicily, Italy) experienced one of the strongest known unrest episodes, related to the intrusion of about 30 × 10 6 m 3 of magma 1 inside the upper crust below the volcano, as evidenced by different volcanological and geophysical information [2][3][4][5][6] . Nevertheless, this major unrest event resulted in a relatively modest, short-lived eruption.
In investigating such anomalous behaviour of the magma supply system, we tested using a crustal stress-meter, i.e., the slope b of the Gutenberg & Richter Frequency-Magnitude Distribution 7 (FMD; logN = a−bM, where N is the cumulative number of seismic events with magnitude above or equal to M and a represents the productivity), otherwise acknowledged as the b-value. It describes the average size distribution of earthquakes and is sensitive to differential stress [8][9][10][11][12][13][14] and material properties [15][16][17] , which in turn are also conditioned by fluid pressure and thermal gradients 18,19 . In particular, we investigated the seismicity and space-time variation of the b-value, considering about 13,700 volcano-tectonic earthquakes, which occurred at Mt. Etna between 2005 and 2019 ( Fig. 1).
This methodology enables investigating large crustal volumes, wherever they are "illuminated" with a sufficient number of earthquakes. We speculate that stress concentration zones along magma pathways may temporarily inhibit the occurrence of large eruptions. Our study suggests that for better outlining and understanding plumbing systems and the volcanoes' state, spatial and temporal analysis of seismicity and b-values could be applied.
Considerations on the physical mechanisms controlling the bvalue. The average size distribution of earthquakes, described by the slope of the FMD (i.e. the b-value) is inversely dependent on the mean magnitude, and thus to some extent on the mean length of seismogenic structures 20,21 . However, it has been confirmed many times in both laboratory 8,10,11,13 and field studies 9,12,14 that it depends not only on the material properties and fractal dimension of fractures and fault systems [15][16][17] but also closely on the local stress conditions. The b-value averages around 1 on a global scale 20,22 , which means a nearly tenfold increase in the number of earthquakes for successively smaller magnitudes, but local conditions may lead to remarkable anomalies.
Pockets of b-values considerably greater than 1, embedded in average valued crust 23 , are generally observed in volcanic and geothermal areas 24 , and an almost single school of thought 16,19,21,[24][25][26][27][28][29][30][31][32][33][34] trace them to the neighbourhood of active magma storage and conduits, corresponding to volumes repeatedly intruded by magmatic fluids, and thus intensely fractured by mechanical and thermal stress 18,19 . Consequently, a high fractal dimension of seismicity characterizes these crustal regions 16,17 , which means that a dense net of relatively smallscale fractures almost completely fills the seismogenic volume, with the preferential occurrence of frequent low-magnitude earthquakes (i.e. high b-values). Slip along fault surfaces is further favoured by the high fluid pressure, which has the effect of locally reducing the normal effective stress 9,26,27,31,34 .
Conversely, values of b considerably lower than 1, which result from a statistically reduced number of small-energy earthquakes, are related to stress concentration zones 8,9 , as asperities along fault planes 14,35 , where a relatively large volume is available for faulting 9,20 .
Temporal variations of the b-value must be interpreted in the framework of possible time-changes in the relevant stress conditions 10,11,13,14,36 . Experimental and numerical results show that increasing differential stress is accompanied by the closing of microcracks and their propagation and coalescence into larger ones, resulting in higher-magnitude earthquakes and thus in decreased b-values 11 . A rupture is likely to connect one subvolume to another in a highly stressed environment because the system contains more energy 8,12 . Conversely, an increase in the bvalue is believed to be the reflection of decreased differential stress, which can derive from the injection of pressurized fluids or the sudden release of seismic energy soon after large earthquakes 14,36 , resulting in a temporary higher proportion of small seismic events.
Etna volcano activity from 2005 to 2019. From 2005 to the beginning of 2008, the volcanic activity at Mt. Etna was confined to the summit area (within the Etna summit caldera, Fig. 1), with occasional powerful lava fountains or strombolian eruptions, and associated modest lava flows 37 (Fig. 2). Between 13 May 2008 and 6 July 2009, a major effusive eruption took place from vents close to the summit craters 38 . During the following year and a half, specifically until the end of 2010, eruptive activity diminished and 40 × 10 6 m 3 of magma is thought to have accumulated inside the volcano (according to the steady-state balancing law proposed by ref. 2 ; Fig. 2a). Summit eruptions vigorously resumed from January 2011 until the first half of 2016, with intense and frequent paroxysmal events that probably ensured the release of the formerly stored magma 2 (Fig. 2a).
Except for a few episodes (e.g. on 19 December 2009), the cumulative number of earthquakes and related seismic strainrelease ( Fig. 2b; seismic energy is computed following the relation reported in ref. 39 ) shows a linear trend from 2005 to the first half of 2016, despite the remarkable volumetric gap of erupted magma between mid-2009 and the end of 2010, also accompanied by intense inflation of the volcano 40 .
Afterward, a major change in the seismic pattern took place at Mt. Etna from the second half of 2016; in fact, exponential growth in the cumulative number of earthquakes can be observed (blackdashed arrow, Fig. 2b). Most of the earthquakes were located down to 8 km below sea level (all depths are relative to sea level), (Fig. 2c). Together, a decrease in erupted volumes 2 and an almost continuous long-term inflation 3 were recorded. During this time, up to the end of 2019, eruptions took place exclusively from the summit craters, through moderate to strong strombolian activity, associated with small lava flows; a further~50 × 10 6 m 3 of magma is believed to have accumulated within the supply system at the end of 2019 2 (yellow-dashed arrow, Fig. 2a).
Between 24 and 27 December 2018, Mt. Etna underwent one of the strongest known unrest phases, evidenced by a marked release of the seismic strain (Fig. 2b), a fast increase in the volcanic tremor, strong deformation 6,41 , and intense degassing 2 . This unprecedented unrest was produced by the intrusion of an estimated~30 × 10 6 m 3 of magma, which aborted its ascent at about sea level 1 . It was only followed by a short-lived eruption (red-dashed line, Fig. 2) that mainly occurred through a fissure opened within the Etna summit caldera, during which onlỹ 3 × 10 6 m 3 of magma were emitted to the surface 1-3 . Geochemical data 42 indicate a long-lasting prelude stage to the eruption begun in 2017 and involved magma-fluid accumulation at depths >7 km, followed by pressure build-up at 2-5 km deep, 6-7 months before the eruption; Etna was extraordinarily overpressurized before the onset of the December 2018 event 42 .

Results
Seismicity and b-values at Etna volcano between 2005 and 2019. To investigate seismicity patterns and stress conditions beneath Mt. Etna, with particular regard to the changes induced by the December 2018 major unrest, we performed a spatial analysis of earthquake foci distribution and related b-values over two time-spans; only the earthquakes that occurred between 2005 and 24 December 2018 (right before the unrest) were included at first, then we considered the entire 2005-2019 catalogue (seismicity and b-value sections and maps in Figs. 3 and 4).
The two models were subtracted to highlight relevant changes after the intrusive episode (percent b-value difference sections and maps in Figs. 3 and 4). If on the one hand, a certain degree of approximation subsists since the two models are mediated over two long and largely overlapping periods, on the other hand, this approach ensures a proper number of earthquakes to calculate reliable b-values over wide crustal volumes 23 .
Earthquakes were accurately relocated using the "tomoDDPS" software 43 and a 3D velocity model derived by the integration of passive 44 and active 45 seismic data, revised and optimized for this work (Patanè, D., personal communication; see "Methods" for details). The resulting earthquake foci distribution within the crustal volume beneath Etna volcano is shown in Fig. 3a, b. The b-values were estimated in those nodes of a 3-dimensional grid spaced 250 m, for which a minimum number of events (N min ) above the magnitude threshold of complete recording (say "magnitude of completeness" 46,47 , Mc), calculated locally for each of the nodes, was found over a maximum search radius of 3 km. To make the results less dependent on an arbitrary a priori choice of the N min , different b estimations were computed using four N min values (i.e. 50, 100, 150, and 200); the mean of such Given the high degree of complexity in computing b-values and estimating the associated uncertainties, only marked spatial and temporal b-differences were ascribed to effective variations in the physical condition of the crust 22,23 . Such differences were statistically evaluated using the p-test by Utsu 48 , which estimates the probability that two samples may come from the same population; differences between b-values can be considered significant whenever they are different at the 95% confidence level (i.e. p ≤ 5%) or higher 23 (see "Methods" for further details).
The following provides insights into different crustal sectors beneath the volcano.
Mid to shallow Etna crustal sector. An alignment of earthquakes, mostly localized above 10 km depth and trending nearly WNW-ESE, was detected in the south-western lower flank of the edifice; it appears to extend the San Gregorio -Aci Trezza fault (SG-AT, Fig. 1a) to the west. This arrangement (along the grey-dashed line, Fig. 1b) bounds a region that occupies the southernmost sector of the study area, where seismicity is rather low (Sector A, Fig. 1b); therefore, it was not possible to calculate b-values therein. Some authors 49 inferred the existence of a regional dextral shear zone that controlled the emplacement and distribution of the oldest Etnean volcanic centres, running along this same belt.
Seismic occurrence noticeably increases on moving to the north; the earthquakes confined within the first 10 km appear to be clustered along preferential NNW to NW directions, up to the Mt. Etna summit area (along the green-dashed lines within Sector B, Fig. 1b). The mapped faults follow the same trends (e.g. the Ragalna, Trecastagni, and the offshore Acireale Lineament continuing landward with the S. Venerina, S. Tecla, and Fiandaca system; R, TC, AL, SV, ST, F, respectively, Fig. 1a). In this sector, values of b less than or around 1 have generally been found down to almost 1.5 km (Fig. 4a, b, f, g); this may be due to local asperities that sometimes produce relatively large earthquakes when triggered by magmatic intrusions (e.g. the 27/10/2002 event, with local magnitude, M L = 4.7, which produced extensive surface faulting along the SV fault 50 , or even the moment magnitude, Mw = 4.9 earthquake that occurred on 26/12/2018 between the F and ST faults, Fig. 1a; see also refs. 4,41 ). At 3-4 km of depth, the b-values locally are up to over 2 ( Fig. 4c, d, h, i), suggesting a considerable degree of interaction between magmatic fluids and faults also at a distance from the central sector of Etna volcano, due to intense crustal fracturing and consequent elevated fluid-phase mobility. Deeper than 4 km, a general reduction of b is observed again (Fig. 4e, l).
The above-discussed seismic clusters do not continue further north of the Etna summit caldera, on the sides of which seismicity is generally between about 2 and 8 km depth and appears to be arranged along preferential E-W to NE-SW directions (above the purple-dashed line, within Sector C, Fig. 1b). On the western side, this seismogenic belt extends from the southern summit area to over 10 km further west ("Western Seismogenic Zone", WSZ); on the other, it is shifted to the north, running for almost 15 km, from immediately east of the summit caldera, along the northern sidewall of the "Valle del Bove" to the town of Riposto ("Eastern Seismogenic Zone", ESZ).
The analysis of this crustal sector highlights the occurrence of a "very high b-value zone" beneath the Etna summit, mainly between 1 and 6 km depth, which likely corresponds to an "intermediate Etna magma storage" (Figs. 3a, c and 4b-e). Considerable high b anomalies at similar depths were found in other volcanoes worldwide and likewise interpreted 24 , because the repeated use of the plumbing systems produces heterogeneous, shatter, and fluid-intruded crustal volumes, characterized by a high fractal dimension of seismicity 16,17,19 , or rather by a larger proportion of small earthquakes and vice versa. Moreover, this interval is consistent with the depth at which most of the volatiles evolve from magma 51 , with the effect of locally increasing fluid pressure, and where open cracks may effectively exist in the host rock 24 . Within this sector, a narrow volume of "relatively low bvalues" (Figs. 3c and 4c, d) stands right beneath the Etna summit caldera, at depths comprised between 2.5 and 4.5 km. Differences in b between the "intermediate Etna magma storage" and the enclosed volume with "relatively low b-values" can be considered statistically significant, since the p-value test by Utsu 48 gives more than 99% confidence (i.e. p ≪ 1%).
After the December 2018 dyke intrusion, a marked decrease in b (down to −70% difference; p ≪ 1%) is observed within the "intermediate Etna magma storage" (Figs. 3b, d and 4g-i; region 2, in Figs. 3e and 4n-p), due to reduced fluid pressure after the eruption; however, at the resolution of our data, the enclosed volume with "relatively low b-values" appears to persist to some extent, since it generally undergoes a much less marked reduction in b (locally around −15%, but also down to −45%; p ≪ 1%), (Figs. 3d and 4h; region 3, Fig. 3e). A considerable increase (up to 50% difference; p ≪ 1%) is observed along the "central plumbing system" (Fig. 3c, d; region 4, Figs. 3e and 4o-q), breaking around the region with "relatively low b-values" (black-dashed arrow, Fig. 3d). Finally, slightly increased b-values (up to around 20% Beneath the summit area, a NNW-trending seismic cluster (within the red-dashed oval, Fig. 1b), mostly confined within 1 km depth, with some events localized down to 5 km deep, marks the shallower portion of the Etna plumbing system. Low bvalues (Figs. 3c, d and 4a, f) are locally related to dyke propagation, which episodically produced earthquakes with M L up to 3.9, within the observed period.
Deep Etna crustal sector and plumbing system. At greater depths, a roughly NE-SW striking seismogenic volume steeply sinks to the northwest, reaching 10-20 km deep in correspondence to the "Parmentelli-Ragalna deep cluster" 52 , (blue-dashed oval, Figs. 1b and 3a) and north of the "Valle del Bove", and over 30 km beneath the town of Bronte (Fig. 1a). It is characterized by several thick clusters of earthquakes with M L up to 4.8, but most frequently around 1.5, and b-values around 1, both before and after the December 2018 intrusion (Fig. 3a, b).
Regardless of whether depth is a universal factor controlling b 53,54 or not 12 , the increase of frictional resistance to failure with confining pressure makes qualitative sense; under these conditions fractures tend to grow larger, producing a systematic lowering of the b-value 11,12 . At a global scale 9,53,54 , this was found to hold true generally down to nearly 15 km, while at greater depths an increase of b has been observed.
Beneath the north-western slope of Mt. Etna, no considerable increase in b-values occurs with increasing depth (Fig. 3a, b), likely due to the presence of an allochthonous seismogenic crustal slab with good mechanical properties, which keeps the b-value low 15,19,35 down to 25 km depth. Considering the eastern Sicily collisional geodynamic setting 55 , such a "NW sinking seismogenic volume" (Fig. 3a, b) should correspond to the "Hyblean foreland carbonate slab", nearly cropping out offshore of the town of Acireale 49 , (Fig. 1a) and steeply subducting north-westward beneath the volcano. Relatively low b-values below 20 km depth may also be favoured by inverse faulting style along NE-trending planes 44 since thrust faulting has shown to globally hold smaller b than that of normal or strike-slip events 12,22,35 .
From 30 km up to around 10 km, a deep "aseismic cusp" (Fig. 3a, b) occurs between the above-discussed seismogenic zones, narrowing upwards towards the central part of the volcano. Seismic tomographic studies report a low-velocity, triaxial ellipsoid region extending beneath the entire volcanic area between 25 and 15 km and passing up to around 10 km to a narrow high p-wave velocity body 56 . Such an aseismic cusp might be an expression of the partially molten, "deeper Etna plumbing system", where deformation mainly occurs aseismically, confined to the north-west by the subducting "Hyblean foreland carbonate slab".
A poorly earthquake-populated volume overlies the deep "aseismic cusp" in the central-southern portion of the volcano, almost up to sea level; it matches with a well-known "High p-wave Velocity Body" (HVB, Fig. 4f-l), generally interpreted as high-density cumulates, fractionated by magma during its ascent, stocked and solidified at depth 56 .
b-value monitoring at Etna volcano. Using seismology to investigate magma movement within the crust is one of the most intriguing long-term challenges of geosciences. From a monitoring perspective, the b-value time analysis (e.g. [57][58][59] ) must be viewed as a proxy of the average crustal stress variations, and could play an important role in tracking magma; nevertheless, several authors 26,35,36,60 confine its temporal variations to a secondary role. This is largely due to several issues that may derive from artifacts linked to changes in network configurations, reporting, or analysing procedures (which here can be almost neglected; see "Methods"), and also because the b-value is believed to primarily depend on stationary properties of the crust 23,26 .
We point out that spatial b mapping is more resolute and accurate than its time variations because by taking into account the whole catalogue and sampling constant, statistically significant N min above Mc around each close-spaced node of a threedimensional grid, distinct b-values can be calculated for relatively small adjacent sub-volumes; in this way, anomalously seismogenic zones, even of limited size, can be distinguished. Conversely, temporal b variations may not be diagnostic if calculated in too narrow regions; as a matter of fact, several earthquakes occurring in neighbouring sectors, characterized by diverse seismic patterns, must be included for reliable computing, with the effect of introducing non-local trends. Such issues can be bypassed to some extent if b time series are calculated within large volumes, which overall show relatively uniform characteristics in the spatial b analysis.
For monitoring purposes, it is more reasonable to investigate the crustal sector of the volcano directly involved in the final stages of magma transfer up to the surface. Therefore  Fig. 5). We used an event-by-event sliding window of 200 earthquakes that occurred within a maximum radius of 3.3 km, with N min above Mc (corrected by a factor of +0.2 47 ) equal to 50; each calculated b was associated at the time of the last earthquake. The most interesting results were achieved by computing the percent differences between each calculated b and the former one, then by plotting the cumulated series of such variations (see "Methods" for details).
Before the end of 2016, the b-value shows an overall regular trend, except for an abrupt fall on 13 May 2008 (Fig. 5a, b); it relates to stronger than normal earthquakes during the opening of the dykes that had drained a major lava effusion until 6 July 2009. The more pronounced drop in TS-1b, compared to TS-2, is linked to the propagation of fractures mainly at shallow depths. During the following year and a half, the remarkable volumetric gap of erupted magma led to a gradual slight rise of the b-value in the time series (Fig. 5b). A similar trend is also observed during the subsequent resumption of summit paroxysmal explosive activity in January 2011, which discontinuously lasted until mid-2016 and is believed to have drained most of the stored magma 2 . Such slight rise in b in the time series can be related to a certain increment of fluid pressure within the plumbing system, which had the effect of diminishing the differential stress, thus favouring the preferential occurrence of frequent, relatively small earthquakes during both the volumetric gap of erupted magma and the subsequent summit paroxysms. However, the observed variations are very small (on the order of 0.1 percent b differences), which suggests that only minor changes in the stress conditions perturbed seismicity patterns during Etna volcanic activity.
The main change instead occurred after the middle of 2016; the relevant decrease in erupted magma volumes (Fig. 2a) is accompanied by the exponential growth of the cumulative number of earthquakes (black-dashed arrow, Fig. 2b) and strongly increased b-values (Fig. 5a, c). The onset of such remarkable growth in b is observed at the beginning of 2017, corresponding very well with the long-lasting magma-fluid pressure build-up that preluded the December 2018 eruption 42 .
The magmatic volatiles injected within the plumbing system and the net of fractures cutting through its surroundings were not efficiently released this time, building up a huge fluid pressure inside the volcano, which in turn increases the b-values.
A very important event is observed on 01 December 2018 (Fig. 5d), when the b-values in TS-2 (i.e. within the "intermediate Etna magma storage") strongly increases, reaching the apex on the 6th of the month, about 19 days before the 24 December major unrest (more than 1.5% difference in 5 days). This trend is likely the effect of an overpressurization of the plumbing system between 1 and 6 km depth, leading to the preferential occurrence of many relatively smaller-than-normal earthquakes shortly before the unrest. During the same time, TS-1b (i.e. within the shallower Etna plumbing system) shows steady values. An abrupt drop of b started on 22 December only in TS-2 but became sharper on the 24th of the month in both TS-1b and TS-2; it marks the onset of the major unrest event, being related to stronger earthquakes due to the dykes' propagation, which at first involved the crustal region between 1 and 6 km depth (TS-2), and then propagated to the shallower sectors (TS-1b), leading to the eruption.
A marked increase in TS-2 on 25 March 2019 is observed, followed by a similar trend in TS-1b on 05 April 2019 (Fig. 5a, c), which may reflect magmatic fluid-pressure build-up within the "intermediate Etna magma storage", and subsequent transfer to the shallower sectors of the plumbing system. Soon after, in May 2019, moderate summit eruptions gradually resumed. At the end of August 2019, almost continuous moderate to strong summit strombolian activity, associated with small lava flows, was observed, but the time series do not show relevant variations. Finally, the marked increase in TS-2 on 30 October 2019, accompanied by almost steady values in TS-1b, is not associated with any relevant variation in eruptive activity.

Discussion
The origin of Etna's magma in the mantle is still the object of debate, especially considering its collisional geodynamic setting 55 , atypical for the emplacement of such a very active basaltic volcano. Nevertheless, from the distribution of earthquake foci and related b-values, we found that the rigid downgoing "Hyblean foreland carbonate slab" establishes a mechanical boundary, limiting to the north-west the "deeper Etna plumbing system", which is highlighted by an aseismic cusp (Fig. 3a). This arrangement favours the conveyance of magma from depth toward the present-day central part of the volcano and may have also conditioned, over time, the emplacement of a central-type stratovolcano. Ascending magma directly pushes on the bounding slab, sometimes triggering earthquakes with M L > 3 in specific sectors, such as the "Parmentelli-Ragalna deep cluster" 52 , (bluedashed oval, Fig. 3a).
Magmatic intrusions are also responsible for high seismicity characterizing the overlying mid-to-shallow Etna crustal volume, as revealed by the available focal mechanism solutions, which show a radial distribution of P-axes, indicating a pressure source located beneath the summit caldera 4,6,61,62 . Nevertheless, earthquake locations are not distributed radially from the volcano centre, but rather organize into distinct, well-demarcated seismic clusters, which can be related to (i) fracturing processes induced by magmatic intrusions and/or (ii) reactivation of pre-existing fault zones during magma uprising. In this latter case, intruding magmatic bodies would load the bedrock, triggering seismicity along faults originally delineated under the effect of tectonic forces, unravelling to some extent the regional structural pattern.
Beneath the summit caldera, specifically between about 1 and 6 km depth, very high b-values mark the depth at which gas exsolution and vesiculation of the ascending magma induce further recurrent fracturing processes, which over time produced a permanent fractured and heated crustal region, characterized by a high fractal dimension of seismicity. Here, a dense network of cracks and conduits overall facilitates the establishment of a persistent "intermediate Etna magma storage" (Figs. 3a, c and 4b-e). If the foregoing is true, then even more so is its central portion expected to show similar characteristics, since it stands right beneath the volcano summit, in a region that has been repeatedly injected by magma throughout Etna history. Instead, a limited volume with "relatively low b-values" has been detected there (Figs. 3c, d and 4c, d, h).
Concerning the origin of the relatively low b-value region, one may think that in an open-conduit condition, like that of Etna volcano, magmatic volatiles may more easily escape along the central axis of the plumbing system, avoiding the building of such a huge fluid pressure to remarkably increase the b-value, as instead observed in the surroundings. A weak point of this interpretation is the narrowness of this relatively low-b volume since it only occurs between 2.5 and 4.5 km depth. Moreover, despite the local slight variations in b after December 2018 (region 3, Fig. 3e), it appears to persist to some extent, being laterally bypassed by the effects of the intrusive episode on the central plumbing system (black-dashed arrow, Fig. 3d), while more marked changes of the b-value occurred in all the neighbouring sectors (regions 1a, 1b, 2 and 4; Figs. 3e and 4m-q). The relatively localized stability of b-values in time might decouple its origin from the episodic occurrence of stronger earthquakes during dyke propagation, or variation of fluid pressure, otherwise, such a big unrest phase would also have noticeably perturbed the local stress conditions and the b-value.
Another possibility is that local high differential stress may have contributed to the occurrence of such zone with "relatively low b-values" within the "intermediate Etna magma storage". Albeit in a different volcanic context, some authors 27 inferred that the fast ascent of magma that pressurizes the system without degassing can locally increase the applied stress beneath volcanoes, reducing b-values. This can hardly be our case, because the exsolution and diffusion of magmatic volatiles within the hosting rocks are testified by the surrounding "very high b-value zone", in which b diminishes notably after the December 2018 eruption (region 2, Figs. 3e and 4n-p), as a consequence of reduced fluid pressure due to intense degassing 2 .
Furthermore, the role of thermal gradients 18,23 in controlling temporal variations of b-values is very difficult to invoke because, aside from heat transport by magmatic volatiles, whose main effect is that of reducing the applied stress and increasing b, thermal conduction in the crust occurs over much larger times 27 .
As a final issue, a tectonic cause is proposed. The volume with "relatively low b-values" is located right where the WSZ and the ESZ are shifted with respect to Etna's centre (Figs. 1b, 3c and 4c, d,  h); in the hypothesis that they represent two en-echelon fault zones, the regional NW-oriented maximum horizontal stress field (SHmax, Fig. 1b) affecting the Etna region 63,64 may induce dextral strike-slip kinematics on them, in turn producing a highly stressed restraining step-over within the plumbing system. In this way, a growing compressive tectonic strain may have gradually created a physical obstacle that became critical for magma transfer up to the surface, leading to the continuous accumulation inside the volcano. It may also explain why no eruptive events capable of efficiently emptying the plumbing system occurred from the second half of 2016 up to the end of 2019, determining a considerable volumetric gap of erupted magma, accompanied by an increase in the number of earthquakes (Fig. 2b, c) and b-values in the time series (Fig. 5a, c). Such relevant changes in the seismicity pattern were not observed instead during the former 2009-2010 gap (Figs. 2b, c and 5a, b); this observation points to a different condition, which seems not to include any anomalous accumulation of magma inside the volcano. As a matter of fact, this latter gap starts at the end of a large effusive eruption, which may have considerably depleted the magma reservoir.
Under the action of NW-trending horizontal stress, maximum extension is predicted to be roughly ENE-WSW oriented 65 , in agreement with the NS to NNW direction of most Mt. Etna dykes 1,4,41,56 . This is also the trend of the seismicity characterizing the shallower portion of the plumbing system (within the reddashed oval, Fig. 1b), suggesting that magmatic intrusions exploit zones of weakness originally shaped by regional tectonics 49 .
According to this latter hypothesis, even though tectonic processes are related to the movement of huge crustal plates, they might locally result in small-scale deformation that can be even faster than previously expected, so rapid as to severely perturb the isostatic uprising of magma within the crust. However, this framework is only conceptual, because, as stated above, most focal mechanisms indicate that seismicity above 10 km depth is mainly controlled by magmatic intrusions 4,6,61,62 . For this to hold, one must consider the possibility that, aside from the volcano dynamics, which episodically results in relatively strong earthquakes, regional tectonics continuously control fault kinematics in the background.
Regardless of the causative source, we observe a major change in seismicity patterns and Etna behaviour after mid-2016, and this condition persisted at least up to the end of the observed period, despite the occurrence of one of the largest known episodes of unrest in December 2018. Increased confining stress within magma pathways may account for such a condition because it can inhibit bubble expansion and then the transition to foamy magma; in this way, the capacity to drive large volumes up to the surface and produce paroxysmal explosive eruptions (lava fountains) can be temporarily limited. Upward migration of magmatic volatiles is also reduced, which favours their spreading within the hosting rocks, increasing b-values in the time series from the beginning of 2017 (Fig. 5). The arisen condition would result in almost continuous summit activity, like that observed from August 2019 until the end of the year, fed by a reduced fresh and foamy magma supply from the deeper plumbing system and the re-mobilization of the liquid residue that fills the shallower sector. After the December 2018 eruption, it should be noted that at shallow depths a moderate rise of the b-value is observed again just northeast of the summit caldera (region 1b, Fig. 4m, n), suggesting a new localized increase in the fluid pressure beneath the volcano.
The time series (Fig. 5) show strongly increased b-values within the "intermediate Etna magma storage", corresponding to an overpressurization of the system starting on 01 December 2018 and culminating on the 06th of the month, about 19 days before the 24 December major unrest event. This latter was preceded by an abrupt fall of b that started 2 days in advance within the magma storage and then propagated to the shallower sectors of the plumbing system, tracking dykes' propagation for depth to the surface. This observation may represent a compelling case that monitoring the b-value could be used not only to investigate magma dynamics, its eruption, and then to characterize the state of volcanoes-as observed at Mt. Etna before, during, and after the 24 December 2018 strong unrest phase -but also for improving eruption forecasting.
Previous studies at Mt. Etna found evidence for increased bvalues from 1981 to the first half of 1987 57-59 , during a long period of intense eruptive activity but often supplied by gasdepleted magma. Spatial variations of b have been investigated on this volcano between 1990 and 1997 28 and also from 1999 to 2001, before the July-August 2001 eruption 29 . Finally, time and space variations of b have been investigated from August 1999 to December 2005 30 . Nevertheless, a former analysis comparable to that here presented is not available, and it cannot be ruled out that similar conditions also occurred in the past, leading to major eruptions (e.g. 1992-1993, 2001, 2002-2003).
The study of the b-value needs to be framed in broader correlation with other monitoring measurements to better define patterns linked to volcanic activity. Nonetheless, our results indicate that the systematic spatial and temporal analysis of this parameter may offer an opportunity to investigate the volcano state and improve the assessment of impending volcanic eruptions, this latter being the biggest concern of the local population and civil protection authorities.

Methods
In this study, we analysed the spatial-temporal distribution and the b-value of the seismicity recorded at Etna volcano from 2005 to 2019 66 by a seismic network of over 30 stations, almost all equipped with broadband three-component seismometers, managed by the "Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo" (Catania, Italy).
Establishing spatial and temporal variation of the b-value. Reliable computation of the b-value must account for possible artifacts, i.e., unreal changes in seismicity patterns, linked to how the observations are bundled and analysed 67 . These problems can be reduced in case a large number of earthquakes is included in the FMD because for a small sampling the log-linear fit of the distribution is degraded by the statistics of small numbers 68 . Unfortunately, further issues on volcanoes derive from the high background seismic noise (e.g. volcanic tremor), which may cause decreases in the detection of the smaller events, making the earthquake catalogues incomplete and thus leading to deviation from linearity at low magnitudes 46 . Furthermore, non-double-couple seismic signals, related to volcanic dynamics (such as explosion quakes or long-period events), may lead to bimodal FMD 16,24 and must be discarded.
Accurate seismic monitoring thus is needed for a proper analysis of the b-value, especially in volcanic areas. The Mt. Etna seismic network, operating since 1989, was greatly enhanced in both the number of stations and hardware at the beginning of 2005, making it one of the best seismically monitored volcanoes worldwide. Network improvement produced remarkable accuracy in hypocentre locations, mitigating the effect of artifacts linked to systematic mislocation and magnitude estimation, and also enhanced the capability to even detect the less energetic earthquakes, lowering the Mc. Furthermore, a unique magnitude scale (i.e. M L ) has been employed to estimate earthquake size; this avoids further systematic problems that may derive from the translation between different scales since the measurement of the b-value strictly rests on the accuracy in the determination of magnitudes 22,23,69 . Finally, continuity and homogeneity of recording are assured, since no important network modification or failure occurred, and no changes in analysing procedures were made.
M L strictly rests below 5 and source depth is shallower than about 30 km throughout the whole seismic catalogue, likely ensuring no breakdown in the selfsimilarity of earthquakes' mechanism 9,12,16,60 .
There is a question whether it is better to use the entire catalogue or its declustered part for mapping b-values in volcanic areas, and several authors 31 found that declustering led to only minor differences accompanied by a lack of information in some volumes. Together, as well as other volcanoes around the world, earthquakes at Mt. Etna mainly occur as "swarms", whereas foreshockmainshock-aftershock sequences are rarely observed (e.g. refs. 70,71 ). Since earthquakes in volcanic swarms are considered as independent, i.e., triggered by mechanical processes that are not controlled by previous events (e.g. ref. 72 ), and there are no sequences in the considered dataset, reliable declustering procedures cannot be applied to the Etna seismic catalogue.
Our analysis comprises the following steps.
Earthquake 3D relocations. About 13,700 earthquakes, with M L in the range 0.2-4.8, were relocated using the 'tomoDDPS' software 43 and a 3D velocity model derived by the integration of passive 44 and active 45 seismic data (the last acquired during the "TOMO-ETNA experiment", FP7 "MED-SUV" 73 ), revised and optimized for this work (Patanè, D., personal communication) to improve accuracy in hypocentre location. The algorithm has the advantage of using a combination of both absolute and differential arrival-time readings between couples of close-spaced earthquakes, having the same velocity profile along the ray path between the source and a common receiver. The velocity paths being almost identical, the difference in arrival times is reasonably attributed only to the spatial offset between the seismic events.
Earthquake 3D relocations produced a better clustering of the earthquakes and reduced the residuals between the observed and theoretical arrival times by about 44%, with an average value of 0.02 s. The final locations are affected by uncertainties of less than 0.7 km in the epicentral coordinates and less than 1 km in depth.
Mapping the b-value. The volume of interest was subdivided into a threedimensional grid spaced 250 m. For each node, the b-value was calculated taking into account the smallest set of closest earthquakes within a maximum search radius of 3 km, such that the number of events with a magnitude above the local Mc was greater than a fixed quantity (N min ). If these conditions were not satisfied, a null value was associated with the local b. The more events a data sample includes, the more robust the estimate of the b-value 23,35,68 . Monte Carlo sampling and bootstrapping of real data, as well as analysis of synthetic catalogues, shows that the error in b is large for small numbers of events and decreases rapidly to~15% for 50 events and less than 10% for more than 150 events, reaching about 3% for 1000 events (numbers represent 1σ and are of the same order as the standard error estimate proposed by Shi and Bolt 74 ). To mitigate the effect of an arbitrary a priori N min choice on the b-value, we considered the mean of the estimated b, derived from four different N min (i.e. 50, 100, 150, and 200); for the bins where less than four b-values are returned, only the available value/values were used. The associated uncertainties of final b values were calculated by using the error propagation formula. Information about the radius for the different bins using the two end members of N min (i.e. 50 and 200) and the mean of the four N min values are shown for the W-E vertical section of Fig. 3 in Supplementary Fig. 1 (see "Supplementary Information").
The estimation of reliable b-values depends critically on the choice of the Mc. In particular, the FMD increasingly deviates from a linear power low for earthquakes with magnitude smaller than Mc, since only a fraction of them can be detected. It follows that underestimating the Mc may result in unrealistic b-values. Conversely, if a too high Mc is considered, a significant number of earthquakes that would fit the linear power-low will be excluded, severely increasing the uncertainties in the estimation of the b-values 23,35,68 and, at the same time, decreasing the resolution for mapping 23 . Woessner and Wiemer 47 reported that, when considering subsets of earthquakes, one can find locally higher Mc than that calculated for the whole dataset, and they discuss the need for identifying such subsets when analysing spatial and temporal variations of seismicity parameters. Accordingly, also in line with several previous studies in volcanic areas (e.g. refs. 21,26,31 ), we considered the local Mc for every minor subset for b-values computing.
Mc was calculated using the maximum-curvature method 46 , which consists of defining the point of the maximum curvature by computing the maximum value of the first derivative of the FMD curve. Mc for the different bins ranges between about 0.7 and 1.9, being 1.1 for the entire catalogue. The estimation of the b-value and its associated uncertainty was performed by using the Tinti and Mulargia 75 maximum-likelihood formulation, which is specific for binned magnitudes and bias-free also for a small dataset. More than 95% of the bins, in both the periods of observation (i.e. from 2005 to 24 December 2018, right before the major unrest event, and from 2005 to 2019) have Mc above the bulk Mc calculated on the entire dataset (i.e. 1.1). For this reason and to also avoid excluding the local contribution of small volcano-tectonic earthquakes that would fit the FMD, we considered it better not to perform any magnitude-based pre-cut to the catalogue, which could be poorly diagnostic for the analysis of single subsets 47 . Following our approach, FMD above Mc for the different bins shows maximum goodness-of-fit with a linear power law (estimated following Wiemer and Wyss 23 ) of 92.8% and 94.2%, and average values of 62.4% and 61.7%, respectively for the two above periods.
To validate our approach, as a common practice when dealing with b-value computation 23 , we statistically evaluated the spatial b anomalies for each of the two analysed periods (i.e. 2005 to 24 December 2018 and 2005 to 2019) and the temporal differences between them by using the p-test by Utsu 48 , which estimates the probability that two samples may come from the same population, explicitly accounting for the number of earthquakes contained in each sample. In particular, the test starts from the condition of two earthquakes' populations to be the same and verifies how true is this condition. The hypothesis of belonging to the same population is more likely if the two sets share a considerable number of earthquakes, making the difference between two b-values less significant. Differences between b-values were considered statistically significant whenever they are different at the 95% confidence level (i.e. p ≤ 5%) or higher 23 .
Constructing b-value time series. The time-variant b was calculated using the same above-reported methodologies, but averaging the results from selected nodes within two crustal volumes of particular interest, each showing relatively uniform characteristics in the spatial analysis of the b-value (region 1b, Fig. 4m, n, TS-1b in In the time series analysis, to highlight the temporal variations of b, we performed a further processing step that consists in computing the percent differences between each calculated b and the former one. To enhance temporal smoothing, the former b-value was averaged with its preceding values (if present) within a 24-h window; the best result was then achieved by plotting the cumulated series of such variations (Fig. 5).

Code availability
The analytical codes for b-value computing were written with the MatLab software, by using formulas reported in the literature. The codes are available from Dr.