Seismic Images of Pressurized Sources and Fluid Migration Driving Uplift at Campi Flegrei Caldera During 2020–2024

1. Introduction

The Campi Flegrei Caldera (CFC), located near Naples, Italy, is one of Europe’s most dynamic and potentially hazardous volcanic systems. The structural architecture of the CFC still remains a topic of debate (e.g., [1] and references therein). However, many researchers propose the existence of a nested caldera structure, primarily shaped by two major caldera-forming eruptions (Figure 1a): the Campanian Ignimbrite eruption around 40 ka (e.g., [2]) and the Neapolitan Yellow Tuff eruption around 15 ka (e.g., [3]). These two eruptions, part of the caldera’s eruptive history, were marked by intense and recurring volcanic activity. After the Neapolitan Yellow Tuff eruption, there were periods of relative calm, interspersed with at least 70 eruptions of moderate size, mostly explosive in nature. There were also a few larger events causing minor caldera collapses, including the Agnano-Monte Spina eruption around 4,500 years ago. The most recent eruption occurred in 1538 AD, leading to the formation of Monte Nuovo.

Since the Roman era, the CFC has experienced a gradual subsidence of about 1–2 cm per year, occasionally interrupted by periods of rapid ground uplift [4]. Since 1950, the caldera has undergone several periods of instability, characterized by alternating phases of uplift and subsidence. The first three events that disturbed the gradual subsidence of the caldera were characterized by similar uplift events in terms of duration and intensity: 1950–1952 (about 74 cm), 1970–1972 (approximately 159 cm), and 1982–1984 (about 178 cm) [5]. The 1970–1972 unrest recorded a minimal seismic activity, whereas the 1982-1984 crisis was accompanied by nearly 16,000 low-magnitude earthquakes, the largest of which had a magnitude of 4.0. Following this, a period of subsidence allowed a partial restoration of ground uplift.

In 2005, ground uplift resumed and in the area of maximum uplift, the Rione Terra (R in Figure 1a), is now 25 cm higher than its peak in 1984 and 129 cm higher than in 2005 (inset in Figure 1b). Seismic activity during this uplift phase has progressively intensified, especially since 2018. Between 2023 and 2024, over 30 earthquakes with magnitudes greater than 3.0 were recorded, including a magnitude 4.4 event on May 20, 2024—the strongest earthquake since 1985. This phase of uplift is also associated with increased degassing in the Solfatara-Pisciarelli fumarole area, which is the site of most seismic activity (Figure 1a).

Numerous geological and geophysical studies, including various tomographic investigations, have been conducted to improve the understanding of the CFC subsurface structure. Aster and Meyer [6] analyzed 228 seismic events recorded by a temporary network from the University of Wisconsin-Madison during the 1982–1984 bradyseismic crisis, deriving both P-wave (Vp) and S-wave (Vs) velocity models of the caldera.

Subsequent analysis mainly utilized datasets from the 1982–1984 seismic activity and the 2001 SERAPIS active seismic experiment [7,8,9,10,11,12,13] important to consider that the results derived from datasets of earthquakes that occurred during the 1982–1984 crisis should be extrapolated to more recent years with caution, as subsurface conditions may have changed since then.

Figure 1. a) A digital elevation model (DEM) of the Campi Flegrei Caldera, illustrating the caldera rims. The outer caldera margin is represented by an orange dotted line, while the inner caldera margin is marked with a yellow dotted line. The morphological boundary of the resurgent dome is delineated by a red dotted line (refer to Figure 2 in [16] for more details). The map also shows seismic stations (red triangles) and earthquake epicenters (white circles), which have been re-located using a 1D velocity model. b) Histogram displaying the number of daily earthquakes (black) from 2020 to June 30, 2024, with grey-shaded areas highlighting two distinct periods during which separate tomographic inversions were conducted for comparison. The inset includes: 1) a histogram of daily earthquakes since 1983 (black), with the count displayed on the left axis, 2) a cumulative plot of released energy (blue line), and 3) the vertical daily displacement at RITE (red line, Rione Terra, Pozzuoli). The values for 2) and 3) are provided on the right axis.

Figure 1. a) A digital elevation model (DEM) of the Campi Flegrei Caldera, illustrating the caldera rims. The outer caldera margin is represented by an orange dotted line, while the inner caldera margin is marked with a yellow dotted line. The morphological boundary of the resurgent dome is delineated by a red dotted line (refer to Figure 2 in [16] for more details). The map also shows seismic stations (red triangles) and earthquake epicenters (white circles), which have been re-located using a 1D velocity model. b) Histogram displaying the number of daily earthquakes (black) from 2020 to June 30, 2024, with grey-shaded areas highlighting two distinct periods during which separate tomographic inversions were conducted for comparison. The inset includes: 1) a histogram of daily earthquakes since 1983 (black), with the count displayed on the left axis, 2) a cumulative plot of released energy (blue line), and 3) the vertical daily displacement at RITE (red line, Rione Terra, Pozzuoli). The values for 2) and 3) are provided on the right axis.

