Introduction

One in ten of the world’s historically active volcanoes are large calderas more than 5 km across1,2. They frequently pass through several episodes of unrest over decades before they erupt1,2. The episodes belong to a single evolutionary sequence, so that the behaviour of each episode depends on the cumulative effect of its predecessors3,4. The non-eruptive episodes can generate mistrust among communities about warnings of eruption4,5. They can also produce structural changes in the crust that render obsolete conventional forecasts of eruption that assume episodes will repeat their previous behaviour.

The Campi Flegrei caldera in southern Italy is a remarkable example. About 12–15 km across, it is the largest active caldera in Europe and extends west from the outskirts of Naples to the Tyrrhenian Sea (Fig. 1). About a third is partially submerged beneath the Bay of Pozzuoli; the remaining two-thirds are home to more than 360,000 people6. The volcano has been restless since 1950. It last erupted in 1538 after an interval of about 3000 years7. Previous intervals have been as short as decades or centuries7, so that a return to eruption after nearly 500 years is a realistic possibility.

Fig. 1: Unrest at Campi Flegrei since 1950.
figure 1

a Ground movement has occurred across the central 80 km2 (blue ellipse) of the caldera (dashed outline). Most of the VT seismicity since 1983 (Fig. 3) has been focussed around the Starza marine terrace and fault system (white curve) between Monte Nuovo and the fumarolic fields of Solfatara and Pisciarelli, which are next to each other. b From the surface downwards, the crust beneath the central region of the caldera consists of predominantly pyroclastic volcanic deposits (primary and reworked) to depths of as much as 2.3 km, overlying marine silts and clays, which become thermally altered (thermo-metamorphosed) at depths of about 2.5–3.0 km and then grade (at unknown depth) into bedrock that contains solidified magmatic intrusions and extends to the top of a melt-rich layer about 8 km below the surface (magenta10). The melt-rich zone feeds magma and gas (notably CO2) to shallower depths. Measured temperatures of about 400 °C in the upper levels of the thermo-metamorphic horizon16,19 favour silica mineralisation and formation of low-permeability horizons (grey band)19. Uplift between 1950 and 1984 is consistent with intrusions of magmatic sills to depths of about 3 km (pale magenta)15. Pore pressures are expected to be hydrostatic above the thermo-metamorphic zone, but to reach higher values (that may approach lithostatic) where magmatic gas accumulates underneath48. Precipitation of Ca-rich minerals favours the formation of a stronger cap rock at depths less than 1.5 km18 (outlined by the blue dashes). Almost all the seismicity during the 1982-1984 unrest (red dots) occurred in the hydrothermal system between the cap rock and thermo-metamorphic horizons. Since 2005 (green dots), seismicity has also occurred through the cap rock in the vicinity of Solfatara and Pisciarelli (black triangle; see also Fig. 3) which may be favouring the escape of CO2 (black arrow).

Four episodes of uplift have raised the coastal town of Pozzuoli, near the centre of greatest movement, by more than 4 m (Fig. 2) and twice triggered evacuations of about 40,000 people8. The first three uplifts each occurred over two-year intervals between 1950 and 1984 at mean rates on the order of 0.1–1 m per year9. The fourth uplift has been continuing since 2005 at a mean rate on the order of 0.01–0.1 m per year—that is ten times longer and more slowly than its predecessors. Here we argue that the change in uplift behaviour has been caused by magmatic gas accumulating about 3 km below the surface and inducing slow rupture in the overlying crust. Once rupture is complete, it will offer magma and gas newly ascending from greater depth an easier access to the surface than previously available. Understanding the structural changes occurring today is thus essential for making pre-eruptive scenarios more reliable.

Fig. 2: Variations in uplift and VT seismicity between 1950 and 2022.
figure 2

a Time series for measured and corrected uplift at Pozzuoli (solid and dashed black curves) and monthly rates of seismicity (columns). After the VT event rate reached a peak in 1984, Campi Flegrei remained almost aseismic during three decades of slow subsidence and uplift. Persistent VT seismicity recommenced in 2015–2016, when the corrected uplift returned to its position in 1984 (magenta arrows). We initially interpreted the decay in VT rate after 1972 as evidence of partial relaxation of the crust4; however, some or all of the events may have been caused by the illegal fishing practice of detonating explosives to drive sea-floor creatures to the surface. b Between 2015 and 2020, VT event rates accelerated to a steady mean value of about 175 events per month. The acceleration began when the corrected uplift had been restored to its 1984 position (magenta arrow). Data from INGV-Sezione di Napoli (Osservatorio Vesuviano)71,72.

Caldera-wide ground movement and local seismicity are key measures of Campi Flegrei’s proximity to bulk failure4 (when the deforming crust breaks across all or most of its thickness). These measures describe changes in the stress applied to the crust by subsurface pressurisation. The volcano is underlain by a melt-rich zone which, with a top about 8 km below the surface10, has intermittently fed eruptions and intrusions of magma to depths of about 3–4 km11,12,13,14,15 (Fig. 1). Across the centre of the caldera, the top 2.0–2.3 km of crust consists of predominantly pyroclastic volcanic deposits, both primary and reworked. These lie on top of marine silts and clays, which become thermally altered (thermo-metamorphosed) at depths of about 2.6 km and then grade into bedrock.

