Hydrological and volcano-related gravity signals at Mt. Somma–Vesuvius from ∼20 yr of time-lapse gravity monitoring: implications for volcano quiescence

SUMMARY

We report on about 20 yr of relative gravity measurements, acquired on Mt. Somma–Vesuvius volcano in order to investigate the hydrological and volcano-tectonic processes controlling the present-day activity of the volcano. The retrieved long-term field of time gravity change (2003–2022) shows a pattern essentially related to the subsidence, which have affected the central part of the volcano, as detected by the permanent GNSS network and InSAR data. After reducing the observations for the effect of vertical deformation, no significant residuals are found, indicating no significant mass accumulation or loss within the volcanic system. In the north-western sector of the study area, at the border of the volcano edifice, however, significant residual positive gravity changes are detected which are associated to ground-water rebound after years of intense exploitation of the aquifers. On the seasonal timescale, we find that stations within the caldera rim are affected by the seasonal hydrological effects, while the gravity stations at the base of the Vesuvius show a less clear correlation. Furthermore, within the caldera rim a multiyear gravity transient is detected with an increase phase lasting about 4 yr followed by a slower decrease phase. Analysis of rain data seem to exclude a hydrological origin, hence, we hypothesize a deeper source related to the geothermal activity, which can be present even if the volcano is in a quiescent state. We infer the depth and volume of the source by inverting the spatial pattern of the gravity field at the peak of the transient. A volume of fluids of 9.5 × 10⁷ m³ with density of 1000 kg m⁻³ at 2.3 km depth is capable to fit reasonably well the observations. To explain the gravity transient, simple synthetic models are produced, that simulate the ascent of fluids from a deep reservoir up to the depth of 2.3 km and a successive diffusion within the carbonate aquifer hosting the geothermal system. The whole process appears to not significantly affect the seismicity rate and the deformation of the volcano. This study demonstrates the importance of a 4-D gravity monitoring of a volcano to understand its complex gravity signals that cover different spatial and temporal scales. Discriminating the different contributions that mix up in the observed gravity changes, in particular those due to hydrologic/anthropogenic activities form those due to the geothermal dynamics, is fundamental for a complete and reliable evaluation of the volcano state.

1. INTRODUCTION

The core goal of time-lapse gravity measurements at active volcanoes is the understanding of processes of underground mass transfer inducing mass/density changes over time, and is aimed at assessing the volcano dynamics and the resulting hazard (e.g. Williams-Jones & Rymer 2002; de Zeeuw-van Dalfsen et al. 2006; Carbone et al. 2017; Johnson et al. 2010; Greco et al. 2022). When time-lapse gravimetry is applied to active volcanic areas, the joint monitoring of ground deformation is pivotal. Indeed, the dynamics of hydrothermal fluids or magma is a coupled process in which both ground deformation, due to pressure variations at the source, and mass variations act jointly (Dzurisin 2006; Battaglia & Hill 2009). A proper separation of the deformation and mass-induced contribution, to the observed gravity variations is the key to a better understanding of the underground processes at work and their nature (e.g. magma vs hydrothermal fluids migration etc., Zerbini et al. 2001; Christiansen et al. 2011). In this framework, the expected contribution from time-lapse high-precision gravimetry is relevant.

During the last 20 yr, the improvement of the gravimetric network at Mt. Somma–Vesuvius (SV) volcano, in Southern Italy, provided the opportunity to obtain long-term series of gravity changes that can give, when jointly analysed with ground deformation (De Martino et al. 2021) and other hydrological data—important insights into the quiescent dynamic of this volcano and the contribution of fluids and water circulation to mass variation and/or redistribution (Berrino et al. 2013, Fig. 1). This goal is crucial especially in high-risk volcanic areas such as the case of SV volcano (Carlino 2021).

