Tracking the 2007–2023 magma-driven unrest at Campi Flegrei caldera (Italy)

To assess eruption probability during volcanic unrest, it is essential to determine the distribution and migration of magma within the shallow plumbing system and track its evolution. Key observables are the depth, size, and connectivity of melt pockets and transfers among them, with the goal of identifying any ascent of shallow magma. This is especially critical for long-lasting unrest beneath densely populated areas like Campi Flegrei (CF), Italy, where surface deformation is widely monitored and inverted to infer source location and size. A major limitation is the uncertain nature of pressure sources (magmatic vs. hydrothermal), especially for intermediate depths (3–4 km) such as those inferred at CF. After the last eruption in 1538 CE, CF experienced three 20th-century unrest episodes (1950–52, 1969–72, 1982–84) with up to ~1.8 m uplift at Pozzuoli. Since 2005, uplift has resumed, reaching ~119 cm by Dec 2023, accompanied by increasing seismicity and intensified degassing at Solfatara-Pisciarelli. Petrological and geochemical data suggest a deep magmatic system at ≥8 km. The research question here is to constrain the location, size evolution, and, crucially, the nature (magmatic vs. hydrothermal) of the shallow deformation source and its relation to the deeper reservoir by merging geodetic/numerical modeling with petrological calculations.

Previous geodetic inversions during CF unrest (including 2007–2023 subperiods) consistently locate the main deformation source deeper than 3 km and often with sill-like geometry. Modeling in homogeneous media tends to yield shallower depths (~2.5–3 km) than heterogeneous models (~3.6–5 km). A deeper reservoir at ≥8 km has been less frequently incorporated, although some studies suggest magma accumulation at >7 km contributing to deformation since at least 2015. The nature of the ~4 km source has been debated (magmatic vs. hydrothermal), with multidisciplinary geophysics/geochemistry broadly agreeing on a magmatic deeper source (~8 km), but the shallow source nature remained unresolved. A recent 4D tomography indicates a 4.5–5 km magmatic body beneath the caldera center since 2019. Historical/archaeological records imply deep reservoir pressurization around the 1538 eruption. Differences among prior studies reflect varying datasets (e.g., Sentinel-1 vs. COSMO-SkyMed), processing, modeling assumptions (homogeneous vs. heterogeneous crust), and number/type of sources.