The volcano has an active hydrothermal system to depths of about 2.5–3.0 km (Fig. 1). Meteoric steam and volcanic gases are persistently released from fumaroles; the best known of these are at Solfatara and Pisciarelli, which are less than 0.6 km apart and about 2.3 km east–north east of Pozzuoli, a coastal town near the centre of the caldera (Fig. 1). The permeability directing fluid flow broadly decreases with depth from the surface deposits to bedrock16,17. However, local combinations of rock composition, pressure and temperature are believed to have favoured the formation of two horizons with lower permeability than the levels immediately underneath: one at depths of less than 1.5 km, caused by the precipitation of Ca-rich minerals in volcanic tuffs18, and another at depths of about 2.5–3.0 km, where silica precipitates in clays and silts near the thermo-metamorphic horizon19 (following the classic Fournier model of hydrothermal systems20). The melt-rich zone below 8 km is persistently degassing, most notably carbon dioxide21,22. A proportion feeds surface emissions, while the remainder remains below the low-permeability horizon at 2.5–3.0 km19, which also coincides with inferred levels of sill emplacement11,12,13,14,15 (Fig. 1).

Historical ground movements have been largest close to Pozzuoli23,24,25. Since Roman times, they have been characterised by gentle subsidence at a mean rate of about 170 cm a century26, interrupted by intervals of net uplift during (1) the seven decades of ongoing unrest, (2) the 100 years before the 1538 eruption26,27, and, possibly, (3) a non-eruptive uplift in c. 600-700 AD26,28 (see Methods).

Uplift and subsidence since 1950 have decayed radially away from Pozzuoli12,13,24,25. The movements have developed in four stages. The first three lifted Pozzuoli by comparable amounts over similar durations9: about 74 cm in 1950–52, 159 cm in 1969–72, and 175 cm in 1982–84. The fourth stage has continued since the end of 1984 for nearly four decades, first with twenty years of subsidence that lowered Pozzuoli by 93 cm (at a mean rate about 2.6 times faster than the background value), followed by eighteen years of uplift that by February 2023 had raised the ground to almost 10 cm above its starting position (Fig. 2).

The uplifts have been accompanied by nearly 31,000 brittle-failure earthquakes, which are commonly described as volcano-tectonic (VT) earthquakes29, most of which have occurred since 1982. Although a small number have reached a magnitude of 4.0, 98% of the events have had magnitudes of 2.5 or less30,31. About 80% had epicentres along the onshore side of the coastline between Monte Nuovo and Solfatara—a distance of about 5 km running just north of Pozzuoli31,32 (Figs. 1, 3). The trend follows the fault system of La Starza (Fig. 1), where marine terraces show that a relative upward displacement of as much as 40 m occurred on the landward side between about 8000 and 4000 years ago33,34,35.

Fig. 3: Distribution of VT events 1982-1985 and 2005-2022.
figure 3

Seismicity has been concentrated in the vicinity of the Starza fault system (Fig. 1). a The 1982-1985 events occurred along the whole system, (b) at depths mostly between 1 and 3 km. c The 2005-2022 events have been restricted to a smaller volume at the eastern end of the system. d Most have occurred at depths less than 1.5 km, apart from an inclined trend 1–3 km beneath Solfatara; the same trend may also have been present in 1982–1984. The black lines show the same distances in each pair of time intervals. Data from INGV-Sezione di Napoli (Osservatorio Vesuviano)72.

The cumulative number of VT events is a proxy measure of co-seismic inelastic strain (see Methods) and has shown repeatable patterns with time and deformation at strato-volcanoes36 and large calderas3,4,36, as well as in scaled rock-physics experiments37. At Campi Flegrei, the number of VT events per metre of uplift increased exponentially across the three rapid uplifts to c. 20,000 events m−1 in 19844. The exponential trend is consistent with the progressive stretching of an elastic-brittle crust3,4,36 and an evolution from elastic (movement without faulting) to quasi-elastic (elastic movement with subordinate faulting36) behaviour. It continues under an increasing maintained differential stress (Fig. 4), which is the difference between the differential stress supplied by a pressure source and lost by faulting. The final approach to rupture occurs at the end of the quasi-elastic regime, after which deformation is accommodated inelastically by fracturing and faulting alone4,36. The rate of stress drop by faulting balances the rate of stress supply, so that the differential stress remains constant and the number of VT events increases in proportion with additional deformation (Fig. 4).

Fig. 4: Regimes of crustal deformation.
figure 4

a Ideal evolution from elastic (white), through quasi-elastic (yellow, A-B) to inelastic (magenta, B-C) deformation before rupture at C. The stress drops after rupture (see Fig. 5b). b The corresponding seismicity-deformation trend follows the difference between actual (solid curve) and linear elastic (dashed line) stress-deformation trends in (a). It increases exponentially (quasi-elastic regime) and then linearly (inelastic regime) with strain36. Since strain is proportional to deformation, similar shapes of trend appear on stress-strain diagrams. Notice that the trend A-B resembles non-linear elastic behaviour; such a description gives a misleading impression that only elastic deformation is taking place, instead of elastic deformation in unbroken rock combined with subordinate rock fracture and slip, which we describe as quasi-elastic36. c Changes in the number of VT events with corrected uplift during Aug 2010–Jan 2015 (squares), Jan 2015–Jan 2020 (triangles), and May 2020–June 2021 (diamonds). The exponential increase began in 2015–2016 (Fig. 2) and became linear in 2020. Best fit trends yield ΣN exp (hcor/49) and dN/dhcor = 120 events per cm (both with r2 = 0.99). d Similar transitions are seen for completeness magnitudes of 0.2 and 0.5; the corresponding pairs of best-fit trends give ΣN exp (hcor/36), dN/dhcor = 29 events per cm and ΣN exp (hcor/37), dN/dhcor = 14 events per cm (all with r2 > 0.98). ΣN denotes the total number of events and hcor is in cm. Supplementary Figs. 1 and 2 analyse the variances of the exponential and linear trends (Supplementary Information). Data from INGV-Sezione di Napoli (Osservatorio Vesuviano)71,72.