The most recent tomographic studies have focused on analyzing velocity changes [14,15], emphasizing time-dependent variations in the velocity structure. The repeated application of three-dimensional (3D) seismic tomography, or 4D tomography (in space and time), can be considered a powerful tool to detect and track changes in the velocity structure under an active volcanic system. In time-lapse tomography, 4D velocity changes are calculated by comparing differences between 3D images generated from data collected over different time intervals. 4D seismic tomography is routinely used in mature fields to monitor the effectiveness of enhanced oil recovery (EOR) techniques and in different tectonic contexts, enhancing our understanding of subsurface dynamics. For example, it can track the movement of injected fluids, such as CO2 or water, and their impact on oil displacement.

In volcanic environments, changes in both Vp and Vs and especially temporal variations in the P-to-S wave velocity ratio (Vp/Vs) have been documented, such as at Etna during its activity between 2002 and 2003 [17]. The authors interpreted the observed lower Vp/Vs ratios as evidence of the intrusion of volatile-rich basaltic magma feeding the eruption. Similarly, significant changes in seismic velocities were detected beneath Mammoth Mountain [18], where the variations were consistent with CO₂ migration into the upper 2 km, with corresponding depletion in surrounding areas, correlating with surface venting zones. Furthermore, in the Geysers geothermal field, between 2008 and 2014, the strengthening of a low-Vp/Vs anomaly after 2010 was attributed to alterations in the water injection regime [19]. Also at the Nevado Del Ruiz, the temporal velocity changes from 2000 to 2016 indicated that the volcanic system underwent changes related to the ascent of new magma, which interacted with preexisting magma bodies, leading to the recent emplacement of a small dome at the bottom of the active crater [20].

Variations in seismic properties are typically linked to changes in stress and deformation, fluids and gases migration, and magma accumulation beneath volcanoes. However, the application of time-repeated tomography can be limited by differences in ray coverage between time periods (Figures S1 and S2), which may result from not always optimal distribution of seismic sources (local earthquakes) and receivers. To address this, careful data selection and inversion procedures are necessary to ensure that the time-resolution of the images is consistent across periods, thereby preventing spatial anomalies from being misinterpreted as temporal changes.

In this work, we present new three-dimensional (3D) Vp, Vs, and Vp/Vs models for Campi Flegrei during the 2020–June 2024 period, which saw a progressive and significant increase in seismicity during the ongoing uplift phase that began in 2005. First, we performed a tomographic inversion for the entire period from 2020 to 2024 and then compared the resulting seismic images with those from Battaglia et al. [11], which used data from the 1982–1984 seismic crisis and the 2001 SERAPIS active experiment. Additionally, in light of the significant increase in both the number and intensity of seismic events starting in 2023, we compared the tomographic images obtained for the 2023–2024 and 2020–2022 periods by analyzing velocity variations. This comparison aimed to identify any changes in velocity that could be linked to the ascent of fluids and/or magma in recent years.

2. Materials and Methods

2.1. Tomographic Inversion for the Vp, Vs and Vp/Vs Structures

It is well known how the accuracy and reliability of results obtained from tomographic inversion depends on several factors, including the inversion code used and the initial velocity model. The choice of starting velocity model is a critical aspect of the inversion process, as it significantly influences the accuracy and reliability of the final results. Selecting an appropriate initial model helps minimize potential biases and artifacts that could arise during the inversion [21], allowing for a more precise depiction of the subsurface structure.

In this study, we utilized various 3D starting models to thoroughly assess their impact on the inversion results. These included the 3D velocity models derived from the 1D velocity model from Vanorio et al. [9], with different Vp/Vs ratios ranging between 1.7 and 1.78, as well as the 3D velocity models obtained in Battaglia et al. [11]. This approach enabled exploring how different initial conditions could influence the final imaging of the subsurface structure, ensuring that the detected features are robust.

Moreover, the choice of the tomographic inversion code plays a pivotal role in accurately imaging the subsurface structure, as different algorithms have distinct approaches to solving the inversion problem and handling data uncertainties [22]. To enhance the reliability of our results, we used and compared two well-established inversion codes of the SimulPS [23] and tomoDD [24] families. The SimulPS14 code [25,26] is particularly effective for generating initial velocity models, thanks to its robust handling of local earthquake datasets and ability to make preliminary adjustments to the velocity structure. Conversely, the tomoDD code excels at refining these models through its double-difference approach, which reduces residual errors and enhances the resolution of subsurface features. Double-difference (DD) algorithms minimize the travel time differences between pairs of earthquakes observed at common stations. This method focuses on relative locations of earthquakes, which improves resolution in clustered seismicity. Moreover, the double-difference approach reduces the sensitivity to errors in the initial velocity model, making it computationally efficient for resolving small-scale features. The tomoDD version used in this study, named tomoDDPS [27], has advanced features for adjusting S-P times to account for different propagation velocities of P and S waves through various materials. Then, tomoDDPS is able to provide a more robust Vp/Vs model than that from the original DD tomography code (tomoDD), obtained simply by dividing Vp by Vs [28].

By employing both codes, we were able to cross-validate the results, assess the influence of the chosen starting models, and achieve a detailed and reliable representation of the Campi Flegrei Caldera’s complex subsurface. This combined approach ensured that our inversion outcomes were not overly dependent on a single method, thereby enhancing the robustness of our findings.

2.2. Dataset and Parameterization of Inversions

