Inflating Source Imaging and Stress/Strain Field Analysis at Campi Flegrei Caldera: The 2009–2013 Unrest Episode

Abstract

In this study, we analyze the 2009–2013 uplift phenomenon at Campi Flegrei (CF) caldera in terms of temporal and spatial variations in the stress/strain field due to the effect of an inflating source. We start by performing a 3D stationary finite element (FE) modeling of X-band COSMO-SkyMed DInSAR and GPS mean velocities to retrieve the geometry and location of the deformation source. The modeling results suggest that the best-fit source is a three-axis oblate spheroid ~3 km deep, which is mostly elongated in the NE–SW direction. Furthermore, we verify the reliability of model results by calculating the total horizontal derivative (THD) of the modeled vertical velocity component; the findings emphasize that the THD maxima overlap with the projection of source boundaries at the surface. Then, we generate a 3D time-dependent FE model, comparing the spatial and temporal distribution of the shear stress and volumetric strain with the seismic swarms beneath the caldera. We found that low values of shear stress are observed corresponding with the shallow hydrothermal system where low-magnitude earthquakes occur, whereas high values of shear stress are found at depths of about 3 km, where high-magnitude earthquakes nucleate. Finally, the volumetric strain analysis highlights that the seismicity occurs mainly at the border between compression and dilatation modeled regions, and some seismic events occur within compression regions.

1. Introduction

The Campi Flegrei (CF) volcanic area is a collapsed caldera with a well-documented and long history of ground inflation and deflation (bradyseismic activity). In the last 15 ka, ground deformations at long (resurgence) and short timescales have affected the inner caldera region [1,2,3,4]. Since 1905, ground-based topographic leveling networks, recently integrated with a global positioning system (GPS) monitoring network and satellite data, have recorded the unrest episodes of the Neapolitan Yellow Tuff (NYT) caldera. The major uplifts of 1.7 and 1.8 m were estimated at Pozzuoli harbor in 1970–1971 and 1982–1984, respectively [5,6,7,8]. After two decades of prevailing subsidence, the CF area experienced uplift periods (e.g., 2005–2006; 2012–2013), alternating with subsidence phases or stationary deformation trends [7,9,10,11]. The pattern of recent uplifts was almost radially symmetric with decreasing deformation velocity values moving from the Pozzuoli harbor towards the borders of the NYT caldera (Figure 1) [12,13,14].

One uplift episode occurred between 2012 and 2013, reaching ~11 cm with a maximum rate of 3 cm/month [15]. This event started on 7 September 2012 with a seismic swarm and the ground deformation observed from 2012 to 2013 has been interpreted by analyzing satellite interferometry and GPS measurements as resulting from a magmatic sill intrusion of 0.0042 km³ at shallow depth (~3.1 km) [10]. This is a major departure from other episodes of caldera deformation that occurred in the last 28 years, almost exclusively driven by shallow hydrothermal processes [13,16].

In the last few decades, many attempts have been made to reconstruct and identify the source(s) of the observed ground deformation of this significant uplift episode [13,14,16,17]. In particular, Amoruso et al. [17] showed the existence of pressurization thin quasi-horizontal ~4 km deep source with a negligible contribution of a shallower source underneath the Solfatara crater; Trasatti et al. [14] inverted COSMO-SkyMed (CSK) data, identifying their best fit model as a north-dipping mixed-mode dislocation source lying at a ~5 km depth.

The purpose of this study is to model and analyze the deformation source responsible for the unrest episode by exploiting the X-band CSK satellite SAR data and GPS measurements of 11 inner NYT caldera stations. We model the pressure change evolution of the source and evaluate the spatial and temporal variations in the acting stress/strain field beneath the caldera. In particular, we follow two main modeling steps: (I) we search for a 3D elastic source model responsible for the observed ground velocity occurring in the 2011–2013 time interval by jointly inverting the mean velocity deformation from DInSAR [15], retrieved from CSK images acquired along ascending and descending orbits, and GPS stations [11]. (II) After the location, geometry, and pressure change of the source underneath the CF caldera are inferred, we model, through a time pressure change function, the displacement time series between 2009 and the end of 2013, recorded by cGPS stations. The final aim of this study is to analyze the modeled spatial and temporal distribution of the shear stress (τ_MAX) and volumetric strain (ε_v), investigating their link with the seismic swarms, which occurred between 2009 and the end of 2013.