Inelastic behaviour immediately precedes bulk failure across the rock being deformed. In a previous analysis of the first three uplifts4, we argued that the crust had reached the end of the quasi-elastic regime in 1984, before the onset of slow subsidence. We predicted that inelastic behaviour would occur when the following phase of slow uplift had restored strain in the crust to its value before subsidence. Subsequent events have confirmed our prediction and show that Campi Flegrei is today deforming in the inelastic regime, following behaviour not previously witnessed during its ongoing unrest.

Results

Evolving regimes of deformation

Campi Flegrei’s unrest represents a perturbation from its equilibrium dynamic state of background subsidence. The amount of perturbation is given by the corrected uplift, hcor, for which background subsidence has been removed from the observed ground movement4,15,26 (see Methods). Following our previous studies4, therefore, we have used hcor to investigate the source of unrest.

Seismicity nearly ceased during the subsidence after 1984 and only re-emerged after the onset of slow ground uplift in 2005. VT events were sporadic until 2015–2016, after which they became more frequent before settling in 2020 to a mean rate of about 175 events per month (with a range from c. 100 to 250; Fig. 2). By April 2022, the ground had returned to its position at the end of the 1982–1984 crisis.

Using the corrected uplift does not alter the general shape of the time series, but does affect the amount and timing of changes in movement (Fig. 2). Thus, hcor was about 77 cm in 1950–52, 164 cm in 1969–72, and 184 cm in 1982–84; it subsequently decreased by 63 cm until 2000, before starting a slow increase. It reached its 1984 position (420 cm) in 2015–2016 (Fig. 2), at which time the number of VT events resumed an exponential increase with corrected uplift (Fig. 4). The trend continued until the first half of 2020 (when hcor was 460 cm). It then changed to a linear relation between the number of VT events and hcor, coinciding with our predicted transition to inelastic behaviour4.

A change from exponential to linear trends is found using all the recorded VT events and subsets with event magnitudes greater than 0.5 and 0.2 (Fig. 4), which correspond to the completeness magnitudes (Mc) for the monitoring networks operating, respectively, before 1985 and since 200038,39. Visual inspection suggests that a linear trend started at a corrected uplift of c. 460 cm for all VT events and for events with Mc ≥ 0.2, although, for Mc ≥ 0.5, such a change cannot confidently be distinguished from an extrapolation of the preceding exponential behaviour. However, a quantitative variance analysis reveals that exponential-linear trends do provide superior fits for all three data sets, with linear behaviour emerging at corrected uplifts of 455–460 cm (Supplementary Figs. 1, 2). The evolving pattern is therefore robust.

The re-emergence of quasi-elastic deformation when hcor had returned to its 1984 level reveals a close connection between corrected uplift and the conditions for faulting: thus, hcor measures the strain imposed on the upper crust by underlying pressure sources and the strain, in turn, determines the onset of faulting and VT seismicity. The connection is consistent with geodetic inversions, which attribute the current slow uplift to increases in pressure at depths of 3–4 km, beneath the hydrothermal system and at the levels inferred (1) for magma intrusions during the uplifts until 198425,40,41,42,43 and (2) for zones near the thermo-metamorphic horizon that are potentially supersaturated with magmatic gas18,19. The connection supports our hypothesis that corrected uplift is an appropriate field measure of strain in the crust. It also indicates that background subsidence is determined by an independent process operating at a location different from that controlling hcor, such as tectonic extension44, contraction of the melt-rich source45 and the compaction of volcanic deposits within the hydrothermal system45; if this were not the case, ground uplift would be the measure more closely connected to the other signals of unrest.

The quasi-elastic trends are described by the relation Total Number of VT events, ΣN exp (hcor/λch), where λch is the extension required for tensile failure36 (see Methods). For deformation between 2015 and 2020, λch is 36 and 37 cm for VT events with magnitudes greater than 0.2 and 0.5, but 49 cm when using all VT events (Fig. 4). The similarity of the first two values suggests that the VT catalogues were complete and that the magnitudes scaled inversely with their frequency; the larger value for all VT events may reflect the incomplete recording of earthquakes with magnitudes below 0.2.

The ratio hcor/λch also gives the ratio of differential stress to tensile strength, Sd/σT, in the deforming crust4,36 (see Methods). For deformation in tension, constraints from fracture mechanics show that Sd/σT (and, hence hcor/λch) cannot exceed 4 before the crust enters the inelastic regime4 (Fig. 5). This maximum value had been attained by the end of the 1982–1984 crisis4 (Fig. 6). We had thus previously argued that the crust in 1984 had reached the final transition before bulk failure4. In light of the caldera’s subsequent behaviour, we now suggest that bulk failure did in fact begin in 1984 and that the current uplift is continuing the process.

Fig. 5: Criteria for bulk failure.
figure 5