The majority of earthquakes in the Campi Flegrei Caldera take place at depths less than 3 km and exhibit a nearly elliptical distribution at the large caldera scale (Figure 1), as shown in the reference catalog of the Istituto Nazionale di Geofisica e Vulcanologia - Osservatorio Vesuviano [29](INGV-OV, https://doi.org/10.13127/gossip). The seismic activity mostly occurs inland, particularly under the Solfatara-Pisciarelli area and to the north of Pozzuoli, while it deepens offshore.

Typically, earthquakes have a duration magnitude of Md ≤ 1. However, since 2020, and particularly from early 2023 onward, there has been a significant rise in the average monthly magnitude. This increase includes over 30 earthquakes with magnitudes exceeding 3.0, the largest being a magnitude 4.4 earthquake recorded on May 20, 2024.

In this study, we considered only located events and the associated P- and S-wave arrival times from the INGV-OV earthquake database [29].

Typical lower magnitude events (Md < 2) present 6–10 P and 2–4 S arrival times. We extracted arrival times from 28 stations of the INGV-OV network (triangles in Figure 1), located within 15 km from the catalog epicentres. Since the data were collected by different operators at various times, the uncertainties in phase picking are not always reported and therefore not well known. Generally, the precision to which arrival times could be measured is ~0.01 s for P phases and ~0.2 s for S phases. Therefore, we decided to perform two different inversions, where we consider a uniform picking uncertainty across the entire dataset: one with a weight of 0 for P-phases and 1 for S-phases, and the other with a weight of 1 for P-phases and 2 for S-phases. From the initial dataset, we selected earthquakes that had at least 6 P and 3 S phase reads, azimuthal gaps smaller than 200, and RMS time residuals smaller than 0.2 s. This selection provided a refined dataset consisting of 4,161 earthquakes from the original 4,710 events of the INGV-OV database, with a final number of 41,272 P-wave readings and 14,683 S-wave readings.

As previously mentioned, we began by applying the SimulPS14 code, using the 1D Vp model from Vanorio et al. [9] as the starting model (hereafter referred as SimulPS-Vanorio). After carefully evaluating the 3D results with Vp/Vs ratios ranging from 1.7 to 1.78, we chose to use a constant Vp/Vs ratio of 1.7, which minimizes the residual. In the second step, the resulting Vp, Vs, and Vp/Vs models were used as starting points for further refinement with the tomoDDPS code (hereafter referred as tomoDDPS-SimulPS-Vanorio scheme). Furthermore, we also performed an inversion with tomoDDPS using the 3D velocity models from Battaglia et al. [11] (hereafter referred as tomoDDPS-Battaglia scheme) and compared the results obtained from these two final velocity models.

In the tomographic inversions using both SimulPS-Vanorio and tomoDDPS-SimulPS-Vanorio to resolve subsurface velocity structures, the inversion parameterized the study area with a grid of nodes spaced 1 km horizontally and 0.5 km vertically, extending to a depth of 8 km. Conversely, the tomoDDPS-Battaglia inversion used a finer grid spacing of 0.5 km in the longitude, latitude, and depth directions, providing a more detailed resolution.

To avoid overfitting and ensure stable inversion results with SimulPS14, we applied appropriate regularization techniques and smoothing constraints, controlling the trade-off between data fit and model smoothness. We chose the more suitable parameterization after some trial inversions, looking at the resulting model, the final data misfit and the resolution matrix. To find the best damping parameter, we ran some single-iteration inversions by varying the damping until the correct balance of misfit reduction and model complexity was found. Then a damping value of 30 for Vp and 25 for Vp/Vs ratio was chosen.

In the tomoDDPS inversion process, there are several important parameters that need to be determined: the weights for each phase, the weighting scheme for the absolute and differential catalog data, and damping. The phase weights used in this study were 1 for the P phase (0 in phase picking) structures and 0.5 (1 in phase picking) for the S phase structures. The weight of the P phase was set higher than that of the S phase because determining the S arrival time is generally more challenging and involves greater uncertainty compared to the P arrival time. For the weightings of the catalog and differential data, we gave a greater weight to the absolute data at the beginning of the inversion to obtain the velocity structure globally and at the end of the inversion a greater weight for differential data was employed to obtain the velocity structure in more detail [24]. In this study, the inversion process initially applied a weighting factor of 100 to the absolute catalog data in comparison to the differential data. Assessment of the optimal damping parameter was made based on the condition number (CND). The CND value was the ratio of the largest to the smallest eigenvalues which is between 80 to 100.

By inverting the 2020–2024 dataset using the tomoDDPS-Battaglia scheme, we processed 41,272 P-phases and 14,683 S-phases, along with a total of 388,166 P and 107,281 S differential times. After 11 iterations, we obtained the P and S velocity structures and determined the earthquake locations. For the unweighted catalog phases, the starting RMS value of 0.094 s decreased to 0.064 s, representing an approximate reduction of 32%. Conversely, for the weighted catalog, the starting RMS was 0.232 s, which was reduced to 0.043 s, resulting in an 82% decrease. Comparable results were achieved with the tomoDD-SimulPS-Vanorio scheme. For the unweighted catalog phases, the starting RMS of 0.082 s decreased to 0.064 s, yielding an approximate reduction of 22%. For the weighted catalog, the starting RMS value of 0.142 s was reduced to 0.043 s, representing a decrease of 70%.

3. Results

3.1. The 2020–2024 Inversion

For the entire dataset spanning 2020–2024, we initially compared the results from tomoDDPS-SimulPS-Vanorio and tomoDDPS-Battaglia. Although some details differ from each other, the main features of the tomographic images obtained are similar. Given the strong similarities between the two velocity models, we focused specifically on the tomoDDPS-Battaglia results. This choice was further justified by the following two reasons: 1) the better accuracy of Battaglia’s models as starting models, particularly in terms of the Vp structure, which benefited from denser sampling due to the SERAPIS active experiment, and 2) the greater flexibility of tomoDDPS in selecting initial velocity models, including the option to start from complex 3D models derived from other studies, which can lead to more accurate final models [27]. Additionally, using Battaglia’s models allowed us to integrate our tomography with external regions where the resolution was lower, enabling to show a more accurate representation of the subsurface structure of the Campi Flegrei Caldera, including offshore areas. However, the checkerboard and resolution tests show that owing to the distribution of the seismic sources and receivers, the model can be estimated correctly in a limited volume (Figures S3 and S4, see Supplementary material for further details). Therefore, our interest will be focused in the central part of CFC.