Figure 1. Campi Flegrei caldera. SRTM DEM of the area with red and yellow dotted lines indicating the Campanian Ignimbrite (CI) and Neapolitan Yellow Tuff (NYT) caldera borders, respectively (modified or adapted from [18]). The black lines show the faults and outer/inner rims [19]. The white squares indicate the locations of the caldera main eruptive centers. The UTM-WGS84 projection zone is 33 N.

Figure 1. Campi Flegrei caldera. SRTM DEM of the area with red and yellow dotted lines indicating the Campanian Ignimbrite (CI) and Neapolitan Yellow Tuff (NYT) caldera borders, respectively (modified or adapted from [18]). The black lines show the faults and outer/inner rims [19]. The white squares indicate the locations of the caldera main eruptive centers. The UTM-WGS84 projection zone is 33 N.

2. Deformation Data and 3D Source Model

We inverted DInSAR and GPS mean velocities deformation to infer the location, geometry, and pressure change of the 2001–2013 active source underneath the CF caldera. The DInSAR data are obtained from CSK images acquired along ascending and descending orbits from February 2011 to December 2013 and are the same data set used in [15]. We used 11 of the 14 available GPS stations [11] (Figure 2b); in particular, we did not consider the data from LICO, QUAR, and FRUL because their deformation velocities are too low to be considered in the modeling procedures (≤2.5 mm/year) and localized outside the NYT caldera border.

Both the ascending and descending LOS mean deformation velocity maps (Figure 2a) record a maximum of about 0.06 m/year at the Pozzuoli harbor, with the velocity pattern accommodated within the NYT caldera boundaries. This pattern is in agreement with the vertical velocity of GPS data (Figure 2b), which is equal to 0.065 ± 0.003 m/year at Pozzuoli harbor (RITE station), and exhibits a progressive decrease in velocity with distance, measured by GPS stations.

2.1. Model Strategy

The 3D finite element (FE) model [20] allows us to use the geophysical information available over the considered area. We employed the structural mechanical module within the COMSOL Multiphysics software package [21], and we performed a stationarity study, in isothermal conditions, by considering the medium as linearly elastic. In addition, we estimated the heterogeneities of elastic parameters by using seismic tomography [22] (Figure 3a).

Several authors [10,14,17,23] have assumed that the effects on the deformation of the hydrothermal system were negligible at shallow depth during the 2011–2013 uplift episode and that the medium surrounding the source exhibited linear elastic behavior. In particular, Troiano et al. [24] demonstrated the validity of this assumption, showing that the ground deformation associated with the uplift phase is nearly linear while the subsiding phase is characterized by a nonlinear viscosity-controlled trend. We built a 3D geometry of the area below the CF caldera, which extends for about 22 km × 16 km along the east and north, respectively, to a depth of 8 km; it also includes the topographic gradients at the CF ground surface (Figure 3). For such a computation, the caldera FE volume is discretized into a mesh of tetrahedral elements, with denser sampling in the source region. The use of a tetrahedral element discretization leads to the generation of a fine mesh that is readily aligned to the source and to the topographic surface. In particular, the entire numerical domain is discretized by considering 217,100 tetrahedral elements (Figure 3a) and by using a quadratic shape function, with the element size ranging from 200 to 800 m. As external boundary conditions, we apply a free surface condition for the upper boundary domain, which corresponds to the caldera topography; the bottom face is held fixed, while rollers are applied at the four lateral faces of the considered numerical domain. For the initial condition, we assume that the stress field is reproduced by applying gravity acceleration under elastic conditions: this condition allows the model to be compacted under the weight of the rock successions (gravity loading) until it reaches a stable equilibrium. This pressure distribution represents the lithostatic stress, starting from an atmospheric value at the surface and linearly increasing with depth (Figure 3b).