a The Mohr-Coloumb-Griffith failure curve (broken curve) shows the combinations of effective normal and shear stress that permit bulk failure73,74,75. The effective normal stress is the difference between normal stress and pore pressure; normal stresses are positive in compression and negative in tension. The axes have been normalised by dividing the shear and effective normal stresses by the tensile strength. The Mohr circles (continuous semi-circles) describe states of stress between the maximum and minimum effective normal stresses (corresponding to where the circles meet the normal-stress axis on the right and left). Failure occurs where the failure envelope meets a Mohr circle. It is tensile when the curves meet along the normal-stress axis. The minimum effective normal stress then equals the crust’s tensile strength, so that the [Effective Normal Stress]/[Tensile Strength] is -1. The diameter of a Mohr circle gives the ratio Sd/σT. The largest circle that can meet the failure envelope at tensile failure73,74,75 has a value of Sd/σT = 4. Conditions for failure when Sd/σT = 2 are shown for comparison. b Internal rupture (white ellipse) and relaxation (A–B–C) explains stress drops without measurable changes in surface and corrected uplift or in conditions in the pressure source (purple). The trend for simple elastic relaxation (A–B–D) would produce a measurable decrease in surface and corrected uplift and require a reduction in deforming pressure. Campi Flegrei’s behaviour is consistent with internal rupture (A–B–C).

Fig. 6: Reconstructed stress history, 1950-2022.
figure 6

a Evolution from elastic to inelastic deformation. Between 1950 and 1984, deformation of Campi Flegrei’s crust evolved from elastic (no or low seismicity until 1969–197023) to quasi-elastic (1970–19844), during which the number of VT events (ΣN) increased exponentially with corrected uplift (hcor: ΣN exp (hcor/λch) for λch = 100 cm). hcor is taken as 0 for zero differential stress in 1950 (hcor50). By late 1984, the ratio hcor/λch (= Sd/σT1) had reached 4.0−4.2, which is the maximum value expected before the transition to inelastic deformation. The transition may have been crossed, after which deformation changed to a slow subsidence and then uplift. Another exponential trend began (see also Fig. 4) after hcor had returned to its 1984 value (dashed line); this changed in 2020 (solid line) to a linear increase characteristic of inelastic deformation (Fig. 4). b The reconstructed stress (Sd) - uplift (hcor) path for 1950–2022. The dates for each stage are A–B, 1950–1984; B–C, 1984; C–D, 1985–2000; D–C (return), 2000–2015; C–E, 2015–2020; and E–F, 2020–2022. The diagonal dashed lines show ideal elastic trends. c The differential stress increased to 4σT1 as deformation evolved from elastic (white), through quasi-elastic (yellow), to inelastic (magenta). Rupturing began at the second transition, triggering a stress drop of 3σT1 (lettering as in b). d The stress drop in 1984 reset the value of hcor to 315 cm for zero differential stress (hcor84, extrapolating backwards the elastic relaxation along C–D). The corrected uplift increased again in 2000 (D), passing from elastic to quasi-elastic behaviour in 2015 (C) and to inelastic behaviour in 2020 (E). The 1984 rupture reduced the bulk tensile strength in the zone of failure to σT2.Inelastic behaviour since 2020 is consistent with Sd = 4σT2 = 1.5σT1.

During the episodes of quasi-elastic uplift until the end of 19844, λch was 100 cm for events with magnitudes greater than the contemporary completeness value of 0.5. The corresponding value of 37 cm for 2015–2020 (Fig. 4) thus shows a drop in λch by a factor of 2.5–3.0 (≈100/37). Using the subscripts “1” and “2” to distinguish the quasi-elastic sequences ending in 1984 and 2020, the ratio λch1/λch2 is (σT1/σT2)(E2/E1) for a constant geometry of deformation, where σT and E denote the crust’s tensile strength and Young’s modulus (see Methods). The modulus represents the elastic resistance to deformation across the deforming crust and decreases with the amount of fracturing, or damage46; hence we expect E1/E2≥1. Rock-physics experiments suggest that E decreases by less than 30% before rupture47, so that 1 < E1/E2 < 1.4. We have no evidence at Campi Flegrei for major changes in damage across the crust as a whole between 1984 and 2020 and so expect E1/E2 ≈ 1. We thus interpret λch1/λch2 as the ratio σT1/σT2 of the crust’s tensile strength before 1984 and in 2015. Hence, our calculated values for λch1 and λch2 suggest that the effective tensile strength is today about 2.5–3.0 times smaller than it was in 1984.

The decrease in strength reflects conditions within the zone of rupture, rather than the whole crust. Tensile resistance drops from its starting value to zero as fractures grow and link together before bulk failure is complete. The time taken for the decrease is typically much smaller than the timescale of observation, so that bulk failure normally appears to be instantaneous. At Campi Flegrei, however, the rupture that began in 1984 appears to have halted before complete failure occurred. Factors favouring incomplete failure are the relaxation of tensile stresses across the incipient rupture and fracture growth being halted against a resistant horizon. In either case, we suggest here that the limitation was overcome in 2020, after which the approach to bulk failure resumed.

