In this section we discuss and analyze closer the strong relationships between deformation and seismicity as emerges in the period 2005-11/2023, and particularly in the period 2019-11/2023, also including a comparison with an analogous relation obtained for the 1982-84 crisis, and discuss the possible evolution of the current bradyseismic crisis under the hypothesis that the caldera will continue to inflate and behave according to the same trends.

Figures 4 and 5 clearly show the coupled behavior of deformation and seismic activity at the different time scales. This coupling is well evident both in the long term (i.e., at the decadal scale) as well as in the short/medium terms (i.e., at the week/month/year scales). At the decadal scale, the accelerating deformation of the system, since 2010 well approximated by a parabolic function, is paired with a super-exponential acceleration of the seismic count (either total or above a given magnitude). The reasons for the quasi-parabolic increase of the uplift are beyond the objectives of this manuscript and the analyses here developed. However, the understanding of the cause of this accelerating process is an interesting future goal and questions on the possibility that the driver of the uplift is itself an accelerating cumulative process, such as, for instance the injection of gas or magma in the upper crust from a deeper source. At the shorter time scales, any temporary increase of the deformation rate is closely mirrored, with an apparent time lag up to 30–45 days, in an analogous increase of the seismic rate, with an amplification effect with time clearly visible since 2021 and even more since the beginning of 2023. Such a time gap between the two variables should be also further investigated possibly extending the analysis to shorter time scales (up to days) and to the consideration of the spatial distribution of faults/fractures and the structural features of the caldera.

Such a strong dependence between deformation and seismic activity is well described by the relationship between vertical uplift of the caldera (here represented through the RITE station) and the cumulative number of earthquakes, total or above a given magnitude. As shown in Fig. 8, such relationship can be well represented by two exponential functions, with exponents increasing in time. Such a dependency resembles the well-known behavior of quasi-elastic inhomogeneous materials and rocks under an increasing differential stress, either compressive or tensile, at both low- and high-pressure conditions35,38,51–54. Patterns of fracturing of brittle rocks in laboratory clearly show that at low stresses fracturing activity is null or very low and the behavior is close to elastic, i.e., there is a quasi-linear relation between stress and strain. At stresses of about half of the values at bulk failure the fracturing activity begins to increase in an exponential way until the strain has typically reached about 90–95% of the failure value. Under such conditions the stress increases in a non-linear way with respect to strain, producing a curve strain-stress that is downward convex. This condition of the material is also characterized by dilatancy, i.e., an increase of the volume due to a quasi-homogeneous formation of small cracks within the rocks. At larger increasing stresses, the fracturing activity typically accelerates very quickly producing the coalescence of fractures in a major fault and eventually leading to bulk failure of the brittle rock. For specific properties of the rock a variety of strain-hardening or strain-softening type behaviors of the material can also be observed before failure. Eventually, when the amount of stress has reached its maximum value, most of the supplied energy is spent in fracture and fault movements and the behavior of the rock becomes inelastic, i.e., the deformation continues under a constant stress. This regime is typically characterized by a proportional increase of seismic events with deformation.

It is important to observe that the nature of the exponential relationship between strain and cumulative seismic count depends on the mechanical properties of the crustal rocks52,54. These in turn depend on lithology as well as on the local conditions at depth such as pressure, temperature, water content as well as on other rock properties such as porosity, compaction and alteration of components and minerals. The underground distribution of these properties and conditions are largely unknown for the crust of CFc although several mechanical properties data exist based on measurements on core samples and outcrops55,56. Data overall show an elastic-brittle behavior of the CFc rocks although remarkable variations of key properties, such as elastic modulus and Poisson ratio, were measured as a function of porosity, temperature and water saturation conditions. It is therefore important to underline that the empirical relationship between deformation and seismic activity as illustrated in Fig. 8, and here interpreted as a progressive evolution of the quasi-elastic response of the shallow crust to an increase of stress, should be simply considered as an empirical relation valid for the superficial crust as a whole, i.e. assuming the crustal rocks as a single inhomogeneous medium with a mean elastic-brittle behavior. It should be noted that this relationship also potentially depends on additional processes which are not considered herein and that likely contribute to the deformation observed such as pore-pressure and thermal effects associated with the presence of fluids of magmatic, hydrothermal and exogenous origin in the rock matrix31,32,57.