Subsequently, to model the observed deformation field, we considered the hydrostatic load applied to the 3-axis spheroidal filled source, also taking into account its internal body force. In order to retrieve the best-fit geometrical source features and pressure change, we performed an optimization procedure based on the Monte Carlo method allowing us to combine the efficiency of the derivative-based algorithm and the robustness of the random search algorithm [25]. We searched for the 3-axis spheroidal source (R_x,y,z), source center (C_x,y,z), source orientation (dip and strike), and pressure change value (

Table 1. Source parameters. R_x,y,z: spheroidal axes; C_x,y,z: source center; dip and strike: source orientation and pressure. Second and third columns: bounding of unknown parameters. Fourth column: estimated values with the related error.

Parameter	Lower Bound	Upper Bound	Estimated Value
Rx (m)	1500	4000	3000 ± 200
Ry (m)	1500	4000	2000 ± 200
Rz (m)	200	1100	580 ± 100
Cx (m)	425,000	427,000	426,330 ± 150
Cy (m)	4,518,000	452,000	4,519,300 ± 150
Cz (m)	−5000	−2000	−2990 ± 100
Dip (°)	−25	25	−4.5 ± 10
Strike (°)	−25	25	5 ± 3
ΔP (MPa)	10	40	24 ± 3

) by exploiting 18,270 and 17,349 DInSAR points for ascending and descending LOS mean velocity maps (Figure 2a), respectively. Simultaneously, we took into account the horizontal and vertical velocity components of 11 GPS stations (Figure 2b). The best-fit solution was selected from 1500 forward models (ascii file named TableS1_OBJ in Supplementary Material) via the Root-Mean-Square Error (RMSE) criterion:

$$\mathrm{RMSE}=w \sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n\left(m_i-d_i\right)}{n}},$$

where n is the number of data points, m denotes the model, d the data, and w represents the inverse-variance weight for each data set. We considered on average DInSAR velocity uncertainty equal to 0.003 m/year [26]; for GPS data, the uncertainty is equal to 0.001 m/year for northward and eastward velocities and 0.003 m/year for vertical velocity [11]. The comparisons between the optimized model and DInSAR and GPS velocities are shown in Figure 4.

2.2. Stationary Model Results

The inferred source is an oblate-like spheroid, located beneath Pozzuoli harbor at a depth of about 3 km b.s.l., ENE–WSW orientated and almost horizontal (about 4.5° west-dipping). Its hydrostatic pressure change is equal to 24 MPa. We specify that this value represents the modeled source pressure change necessary to overcome the lithostatic stress and allows the upper FE domain crust to instantaneously deform. All the search bounds and the estimated source parameters are reported in Table 1.

The estimated parameter errors are acceptable in a very narrow range of 10%, except for the dip angle that, although characterized by higher uncertainty, does not significantly change the spatial orientation of the inferred source, which appears to be sub-horizontal.

The modeled ascending and descending DInSAR LOS mean deformation velocity, the horizontal and vertical velocity components of GPS, and the associated residual maps are depicted in Figure 4a–c, respectively. The RMSE values for the LOS components are 0.0034 m/year for ascending and 0.0025 m/year for descending maps; for GPS, velocities are 0.0049 m/year for horizontal and 0.0013 m/year for vertical components. These results indicate a good fit between the data and the model for both LOS and GPS mean velocity components. Predominantly in the ACAE GPS station, a difference between the modeled and measured horizontal and vertical components is evident; it could be associated with the hydrothermal shallow source beneath the Solfatara crater, as proposed by Amoruso et al. [17].

In addition, in Figure 5a, we compare the measured and modeled DInSAR LOS mean velocities along an E–W profile (with 4520 km N, UTM coordinate), crossing the maximum deforming area, where a good fit is observed. As a further comparison, we show in Figure 5b the vertical section of the Vp seismic tomography [22] along the above-mentioned profile and the retrieved active spheroidal source; this highlights that the source is located just below the low-velocity zone of Vp < 2.8 km/s, which corresponds to a fractured region filled by hydrothermal fluids. This finding may support the hypothesis that the observed ground deformation is ascribed to a magmatic-related source rather than to the action of shallow pressured hydrothermal fluids; this is in agreement with previous analysis [27], in which the authors demonstrate that magmatic volatiles released by decompressing magmas at a critical degassing pressure (CDP) can drive volcanic unrest.

Finally, to fully interpret the achieved results, we computed the total horizontal derivative (THD) of the vertical component of the modeled velocity field (see THD_mod in Figure 6). We recall that the THD technique is an edge detection filter commonly employed for analyzing potential field data [28,29] and is based on analyzing the maxima of the horizontal gradient magnitude, as computed from the first-order x- and y-derivatives of the field. Accordingly, we used the THD technique to perform a basic detection of the deformation source boundaries. The results, reported in Figure 6, show high values of THD that overlap with the projection of oblate-like spheroid source boundaries at the surface, both for the modeled (THD_mod) and measured (THD_obs) fields already proposed by Pepe et al. [15]. The goodness of fit between THD_mod and THD_obs is also highlighted along an NNW–SSE profile, as shown in Figure 6b.