We have reconstructed the sequence towards bulk failure using four constraints from the trends between VT seismicity and corrected uplift: (1) Sd = 0 at the start of unrest in 19504 so that, for 1950–1984, [Change in hcor]/λch1 = [Change in Sd]/σT1 = Sd/σT1 ≈ 4; (2) by 2015, Sd had regained the value it had at the end of 1984 so that, for 2015–2020, [Change in hcor]/λch2 = [Change in Sd]/σT2 = 1.35; (3) comparing the quasi-elastic trends, λch1/λch2 = σT1/σT2 = 2.7; and (4) in 2020, the VT-deformation trend changed from exponential to linear. Utilising these constraints, we have reasoned as follows:

  1. 1.

    Inelastic behaviour since 2020 (Constraint 4) indicates deformation under a constant differential stress (Sd2).

  2. 2.

    Working backwards in time (Fig. 6), Sd2 was reached following a quasi-elastic increase in stress by 1.35σT2 in 2015–2020 (Constraint 2). This, in turn, followed the 1984–2015 subsidence and uplift of hcor by 63 cm; because the movement was almost aseismic, we associate it with an elastic relaxation and subsequent re-straining of the crust. Before 1984, the elastic modulus allowed the differential stress to increase at 1σT1 per metre, which is consistent with E on the order of 10 GPa (see Methods). Assuming that the modulus has remained unchanged, the variation of 63 cm corresponds to a change in stress equivalent to 0.63σT1.

  3. 3.

    The crust has twice reached the transition between quasi-elastic and inelastic behaviour—in 1984 and 2020 (Constraints 1 and 4). The 1984 transition occurred when Sd1 ≈ 4σT1. By analogy, we assume that Sd2 ≈ 4σT2 in 2020. The tensile strength in the zone of rupture was 2.7 times smaller than in 1984 (Constraint 3). Hence Sd2 = 4σT1/2.7 = 1.5σT1, so that the stress at the start of elastic relaxation in late 1984 (and after its recovery in 2015) was 1.5σT1 − 1.35σT2 = 1.5σT1 − 1.35(σT1/2.7) = σT1.

  4. 4.

    We conclude that the applied stress dropped by 3σT1 (from 4σT1 to σT1) in 1984. Stress drops are accommodated by relaxing the crust. The drop in 1984 occurred without measurable changes in either surface or corrected uplift. This can readily be explained by relaxation of the crust around the opening of a new internal rupture (Fig. 5), which has since been reactivated and is now extending within the zone of VT seismicity at the eastern end of the Starza fault system (Fig. 6).

  5. 5.

    The partial relaxation in 1984 means that the uplift for zero differential stress is no longer its value in 1950. Starting after the 1984 stress drop (to Sd = σT1), backward extrapolation along the elastic trend gives an effective value of 315 cm for hcor when Sd is zero (Fig. 6). As a result, the proxy measurement for deformation is currently hcor – 315 (in cm).

In 1982–1984, most of the VT hypocentres along La Starza occurred within a volume of c. 25 km3 at depths to 3 km (Fig. 3). During the slow uplift since 2005, they have been restricted to c. 5 km3 towards the eastern end of the system, mostly at shallower depths of 2 km or less; an exception is a narrow, inclined zone of hypocentres that extends to 3 km beneath Solfatara and may mark a pathway for gas escape22. Overall, however, fracturing has become concentrated in a volume as much as five times smaller than its original value.

The concentration around a fault system and the general shallowing of VT seismicity (Fig. 3) suggests that rupture began in 1984 at depths of about 3 km and is now propagating towards the surface. It is thus consistent with the breaking of a sealed fault (or adjacent rock) as the crust is bent upwards around the periphery of a pressure source beneath Campi Flegrei’s hydrothermal system (Fig. 7).

Fig. 7: Mechanisms of slow subsidence and uplift, 1984-2022.
figure 7

a Campi Flegrei’s long-term melt-rich zone feeds shallow intrusions (magenta), while degassing (arrows) maintains crust rich in CO2 (green) below the low-permeability base (brown) of the hydrothermal system (blue). The surface cap for the hydrothermal system (olive)18 is overlain by the uppermost volcanic deposits (pink). Fractures beneath Solfatara and Pisciarelli (black dashed line leading to black triangle) allow the leakage of hydrothermal steam and CO2. b Changes with time in the CO2/H2O ratio of gases collected from Solfatara-Pisciarelli (1984-2022)22,72. Both. 1984–1986. Repeated intrusions of sills (1950–1984) leads to partial rupturing through the low-permeability zone in 1984 (double black lines) and an increase in the escape of trapped CO2 towards Solfatara-Pisciarelli (a, left; b, pale yellow). 1986–2000. CO2 pressure decays during withdrawal (a, middle, white), continuing slow subsidence and promoting a decrease in the CO2/H2O ratio. Subsidence closes fractures and enhances resealing by mineralisation. Young seals are intermittently broken to allow temporary surges in CO2 escape (b, pale green). 2000-2015. The pressure drop increases the gas pressure gradient above the melt-rich zone and, hence, also the flux of ascending CO2 (a, middle, arrows). CO2 accumulates as the low-permeability horizon re-heals. Slow uplift begins with almost no seismicity until 2015. The CO2 flux to the surface increases as the local pressure gradient becomes larger (b, blue). 2015-2022. The crust stretches beyond its value in 1984, triggering the return of persistent VT seismicity at shallow depth and continuing the approach to completing crustal rupture (a, double black lines) with a concentration of seismicity at shallow depth (a, b, orange).

Progressive rupture offers new insights into the changing style of deformation after 1984 (Fig. 1). The abrupt and permanent uplifts before 1984 are consistent with intrusions of magmatic sills4. The slower and reversible movements after 1984 are more readily explained by variations in the flow of fluids through the crust4,48,49, as recorded by changes in the compositions of gases from fumaroles within the Solfatara cone and outside its perimeter, at Pisciarelli49.