Figure 2. Maps showing the final solutions for Vp (a), Vs (b) and Vp/Vs (c) models. Relocated 3D earthquakes (black circles) occurring within ±0.25 km from each layer are depicted. S indicates the Solfatara area. The blue continuous lines delineate the well-resolved areas. Outside of these regions, the velocity models from Battaglia et al. [11] are included.

Figure 2. Maps showing the final solutions for Vp (a), Vs (b) and Vp/Vs (c) models. Relocated 3D earthquakes (black circles) occurring within ±0.25 km from each layer are depicted. S indicates the Solfatara area. The blue continuous lines delineate the well-resolved areas. Outside of these regions, the velocity models from Battaglia et al. [11] are included.

Figure 2a-c show the Vp, Vs and Vp/Vs tomographic images and 3D earthquake locations obtained for the period 2020–2024 at various horizontal depths.

According to synthetic tests and resolution criteria (Figures S3 and S4, see Supplementary material for further details), the 3D Vp, Vs and Vp/Vs structures are well-resolved in the central part of the study area down to a depth of about 4.0 km. The outline of the inferred region with reliable results is marked by a thick solid blue line (Figure 2a-c). In regions outside the area marked by the blue line, where our model resolution is poor, the displayed velocity models are those from Battaglia et al. (2005). To allow a direct comparison with the tomographic images of Vp models presented in Vanorio et al. [9] and Battaglia et al. [11], we maintained the color scales used in their figures, with only slight modifications.

The tomographic images shown in Figure 2a–c and Figure 3a–c, and the 3D representation of the Vp/Vs ratio in Figure 4 were obtained through 3D interpolation using Voxler software (Golden Software) with the Inverse Distance to a Power method.

The 3D approach, compared to 2D layer interpolation, typically presented in some tomographic studies, generally provides a more accurate and detailed representation of subsurface structures. Indeed, by accounting for variations in all three spatial dimensions, 3D interpolation can better capture the complexity of the subsurface, including vertical changes and subtle structural features that may be missed or not visible by 2D interpolation methods.

Key features of our Vp model include:

•

a wide relative low P-velocity anomaly (2.0 < Vp ≤ 4.0 km/s) in the central part of the caldera, mainly extending inland and down from approximately 0.5 km to 3.0 depth (L1 in Figure 2a). Between 2.5 and 3 km depth, this anomaly takes on an almost L-shaped aspect.

•

a relative low Vp (4.0 < Vp ≤ 5.0 km/s) region located at a depth of 3.0 and 4.0 km (L2 in Figure 2a), showing an almost elliptical shape and encircling a small high Vp anomaly (5.5 < Vp ≤ 6.0 km/s) at 4.0 km depth. This elliptical anomaly appears to delimit the morphological boundary of the internal resurgent dome [16], within which deeper earthquakes are mainly concentrated.

•

some portions of the arc-like high P-velocity anomaly (3.0 < Vp ≤ 5.0 km/s), previously identified in earlier tomographic studies [7,8,9,11] (H1 in Figure 2a). This anomaly, attributed to the buried rim of the Campi Flegrei outer caldera, is primarily located along the southern and western borders of the Gulf of Pozzuoli; it extends from about 0.5 km to approximately 2 km in depth, with seismic velocities increasing with depth. It has been interpreted as consolidated lavas and/or tuffs and interbedded lava sequences that signify the buried rim of the caldera [7].

Regarding the Vs structure, the most notable features of the model include:

•

a wide relative high-Vs anomaly (2.5 < Vs ≤ 3.0) located inland at about 2 km depth, between Pozzuoli and Campi Flegrei, coinciding with the region where a great part of earthquakes occurred (HS in Figure 2b).

•

a relatively low-Vs anomaly (2.0 < Vs ≤ 3.0 km/s), located at a depth of 3 to 4 km (marked as LS in Figure 2b and Figure 3c), is primarily situated offshore and south of Pozzuoli.

▪

Analyzing the Vp/Vs structure, we obtain slightly more conservative absolute values of Vp/Vs with respect to the previous studies. They range from 1.6 to 2.5, suggesting that direct inversion of Vp/Vs from S-P data is probably more stable than deriving Vp/Vs from separate Vp and Vs models.