3. Stress and Strain Time-Dependent Study

In order to simulate the stress and strain field evolution generated by the retrieved active source, we performed a time-dependent modeling in elastic and isothermal conditions by exploiting the high temporal rate of cGPS measurements between 2009 and the end of 2013 [11]. In detail, assuming that the oblate-like spheroid source shape had not changed, we first extracted the time-dependent function from the GPS RITE signal by normalizing its vertical component (see Figure 7a), and then we multiplied this function by an unknown scalar value of pressure change (ΔP). To determine this unknown value, a linear method was used, estimating a ΔP equal to 78 MPa, and the retrieved ΔP time function that reproduces the observed cGPS measurements is reported in the upper panel of Figure 7b. The FE temporal analysis was performed by considering a time step equal to 0.1 year.

Subsequently, in order to estimate the ΔV time history, we calculated the source volume increase by considering the time-dependent FE modeled volumetric strain field. More specifically, the volumetric strain (ε_v) is given by ε_v = ε_x + ε_y + ε_z = ΔV/V, where V is the original volume, ΔV is the change in volume due to the deformation, and ε_i with (i = x, y, z) are the principal strain components. The bottom panel of Figure 7b shows the retrieved ΔV curve that slightly grows up to June 2012 (2012.5), then accelerates between June 2012 (2012.5) and March 2013 (2013.2) time interval, after which a plateau tendency is observed until the beginning of 2014. The total source volume variation is about 6∙10⁶ m³, with a peak of 6.5∙10⁶ m³ reached during the summer of 2013, when the maximum ground deformation is observed. These values are in agreement with those from other studies based on the same data set (i.e., ΔV_PTE in [17]).

In Figure 8, we compared the selected cGPS displacement time series of the three components (vertical, east–west, and north–south), plotted with triangles, with the time-dependent FE modeled components represented by the continuous lines. The results are generally good in the whole part of the NYT inner caldera, although some residuals are observed for the north–south component (blue triangles and lines). This discrepancy is particularly evident at the ACAE GPS station, as already highlighted in the previous stationary model results (Figure 4c).

Starting from the retrieved time-dependent FE model results, we evaluated the spatial and temporal three-dimensional distribution of the stress and strain fields generated by the active source to evaluate their possible link with the seismic swarms that occurred between 2009 and the end of 2013. In relation to the stress, the theory states that failure can occur when the maximum shear stress is equal to the maximum simple tension, and it is represented by the difference between the principal stresses (τ_MAX):

$${\mathsf{\tau }}_{\mathrm{MAX}}=\left(\frac{{\mathsf{\sigma }*}_1-{\mathsf{\sigma }*}_3}{2}\right)$$

In this case, the shear stress was calculated by removing the lithostatic stress (σ* = σ_tot − σ_lith). Figure 9 shows a time-lapse of the modeled shear stress and ground deformation at the end of 2013 (see also Video S1 in the Supplementary Material for the movie of the deformation and stress fields between 2009 and 2013).

The evolution of shear stress is shown in Figure 10a. It shows that the shear stress release, due to the source effect, is characterized by slightly alternating phases of low and high values between 2009 and the end of 2013. By comparing the temporal evolution of shear stress with the observed seismicity (M > 0) underneath the CF caldera, it seems that the shallow seismic clusters (blue and light blue filled circles) occur before the peaks, showing no link with the modeled shear stress, while the deeper seismic cluster (black filled circles) appears at the beginning of the first significant increase in the shear stress curve (see 2012.7 in Figure 10a). The 3D shear stress distribution allows us to discriminate the different causes of the recorded seismic clusters. Indeed, the 2010.3 and 2012.4 shallow seismic clusters (between 1 and 2 km of depth) may be interpretable as the effect induced by the fluid in the hydrothermal reservoir [16], and they are located between the Solfatara and Agnano craters. There is no evidence of these events in the modeled shear stress. Conversely, the August 2012 (2012.7) seismic cluster, characterized by deeper hypocenters (deeper than 2.5 km), is mainly associated with the stress induced by the modeled source; this last event, characterized by high earthquake magnitude, is located in an area with high shear stress. In the map shown in Figure 10d, the main spatial correspondence between the modeled stress and seismic cluster is evident at the NW of the source, where a negative value of τ_max is found. Although the August 2012 (2012.7) cluster is not directly generated by the source, we interpret it as the effect of stress concentration along preexisting structures or old volcanic chimneys that feed the monogenic cones, as in the case of the Gauro and Cigliano [30].