The CO2/H2O ratio, in particular, has been used to measure the relative contributions from magmatic (CO2) and hydrothermal (H2O) sources49. Following a general decrease between 1986 and 2000, its value has since shown small fluctuations about a sustained increase (Fig. 7). The changes mimic the mean and fluctuating patterns of ground movement, but their relations at different timescales are contradictory. Thus the short-term fluctuations in CO2/H2O occurred after those in ground uplift (within months during both mean subsidence and uplift)22, whereas the sustained increase in CO2/H2O in 2000 occurred 48-60 months before the start of slow ground uplift (Fig. 7). Apparently, therefore, the same mechanism cannot have controlled the connections between uplift and gas chemistry on short and long timescales.

The contradiction disappears when we use corrected uplift, and not ground uplift, as the measure of pressure changes causing unrest. As was found for the ground uplift, fluctuations in hcor occurred months ahead of those in CO2/H2O (Fig. 7). However, hcor also began to increase in 2000, before the sustained increase in CO2/H2O (Fig. 7). The chronological consistency at all timescales favours a common mechanism for coupling the changes in deformation and gas composition since 1984; the connection with corrected uplift suggests a control from pressure changes at depths below the hydrothermal system; and the short time differences of months between changes in both signals suggests that the pressure changes are related to the same mechanism for accumulating CO2.

A popular explanation for the connection between uplift and gas composition is that injections of magmatic gas mix with and heat hydrothermal fluids and, by raising pore-fluid pressure, drive caldera-wide uplift and the escape of fluids temporarily enriched in CO250,51. Implicit assumptions are that pressure changes are concentrated within the hydrothermal system and that uplift is a consequence of gas injection. Our analysis suggests that the source of pressure changes instead lies below the hydrothermal system, in which case the fluid-mixing hypothesis does not apply. We propose instead that rupturing of the crust began in 1984 in response to the strain accumulated in the crust by the emplacement of sills during the three episodes of rapid uplift between 1950 and-19844. The initial rupture broke open a new connection between the zone of CO2 storage beneath the hydrothermal system and Solfatara-Pisciarelli at the surface (Figs. 1, 7); in contrast to previous interpretations, therefore, we view the changes in gas flux to be a consequence of uplift. The source pressure decreased as CO2 escaped, triggering subsidence in addition to the background value from the end of 1984. The rate of subsidence gradually slowed as the source pressure decreased (Fig. 7). The CO2/H2O ratio at first increased in 1984-1986 with the influx of CO2 through the new connection; it then followed a long-term decrease as the rate of gas supply declined (Fig. 7).

From Darcy’s Law, the rate of permeable flow is proportional to the pressure gradient driving fluid motion. While the withdrawal of CO2 from the storage zone reduced the gradient of gas pressure to the surface, it simultaneously increased the gradient driving the ascent of CO2 from depth, where fluid pressure remained unchanged. By 2000, the rate of arrival of new gas exceeded the rate of withdrawal through Solfatara-Pisciarelli. The pressure beneath the hydrothermal system increased, initiating the slow uplift and increase in CO2/H2O (Fig. 7).

Our reconstruction suggests that the geophysical and geochemical signals at Campi Flegrei since 1950 belong to a cascading chain of unrest, as intrusions of magma until 1984 stretched the crust sufficiently to perturb the subsurface flow of gas. The depth and location of changing gas pressure overlap with those inferred for shallow magma intrusions. This is consistent with the geometric similarity of caldera-wide uplift before 1984 and since 2000, independent of the pressurising mechanism12,25. The overlap may not be a coincidence. Low-permeability horizons, such as magmatic sills and the mineralised claystones and siltstones19, are natural locations for trapping underlying fluids and allowing pore pressures to approach lithostatic values20,52,53 (Fig. 1). Near-lithostatic pore pressures in turn favour the arrest of ascending magma and its emplacement as a sill54. It is thus easy to imagine a self-reinforcing process in which low-permeability horizons favour high pore pressures, high pore pressures favour sill emplacement and sill emplacement favours high pore pressures still further.

The close temporal associations between corrected uplift and (a) the values of CO2/H2O and (b) the emergence of persistent seismicity support the use of hcor, and not surface uplift, as a proxy for the strain caused by unrest. The difference between the two measures of movement is too small to significantly affect interpretations of geodetic changes during short, rapid uplifts15. Since 1950, however, the total increase in corrected uplift is about one-third greater than that for ground uplift (Fig. 2), indicating a closer approach to crustal rupture than is evident from surface movements alone.

A pressure control beneath the hydrothermal system removes the need to invoke fluid mixing within the hydrothermal system itself. We cannot exclude the possibility that some gas may also filter into the rest of the hydrothermal system, but suggest that the amount is too small to cause significant caldera-wide deformation; hence the primary source of pressure change remains below the hydrothermal system. The observation that hcor is higher in 2022 than its position in 1984 suggests that either the CO2 pore pressure is greater than before or that the zone of gas storage has become larger.