The SV volcanic complex has been investigated by different geophysical and geochemical methods (Zollo et al. 1996; Berrino et al. 1998; De Natale et al. 2001; Manzella et al. 2004; Caliro et al. 2011; Del Pezzo et al. 2013; D'Auria et al. 2014; Linde et al. 2017) providing important insights in understanding its structure and dynamics. The volcano, located east to the city of Naples, is formed by an ancient small caldera (the Mt. Somma) and a youngest volcanic edifice (the Vesuvius) reaching a maximum altitude of 1281 m a.s.l. (Figs 1a and b, Cioni et al. 1999). It is structurally located along a main faults system trending NE-SW (Fig. 1a). The oldest products of the SV volcano have been encountered in the Trecase well east to the crater (see Fig. 1b for location), and date back around 350 ka (Brocchini et al. 2001), while the last volcanic activity took place in 1944. The volcanic structure, inferred form the Trecase well stratigraphy, mainly consists of (i) a top layer of lavas interbedded with pyroclastic rocks and (ii) an intermediate to bottom layer mostly consisting of marine, sands and clays deposits interbedded with volcanoclastic deposits. This structure overlies the Mesozoic carbonate Apennine platform that was intersected by the drill at about 1800 m below the volcano (Cioni et al. 1999; Brocchini et al. 2001).

The present quiescence state of the volcano is characterized by a quite steady (hundreds of events per years) and low-energy seismic activity (maximum duration magnitude M_d = 3.6) which is mainly located along the main crater axis down to a depth of 4 km b.s.l. (D'Auria et al. 2014). A minor cluster of earthquakes is located between 5 and 8 km b.s.l. The bottom of the latter level roughly corresponds to a rheological transition, possibly throughout a zone of partial melting of rocks (Zollo et al. 1996; Carlino 2018). The stress generating the shallower seismicity has been correlated to the loading of the volcano which slowly deforms the deeper clays strata having different rheology. The response of the volcano edifice is sagging or spreading with fracturing in the shallow part (Borgia et al. 2005; De Matteo et al. 2022). The deeper seismicity is the result of the interaction between the regional stress field and the dynamic of the hydrothermal system of the volcano (Chiodini et al. 2001; D'Auria et al. 2014). The NeVoCGNSS and DInSAR data of the last tens of years at SV show that the central part of the volcano has undergone to a slight subsidence, with the maximum of about 1–2 cm yr⁻¹ along the crater axis (Lanari et al. 2002). In particular, the data provided by the NeVoCGNSS network of INGV (Osservatorio Vesuviano) show a clear difference in the subsidence rate, being greater at stations located at the base of the Vesuvius crater and within the caldera structure, and almost negligible at the stations located outside the caldera rim (De Martino et al. 2021).

Currently, a slight fumarolic activity with temperature around 60 °C occurs just inside the Vesuvius crater. The fumarolic fluids have a typical hydrothermal signature with major content of H₂O and a rich CO₂ gas phase (Chiodini et al. 2001). The CO₂ phase is produced mainly through metamorphic reactions involving marine carbonates (although a small contribute of magmatic CO₂ cannot be ruled out). The high-temperature hydrothermal system (∼ 450 °C) feeding the fumarolic activity is located along the volcano conduit at a depth of about 2 km (Chiodini et al. 2001). Furthermore, chemical analyses of fluids at SV indicate that the geochemical behaviour of the volcano aquifer is strictly controlled by variations in the input of deep-seated volcanic gases (Federico et al. 2013). The hydrogeology studies of the SV (Corniello et al. 1990; Celico et al. 1998) have shown that the volcano is characterized by a shallow and high permeable aquifer, hosted in fractured lavas and coarse-grained pyroclastic deposits, and a deeper one whose circulation takes place within the fractured carbonate basement at depth > 1.8 km, fed by the slow seepage waters from the Lattari and Sarno Mountains (Fig. 1a). Furthermore, in the southern and western aquifers the contributions of deep fluids appear more evident, while the eastern and northern sectors are dominated by shallow meteoric fluids (Federico et al. 2013).

Time-variable gravity survey at SV volcano have been performed since 1982, even if at a limited number of stations, while, since 2003, a most stable gravity monitoring network was established (Berrino et al. 2013). Analyses of gravity variations were performed by Berrino et al. (1997) and Berrino et al. (2013). In particular, Berrino et al. (2013) analysed both long- and short-term gravity variations inferring that only in the former case significant changes occurred. The authors took into account the gravity time series from 1982 to 2012 and considered negligible the contribute of the volcano subsidence to the gravity changes. They concluded that, during the studied period, a considerable redistribution of mass occurred at SV, essentially by density change at constant volume and possibly due to fluids migration in unsaturated porous medium (Berrino et al. 2013). Nevertheless, when observing the long-term gravity changes at slowly deforming volcanoes, the ground deformations need to be taken into account in order to separate the loading and/or compaction components from those due to the mass changes (Battaglia et al. 2008). Furthermore, long-term residual gravity changes and the short-term variations must be carefully analysed in the light of both natural hydrological cycles and possible anthropic activity (e.g. aquifers exploitation and water withdrawal). The latter activity is of particular importance in highly urbanized areas such that of SV volcano.

In the present work, we jointly analyse the long-term data of time-lapse gravity surveys and ground deformation (NeVoCGNSS, De Martino et al. 2021) at SV from 2003 to 2022. The long-term variation of gravity at each station is used to infer the corresponding mean gravity rate of change and its spatial field. Hereinafter, through ‘long-term’ we refer to the linear trend extracted from fitting the observations with a linear model. The long-term pattern of ground deformation, inferred from NeVoCGNSS network, is used to obtain the residual gravity field corrected for the effect of the observed vertical ground displacement through the free-air gradient (FAG). We analyse the result to obtain the long-term gravity imprinting of SV volcano and its relationship with internal volcanic and non-volcanic processes. We also analyse the short-term gravity changes to evaluate possible influence of the seasonal hydrological component. A fairly clear picture is retrieved by relating long-term gravity variation to volcano subsidence and to anthropic activity in its western sector. A detailed analysis of the gravity time-series reveals the presence of a multiyear gravity transient affecting the volcano edifice. To explain such gravity variations, we suggest a process of deep recharge of fluids, responsible of water migration (possibly a mixing of hydrothermal and meteoric origin) along the central sector of the volcano, and its subsequent radial diffusion within the shallower aquifer units.

2. GEODETIC DATA SETS AND PROCESSING

The measured variations are processed by taking into account the contribution of vertical ground deformations due to the slow subsidence, which the volcano manifested during the study period, to deduce the long-term residual gravity due to mass redistribution processes at depth. Since 2003 the INGV (Osservatorio Vesuviano) has performed repeated time-lapse gravity measurements on a network of more than 30 gravity stations in the Vesuvius area, using a relative D model LaCoste&Romberg gravimeter (LCR-D85, Berrino et al. 2013). The campaigns are performed at least once per year, with 2-yr gap during the COVID lockdown period, but until 2012, even more than one survey per year had been acquired. All the surveys are referred to an absolute reference station located outside the volcanic area, in the city centre of Naples (Largo San Marcellino, NArif in Fig. 1), which is considered a stable site (Berrino 1995). This station has been periodically surveyed by absolute gravity measurements to check its long-term stability (Berrino et al. 2013). Absolute measurements were also performed at the Osservatorio Vesuviano station (GVES20) from 1994 to 2010 and fitted well with the gravity variation from time-lapse gravity data (Berrino et al. 2013).

The gravity network has changed over the years, in particular the number of the stations has increased up to 36 (blue dots in Fig. 1b). Most of the stations (30) have a time-span of observations of more than 12 yr, while other six stations were established in 2014, in order to densify the network in the area around the crater. The network covers an area of about 200 km² centred on the volcanic edifice; the average distance between nearby stations is about 5 km in the most external areas, while, approaching the crater, the station-to-station distance decreases to 1 km. Unfortunately, GNSS measurements are not performed at each gravity station during the gravity campaigns. In order to monitor the height changes, we rely on the data provided by eight stations in the NeVoCGNSS network, supplemented by a few DinSAR points close to the volcano summit (Crosetto et al. 2020; De Martino et al. 2021). The network of continuous GNSS stations covers quite well the study area (yellow crosses in Fig. 1b). The interpolation of these data allows us to obtain a reliable reference ground deformation rate for each gravimetric station.

Time-lapse gravity measurements have been performed using a standard procedure that ensures data quality suitable for the detection of μGal-level gravity changes (Berrino et al. 2013). Gravity measurements are corrected for the Earth-tides and atmospheric pressure change effects. The Earth-tides effects are removed by synthesizing the Tamura's gravity potential catalogue (Tamura 1987). For a more accurate reduction of the tidal effects in the next future, we aim to use synthetic tides retrieved from local gravity records (Riccardi et al. 2023). At least three gravimeter readings are taken at each measurement point of the network. To get a reliable model of the daily instrumental drift, measurements are taken in the same site at the beginning and at the end of each survey day. Gravity measurements are collected according to occupation schemes of sites in concatenated closed loops to be able to produce a network adjustment in a least-squares sense (Hector & Hinderer 2016). The gravity observations, corrected for the above-mentioned effects, are then used to obtain the gravity differences with respect to the reference station in Naples. Finally, the so called ‘double differences’, namely the difference between the values (differences) obtained from successive gravity surveys, provide the time-lapse gravity variations. The procedure generally provides an error of few μGal for the time gravity changes and closure error on the whole network around ± 10 μGal.

For this study, we first inspect each time-series to perform a quality check and remove gross errors; through this procedure we excluded six stations from the original 30 that do not meet quality requirements (i.e. short-time span and low temporal resolution). We then apply different analysis techniques (linear trend extraction and low-pass filtering of the time-series) to highlight gravity variations at different spatial and temporal scales. Continuous GNSS data from NeVoCGNSS network are jointly analysed to account for the gravity changes induced by the long-term vertical movements. Forward and inverse modelling schemes are implemented to interpret more quantitatively the gravity observations.

3. RESULTS

Time-series of gravity variations from nine selected sites allow the characterization of the gravity field evolution and long-term trends in the different sectors of the study area (Fig. 2). The time-series of the vertical component provided by the NeVoCGNSS network are used to decouple the effects of vertical deformations from Newtonian variations (far right-hand column in Fig. 2).

Generally, the time variations are characterized by large amplitude oscillations with seasonal and/or multiyear periods superposed on long-period (LP) trends; in most of the stations the trends show a systematic increase of gravity. The multiyear oscillations have amplitudes of about 50 μGal, while the general rates are smaller and lie in the –1 to +8 μGal yr⁻¹ range (Table 1). Hereinafter, we analyse the gravity variations occurring at the SV considering three different temporal scales: short-term (seasonal), long-term gravity rate and multiyear transients. We mainly focus on the multiyear transients whose characteristics deserve a more detailed analysis.

3.1 Short-term gravity variations

The short-term gravity variations, where by short-term we mean differences between successive campaigns, are expected to be mostly influenced by the seasonal hydrological component and, in particular, by the local discharge/recharge process occurring in the aquifers in the Vesuvius area.

Most of the seasonal variations likely occur within the shallower aquifer (Celico et al. 1998), where the amplitude of the oscillations is in the range 0.5–3.5 m (Celico et al. 1998). The largest variations (3.5 m) occur within the plains surrounding the volcanic edifice, which collect both the meteoric waters drained by the Vesuvius and the waters provided by the carbonate aquifers of the Sarno and Lattari Mountains (Celico et al. 1998).

Considering such values for the water level variations, we should observe gravity differences in the range of 5–50 μGal between successive acquisitions, assuming an effective porosity of ∼30 per cent (for an equivalent sand deposit, Schön 2015) and a Bouguer slab approximation of the gravity effects (about 12 μGal m⁻¹). The recorded gravity variations generally lie within this interval (Figs 3a and b) but no relevant gravity changes are observed at the stations located within the caldera rim and those in the plains.

The stations within the caldera rim show coherent gravity trends (Fig. 3a), while a lack of coherence characterizes the gravity time-series from the stations at the base of volcano (Fig. 3b). The stations within the caldera rim show a quite good agreement with seasonal recharge, as shown Fig. 3(c), where the detrended cumulative rain is reported for three rain gauge stations located around the volcano (see Fig. 1b for location).

The short-term gravity changes display complex spatial patterns but, in many cases, one can distinguish the variations within the intra-caldera area from those from the surrounding plains (Figs 3d and e). The intra-caldera stations show gravity changes that seem mostly in accordance with periods of excess/deficit of rain; for instance, the snapshot of Fig. 3(d) shows a gravity decrease in the stations within the crater for the period 2004 May–November that agrees with the decrease in the detrended cumulative rain. A similar situation occurs also for the period 2008 December–2009 May, but in this case, we observe an increase of gravity associated to a period of rather intense recharge. For the stations located in the plain the spatial patterns are generally less coherent with alternating wet/dry periods (see detrended rain plot), and this may be the effect of the superposition of gravity signals of other non-local hydrologic contributions (i.e. far-field recharge from the carbonate aquifers).

3.2 Long-term gravity rate and volcano subsidence

The long-term gravity rates depict a clear spatial pattern, mostly coherent with the vertical rates of subsidence (Figs 4a and b) obtained by the interpolation of the GNSS/InSAR vertical rates (Lanari et al. 2002; De Matteo et al. 2022). The largest subsidence is observed in the crater area (2 cm yr⁻¹) and then it radially decays down to less than 1 mm yr⁻¹ at about 6–8 km distance from the volcano centre, in correspondence of the GNSS stations TERZ, PRET, ONPI and SANA (Lanari et al. 2002).

The subsidence causes a gravity increase due to the combined effect of height change and mass change of the topography; such effect is referred in literature as deformation-induced topographic effects (DITE, Vajda et al. 2019). The optimal approach to take into account this effect would be to measure the gravity gradient on the field and use it in combination with a modelling of the topographic effect to reduce the observations, but up to now all these data, that require a dedicated study, are not available for the SV. However, as suggested by Vajda et al. (2019), in case of prominent and rugged topography, as that of our study area, the FAG may be a first order approximation of DITE.

Hence, we perform the ‘residualization’ with the normal FAG (−310 μGal m⁻¹) in order to test, just in a first approximation, the compatibility of the observed gravity spatial variations with the subsidence pattern.

The residual gravity field depicts an overall decrease of the gravity rates towards the crater area (Fig. 4d); the central-northern sector of the volcano is characterized by slightly positive values, while in the southern sector of the volcano the values are close to zero. These non-null values can be considered poorly significant since they fall within the calculated standard error (Fig. 4e) of the fit. From these considerations, at least in the central part of the volcano and within the structural boundary of the Mt. Somma caldera, the long-term gravity increase observed at several stations can be ascribed to the volcano subsidence process, with negligible contribution from subsurface mass redistributions.

Moreover, we find some relevant gravity rate residuals at several stations located WNW of the volcano edifice (GVES14), far from the subsiding area. These gravity variations are significant and are above the standard error of the observations (Fig. 4e). To further check the reliability of the procedure to retrieve the gravity residuals the aforementioned analysis is performed using the DInSAR observations provided by the European Ground Motion Service (EGMS; Crosetto et al. 2020). The results are in good agreement with the GNSS outcomes (see Supporting Information, in particular Figs S1 and S2).

3.3 Hydrological component of residual gravity trends

Along the WNW direction the residual gravity rate increases up to the area of San Sebastiano (GVES14 station), 5 km NW to the crater, with a maximum value greater than 6 μGal yr⁻¹. The area affected by such anomalous increase of gravity has been interested by water withdrawal, for public use, from many wells (at least 180) WNW of the SV (Fig. 1; Lufrano and San Sebatiano areas). The aquifer was intensively exploited from about 1940 to 1990, causing a significant lowering of the water table. The available data from Lufrano wells field indicate that the peak of withdrawal was attained in 1989 when about 90 × 10⁶ m³ yr⁻¹ of water were removed (Corniello et al. 2003). In the period between 1979 and 1989, the water table dropped down by about 12 m in the NW sector of SV, in the Lufrano area, and by about 4–8 m closer to the SV at San Sebastiano. The water level started to slowly raise again after 1989, as a consequence of the progressive reduction of water withdrawal from Lufrano and San Sebastiano wells; an almost complete stopping of the pumping was reached around the 2000 (Corniello et al. 2003). Accordingly, the wells data indicate a general rising of the water table, between 1989 and 2003, of the same amount of the lowering occurred between 1979 and 1989, at Lufrano, and of about 2–5 m in the San Sebastiano area (Allocca & Celico 2008). The different recovering of the water level can be ascribed to the different hydraulic transmissivity of the aquifer (Celico et al. 1998), even if it must be considered the uncertainty in the assessment of the operating wells. Later studies evidenced the occurrence of groundwater rebound and flooding phenomena in the sector between Naples and SV, indicating that the rising of the water table of this area was still active at least until 2015 (Allocca et al. 2022). We think that the San Sebastiano gravity station (GVES14), which is very close to the wells field, is very sensitive to the local water level variation that until 2018 was probably still ongoing. In fact, assuming a drawdown in the San Sebastiano area of 8 m in the period 1979–1989 and a successive rebound of 2 m in the period 1989–2003, a complete recovering of the drawdown requires 6 m of rise in the period 2003–2022. A 6 m rebound in 19 yr implies a water rate of about 30 cm yr⁻¹. The gravity rate of 6.5 μGal yr⁻¹ corresponds to a higher water level rise of about 50 cm yr⁻¹, assuming a realistic porosity of 30 per cent. In any case, we have to consider that gravity data at San Sebastiano (GVES14) are available until 2018, so we have no data in the past 5 yr and we do not know if there was a change in the gravity trend, maybe to lower values. Considering the gravity data at GVES14 until 2018, about 7,5 m of water rebound are expected, a value that is not too far from the actual 6 m that we expect from the above considerations. Similar considerations can be made for the Sant'Anastasia (GVES11) and Volla (GVES13) stations (for Location see Fig. 1b); in this case the positive rates (cf. Fig. 2) are lower (respectively, 1.3 and 2.9 μGal yr⁻¹), suggesting conditions closer to the hydraulic equilibrium with respect to San Sebastiano (water level rates of about 10 and 20 cm yr⁻¹, assuming porosity of 30 per cent).

3.4 Multiyear analysis of gravity variations

In order to detect multiyear transients, we analyse the time-series after removing the long-term gravity rates and after applying a low-pass filter to reduce the effect of the seasonal oscillations. The original observations are oversampled to daily temporal resolution and then low-pass filtered, with cut-off period of about 3 yr, similarly to what done by Greco et al. (2022), for analysing the absolute gravity time-series acquired at Mt. Etna.

The filtered time-series of the gravity variations (Figs 5a and b) from stations located within the caldera rim of Mt. Somma and from stations located at the base of volcano show different trends. In particular, in the former case, we notice that all the stations in Fig. 5(a) record an average gravity increase of 30 μGal from 2008 to 2011 (Figs 5a and b; thick black line shows the stacked average). Later on, the average gravity decreases, going back to the pre-increase value in 2015. On the contrary, the stations located at the base of the volcano (Fig. 5b) show less coherent gravity changes during the whole period, with variations mostly within 10 μGal.

The corresponding spatial patterns related to the phases of ‘rising’ and ‘decreasing’ of gravity within the volcano are shown in the Figs 6(a) (2007–2011) and (b) (2011–2015). In both phases, the spatial patterns depict large-scale gravity changes with maximal amplitudes > 50 μGal and half-wavelengths of about 10 km centred on the volcano edifice.

Such observed transient can hardly be explained by the multiyear variability of the subsurface hydrology. This is supported by both the peculiar spatial pattern observed, which does not involve any relevant gravity change in the plains surrounding the volcano, and also by the time-series of the de-trended cumulative rain (Fig. 5d), showing no correlation with the gravity trends at these long-term temporal scales. More comparisons with hydrological models are reported in Supporting Information (Fig. S3). A possible mechanism to explain the observed gravity changes might involve the circulation of hydrothermal fluids at depth and their interaction with the water in the carbonate aquifers located below the SV (Celico et al. 1998; Chiodini et al. 2001).

We figure out the mass causing the gravity changes in the central part of the volcano, by first applying the Gauss theorem, for the period 2011–2007 and 2015–2011 and we find similar masses for the increase and decrease phases (+5.7 × 10¹⁰ and −4.0 × 10¹⁰ kg). Assuming a density of the source of 1000 kg m⁻³, we obtain volumes of + 5.7 × 10⁷ and −4.0 × 10⁷ m³, respectively. The volumes are clearly biased by the choice fluid density and porosity of the fluid-hosting medium: considering the quiescent state of the volcano since 1944, as evident from all the other geophysical and geochemical measurements which exclude shallow magmatic activity (Caliro et al. 2011; D'Auria et al. 2014) we can consider 1000 kg m⁻³ as an upper limit for fluid's density. Considering a realistic porosity of 0.3 has the effect to further enlarge the volume required to host the fluids by a factor of 3. As a consequence we obtain volumes of 1.5 × 10⁸ and −1.2 × 10⁸ m³ for the rising and decreasing phases, respectively.

4. INVERSION AND MODELLING OF MULTIYEAR GRAVITY VARIATIONS

To get more insights into this process we apply a more sophisticated method in which we invert the gravity differences observed between 2007–2011 in terms of horizontal location, depth and volume of a spherical source filled with a fluid with density 1000 kg m⁻³. The spherical source is chosen for its simplicity given that it requires the minimum number of parameters (4) to be inverted compared to other commonly employed sources in gravimetry (i.e. ellipsoid; Nikkhoo & Rivalta 2023); this is an important aspect to ensure a robust convergence to a solution, considering also the limited number of gravity observations available. Moreover, the relatively simple spatial pattern of the observed gravity variations, with a strong radial symmetry, further justifies the choice of a rather simple geometry of the source, as the spherical one.

For the computations we do not take into account the effect of porosity; in any case, as outlined before, including porosity affects the volume estimations in the same way as changing the average density of the fluids. Hence, for instance, assuming a porosity of 0.3 just increases the volume by a factor of roughly 3.

The inversion is performed by minimizing an objective function through the Nelder–Mead simplex method, which is an iterative derivative-free method apt for the optimization of nonlinear functions of multiple variables (Nelder & Mead 1965). We use the Matlab implementation given by the fminsearch function; further details on the algorithm can be found in Lagarias et al. (1998) and on the Matlab manual (https://it.mathworks.com/help/optim/ug/fminsearch.html). The cost function to be minimized is defined as follows:

where |$g_{\mathrm{ obs}}^i$| and |$g_{\mathrm{ mod}}^i$| are the observed gravity and modelled gravity at the ith point, respectively and N is the number of observations. |$g_{\mathrm{ mod}}^i( {\boldsymbol{x}} )$| and |$f( {\boldsymbol{x}} )$| are functions of the parameters vector |$( {\boldsymbol{x}} )$| to be inverted. For the inversion we set as stopping criteria two tolerance values between successive iterations: the first considers the change of the objective function and the second one the change of the inverted parameters. These tolerances represent relative bounds, i.e. the tolerance on the change of the objective function is defined as |$| {f( {{{\boldsymbol{x}}}_{{\boldsymbol{k}} + 1}} ) - f( {{{\boldsymbol{x}}}_{\boldsymbol{k}}} )} |/( {1 + | {f( {{{\boldsymbol{x}}}_{\boldsymbol{k}}} )} |} )$| and the change on the inverted parameters is defined as: |$| {| {{{\boldsymbol{x}}}_{{\boldsymbol{k}} + 1} - {{\boldsymbol{x}}}_{\boldsymbol{k}}} |} |/( {1 + | {| {{{\boldsymbol{x}}}_{\boldsymbol{k}}} |} |} )$|⁠, where, || denotes the L²-norm of the vector parameters. As it is evident, both the tolerances are dimensionless quantities and are both set to be 10⁻⁴.

The model that better explains the data are a sphere of radius 300 m (± 20 m), volume of 9.5 × 10⁷ m³ located at a depth of about 2400 m (± 570 m) below the sea level. The horizontal position corresponds roughly to the crater axis (North 451 873 m (± 310 m); East 4 519 660 m (± 420 m)). The uncertainty on parameters, shown in parentheses, is calculated as the standard deviation of 5000 inversions, where input data are randomly perturbed by Gaussian noise with standard deviation equal to that of the residuals.

The quality of the inversion process can be inspected by the maps and plots of Fig. 7, which illustrate the comparison between observations (Fig. 7a) and model outcomes (Fig. 7b), together with the residuals (Fig. 7c). The plot in Fig. 7(d) shows, as a function of the distance from the crater, the fit between observations (red asterisks) and modelled gravity changes (black circles). We see that, overall, we can reproduce quite well the long wavelength gravity change observed at the stations far from the crater area. Residuals of the inversion are mostly within the ± 20 μGal range (Fig. 7d) with the largest values for the stations in the crater area. The inversion process reduces the energy of the observations by more than 60 per cent, so that the RMS reduces from more than 30 μGal to about 10 μGal.

The source volume obtained through the inversion is in the same order of magnitude as that obtained by applying the Gauss theorem. Moreover, the depth of this source is compatible with the depth of the hydrothermal system of SV (Chiodini et al. 2001; Del Pezzo et al. 2013). To prove the robustness of our results, we perform further inversions, testing different geometries (prismatic body), using the raw observations and excluding/including some points (see Supporting Information: Figs S4–S7 and Table S1). All inversions are quite consistent with the one discussed above, and they all suggest the presence of a rather compact fluid mass at about 2–3 km within the volcano conduit. The actual observations seem not to be able to distinguish between different geometries of the source, hence a spherical or prismatic source provide equally satisfying fits to the observations.

The inverted mass represents a snapshot of the fluid position at the acme of the transient (i.e. end of 2011); however, we need to further justify where/when this mass originates and where it goes in the sequent phase associated to the gravity decrease.

We hypothesize a scenario that involves the rising of fluids from depth, that is then further validated with the aid of some synthetic models, to explain the gravity transient. We provide a possible mechanism of fluid rising compatible with the quiescent state of the volcano (see discussions section for further details). The scenario involves the following three phases:

• Gravity increasing phase: first, in 2007, a volume of about 9.5 × 10⁷ m³ of hydrothermal fluids (mainly water) rises from the deeper portion of the volcano towards the shallow hydrothermal system.

• Peak of the gravity: the fluids reach in 2011–2012 the shallow hydrothermal system, located at a depth of about 2300 m b.s.l within the volcano conduit. In this phase, the maximum gravity variation with respect 2007 is observed; the fluids are mostly stationing at this depth for about 1 yr, since gravity is not significantly varying.

• Gravity decreasing phase: the fluids are radially diffused mainly through the fractures of the carbonate aquifer, so the mass now is spread over a wider area with the effect of decreasing the observed gravity, especially at the stations at higher altitudes.

To verify more quantitatively the compatibility of this process with the gravity observations we implement some simple forward models that simulate the gravity response due to the ascent of fluids (i.e. to explain the 2007–2011 gravity increase) and due to the sequent fluids diffusion phase (i.e. 2011–2015 gravity decrease).

The model outcomes are compared with the corresponding gravity changes observed at different times and at different radial distances from the crater axis (plots in Figs 8c and d).

As for the ascent phase (2007–2011), we model the gravity effect of a sphere with the same volume of the inverted one (9.5 × 10⁷ m³), but located at progressively shallower depths (2300, 3500, 5000 and 9000 m). The results of the models are reported in Fig. 8(c) with different black lines (solid/dashed/dotted–dashed/dotted; corresponding depths are also indicated) together with the gravity observations, that are shown with colour codes proportional to the time. The observed gravity changes are all calculated with respect to the 2007 (begin of the transient). We see that, for instance, a source located at 5000 m depth would be compatible with gravity changes observed between the 2007 and 2009, while a source at 3500 m with the observations between 2007 and 2010. The beginning of the transient, on the other hand, indicates a quite deep source (9000 m). We remark that in these calculations we neglect the gravity contribution at surface due to the mass loss originated from the removal of the fluids at depth. This can be a reasonable approximation, assuming a depth's source > 9000 m; in this case, in fact, the gravity change at surface is expected to be below 10 μGal in amplitude (see the solid black line in Fig. 8c as a reference). Including such effect would slightly offset the simulated curves, but it would not change substantially the results and consequently the interpretation.

Regarding the diffusion phase the process is simulated by assuming an isotropic radial diffusion of the fluid into the carbonate aquifer. The diffusion causes the fluid mass to be distributed into larger areas; a reasonable approximation to model the process (and its gravity effect) is to spread the fluid volume into progressively larger cylinders (Fig. 8a, coloured circles) and then calculate the corresponding gravity changes at the topographic surface exploiting the prisms approximation (Nagy et al. 2000; implementation adopted from Uieda et al. 2016).

In particular, the calculation is performed as follows:

• we start defining different radii of the cylinders (1000, 2500, 5500 and 7500 m) and then calculating the corresponding thicknesses required to fit the total volume of fluids of the ascent phase;

• the area of the cylinder is then discretized with several prisms with base area of 100 m × 100 m; at each prism we assign a constant thickness (equal to that calculated at point 1) and density of 1000 kg m⁻³;

• the gravity effect is finally calculated at the different stations of the network, taking into consideration their elevation.

The change of the geometry, from a simple spherical body to a cylinder, is justified by the different nature of the process. In any case, as shown in the previous paragraphs, the actual observations do not allow to distinguish between different geometries of the source during the rising phase; inversions only suggest the presence of a relatively compact body at 2000–3000 m depth. Fig. 8(d) shows the gravity effects of the various cylinders along the profile with black lines with different line styles; the same style code is used in the map of Fig. 8(a) to depict the corresponding radial distances from the source, that can be interpreted as a front of the fluid diffusion.

As expected, the greater the radius of the cylinder, the lower the amplitude of gravity near the volcanic axes (Fig. 8d); in particular, we note that in 2011–2012 the diffusion front seems to be limited within 1–2.5 km from the crater, while in 2014 and 2015 the small amplitude gravity variations suggest that the front reaches a distance of 5–7 km away from the volcano axis.

A summary of the process occurring in the volcano is shown in the sketch of Fig. 8(b).

5. DISCUSSIONS

Time-lapse gravimetry in the SV area allows monitoring mass variation at different temporal scales relevant for understanding the volcano dynamics.

Our analysis of the long-term gravity rates testifies that the dominant process, observed within the caldera rim, is related to the subsidence of the volcano edifice (Borgia et al. 2005). In fact, the residual gravity after the correction of the subsidence effect is reduced and close to zero, with most of the residual variations within the estimated error on the rates. Such low residual field, within the caldera rim, further demonstrates that no significant mass accumulation or loss—related to both volcano and hydrological processes—occurred in the considered time-span of nearly 20 yr (2003–2022). The residual field should be in any case taken with a bit of caution since it is a first-order approximation, obtained with a standard gradient (Rundle 1982). Better refinements could be obtained with a campaign of vertical gradient measurements in the Vesuvius area, similarly to what already done at Mt. Etna or at Tenerife (Vajda et al. 2019).

On the other hand, we detect a multiyear gravity transient lasting at least 7 yr, mostly affecting the central sector of the volcano. The transient is constituted by a gravity increase starting at the end of 2007, with a peak reached in 2011 (up to 50μGal), followed by a decreasing phase lasting until 2014–2015.

This transient seems not to be related to LP variations in the surface hydrology and it is likely associated to volcano processes. The inversion of gravity data, in terms of a spherical source, provides a depth compatible with the hydrothermal system at the Vesuvius and a volume of 9.5 × 10⁷ m³, assuming a standard density of 1000 kg m⁻³ (Chiodini et al. 2001; Del Pezzo et al. 2013). Lowering the density value (i.e. more gas and vapour phases present) results in larger volumes; discriminating between all these possible equivalent sources would require a thermodynamic study on the physical–chemical state of the fluids, which is beyond the scope of this study. Porosity is another factor that introduces a bias in the volume estimation: including a realistic porosity of 0.2 or 0.3 increases the required volume of fluids by factor of 5 or 3, respectively. We remark that, for our simplified inversion scheme (a spherical source) all these models generate an equivalent gravity effect and we cannot uniquely distinguish any of them.

We suggest that the obtained volume represents the recharge rising from a deeper degassing source (probably located at 8000–9000 m, Zollo et al. 1996) up to the level of the hydrothermal system, located at about 2000 m depth. We propose that the rising of fluids is favoured by the tensile regime and the opening fractures, according to the deformation style of Vesuvius and its spreading (Merle & Borgia 1996; Borgia et al. 2005; D'Auria et al. 2014; De Matteo et al. 2022). The spreading, tends to reduce the loading stress in the central part of the volcano, a process that is accompanied by a constant low-energy seismicity and a slow depressurization of the volcanic system (Borgia et al. 2005). These processes generate favourable conditions for fluids migration in open cracks in the upper part of the volcano. Furthermore, the fluids storage and diffusion at around 2 km depth, may be also favoured by the presence of a lithological discontinuity between the carbonate basement and the sedimentary and volcanic layers.

From our model we can also roughly estimate the rising velocity of the plume fluids, that seems to be varying between the periods 2007–2009 and 2009–2011. In particular, between 2007–2009 the rising velocity would be around 7 × 10⁻⁵ m s⁻¹, being about 4 km the travel path of the fluids, while in 2009–2011 the velocity would reduce to 5 × 10⁻⁵ m s⁻¹, as the travel path is less than 3 km. The decrease of velocity could be explained as a reduction of the thermal buoyancy effect, as the fluids, during the ascent, release the heat and increase their density. It appears that the rising process has occurred in a quasi-steady state, without significantly perturbing the stress field of the volcano since no significant variations in the seismicity rate or change in the deformation pattern are observed (De Martino et al. 2021). This may suggest that the background seismicity of SV is mainly controlled by the lithostatic loading of the volcano, rather than by the fluids circulation (De Natale et al. 2000; D'Auria et al. 2014). Possibly, only faster transients of pressurized fluids from depth can generate an increase of seismicity, as occurred during the 1995, 1996 and 1999 seismic crises (Bianco et al. 1999; Ventura & Vilardo 1999; Caliro et al. 2011).

We also suggest that, once the fluids arrive at the level of the SV hydrothermal system, they diffuse throughout the surrounding carbonate aquifer. We demonstrate, through simple synthetic models, that this scenario is compatible with the observed gravity variations.

In our models we assume that, once they reach the minimum depth of 2.3 km, the newly arrived hydrothermal fluids are only diffused horizontally throughout the carbonate aquifer, without a significant vertical migration towards the surface. This assumption seems to be supported by the lack of relevant variations in the geochemical indicators as well as in the other geophysical parameters monitored (i.e. strain and tilt, Bollettini Sorveglianza Vesuvio OV-INGV 2023).

The analysis of the gravity time-series also revealed other signals of non-volcanic origin, which seem to be mostly related to hydrology. In particular, we detected short-term components related to seasonal variations and long-term rates, mostly reflecting changes in the aquifer storage units.

As for the short-term gravity variations, it appears that the local recharge may play a non-negligible role in the stations within the crater rim where rains are primarily collected by the ephemeral aquifers of the volcano. No other hydrologic contributions from nearby aquifers are expected here and, indeed, the observed gravity changes result to be quite well correlated with the rain time-series. On the other hand, in the areas of the plains the seasonal hydrology may be more influenced by the contribution of the deeper carbonate aquifer, that has a completely different mechanism of recharge, mainly influenced by the slow water flows from the recharge area of Sarno and Lattari Mountains (Celico et al. 1998; cf. Fig. 1a). This far-field hydrologic contribution may cause a further modulation of the gravity signal causing the loss of correlation with the local hydrology.

For what concerns the long-term hydrological components, the gravity increase observed in three stations located in the plains WNW to the volcano, seems to confirm a long-term rising of the water table as an effect of the reduction of the water withdrawal from wells field installed for civil purposes. The gravity rates would imply an equivalent water level rebound, with rates ranging from about 10 to 50 cm yr⁻¹. We have to remark that these rates cannot be quantitatively checked up to now, since hydrologic time-series for the period of interest are lacking. However, at least part of the observed water equivalent rates from gravimetry may be compatible with the estimates and observations of Corniello et al. (2003).

6. CONCLUSIONS

Our results show that gravity changes at quiescent volcanoes can shed lights on different processes, at various spatial and temporal scales that can be related to volcanic and non-volcanic processes. At Vesuvius, time-lapse gravimetry allows to track volcanological, hydrological and hydrothermal processes occurring at depth, thus reducing the inherent ambiguity on the interpretation of the surface manifestations. The joint analysis of temporal and spatial gravity variations allows to better assess the present-day volcano state. The gravity rates confirm the slow subsidence detected by GNSS and DinSAR data, which affects the central sector of the volcano. Our analysis unveils the presence of a multiyear transient, presently unknown, that we relate to mass movements of fluids (possibly a mixing of magmatic and non-magmatic origin) coming from a deep source (8–9 km depth). In our interpretation the transient would be part of a long-term process experienced by the volcanic system during its quiescence and possibly associated to the spreading of the shallower portion of the volcano.

Finally, our analysis evidences the superposition on the volcano-related gravity changes of a non-negligible hydrologic contribution, that is both seasonal and long-term. The lack of hydrologic data and models prevents quantitative considerations. Improving the knowledge of this component requires that discrete time-lapse observations are complemented by continuous gravity recordings with one or more superconducting gravimeters. This would make it possible to significantly reduce temporal aliasing related to the undersampling of temporal variations in gravity, as a consequence of the low repetition frequency of gravimetric surveys. The state of the art (e.g. Carbone et al. 2017; Portier et al. 2018) acknowledges that a hybrid gravimetric approach, pursued through the joint use of relative measurements on networks, continuous and absolute measurements, would improve the efficiency of gravimetric monitoring.

We remark that evaluating such different contributions in the gravity field is relevant in understanding the processes occurring during the quiescent state of the volcano and it has fundamental implications in recognizing significant variations related to the possible renewal of volcanic activity.

SUPPORTING INFORMATION

Figure S1. First Principal component of the EOF decomposition of the 4-D DInSAR database. (a) Eigenvector of the first mode (blue) and its linear fit (black dashed line). (b) Eigenvalues of the first mode, GNSS sites of the NeVoCGNSS (black dots).

Figure S2. DInSAR subsidence model and effect on the gravity residuals. (a) Subsidence pattern retrieved by EOF analysis of DInSAR observations. White dashed line: trace of profile shown in plots (c) and (d). Red squares: gravity stations. (b) Residual gravity rates after reduction of vertical displacement from DInSAR data. (c) Subsidence velocity profile (black dashed line) extracted across the Vesuvius (topography shown thorough the solid black line). Red dots: DInSAR vertical rates observed inside a buffer area of 1 km around the profile. (d) Gravity rates due to the subsidence (red dashed line), calculated assuming the normal FAG. Asterisks with vertical bars show the gravity observations.

Figure S3. Hydrologic time-series in the study area according to global hydrologic databases (GLDAS model) and local hydrologic observations (rain gauge). Top graph: moisture variations in the first 2 m of soil according to the GLDAS model (red) and seasonal variations retrieved by fitting annual and semi-annual components to the data (black line). Central and bottom plots: comparison between monthly rain of the GLDAS model and of the Pompei rain gauge (cyan vertical bars). Red lines report the detrended cumulative rain. Vertical grey dashed lines mark the timings of beginning, peak and end of the gravity transient observed by the stations within the caldera rim.

Figure S4. Spatial patterns of gravity changes associated to the (a) ‘increase’ and (b) ‘decrease’ phases of the transient.

Figure S5. Performance of the inversion 1 in Table S1. Top: fitting of observations and residuals; see the legend and the text for details. Plot with cyan line shows RMS error as function of inversion iterations. Histograms show frequency distribution of the input and model gravity values and the residuals after inversion.

Figure S6. Performance of the inversion 2 in Table S1. Top: fitting of observations and residuals; see the legend and the text for details. Plot with cyan line shows RMS error as function of inversion iterations. Histograms show frequency distribution of the input and model gravity values and the residuals after inversion.

Figure S7. Performance of the inversion 3 in Table S1. Top: fitting of observations and residuals; see the legend and the text for details. Plot with cyan line shows RMS error as function of inversion iterations. Histograms show frequency distribution of the input and model gravity values and the residuals after inversion.

Figure S8. Performance of the inversion 4 in Table S1. Top: fitting of observations and residuals; see the legend and the text for details. Plot with cyan line shows RMS error as function of inversion iterations. Histograms show frequency distribution of the input and model gravity values as well as residuals after inversion.

Figure S9. Performance of the inversion 5 in Table S1. Top: fitting of observations and residuals; see the legend and the text for details. Plot with cyan line shows RMS error as function of inversion iterations. Histograms show frequency distribution of the input and model gravity values and of the residuals after inversion.

Table S1. Sensitivity of inversion parameters for different inversion schemes and gravity inputs.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

This research was partially funded by the Agreement between Istituto Nazionale di Geofisica e Vulcanologia and the Italian Presidenza del Consiglio dei Ministri, Dipartimento della Protezione Civile (DPC). We thank Giovanna Berrino for providing us the gravity data from 2003 to 2020. The authors are grateful to Daniele Carbone and an anonymous reviewer for their comments and suggestions that greatly improved the paper.

DATA AVAILABILITY

Gravity data are available upon request to the corresponding author

References