Data: The study uses GNSS (21 inland stations; 2 seafloor stations in the Gulf of Pozzuoli) and multi-temporal InSAR from ENVISAT (2007–2011) and COSMO-SkyMed (2011–2023), validated against independent geodetic data (RMSE < 5 mm/yr). Based on velocity changes in RITE time series, seven intervals were modeled: T1 (2007–2011), T2 (2011–2013), T3 (2013–2015), T4 (2015–2017), T5 (2017–2020), T6 (2020–2022), T7 (2022–2023). Seafloor GNSS displacements were included in T5 (2017–2020). InSAR data were spatially downsampled (300 m within ~5 km radius of maximum deformation; 600 m outside). Finite Element Model (FEM): The 3D elastic domain (140×140×60 km) comprises ~150,000 8-node brick elements, with finest resolution (400 m) in the central 15×15×11 km. Free surface is flat given mild topography. The elastic structure is heterogeneous to 5 km based on integrated Vp and Vp/Vs tomographies (active and passive), converted to rigidity (μ) and Poisson’s ratio (ν). Below 5 km, a homogeneous half-space is used (μ=20 GPa, ν=0.25). Source representation: A two-source plumbing system is modeled: (1) a shallower source represented implicitly as a point-source equivalent moment tensor (only diagonal components used; U11 solutions), located horizontally ~600 m S of RITE and allowed to vary in depth per period; (2) a deeper source at 8 km depth represented as a 10×10 km tabular crack under uniform overpressure/underpressure. Potential shallow element-sources are a grid of 1000 cubic elements (400 m side) centered beneath Pozzuoli. Green’s functions were computed by applying unitary dipoles/double couples to each element face (6000 displacement solutions sampled at 4537 surface nodes). The deep source’s surface displacement response for 1 MPa was precomputed; its contribution is scaled in inversion. Mechanical interactions between sources are neglected. Inversion: A two-step framework per period: (i) global optimization (Neighbourhood Algorithm) sampling 51,000 models over 50 iterations, minimizing a combined misfit to GNSS (weight 1/3) and InSAR (ascending/descending) velocities; (ii) Bayesian appraisal (MCMC ~40 chains × 3000 iterations) to obtain parameter marginal distributions and mean models. Free parameters: shallow source depth and three diagonal stress scale factors; deep-source pressure scale factor. Source interpretation and volume change: Principal stress ratios from inversions are mapped to triaxial ellipsoid aspect ratios (b/a, c/a) via grid search using analytical relationships, yielding ΔP/σ3 as a function of shape. Volume change is related to overpressure and elastic parameters. The shallow source is treated as a point ellipsoidal cavity; absolute volumes and volume changes require an initial volume estimate. Based on similarity to 1982–84 unrest and cumulative uplift since 1950, the initial shallow source volume was set to 0.1 km³; subsequent volumes scale from this baseline. The deep-source volume change is obtained by integrating relative node displacements on the two faces of the tabular crack; the unit pressure response gives ΔVdeep=20.8×10^6 m³ per MPa, scaled by the inverted pressure. Strain-field analysis: Forward FEM simulations were run incrementally to the end of each period using the retrieved sources to compute the equivalent strain ε* (J2 norm of strain) at selected depths (700 m, 1400 m) and compare with seismicity distributions. Petrological calculations: Equations of state for H2O-rich magmatic fluids were used to compute molar volumes at 210 MPa (~8 km) and 100 MPa (~4 km). Several scenarios were tested to explain the observed shallow inflation and deep deflation: A) only fluids rising and expanding; B) coupled magma+fluids transfer equal to deep deflation; C) second boiling (crystallization-driven exsolution) at 4 km; D) magma ascent from 8 km releasing excess volatiles during decompression that accumulate at/near 4 km; E) direct magma plus volatiles transfer to 4 km equivalent to the shallow inflation, with deeper reservoir recharge from below. Thermal diffusion timescales for cooling were estimated for spherical and sheet geometries to evaluate scenario C plausibility.

- Surface deformation (2007–2023) is explained by two sources: a shallower source that shallowed and widened over time, and a deeper source at 8 km that underwent persistent but limited deflation. - Shallow source depth evolution: 5.9±0.6 km (T1, 2007–2011) → 5.5±0.2 km (T2) → 5.1±0.6 km (T3) → 3.9±0.2 km (T4–T7, 2015–2023). Depth is least constrained in T1 and T3 due to low signal-to-noise. - Shallow source geometry varies through time (triaxial ellipsoids, sill-like in T2), with major axis generally E–W. From 2015 onward, the source behaves as a thick oblate/spheroidal body with volumetric expansion at constant depth (~3.9 km). - Shallow source volume increased from ~0.10 km³ (2007) to ~0.16 km³ (2023), i.e., net ΔV > 60×10^6 m³ over 16 years. Periodic volume change rates increased from ~1.4×10^6 m³/yr (T1) to ~8–9×10^6 m³/yr (T7). - Deep source (8 km) shows persistent deflation with volume-change rates of order −10^6 m³/yr in some intervals; the largest single-interval change was −1.3×10^6 m³ over 1.25 years (T7). Cumulative deep deflation is ~−5×10^6 m³ (2007–2023), an order of magnitude smaller than shallow inflation. - Seismicity (mostly <3 km depth) forms an annulus centered near Pozzuoli, concentrated beneath/NE of Solfatara. Modeled equivalent strain fields reproduce annular strain maxima at 700–1400 m depths, spatially consistent with seismicity clusters, linking long-term pressurization to seismicity patterns. - Petrological scenarios indicate that neither transfer/expansion of volatiles alone (A) nor coupled transfer of magma+fluids equal to the deep deflation (B) can reproduce the shallow inflation (volumetric ratios too small by factors >5 to >30). Cooling-induced exsolution at 4 km (C) requires unrealistic volumes and timescales (≥182–3000 years) and would not match the observed increasing inflation trend. - The only plausible mechanisms are: D) ascent of magma from ≥8 km releasing excess fluids during decompression, or E) direct transfer of magma+volatiles to ~4 km, with the deeper reservoir recharged from below. Calculations suggest 0.06–0.22 km³ of magma must have risen from ≥8 km to explain the shallow inflation, with the smaller deep deflation implying replenishment from deeper levels to maintain melt fraction. - Overall, magma ascent to depths shallower than 8 km is identified as the ultimate driver of the ongoing unrest, with the shallow source at ~3.9–5.9 km evolving in depth, shape, and volume across 2007–2023.