A further implication of the time-dependent three-dimensional analysis concerns the volumetric strain field [31,32,33]. Where a material is compressed, the changes in strain and length are negative; conversely, where a material is dilated, the changes in strain and length are positive. In Figure 11, we show the strain fields (volumetric strain) related to the source effect only and the seismicity distribution between 2009 and November 2012 (2012.9) with M > 0. A clear correspondence does not seem to emerge from this analysis in any of the horizontal maps (slices at z = −2.5 km) of Figure 11b–e between the occurrence of seismicity and the areas with different strain anomalies (positive or negative) induced by the modeled source. On the contrary, a good correspondence seems to exist, at depth of about −2 km along the sections, between the compression and dilatation regions and the considered earthquakes; we observe that the seismicity occurs mainly at the border between compression and dilatation regions and some events fall within the compression regions.

4. Discussion

In this study, we have analyzed the 2009–2013 uplift at the CF caldera in terms of temporal and spatial variations in the stress/strain field due to an inflating source. We started by searching for a spheroidal source via stationary FE optimization modeling of CSK DInSAR LOS and GPS mean velocities. To reproduce the structural and lithological complexity conditions of the crust beneath the CF caldera region, we built a 3D FE model by considering a realistic distribution of rock heterogeneities derived from seismic tomography. This allowed us to evaluate the impact both of the source geometric anisotropy (as the three-axis best-fit spheroid) and the mechanical heterogeneity of the surrounding rocks on the distribution of the shear stress from 2009 to the end of 2013. The optimization procedure suggests that the deformation best source is a three-axis spheroid compatible with a mostly elongated NE–SW ellipsoid (Figure 4 and Figure 5, Table 1). Its shape, an oblate-like spheroid, and depth (~3 km) are similar to what was found by D’Auria et al. [10], even if these authors found a narrower source. On the other hand, Amoruso et al. [17] and Trasatti et al. [14] found a pressurization sill at ~4 km deep and a north-dipping mixed-mode dislocation source lying at a ~5 km depth. Our retrieved source is located just below the low-velocity zone of Vp < 2.8 km/s (Figure 5), which corresponds to a fractured region filled by hydrothermal fluids [22,34]. This evidence suggests that the observed ground deformation may be ascribed to a magmatic-related source rather than to the action of shallow pressurized hydrothermal fluids. The hypothesis on the source nature is also supported by the seismic attenuation imaging performed by De Siena et al. in 2010 [34] and by fluid-dynamic analysis performed by D’Auria et al. [10].

To verify the reliability of the FE model results, we also performed a THD analysis on the modeled deformation vertical component at the surface (Figure 6a). The THD maxima overlaps with the projection of source boundaries at the surface, although we found some differences between the THD_obs and THD_mod along the selected profile (Figure 6b); this discrepancy can be attributed to the effects of pre-existing caldera structures (that have disjointed the intracaldera volcanic deposit), which were not considered in the proposed FE modeling setup.

To investigate any mechanical effects of the retrieved source on the surrounding crust, we developed a 3D stress/strain time-dependent FE model of the area below the CF caldera during the 2009–2013 unrest episode by comparing the retrieved results with the occurrence of seismicity (Figure 10 and Figure 11). The spatial distribution of the modeled shear stress reflects the different origins of the recorded seismic clusters. We observe that these clusters are substantially localized above the magmatic-related source and along its external edges. In particular, the clusters occurring up to June 2012 (2012.5) are concentrated in the region between the Solfatara and Agnano craters. In this region, no significant values of the shear stresses directly associable to the retrieved source have been found (see cross-sections of Figure 10b,c). According to several studies [16,35,36], these observed clusters could represent the effect of the active hydrothermal system. This interpretation is also supported by the ambient noise tomography and radar interferometry study proposed by Pepe et al. [15], which recognized the existence of a subvolcanic structure at a depth of 0.8–1 km. In particular, the rising of the pressurized hydrothermal fluids along this structure could justify the location of the seismic shallow cluster, the high rate of observed ground deformation, and the modeled shear stress (τ_max of about −0.3 MPa in Figure 10b). However, the August 2012 (2012.7) seismicity affects two distinctive areas: the first is located along the coastline, between Bagnoli and Pozzuoli harbor, with depth ranging between 2.5 and 3 km; the second is clustered in an area south of the Cigliano and Gauro craters, where a relevant negative value of shear stress (τ_max of about −1 MPa) is observed (see map in Figure 10d). We interpret this last seismic cluster as the effects of the modeled stress source along a pre-existing structure or volcanic chimney that led to the formation of the Cigliano and Gauro eruptive centers. Finally, high positive values of shear stress are found at depths of about 3 km surrounding the source, where high-magnitude earthquakes occur.