In 2020, the crust entered the inelastic regime of deformation. Inelastic behaviour is a precursor to further rupture36. Bulk failure, though, is not guaranteed and at least three outcomes can be envisaged:

  1. (1)

    The crust establishes a new steady fluid flow. The corrected ground movement slows to a halt and Campi Flegrei returns to background subsidence.

  2. (2)

    The crust establishes an oscillating fluid flow, controlled by competition between cracks being sealed by mineralisation and opened by accumulating pressure. The corrected ground movement switches between slow uplift and subsidence, as has been seen, for example, at the Yellowstone caldera in the USA for at least a century53.

  3. (3)

    Uplift continues until the upper crust breaks completely. The fracture system near Sofatara-Pisciarelli is a preferred location for rupture. Assuming that the pathway will tap only pressurised fluids, the consequences may range from a general increase in gas flux at the surface to concentrated releases of gas pressure as phreatic explosions55.

In each scenario, the shallow crust at Campi Flegrei will have been partially or completely ruptured. It will thus be structurally weaker than it has been at any time during its unrest since 1950 and so, compared with earlier episodes of unrest, it will be easier for new magma to find a route to the surface. As a result, a magmatic eruption may be preceded by weaker precursory signals than seen to date—for example, smaller rates of VT seismicity or amounts of ground uplift. Such behaviour would resemble that observed at Rabaul caldera, in Papua New Guinea, when its 1994 eruption was preceded by VT event rates ten times smaller than had been recorded during a two-year VT seismic crisis the decade before3,56,57.

Rabaul’s example shows that analogy with past behaviour cannot be relied upon to forecast eruptions from large calderas. Our view that major structural change is occurring in Campi Flegrei’s crust leads to the same conclusion. The sequence from elastic deformation to rupture has also been identified before eruptions at stratovolcanoes reawakening after long repose36. In such cases, it is not unusual for the sequence to have been completed within the months of recognised precursory unrest. The sequence thus parallels increasing levels of alert. Similar parallels are not evident at large calderas, because much longer time intervals are needed to complete the whole sequence. We speculate that the experience of short sequences at stratovolcanoes may be encouraging a subconscious delay in recognising when large calderas have reached a similar state before rupture. Pressing goals, therefore, are (1) for large calderas in general, to enhance current pre-eruptive scenarios58,59 with new criteria that explicitly account for structural changes in the upper crust, and (2) for Campi Flegrei in particular, to reinforce detailed geodetic and gas monitoring to investigate the potential for rupture, especially along the whole Starza fault system.

Methods

Relative location of earthquakes

We obtained relative double-difference locations for earthquakes in 1982–1984 and 2005–2022 using the HypoDD algorithm60 to analyse P-wave and S-wave arrival times. To account for changes in the number and type of seismic monitoring stations at Campi Flegrei, we examined the two time intervals separately, but maintained the same inversion scheme. We used the 12-layer, one-dimensional seismic velocity profile derived from P-wave velocities61 with a value of 1.8 for the ratio vP/vS of P-wave to S-wave velocities (Table 1).

Table 1 One-dimensional seismic velocity model used for relative earthquake location.

From pre-processing, we obtained 8142 and 5037 linked events for the 1982–84 and 2005–2022, respectively. Our loose inversion scheme required three sets of iterations with the Least SQuaRe (LSQR) algorithm, which is suitable for large datasets, for at least 8 observations per pair of events, a maximum distance of 2 km between linked events, and a damping parameter of 200. Processing with HypoDD yielded (1) for 1982–1984, relative locations for 3062 events, median residual time of 0.09 s (range 0.05–0.4 s), and median estimated errors for the three spatial components of 77 m, 79 m, 100 m, and (2) for 2005–2022, relative locations for 2522 events, median residual time of 0.04 s (range 0.02–0.1 s), and median estimated errors of 43 m, 44 m, 42 m.

VT earthquakes and seismic moment

The number of VT events is a proxy measure for the cube root of seismic moment (Fig. 8), which, in turn, is proportional to the amount of co-seismic dislocation and faulting (i.e., the inelastic component of deformation). The seismic moment M0 = μWLD, where W and L are the width and length of a fault, and D and μ are the dislocation and shear modulus. For a constant stress drop, the aspect ratio (W/L) of a fault and the strain (D/L) induced by dislocation can be assumed constant62, so that M0 is proportional to D3. Hence, the co-seismic dislocation DM01/3. Figure 8 shows that the number of VT events increases in proportion with M01/3, implying that it is also proportional to co-seismic dislocation.

Fig. 8: VT seismicity as a proxy measure for co-seismic movement.
figure 8

The amount of co-seismic dislocation associated with a VT earthquake (co-seismic movement) is proportional to the cube root of its seismic moment (M01/3; see Methods). The total number of VT earthquakes is in turn proportional to the cumulative value of M01/3. The examples show the relations for all magnitudes and for magnitudes greater than 0.2; the gradients of the corresponding best-fit lines are 2.13 × 10−6 and 1.05 × 10−6 events (dyn cm)-1/3 (both with r2 > 0.99).

M0 in dyne cm is related to the moment magnitude MW by

$${{{{{\rm{Log}}}}}}\,{M}_{0}=(3/2){M}_{W}+16.1$$
(1)

At Campi Flegrei, Mw can be calculated from the duration magnitude Md, which is provided routinely by the INGV-Vesuvius Observatory (with an a priori uncertainty of ±0.3) using63

$${M}_{W}=0.61({{\pm }}0.04)+0.82({{\pm }}0.03){M}_{d}$$
(2)

Background subsidence at Pozzuoli, Campi Flegrei

An approximately constant rate of background subsidence at Pozzuoli for at least 2400 years has been identified from the changing position of the floor of the Roman market place, Serapeo (or Serapis), which was built close to the port in 200 BC26. The position of the floor has been recorded by reports of restoration since Roman times, contemporary paintings and, since the early 1800s, by direct geodetic surveys26. These measurements have been supplemented by archaeological and palaeontological data from the historic centre of Rione Terra, less than 500 m to the south28.