The study resolves the debated nature of the shallow deformation source at Campi Flegrei by integrating geodetic inversions with petrological mass-balance. The observed progressive shallowing and inflation of the ~central shallow source, together with the smaller-amplitude deflation at 8 km, cannot be explained by hydrothermal fluid migration or in situ second boiling alone. Instead, scenarios consistent with both datasets require direct involvement of magma rising from ≥8 km: either degassing during ascent with fluids accumulating at ~4 km (scenario D) or injection of magma with volatiles into the shallow source (scenario E). The deep-source deflation being much smaller than the shallow inflation necessitates recharge of the deeper reservoir from beneath, maintaining a high melt fraction. The strain modeling reproduces the annular seismicity pattern and its concentration near Solfatara, reinforcing the causal link between long-term source pressurization and seismicity, and highlighting the interaction between hydrothermal (≤3 km) and magmatic (≥4 km) systems. These findings constrain the plumbing system’s evolution and provide actionable insights for monitoring and hazard assessment by indicating that shallow inflation reflects magma ascent rather than purely hydrothermal processes.

- The 2007–2023 unrest at Campi Flegrei is driven by a two-source system: a shallow source that migrated from ~5.9 to ~3.9 km and inflated substantially, and a deeper magmatic reservoir at 8 km that deflated modestly. - Petrological modeling requires 0.06–0.22 km³ of magma to have ascended from ≥8 km to explain the shallow inflation, implying ongoing magma involvement and recharge from deeper levels. - The modeled strain fields align with observed seismicity, demonstrating a direct link between deformation, source pressurization, and shallow seismicity distribution, particularly near Solfatara. - Methodologically, combining heterogeneous-elastic FEM geodetic inversions with petrological simulations enables improved tracking of plumbing-system evolution and better discrimination of magmatic vs. hydrothermal drivers. - Future work should incorporate viscoelastic/thermo-poroelastic rheologies, better constrain shallow source geometry and orientation (including off-diagonal moment components), and leverage expanded offshore geodesy and continuous petrological monitoring to refine mass-transfer estimates and stalling depths. The framework can be applied in near real time to improve eruption forecasting at restless calderas.

- Mechanical assumptions: purely elastic modeling with heterogeneous structure only down to 5 km; homogeneous medium below may bias deep-source stress/displacement fields. Viscoelastic and inelastic processes over decadal timescales (10^18–10^19 Pa·s) could reduce inferred overpressures and volume changes and modify temporal patterns; poroelastic effects not modeled. - Source parameterization: shallow source represented as a point ellipsoid using only diagonal moment tensor components (strike/dip unresolved), and its horizontal position fixed; mechanical interactions between shallow and deep sources neglected. - Data coverage: limited offshore stations (two seafloor GNSS used) and submerged caldera hamper resolution of source orientation; lower signal-to-noise in T1 and T3 increases uncertainty on depth and shape. - Volume baseline: shallow-source absolute volumes depend on an assumed initial volume of 0.1 km³ (scaled from historical unrest), making absolute volumes minimum estimates; compressibility of exsolved fluids implies geodetic ΔV are lower bounds. - Petrological simplifications: first-order EOS-based calculations (focus on H2O) and end-member thermal models; actual volatile mixtures, melt compositions, and transient connectivity may alter precise mass/volume budgets; stalling depth of rising magma remains non-unique.