In respect of the strain field analysis (Figure 11), we observe that the rupture phenomena generating the seismicity seem to predominantly occur between compression (negative values) and dilatation (positive values) areas, emphasizing the role of the magmatic-related source in the nucleation of earthquakes during the uplift trend. This behavior has also been observed in other volcanic contexts [37]. In any case, the analysis of this finding is beyond the scope of this paper, but it may be a suitable subject to be investigated in the future by studying, for example, the trend and distribution of the following seismicity to verify if this correspondence is also confirmed in the last seven years.

Finally, we have summarized, based on the achieved modeling results, the Campi Flegrei volcanic system in a conceptual model. Specifically, from west to east, the uppermost crust beneath the caldera down to an 8 km depth is illustrated in Figure 12. We interpret that the retrieved intermediate magmatic-related source (deformation source #2 in Figure 12) acts as an interconnection between the deeper magmatic reservoir (thermal source #1 in Figure 12), already proposed by Castaldo et al. [38] by exploiting the thermal data, and the hydrothermal shallow system (hydrothermal source #3 in Figure 12). Furthermore, the existence of a magmatic-related source at about 3 km of depth is also confirmed by Chiodini et al. [27], where the authors demonstrate that magmatic volatiles released by decompressing magmas at a critical degassing pressure (CDP) can drive volcanic unrest towards a critical state, indicating that the magma could approach the CDP at Campi Flegrei.

On the other hand, hot fluids rising from the intermediate source through the highly fractured crustal zones feed the shallow hydrothermal source.

This shallow structure of the CF caldera is also in accordance with what was proposed by Kilburn et al. [23].

5. Conclusions

In this work, we created a physically based FE model of a volcanic system beneath the CF caldera following two steps: (i) modeling of LOS DInSAR and GPS mean velocities between 2011 and 2013 in order to retrieve the best geometrical features of the active source and (ii) analyzing the evolution of the stress/strain field distribution between 2009 and the end of 2013 by exploiting the high temporal rate of cGPS measurements. Our modeling approach allowed investigation of the possible mechanical effects of the retrieved source on the surrounding crust. The main findings are as follows:

• The DInSAR and GPS observed deformation velocity patterns can be explained by the action of an oblate-like spheroid source, slightly ENE–WSW oriented and represented by a filled hydrostatic pressurized source at ~3 km deep.
• The maxima of the calculated THD correspond to the projection at the surface of source boundaries, emphasizing the main incidence of the source on the observed deformation, although some discrepancies have been identified at the western region of the caldera where the influence of shallow structures probably becomes predominant in modulating the deformation pattern.
• The spatial distribution of the modeled shear stress during the 2009–2013 unrest episode reflects the different origins of the recorded seismic clusters. Indeed, low values of modeled shear stress are related to the shallow hydrothermal system, where low-magnitude earthquakes occur (0 < M < 1), while high values of shear stress are found at depths of about 3 km, where high-magnitude earthquakes (M > 1) nucleate.
• The volumetric strain analysis highlights that the seismicity occurs mainly at the border between the compression (negative values) and dilatation (positive values) regions, and some earthquakes fall within compression regions.
• The retrieved source (3 km deep), interpreted as magmatic-related, lies between the low Vp tomographic velocity anomaly (hydrothermal system) and the deep thermal anomalous region. It acts as an interconnection between the deeper magmatic reservoir (at about a 5 km depth) and the shallower hydrothermal system (0–2 km depth). Accordingly, the deformation source possibly accumulates magma batches from a deeper magma reservoir.

Finally, we remark that our FE approach allows us to distinguish the induced stress/strain effects attributable to the action of the intermediate magmatic-related source from those attributable to the action of the hydrothermal source. This result is relevant particularly for highly structured resurgent calderas, such as the CF, in which a deep magmatic system, an intermediate region, and a shallow hydrothermal system establish the complexity of the volcano plumbing system.