Re-pressurized magma at Mt. Etna, Italy, may feed eruptions for years

Identifying and monitoring the presence of pressurized magma beneath volcanoes allows for improved understanding of internal dynamics and prediction of eruptions. Here we show with time-repeated tomography clear evidence that fresh melts accumulate since 2019 in three reservoirs located at different depths in the central feeding system. In these three volumes, we observe a significant reduction of seismic wave velocity, an anomaly that has endured for almost two years. Reservoir re-pressurization induced seismicity clusters around the pressurized volumes within high fluid pressure compartments. This indicated a sharp change in volcano behavior, with re-pressurization of the central system replacing two-decade-long, flank collapse-dominated dynamics. The volume where the velocities are altered is remarkable in size, suggesting the injection of new melt, and that erupted lava represents only a small percentage. Our findings suggest that ongoing volcanic recrudescence can persist. Melt accumulating in multiple reservoirs beneath Mt Etna could feed eruptions for decades to come, according to detailed seismic tomography modelling.

D eciphering how volcanoes work is one of the most challenging issues of Earth physics and it is a prerequisite for hazard assessment. Many volcanoes erupt abruptly while others only after some period of unrest, but the unrest duration and the entity of precursory signals are still poorly understood. One major challenge is to forecast how long an eruption will last, how big will be its size, i.e., what will be the volume of lava that is going to erupt. Answering such challenging questions requires a deep knowledge of the magmatic system, good examples and multidisciplinary data. For many decades, our constant focus on Mt. Etna yielded great steps in the understanding of volcanic processes, taking advantage of the intense and disparate activity well monitored by modern multidisciplinary systems 1,2 . The volcano is located at the leading edge of the accretionary prism along the Apennines subduction system in southern Italy, in particular above the tear between the SE-ward faster retreating slab of the Ionian oceanic lithosphere relative to the continental Sicily 3 . Its hybrid origin emphasizes a variety of processes and complex dynamics raised by the contamination of a pristine low viscosity magma sourced by the asthenosphere above a vertical slab window [3][4][5] .
Although the volcano is continuously active, episodes of inflation and increase of eruptions are documented in recent and historic times 6,7 . Since December 2020, volcanic activity has resumed sharply with some paroxysmal events, after more than 1 year of long steady activity starting in September 2019. The present violent eruptive episode with spectacular lava fountains started in mid-February 2021. This sharp reactivation, still ongoing at the time we are writing, follows an overall inflation trend observed by GNSS data and a general increase in seismicity since summer 2020, suggesting a resumption of pressurization in the magmatic system 8 . Petrological analysis shows that the chemical composition of the magma that erupted from the Southeast Crater from December 2020 to February 2021 episodes is one of the most "primitive" and rich in gas of the last 20 years 9 .
The seismicity rate increase was pronounced, with shallow (z ≤ 3 km) clusters in the central area and inter-medial (3 km < z ≤ 12 km) and deep clusters (z ≥ 12 km) along the western and north-eastern flanks (Fig. 1). The occurrence of peripheral deep earthquakes is thought to be a marker of the volcano recharging, heralding eruptive phenomena [10][11][12][13] . It is then not surprising that in a short time, after an Ml 3.1 earthquake on February 11 2021, beneath the summit area, a series of paroxysmal events and lava fountains started at the Southeast Crater. Since mid-February 2021 to April 2021, the number of strong lava fountain episodes constantly increased, accompanied by a sharp increment in volcanic tremor amplitude and infrasonic activity, but with scarce earthquakes and moderate tiltmetric variations 9 . Activity in this period is referred as the first 2021 lava fountain cycle.
In this study, we get for the first time near real-time imaging of the magmatic system with time-lapse local earthquake tomography. We compute changes in P-and S-wave velocity between the present and older periods to trace and locate the new fresh accumulation of melts within the crustal reservoirs. The observed transient anomalies help in understanding the ongoing dynamics and quantifying the percentage of melt recently added to the system, useful to forecast the evolution of the unrest.