▪

Horizontal and vertical cross-sections of the Vp/Vs ratio (Figure 2c and Figure 3c) reveal:

•

two high Vp/Vs (≥ 1.9) volumes (H2 in Figure 2c and Figure 3c) are located in the shallower central part of the caldera, at depths ranging from approximately 0 to 1.5 km, surrounded by regions with normal to low Vp/Vs. The larger of these volumes lies beneath the inland area of Pozzuoli and partially extends offshore.

•

an extended, thick high Vp/Vs layer (H3 in Figure 2c and Figure 3c), with values ranging from 1.9 to 2.5, located beneath a large portion of the caldera at depths of approximately 1.0 to 2.5 km. In the caldera’s center, this high Vp/Vs layer becomes discontinuous.

•

a volume with a low Vp/Vs ratio with values ranging from 1.6 to 1.66 that rises from depth in the central part of the caldera and reaches approximately 1.5 km b.s.l. (L3 in Figure 2c and Figure 3c). This structure widens as it extends downward merging into an almost continuous layer at 3.5 km, that likely extends beyond 4.0 km depth.

The 3D representation of these low- and high-Vp/Vs volumes (Figure 4) in the central part of the caldera provides a clearer view of their shape and distribution at depth . Figure 4 displays the isovalues for both low Vp/Vs (≤ 1.66) and the shallower high Vp/Vs (≥ 2.1) structures. The low Vp/Vs structure appears as a complex, branching volume extending from approximately 1.5 km to 4 km depth. The two high Vp/Vs volumes in the central part of the caldera, situated between 0.2 km and 1.5 km depth, are positioned just above the branched low Vp/Vs region.

Figure 4. 3D representation in the central part of CFC of high and low Vp/Vs volumes based on their isovalues. The isovalues ≥ 2.1 are related to the higher Vp/Vs volumes, whereas isovalues ≤ 1.66 are associated the lower Vp/Vs volumes. 3D hypocentral earthquake locations (red circles) are also shown.

Figure 4. 3D representation in the central part of CFC of high and low Vp/Vs volumes based on their isovalues. The isovalues ≥ 2.1 are related to the higher Vp/Vs volumes, whereas isovalues ≤ 1.66 are associated the lower Vp/Vs volumes. 3D hypocentral earthquake locations (red circles) are also shown.

It is important to note that Figure 4 does not display the very shallow low Vp/Vs anomalies at depths ≤ 0.5 km, nor does it include the thick layer with high Vp/Vs values (between 1.0 km and 2.5 km depth) lacking in the central part of the calder and surrounding the two central high Vp/Vs volumes.

We believe that the branching low Vp/Vs volumes should represent regions where fluid and gas phases ascend toward the surface, a portion of which contributing to the Solfatara-Pisciarelli shallow hydrothermal system (Figure 3c). These findings are consistent with the interpretations of Vanorio et al. [9] and Battaglia et al. [11], who associate the low Vp/Vs anomaly observed beneath the caldera center, between 2 km and 4 km depth, and below 3.5 km in the surrounding areas, with overpressured water and/or gas-bearing rocks (thermo-metamorphic, Mesozoic carbonate and crystalline rocks; see Figure 3 in [30]) at supercritical conditions.

In contrast, the high Vp/Vs thick layer between ca. 1.0 and 2.5 km, and the two shallower high Vp/Vs anomalies at around 0-1.5 km depth in the central part of the caldera, likely indicates the presence of brine-saturated rocks pre- and post-caldera volcano and marine deposits (Figure 3 in [30]), possibly resulting from steam condensation due to lower temperatures measured at shallow depths [31].

A more detailed discussion on these observations is provided in the paragraph 4.0.

Temporal Variation of Vp, Vs and Vp/Vs (2023–2024 vs. 2020–2022)

4D tomographic inversion aims to detect changes in subsurface properties over time by comparing seismic data from different periods. By consistently using the same starting velocity model across multiple inversions, the differences observed in the final velocity models are more likely to reflect true temporal variations in the subsurface rather than artifacts introduced by variations in the initial conditions. This approach reduces a lot but does not zero the bias that might arise if different starting models are used, as variations due to the choice of different initial models could overshadow or distort the actual changes occurring in the subsurface [21]. Thus, maintaining a consistent starting velocity model improves the stability and comparability of the inversion results, ensuring that the detected changes likely can be primarily due to physical changes in the subsurface.

Given the significant increase in seismicity observed since 2023 (Figure 1b), we divide the 2020–2024 dataset into two sub-periods and compared the resulting tomoDDPS-Battaglia models obtained separately for the 2020–2022 and 2023–2024 to identify potential velocity changes. To verify the resolutions of the two seismic tomographic models, we also conducted resolution tests for each period as for 2020–2024 (Figures S5–S9, see Supplementary material).

First, we compared the images obtained from the two separate inversions and even though some details differ the main features of the resulting tomographic images remain consistent. In Figure 5 a,b, the temporal variations in velocity, obtained by subtracting the velocity values of the 2023–2024 model from those of the 2020–2022 model, are shown.

The observed differences highlight the evolution of velocity structures during 2023-2024, revealing areas where seismic velocities have increased or decreased. Specifically, both Vp and Vs show a broad, ring-shaped decrease, generally ranging from -2.5% to -10%. This decrease is particularly prominent at depths between approximately 2 km and 4 km, becoming more evident below 3 km, with some zones inside the ring exhibiting reductions also exceeding -10%.