The interpreted quasi-elastic behavior of the superficial crust appears, however, largely consistent with previous analyses and interpretations of the ongoing unrest26,34–36,40. According to them, the presence of low permeability levels at depths ranging between 1.5 and 3 km could favor the pressurization of the underlying layers due to the accumulation of magmatic gasses, fluids or magma rising from deeper regions. Detailed analyses of the distribution and mechanisms of the ongoing local seismicity also appear to be consistent with a stress concentration caused by the uplift of the central part of the caldera as well as with a possible effect of fluid-driven pore-pressure increase at the reactivated faults involved18,40,41. Similarly, the existence of the observed oscillations with variable periods in the geophysical trends investigated could be interpreted as a result of the unsteady fracturing of the crust under the increasing stress due to the pressurization generated by the arrival of deep fluids or magma batches and their complex interplay with the superficial hydrothermal system20,26,31,32,54.

Nevertheless, a significant difference between the analysis here presented and that of [36] is the fact that, based on the relation between deformation and number of earthquakes, the latter study claims the reaching of the inelastic regime of the CFc crust since May 2020 (and up to June 2021, corresponding to the end of the period analyzed), due to the establishment of a linear relationship between the two geophysical variables. We hypothesize here that such inconsistency is due to the different time period (about 2.5 years shorter than that analyzed here) over which the two analyses have been made. As a consequence, the present study, although generally consistent with the main conclusion of [35, 36] which suggests a progressive weakening of the mechanical properties of the shallow crust of CFc, gives evidence of a stronger sensitivity of the seismic activity on the uplift of the caldera, as will be better described and quantified later on in this section. It is worth underlining that this results is independent from the existing interpretations in terms of rheological behavior of the shallow crust.

To gain further insights on the recent dynamics of the caldera, the relationship between deformation and seismic count for the unrest crisis of 1982-84 was also analyzed. This is illustrated in Fig. 9. Figure 9a shows the cumulative total number of earthquakes versus the ground uplift as measured at Benchmark 25A of the leveling network (close to Pozzuoli, see Fig. 1-SI), whereas Fig. 9b shows the same relationship by considering only earthquakes with Md > 0.5 and the vertical uplift as measured by the tide-gauge at the Pozzuoli harbor (close to the center of the CFc, see again Fig. 1-SI). Despite the lower quality and amount of data recorded during this past crisis with respect to the current one, in both plots the correlation between the two geophysical variables results quasi-linear, R-squared > 99%, and not exponential. The same relation was obtained by considering different thresholds of magnitude. This is also different from the results of [35, 36] who estimated an exponential trend over the period 1982-84 by apparently using a lower number of data points.

The interpretation of such different deformation-number of earthquakes relationship with respect to the current crisis (see Fig. 8) is not trivial or obvious. However, as mentioned above, a quasi-linear relation between deformation and seismic count is typically indicative of an inelastic mechanical behavior of rocks in which most of deformation is accommodated by the propagation and slip of existing fractures and faults. This process appears consistent with most interpretations of the 1982-84 crisis as caused by an intrusion of magma with the formation of a thin sill at about 3–4 km depth14,22–25. We also note that the 1982-84 crisis occurred on a time period much shorter and therefore with deformation and seismic rates about one order of magnitude larger than those characterizing the present crisis. Such a different relationship between deformation and seismic activity appears as further evidence of the different forcing source and/or mechanism associated to the 1982-84 crisis with respect to the present one.

Finally, the ongoing trends of deformation and seismic activity, and particularly their so far robust, empirical exponential relationship as illustrated in Fig. 8, can also be used to gain some insights on the potential evolution of the bradyseismic crisis in case of a further uplift of the caldera with the same trends. From Fig. 8 it is evident how such exponential relationship takes into account and explains the progressively intensifying seismic activity recorded in the last few years due to the continuous, although oscillating, increase of the caldera uplift. Similarly, under the assumption that the same trend of uplift will continue in the future, we could expect a corresponding exponential increase of the number of earthquakes. According to Fig. 8b, in particular, it emerges that, for instance, a further uplift increase of just 20 cm (i.e., from 1.2 to 1.4 m at the RITE GNSS station) could generate, with the current exponential relationship, a number of Md > 0.5 events larger than that recorded since 2000 (i.e., about 2,000 more events vs the about 1,800 recorded since year 2000, see dashed portion of the exponential curve of Fig. 8b). A similar consideration can approximately be done also for the total number of earthquakes as shown in Fig. 8a.