As shown in Fig. 9, the intervals of most data (before 600 and after 1430) can be linked by linear trends with similar mean rates of subsidence (c. 150–200 cm per century) and mean rates of uplift an order-of-magnitude faster (c. 1200–1500 cm per century)26. Although the trends between 750 and 1430 are not well constrained, they have the value of simplicity by assuming that their behaviour lies within the ranges provided by the better-constrained trends. The observation that vertical movements have been dominated by subsidence at similar mean rates suggests that this behaviour is a typical background state controlled by a common mechanism. When background subsidence is removed, the corrected movement shows how Campi Flegrei has been occasionally interrupted in historic times by net uplifts of as much as 16 m at Pozzuoli (Fig. 9); as described in the main text, the latest of these began in 1969 and continues today (2022).

Fig. 9: Vertical ground movement at the Serapeo in Pozzuoli since Roman times26.
figure 9

a Ground movement shows a saw-tooth pattern with time (grey dashes), consistent with (b) a background rate of subsidence of c.150–200 cm per century, interrupted by (c) intervals of uplift at rates about 10 times faster. The step-like pattern of uplift illustrates how Pozzuoli has been raised in the past two millennia by more than 35 m above the level it would have had if only subsidence had occurred.

Stress, strain and VT event rates

For quasi-elastic behaviour, the parent relation between the total number of VT events (ΣN) and applied differential stress (Sd) is3,4

$$\Sigma N={(\Sigma N)}_{0}\,{{{{{\rm{exp }}}}}}[({S}_{d}{-}{S}_{0})/{{\sigma }}_{T}]$$
(3)

where (ΣN)0 and S0 are the number of VT events and the differential stress at the start of the quasi-elastic regime when failure begins at the transition from quasi-elastic to inelastic beghaviour (position “a” in Fig. 4a) and σT is tensile strength. The differential stress for bulk failure depends on the geometry of the deforming system and can take values of up to 4σT (Fig. 5).

For tension and extension in the quasi-elastic regime, SdEε f(ν), the product of Young’s modulus of Elasticity E, bulk strain ε and a function of Poisson’s ratio ν (whose specific form depends on the shape of deformation; e.g., for a thin sheet with similar width and length Sd is (4/5) in pure tension (uniaxial strain)64, but (1 + ν)/(1 − ν2) = (4/3) for low-angle bending65 and an ideal value 0.25 for Poisson’s ratio64). From standard scaling64,66, bulk strain is Chcor/R, where R is the starting length of the feature to be deformed (such as the radius of a pressurising source) and C describes the geometry of deformation. Using the tilde symbol “~” to denote “order of magnitude”, C ~ 1 for the small deformations of ~10-3 or less involved at Campi Flegrei (hcor on the order of metres or less and R ~ 103 m).

Sd can thus be approximated to CE(hcor/R) f(ν) and, when substituted in Eq. (3) yields

$$\Sigma N={(\Sigma N)}_{0}\,{{{{{\rm{exp }}}}}}[({h}_{{{{{{{\rm{cor}}}}}}}}{-}{h}_{0})/{{\lambda }}_{{{{{{{\rm{ch}}}}}}}}]$$
(4)

where h0 is the corrected uplift at the start of quasi-elastic behaviour and λch = (σT/E)(R/C)(1/f(ν)), which describes the increase in the length of the crust along the direction of tension at the start of inelastic deformation (when the tensile stress equals the crust’s tensile strength). For a constant geometry of deformation (R/C)(1/f(ν)) remains unchanged. Hence, comparing values of λch at stages “1” and “2” the value of λch1/λch2 = (σT1/σT2)(E2/E1) and so depends only on the ratios of tensile strength and of Young’s modulus. Equation (4) describes the quasi-elastic trends shown in Figs. 5, 6.

For common rocks breaking in tension67,68, σT/E ~ 10-4–10-3. At Campi Flegrei, λch = 1 m in the quasi-elastic regime (Fig. 6) and so, for R = 2–2.75 km13,14, λch/R indicates a failure strain ~10-4, in agreement with experiment. The direction of tension runs approximately parallel with the surface. The amount of corrected uplift required to develop sufficient tension reflects the geometry of deformation. At Campi Flegrei, deformation is inferred to have started in 1950 (Sd = 0) and failure to have started in 1984 after a corrected upift of about 4 m (Fig. 6), yielding hcorr/λch = ΔSd/σT = Sd/σT = 4, consistent with values for bulk failure in tension (Fig. 5).

A second consistency check is available by setting Sd ~  ~Ehcor/R (see above), or E ~ SdR/hcor. By the end of 1984, Sd ≈ 4σT (Fig. 6). Hence, for R/hcor ~ 103 and σT ~ 106–107 Pa67,68, E ~ 10 GPa. The crust below Campi Flegrei consists mainly of marine sediments at depths of c. 2–3 km and of carbonate bedrock and igneous intrusions at deeper levels10,16 (Fig. 1). Rock physics experiments yield values for the elastic moduli of 5–7 GPa for sandstones69,70 and 20–60 GPa for limestones, trachytes, dolerites and granites67,69,70,70. A magnitude of 10 GPa is thus consistent with the moduli expected for rocks surrounding a pressure source at a depth of c. 3 km.