For Vp, the region of decreased velocity surrounds, at depths of 3–4 km, a well-defined volume oriented roughly NW-SE, where seismic velocities increased by +2.5% to +10% (U1 in Figure 5). This increase is located between Pozzuoli and the adjacent offshore region, with the maximum increment (+10%) occurring between 3.5 and 4.0 km depth.

In contrast, the zone with increased Vs (+2.5% to +10%) (U2 in Figure 5) lies mainly at a depth of 2–3 km, with the maximum increment (+10%) occurring at ca. 2 km, primarily concentrated near and beneath the Solfatara area.

Considering the Vp/Vs ratio, the most significant variation (ranging from -2.5% to -10%) is observed between ca. 1.0 km and 3.0 km depth, with a maximum reduction of -10% occurring at approximately 2.0 km depth, almost in correspondence with the Vs increase (D1 in Figure 5). This observation might align well with seismic changes expected during the depletion of liquid pore fluids and their replacement by steam, a process observed to be progressive and extreme in other geothermal systems such as The Geysers (e.g., [32]). Instead, the observed deeper patterns of Vp and Vs increase at 3–4 km may suggest the presence of an inflating source at depth, due to a small intrusion of magma enriched in supercritical fluids. This intrusion likely influences the surrounding crustal structure through stress changes and fluid migration that reduce seismic velocities..

A more detailed discussion on these observations is provided in the paragraph 4.0.

4. Discussion

In this study, seismic imaging of the shallow volcanic structure beneath the central part of the Campi Flegrei Caldera (CFC) for the period 2020–2024 yields more refined Vp, Vs, and Vp/Vs models compared to previous tomographic investigations, reaching depths of approximately 4 km. The following discussion integrates these findings with recent geochemical and geodetic studies and compares the observed velocity variations in 2023–2024 to assess whether the ongoing unrest—accelerating since 2023—is primarily driven by magma intrusion or hydrothermal processes, a key challenge in understanding the current dynamics of Campi Flegrei. The volcanological literature on CFC presents contrasting interpretations, with several studies attributing the unrest to shallow magma intrusions and many others to hydrothermal perturbations associated with the injection of deep magmatic gases (e.g., [33]).

The Vp and Vs velocity models presented here enhance our understanding of the subsurface structures (Figure 2a,b and Figure 3a,b), though they do not provide clear evidence of low-velocity perturbations indicative of a distinct magma intrusion. The only exception is a relatively not well-confined low-Vs anomaly, located at a depth of 3–4 km (marked as LS in Figure 2b and Figure 3c), offshore and south of Pozzuoli. Differently, the Vp/Vs ratio has provided additional constraints for the identification of the causes of some observed seismic velocity variations (Figure 2c, Figure 3c and Figure 4).

It is well known that seismic velocities generally increase with depth due to the progressive compaction and densification of rocks under increasing overburden pressures. They are also influenced by multiple factors—including rock composition, porosity, fracturing, minerals, fluid saturation, temperature, and pressure—that can alter elastic properties and complicate the interpretation of seismic anomalies. Additionally, it is important to consider that ray-based seismic imaging methods may underestimate low-velocity anomalies and melt fractions, because seismic waves preferentially travel through high-velocity regions, potentially bypassing low-velocity zones [34]. This limitation is particularly significant when dealing with small magma intrusions (<1 km³) or when melt is distributed in thin, interconnected pockets rather than forming a large, coherent body [35].

The analysis of the compressional-to-shear wave velocity ratio (Vp/Vs) has proven particularly useful for constraining the distribution of fluids in volcanic and geothermal environments. In saturated or partially saturated rocks, fluid content and phase significantly influence P-wave velocities while having little effect on S-wave velocities [9]. Fluid phase transitions alter fluid compressibility and, consequently, the bulk modulus [36,37], while shear moduli remain largely unchanged. As a result, low Vp/Vs ratios are generally associated with gas-bearing rocks, which exhibit high fluid compressibility, whereas high Vp/Vs ratios indicate liquid-bearing rocks [9].

One of the most significant features of our tomographic results is the 3D reconstruction of the low Vp/Vs volumes (<1.66) beneath the central caldera (Figure 4). This 3D body appears as a broad, interconnected network of gas-bearing rocks extending from depths of 3.5–4.0 km up to approximately 1.5 km below the fumarolic areas of Campi Flegrei (Solfatara-Pisciarelli). We infer that this complex structure may represent the degassing pathway linking a deeper source to the surface fumarolic areas [38].

Another notable feature is the thick, extensive high Vp/Vs layer (values ranging from 1.9 to 2.5) located between 1.0 and 2.5 km depth beneath much of the caldera (Figure 2c and Figure 3c). This layer is attributed to liquid-bearing rocks and brine-saturated zones resulting from steam condensation. In the central caldera, this high Vp/Vs layer becomes discontinuous, and at shallower depths (0–1.5 km) two distinct high Vp/Vs zones emerge. The larger of these is situated beneath the inland area of Pozzuoli and extends partially offshore [11].