Based on the analysis of the distribution of the number of earthquakes as a function of the magnitude Md, it is also possible to define the best fitting parameters of the exponential law (Gutenberg-Richter) describing this distribution. As shown in Fig. 8-SI, the fitting of the distribution is very good (R-squared > 99%) and the corresponding *b*-value over the investigated 2000–2023 period is equal to 0.87, by assuming a Mc equal to 0.5 (a very similar value was obtained with Mc=0.2 but given the discretization of the earthquake magnitude a value of Mc=0.5 was preferred). A slightly lower *b*-value of 0.83 is obtained for the most recent period 2021–2023, whereas, over the period 2007–2020, we obtained a value of 0.95, equal to that obtained by [45] with a different approach over the same period. We note also that our estimates of the *b*-value are ca. 5–10% lower than those obtained by previous studies over different previous periods19,45,58 evidencing a progressive decrease with time. This is also consistent with recent estimates of the *b*-value of the different volumes of the caldera showing a close relationship between this parameter and the physical and mechanical properties of the rocks and indicating lower *b*-values just at depths of 2.5-3 km where a caprock with brittle-fragile behavior is located59. All this represents further evidence of the overall progressive stress increase and medium weakening of the shallow crust of CFc as clearly observed in many experimental and volcanic systems19,54,60.

Based on the above-described law of magnitude distribution and the expected number of events of magnitude Md > 0.5 as a function of the ground deformation reached by the caldera, it is also possible to estimate the probability of occurrence of an earthquake above a given magnitude associated to a given future uplift of the caldera. Here we just mention that, based on the above-described data and hypotheses, just a further modest uplift of about 20 cm would be associated with a significant number of expected earthquakes with potential impact in the CFc area. Table 1-SI details the expected number of earthquakes above various magnitudes in a population of 2,000 events with Md > 0.5 (corresponding to a further uplift of about 20 cm), as a function of the *b-*value and its associated uncertainty. For instance, by using a *b-*value equal to 0.87, the magnitude-frequency law would expect 0.82 events with Md ≥ 4.5 on average, and 0.30 events with Md ≥ 5, simply by applying the above-described statistical basis. A *b*-value equal to 0.83 would produce an increase of ca. 50% expected large events, whereas a *b*-value equal to 0.95, a decrease of 50% (see Table 1-SI). These estimates, although based on an extrapolation of the reconstructed distribution of magnitudes to values larger than those so far observed, are qualitatively consistent with analyses of potential seismic activity at CFc based on different data and hypotheses40,61.

It should be noted, however, that due to the above-described oscillations in time of the rate of deformation, it is not possible to precisely forecast the expected seismic activity as a function of time. For instance, by hypothesizing future deformation rates in the range of values recorded in the last 6 months or 1 year, short-time forecasts of seismic activity, e.g. on a monthly basis, are affected by an uncertainty significantly larger than that associated with medium- and long-term forecasts which are able to filter out oscillations with shorter periods.

Nevertheless, the above-described super-exponential decadal trend of seismic activity, also given the time-increasing exponents describing the relationship between deformation and seismic count (see Fig. 8), remains suited to be investigated through a probabilistic application of the Failure Forecast Method (FFM), similarly to what done by [39] for the period 2000–2020, in order to explore the time scale of the evolution of the crisis in case of continuation of the observed trends.

In particular, in this new application of FFM the exponent alpha of the seismic rate has been varied between 1.2 and 239,62. The application of the method to the new datasets confirms the main outcomes of the previous analysis providing median values of the waiting times (i.e. the time to be waited to reach the critical time computed since 12/2023) from 0.5 to about 2.5 years with 95%ile values up to about 8 years (see Fig. 9-SI for more details). It is worth noting that these estimates are remarkably dependent on the period of time used in the regression of the inverse rates, in this case we tested 3 and 10 years. As a consequence, these estimates cannot be considered reliable forecasts of the occurrence of specific volcanic events given the large number of false alarms typically experienced in the application of FFM to real crises. Nevertheless, the analysis is of some interest since it indicates, over the parameter space explored, a relatively fast evolution of the crisis (up to few years as median values) in case of continuation of the observed trends.

The bradyseismic scenario above described represents the simple temporal extension of the trends currently observed. Nevertheless, it is uncertain and it represents just one of the plausible evolutions of the ongoing bradyseismic crisis. Based on observation of similar unrest crises at other calderas35,63, the interpreted quasi-elastic behavior of the shallow crust could be interrupted at any time, either by evolving into an inelastic regime of the crust or simply due to the occurrence of superficial magmatic intrusions or events, such as phreatic explosions, triggered by some of the complex, and still largely unknown, processes governing the magmatic and hydrothermal systems of CFc36,37,64. Similarly, the caldera uplift could, abruptly or more slowly, stop, the seismic activity decrease and the caldera enter into a new period of subsidence36. Unfortunately, based on our present understanding, the future evolution of the system is largely controlled by the internal dynamics of the volcanic system of CFc that, at this time, we are unable to observe and forecast accurately.