Results
In time-lapse tomography, 4D velocity changes (in space and time) are computed by differences between 3D images obtained by data that span different time intervals [14][15][16][17] . Applications in many geo-sites and different tectonic processes give strength to this pseudo-4D structural imaging 15,18,19 . The variations of seismic properties are related to changes in the stress and deformation, migration of fluids, and melts accumulation under volcanoes 14,17,19 . Limits in the application of time-repeated tomography derives from differences in ray sampling between distinct time periods, a consequence of the not optimal distribution of sources (local earthquakes) and receivers. Peculiar care in data selection and inversion is usually taken to ensure that the time-resolution of images is similar among periods, avoiding that spatial anomalies could map into spurious time anomalies 15,19 .
Here, we get inference on the present setting of the Mt. Etna magmatic system from differences in body wave velocity between the static 3D image of the volcano 20 and data in the present period (i.e., December 2019-February 2021). To enforce the modelling, we use static background models as the starting models for the inversion.
The distribution of P-and S-wave residuals for earthquakes located with the static 3D model gives a rough, blind, indication that data possibly contain significant velocity changes in time (Fig. 2). The asymmetric trend of residuals toward positive values indicates that the data need a slower model than the static 3D model. Selecting peculiar events and stations, we observe that ray paths that travel beneath the central part of the volcano structure accumulate positive residuals with respect to rays that travel outside the volcano ( Supplementary Fig. 1). The tomography successfully catches this pattern (see Method and data for details). A clear drop in residuals is observed before the sharp increase, suggesting that rays travel through faster anomalies. A similar drop is also observed at the end of the period, concomitant with a renewal of deep seismicity, with rays that sample deep fast anomalies.
The most salient feature of the Mt. Etna deep structure is a high-velocity body interpreted as a huge magma-mush formed during volcanic history 10,12,14,[20][21][22][23][24] . This body, extending from 3 to 15 km depth, is the main feature of any deep images of the volcano, as well as for this model. We catched several velocity changes from the background model (Fig. 3). Here we present and discuss the most intriguing features, relevant for the definition of the magmatic system: A broad volume with a strong Vp reduction, centred at 6 km depth (b.s.l.) and extending between 3 and 9 km depth is found beneath the central area (R-3). The Vp anomaly is well resolved, but the Vs resolution in the same region is very poor (Supplementary Fig. 2), although the potential ray coverage should be adequate based on earthquakes and station distribution. We explain this poor resolution as due to the difficulty in reading the phase onset on highly attenuated S-waves that pass through the anomalous body (see Fig. 4). These are indirect evidence for melt accumulation within the central feeding system. Seismicity is mostly located around the low Vp deep volume (Fig. 3).
Two small, central volumes with clear Vp reduction are centred beneath the summit craters at −1 (R-2) and 1 (R-1) km depth (b.s.l.). These anomalies are resolved and reasonably indicate fresh melt addition in the shallow reservoirs. A clear on/off of shallow seismicity at −1 km depth around R-2 is observed (Fig. 2d).
The three anomalous Vp volumes are different in size but they are all well resolved (according to the Resolution matrix analysis, see Data and Method section). R-3 is deep and elongated suggesting that the principal magma accumulation is vertically distributed at the edge of the high-velocity body. The volume of R-2 and R-1 is smaller, but they emerge as a permanent residence for magma feeding lava fountains and explosive activity of the summit craters. To verify if eventual bias in the velocity change is introduced by an uneven distribution of events, we have inverted the dataset following a conservative approach designed to minimize such changes, in which the static 3D model is computed with the entire dataset to absorb as much heterogeneities as possible 25 . We observe that the transient low-velocity anomalies are similar to those computed with the first approach (Fig. 5).
To further assess the reliability of velocity changes, we have computed a synthetic test where a transient low Vp anomaly is introduced only in the last period. The recovery of the synthetic feature is quite reasonable, although the procedure partially under-sizes the amplitude of the recovered anomaly (Fig. 6).