Geothermal drilling at Mofete and San Vito (Figure 1; [31,39]) confirms that volcanic sediments infiltrated by meteoric water form a porous, fluid-saturated matrix. These wells reveal a multi-reservoir geothermal system extending from 550 to 2700 m depth, characterized by high-salinity fluids, temperatures of 250–390 °C, and distinct fluid origins—seawater-derived at depths <2000 m and meteoric-magmatic below.

These high temperatures of approximately 400°C favor the decarbonation of deep limestone rocks, releasing additional CO₂ and contributing to the intense fumarolic activity observed at the surface in areas such as Solfatara and Pisciarelli.

Deep hot brines ascend through faults and fractures, mixing with shallow groundwater [40], while volcano-tectonic discontinuities facilitate the upward migration of fluids, including deep CO₂-rich mineralized fluids. Exceptionally high geothermal gradients measured in these wells [39,41] suggest a direct link between geothermal anomalies and hydrothermal activity. The interaction between cooler water and a deeper, hotter gas/steam reservoir may trigger rapid phase transitions (flashing), generating steam.

In the sketch model of Figure 6, we overlay the temperature distribution from the 3D Petrillo et al.’s model (Figure 10 in [38]) onto the W-E image of Vp/Vs structure. It is possible to observe a strong correlation with the boundaries of both high and low Vp/Vs regions, particularly in the western and central caldera. Additionally, the sketch model reports the most significant features obtained from comparing P- and S-wave velocity and Vp/Vs ratio (Figure 6) between the 2023–2024 and 2020–2022 periods.

The observed velocity variations at 3–4 km depth beneath the Pozzuoli area and its offshore (a ca. NW-SE elongated volume) indicate a significant increase in Vp (+2.5% to +10%) along with a minor increase in Vs, and a slight decrease in the Vp/Vs ratio. Extensive interdisciplinary research supports that tectonic structures oriented NW–SE and WNW–ESE may significantly influence recent unrest at Campi Flegrei (e.g., [33] and reference therein).

Two possible scenarios may explain these velocity variations:

1) A small magma intrusion enriched in supercritical fluids. In this case, the intrusion would stiffen the surrounding rock, leading to a significant increase in Vp, a modest increase in Vs, and a slight decrease in Vp/Vs. This behavior aligns with observations that high-pressure, high-density supercritical fluids reduce pore space and enhance rock stiffness [9,35].

2) The accumulation of pressurized fluids or high-density, mineral-saturated brines. These fluids similarly stiffen the rock matrix, increasing Vp, though their geochemical signatures differ from those associated with magma intrusions [42].

In both scenarios, the pressurized source influences the surrounding crust through stress changes and fluid migration, creating favorable conditions for rock fracturing and faulting. This is evidenced by the surrounding zones of Vp decrease (-2.5% to -10%), which can be linked to fluid-induced weakening and microfracturing in the adjacent rocks

At shallower depths (2–3 km), a significant increase (+10%) in Vs is observed, encircled also in this case by broader zones of Vs reduction (-2.5% to -10%). A decrease in Vs is typically associated with microfracturing, as the formation of microcracks lowers the shear modulus and, consequently, S-wave velocity. Additionally, fluid saturation reduces effective stress and acts as a lubricant within the rock matrix, further decreasing the shear modulus and Vs [42,43,44].

These observations support the interpretation that regions of reduced Vp and Vs during the 2023–2024 period are more damaged and less rigid, likely due to increased fracturing and fluid infiltration. Additionally, at shallower depths (2–3 km) beneath the Solfatara-Piscirelli area, the region of Vs increase coincides with a decrease in the Vp/Vs ratio. This may indicate the presence of relatively intact, stiff rock in which steam has gradually replaced liquid water under high-pressure and high-temperature conditions, thereby enhancing surface emissions. The resulting steam condensates can form brines and migrate upward along fracture networks, potentially explaining the intensification of fumarolic activity observed at the surface. Such seismic changes are similar to those documented in other geothermal systems, such as The Geysers [32], where the progressive replacement of liquid pore fluids by steam is associated with marked seismic anomalies.

Overall, these findings underscore the complex interplay between fluid migration, phase transitions, and rock mechanics in geothermal and volcanic settings, offering critical insights into the subsurface processes driving unrest at Campi Flegrei. Recent studies, focused on the ongoing unrest, offer varied insights.

Recently, Caliro et al. [45] refined some previous geochemical models [46,47,48] by integrating hydrothermal and magmatic contributions. They propose that the increased H₂S emissions observed from 2018 to 2022 may result from decompression-driven degassing of mafic magma at ~6 km depth, coupled with sulfur remobilization within the hydrothermal system. However, they caution that geochemical data alone cannot precisely constrain the pressure or depth of the magma source.

Anyway, nearly all recent geochemical studies consistently indicate that the hydrothermal system is responding to a deep pressurization source, with magmatic fluids migrating upward and reaching the hydrothermal system. This process induces substantial geochemical and geophysical modifications, including alterations in the gas compositions of fumarolic emissions and progressive changes in the mechanical properties of the crust as detected by seismic and geodetic measurements.

Ground deformation studies [49,50,51] consistently identify a primary magmatic pressure source between 3 and 4 km depth. Tizzani et al. [49] use multi-platform, multi-frequency InSAR data (Cosmo-SkyMed and Sentinel-1) from 2011–2022 to identify an inflating source at 3–4 km depth, interpreted as a pressurized magmatic intrusion. Giudicepietro et al. [50] analyze GNSS and Multi-Temporal DInSAR data from 2015 to 2023, highlighting a broad, bell-shaped deformation affecting the entire caldera (evident since 2021) interpreted as either a sill intrusion or a magmatic fluid accumulation zone inflating at ~3.8 km depth. Both Tizzani et al. [49] and Giudicepietro et al. [50] also identified a very shallow secondary source at 400-500 m of depth located beneath the Campi Flegrei that may be feeding fumaroles in the area. Astort et al. [51] analyze surface deformation at CF from 2007 to 2023, suggest the action of both a shallow and a deep source. They identify a deformation source that progressively widens and migrates upward from 5.9 km to 3.9 km depth, with inflation linked to the ascent of 0.06–0.22 km³ of magma from the deeper source (≥ 8 km), which undergoes minor deflation.

Therefore, the tomographic results presented here well agree with ground deformation studies and modeled source depths. However, even though the 2023–2024 temporal variations in Vp, Vs, and Vp/Vs allowed us to identify a deeper primary pressurized source at 3–4 km depth beneath Pozzuoli (and in its offshore area) and a secondary pressurized source at approximately 2 km depth—likely feeding the Solfatara-Piscirelli fumarolic field—we do not have clear evidence of a magma intrusion. Indeed, if a deeper magma intrusion enriched in supercritical fluids were present, the expected low Vp/Vs anomaly is indistinguishable from that of the surrounding rock, which already exhibits values below 1.66.

5. Conclusions

Distinguishing between magma intrusions and hydrothermal effects remains a critical challenge for understanding unrest at Campi Flegrei. Past interpretations have alternately attributed unrest to shallow magma intrusions or to thermal-fluid disturbances in the shallow geothermal system. Our analysis of seismic data from 2020–June 2024 produced updated 3D models of Vp, Vs, and Vp/Vs down to approximately 4 km depth, consistent with previous tomographic studies. Although these models did not reveal clear low-velocity anomalies typical of magma intrusions—except for a localized, poorly confined low-Vs anomaly at 3–4 km depth offshore and south of Pozzuoli—the comparison of tomographic images from 2023–2024 and 2020–2022 shows significant velocity variations. In particular, we observed a notable increase in Vp at 3–4 km depth beneath Pozzuoli and its offshore area, accompanied by a substantial increase in Vs at shallower depths (2–3 km) and a decrease in the Vp/Vs ratio at around 2 km beneath the Solfatara-Piscirelli area. Both the zones of increased Vp and Vs are surrounded by broad regions of velocity decrease.

One plausible interpretation is that the Vp increase at around 3–4 km depth reflects the effect of a deeper, roughly NW–SE-elongated pressurized source that stiffens the rock— either a magma intrusion enriched in supercritical fluids or an accumulation of pressurized, high-density fluids. This stiffening results in higher P-wave velocities relative to the surrounding, more porous, fluid-saturated rock.

These findings corroborate recent ground deformation analyses, which suggest the presence of a pressurized source between 3 and 4 km depth. In particular, Astort et al. [51] indicate that the inflation source at 3.9 km depth is linked to the ascent of 0.06–0.22 km³ of magma. This represents a very small volume (< 1 km³) that seismic tomography might fail to detect due to its resolution limitations at present, being influenced by factors such as the spacing between seismic sources and receivers, as well as the adequacy of ray coverage. Moreover, if the magma intrusion is enriched in supercritical fluids, the expected low Vp/Vs anomaly may be indistinguishable from that of the volume located at approximately 1.5 km b.s.l. in the central part of CFC, which merges into an almost continuous layer at 3.5–4.0 km depth (Figure 4) and exhibits a low Vp/Vs ratio (1.6–1.66).

Instead, the upward migration of fluids from this deeper source may cause the observed increase in Vs at around 2 km depth, indicating a separate and pressurized process within the shallow hydrothermal system. For instance, the shallower zone—characterized by increased Vs and a decreased Vp/Vs ratio—may indicate an area where, under high-pressure and high-temperature conditions, steam has gradually replaced liquid water, thereby contributing to the active feeding system of the Solfatara-Piscirelli fumarolic field. Meanwhile, the surrounding rock appears weakened by microfracturing, likely promoting increased seismicity.

These depth-related differences likely reflect the interplay between deep pressurization effects and shallow fluid phase transitions, which together produce distinct seismic signatures in Vp and Vs. Overall, our findings confirm recent ground deformation studies and modeled source depths associated with the ongoing unrest and underscore the complex interaction between deep magmatic processes and shallow hydrothermal dynamics. Fully resolving the origin of this pressurized source—and understanding its impact on the geothermal system—requires a robust, integrated approach that combines geophysical, geodetic, geochemical, and petrological data into a unified model for a more precise interpretation of these subsurface processes. Furthermore, we believe that additional multidisciplinary monitoring efforts and targeted investigations—such as a new active seismic experiment—are necessary. Considering that the unrest is still ongoing, the scientific community, local authorities and residents, are increasingly concerned, as a new eruptive episode could necessitate the evacuation of several hundred thousand people. Moreover, even without an eruption, the observed increase in seismicity could have significant consequences for such a densely populated area.