Discussion
Forecasting the evolution of volcanic activity is very compelling and controversial, even for well-monitored volcanoes like Mt. Etna. Short time scenarios are always based on time series of geochemical and geophysical data and changes in monitored parameters, but the information on how much melt is accumulated and at which depth it resides is often unknown or not considered. Earthquakes and low-frequency signals as volcanic tremor indirectly trace the ascent of magma within the crust and in local conduits close to the surface [26][27][28][29] , but estimates on the volumes involved are lacking.
Decades of studies revealed a central high-velocity magma mush beneath Mt. Etna, formed by long-term activity 21,30,31 . The mush incorporates a gabbroic intrusive counterpart of nonerupted lavas 32   of early 2000 was the first successfully observed and well documented by high-quality data from INGV-OE monitoring systems. These paroxysmal episodes followed a period of pressurization of the crustal reservoir, inferred by seismicity distribution and local stress 10,13,41,42 . Also, in that case, we missed indications on melt volumes that were large for having fed a prolonged activity. For all the following intrusive episodes, the identification of fresh melts in the shallow reservoir remained eluded by geophysics. Velocity changes from time-lapse tomography defined only local variation in volumes around shallow intruding dykes 12,14,20,23 .
After a period without important paroxysms, the sharp resumption of volcanic activity in mid-February 2021 represents the first opportunity to investigate the broad dynamics of the volcano. The reappraisal follows an increase of seismic release started in December 2020 where elongated deep (z > 12 km) and shallow (z < 3 km) clusters occurred in concomitance.
The need for a quantitative estimate of magma volumes potentially in play can be answered by time-lapse imaging of the volcano. The abundant and diffuse seismicity is a good prerequisite for its feasibility. We observe clear transient velocity anomalies associable with melt accumulation within three reservoirs located at different depths beneath the central portion of the volcano. The main addition is in a deep reservoir (R-3), broadly elongated between 4 and 9 km depth at the northern edge of the high-velocity magma mush (Fig. 7). This is the principal volume of storage where primitive melts ascending directly from the mantle source contaminate shallow residing magma. R-3 is connected with two shallower and smaller reservoirs (R-2 and R-1) that directly feed the summit craters activity.
Based on the distribution of velocity residuals vs time (Fig. 2a,  b), we hypothesize that this new accumulation started in mid-September 2019, in concomitance with the resumption of the summit activity after the end of the 2018 flank eruptions, with a recent acceleration in December 2020. The portion of melt added in 2 years is indeed significant to be traceable with local earthquake tomography, whose resolution is in the order of kilometres. Although the size of the transient anomalies is broadened by the alteration of temperature and pressure due to the injection, we can infer that fresh melts intrudes and pervades wide volumes of crustal reservoirs. Melt percentage can be extrapolated from velocity perturbations, following the approach used for the Kilauea 43 . By considering P-wave velocity differences and assuming a simplified composition compatible with the petrological characteristics of products and a mid-crustal depth (6-9 km    b.s.l.), we could expect that the −6% of Vp anomaly is compatible with a melt fraction of about 4% (i.e., about 1.4% for 1% of melt, see 44 ). This melt fraction has been added during the past 2 years to the central reservoir (R-3), reasonably after the start of pressurization in mid-September 2019 (see Fig. 2a, b). Such melt fraction could be a lower limit since the linearized tomographic approach tends to under-size the amplitude of recoverable anomalies as indicated by the test.
Our observations permit the development of a conceptual model of volcano dynamics. Very deep crustal earthquakes (z ≥ 12 km) around the volcano are the first signal of new magma ascent from the mantle source (Figs. 1, 2d). P-wave residuals in this period, during which shallow events are minor, show a sharp decrease, suggesting that rays from such deep events travel through the entire crustal high-velocity anomaly. The fluid upraise is followed by the accumulation of melts and pressurization within the main crustal reservoir (R-3). We resolved here the first tomographic evidence for the fresh melt accumulation within the reservoir. Seismicity responds to such broad re-pressurization, clustering within high fluid pressure volumes, as also documented for the late 1990 episode 41 . The increase of seismicity around the volume (Figs. 3 and 6) follows the repressurization of the reservoir, which is internally aseismic. Episodes of magma upraise from the main to the shallow reservoirs (R-2) are documented by rapid on/off of shallow seismicity (z ≤ 3 km, see Fig. 2d). Seismicity traces that deep and shallow melt injection within the reservoirs are almost contemporaneous. Magma residing within R-2 directly feeds the central volcanic activity similar to that active today. At present, the amount of erupted lava is only a small portion of the accumulated melts. We find that the overall repressurization of the main reservoir is indeed very conspicuous and could feed eruptions for years in the future.

Methods
The dataset used consists of local earthquakes selected by the INGV-OE (Istituto Nazionale di Geofisica e Vulcanologia-Osservatorio Etneo) database 45 From the whole dataset, we extracted 3140 earthquakes having at least 6 P and 2 S arrival times (with an average of 15 arrival times), RMS <0.4 s, formal hypocentral errors <2.0 km, azimuthal gap less than <200°(98% <180°) and a minimum epicentral distance of 10 km from the hypocenter to the closest station (Supplementary Fig. 3). The 1D location quality is summarized by histograms in Supplementary Fig. 4. From the selected earthquakes, we derived a total of 36,619 P and 11,789 S arrival times recorded at a maximum of 40 seismic stations deployed within 40 km from the volcanic summit area. Tomography inversions have been computed by using the Simulps14 code 47 , an iterative damped least-squares algorithm that solves for Vp and Vp/Vs in a 3-D grid with velocity values linearly interpolated through the medium and hypocentral/velocity parameters estimates progressively updated after parameter separation. As a starting model, we used the pre-2019 tomographic model 20 obtained by inverting a most comprehensive earthquake dataset collected at Mt. Etna from 2005 to January 15-2019 spanning a 20-year long period. Earthquakes selected for the 2019-2021 period have locations that spread diffusely beneath the volcano, ensuring a ray sampling similar to that of the overall dataset. To minimize eventual artefacts due to difference in ray sampling between different periods, we have also inverted the period using the approach described in 25 . A starting model is first obtained by inverting the entire dataset using conservative damping, to absorb as much signal as possible from the data into the static 3D model. Then, data from the sub-period are inverted using the static model as reference and with higher damping only to find the most robust transient signals.
The geometric characteristics of the grid mesh were the same as in 20 . Nodes are 2 km spaced in horizontal directions while layers are spaced 1 km from 2 km above the sea level down to 6 km depth. In addition, six deep layers are placed in the lower crust at 9, 12, 15, 18, 25 and 30 km to include in the inversion the deep seismicity occurring mainly along the peripheral rim of the volcano.
We set damping values of 80 and 120 for Vp and Vp/Vs models, based on the analysis of trade-off curves. After five iterations, we obtained a final rms of 0.13 s with a variance improvement of 34%.
To assess the reliability of the tomographic model, we compared, for each node, ray sampling and the full resolution matrix, following the method described by 48 . We computed the Spread Function (SF) to quantify the compactness of each parameter averaging vectors and the Derivative Weight Sum (DWS) to quantify the ray density around each node. The two-dimensional plot reporting the SF and DWS values for each node defines a sort of L-shaped trend, with DWS decreasing as SF gradually increases. Following 48 , the upper threshold of SF that ensures a well-resolved node and negligible smearing effects should correspond to the kink of the observed L-shaped trend. Similarly to 20 , SF = 3.0 is the limit value to define the well-resolved region of the model, decided visually by inspecting the averaging vectors of the matrix. For tomographic changes, we retain differences only for nodes where the threshold of SF ≤ 3 is matched both in the starting and in the final model.

Code availability
The seismic tomography software code is free and a version can be downloaded at http:// faldersons.net/Software/Simulps Received: 30 April 2021; Accepted: 14 September 2021;