Pseudo-prospective earthquake forecasting experiment in Italy based on temporal variation of the b-value of the Gutenberg–Richter law

SUMMARY

We present the BVAL method, designed to forecast potentially damaging earthquakes (M_w ≥ 5.0) in Italy based on temporal variations of the b-value of the Gutenberg–Richter frequency–magnitude distribution. The b-value is used as an indicator of stress within the Earth's crust, with lower b-values associated with higher stress levels and an increased likelihood of significant seismic events. This method issues alarms when the b-value falls below a critical threshold. It is optimized using the HOmogenized instRUmental Seismic catalogue data from 1990 to 2004 and validated pseudo-prospectively using data from 2005 to 2022. Our analysis uses the recently developed b-positive (b⁺) method to compute the b-value from magnitude differences, providing resilience against data incompleteness. We compare the performance of the BVAL method with two established models: the Epidemic Type Aftershock Sequence (ETAS) model, which forecasts earthquake rates based on the epidemic principle that each shock triggers subsequent shocks, and the FORE model, which relies on the occurrence of strong foreshocks. Additionally, we evaluate two ensemble models that combine BVAL and FORE through additive (EADD) and multiplicative (EMUL) strategies to balance false alarms and missed events. The EADD model declares an alarm when either BVAL or FORE signals it, while the EMUL model triggers alarms only when both methods agree. We assess the predictive efficiency of these models using the area skill score, derived from Molchan diagrams, which plot the miss rate against the fraction of space-time occupied by alarms. Our results demonstrate that BVAL is less effective than FORE and ETAS at high space-time fractions, but it outperforms ETAS at low fractions (⁠|$\tau \ $| < 2–4 per cent), indicating its potential utility in scenarios where minimizing false alarms is critical. This comprehensive comparison highlights the strengths and limitations of each method, suggesting that integrating multiple forecasting strategies can enhance the reliability of earthquake preparedness and response efforts in Italy.

INTRODUCTION

The slope (b-value) of the frequency–magnitude distribution (FMD) of earthquakes (Gutenberg & Richter 1944) is indicated by some studies (Scholz 1968, 2015; Amitrano 2003) as a proxy of the state of stress within the Earth's crust and then as an indicator of the state of preparation of a future strong earthquake (Gulia & Wiemer 2010, 2019, 2021; Gulia et al. 2018, 2024). The ratio behind such studies is that low b-values are likely to be associated with high stress levels and then with a larger probability of occurrence of strong shocks, and high b-values with low stress levels and a lower probability.

In the last two decades, several works studied the space-time variations of b-value before moderate and strong earthquakes in different parts of the world. Among the others, we can mention Murru et al. (2004) for the 1997 central Italy seismic sequence, Nuannin et al. (2005) for the great 2004 Sumatra earthquake, Tsukakoshi & Shimazaki (2008) for the 2003 Northern Miyagi earthquake in Japan, Chan et al. (2012, 2023) for large earthquakes in Taiwan, Varotsos et al. (2012) before large earthquakes in Greece, El-Isa & Eaton (2014) for Central California and Andaman-Sumatran regions, Montuori et al. (2016) and García-Hernández et al. (2021) for the 2016 Central Italy sequence, Lu (2017) for New Zealand, Öztürk (2020) for Central Anatolia, DeSalvio and Rudolph (2021) for western Nord America and Pastoressa et al. (2023) for Central Italy. In such works, there is a general consensus that strong earthquakes are preceded by low b-values but even the opposite is found in a few cases (e. g. Parsons 2007; Sharma et al. 2023).

In this work, we develop a method to forecast potentially damaging earthquakes (M_w ≥ 5.0) in Italy based on the temporal behaviour of the b-value (called BVAL). Our algorithm issues an alarm every time the b-value goes below a given threshold coming from a higher value. The duration of the alarm is made to vary between a fraction of second to the total duration of the experiment.

We take advantage of the recently developed b-positive (or b⁺) method (van der Elst 2021) to compute the b-value from the distribution of magnitude differences which was proven to be more resilient to incompleteness than methods based on the distribution of simple magnitudes (Tinti & Gasperini 2024).

We optimize the forecasting algorithm based on the analysis of the HOmogenized instRUmental Seismic (HORUS) catalogue (Lolli et al. 2020) from 1990 to 2004 (learning period). Using the same catalogue, we test the algorithm pseudo-prospectively on a successive interval from 2005 to 2022 (testing period) and compare its performance with that of other two forecasting methods already applied in Italy: the Epidemic Type Aftershock Sequence (ETAS; Ogata 1988; Murru et al. 2009; Lombardi & Marzocchi 2010; Biondini et al. 2023; Mizrahi et al. 2023), the FORE (Gasperini et al. 2021; Biondini & Gasperini 2023), and two ensemble models obtained, respectively, as addition and multiplication (Bayona et al. 2022) of BVAL and FORE methods.

In the additive model (EADD), an alarm is triggered if either of the two individual methods detects the precursor signal. Consequently, the prediction window of this ensemble model is the union of the prediction windows from both methods. This design increases the sensitivity of the model, aiming to minimize missed detections of the precursor signal. However, it may result in a higher number of false alarms due to its broader criteria for issuing alerts.

In contrast, the multiplicative model (EMUL) activates an alarm only when both simple methods independently detect the precursor signal. Therefore, the prediction window for this ensemble model is the intersection of the prediction windows of the two methods. This approach makes the EMUL model more prudential, focusing on reducing false alarms by requiring a higher level of confirmation before triggering an alert.

The choice between these models depends on the user's priority: the EMUL model seeks to limit false alarms by demanding agreement from both methods, while the EADD model aims to enhance detection of the precursor signal, even at the risk of generating more false alarms.

See in Table 1, a summary of various forecasting models considered in this work.

ANALYSIS REGION AND EARTHQUAKE CATALOGUE

As application region, we consider the overlap of two tessellations (⁠|${R}_1$| and |${R}_2$|⁠) of the Italian territory spanning from 7° E to 19° E in longitude and from 36° N to 47° N in latitude. Both tessellations consist of square cells with side length of |$L\ = \ 30\sqrt 2 $| km. The tessellation |${R}_2$| is shifted by L/2 both in latitude and in longitude with respect to |${R}_1$| (Fig. 1). The purpose of this dual tessellation is to enhance the spatial alignment between the forecasting features detected by different algorithms with the actual earthquakes. This approach almost doubles the portion of space covered by the alarms, thus reducing the efficiency of the various algorithms, but avoids failing a correct prediction because the earthquake to be predicted, although close, is located slightly outside the alarmed cell as in the sketch of Fig. 2, in which the precursor seismicity (blue circles) is able to predict the target earthquake (red star) using the red grid but not the black grid.

The tessellation |${R}_1$| was developed starting from the initial point |${x}_1$| = (7° E, 47° N) and moving the centre from |${x}_1$| in steps of L km for each cell, until reaching |${x}_2$| = (19° E, 36° N). The centres of the cells in tessellation |${R}_2$| are shifted relative to |${R}_1$| by −L/2 km in longitude and +L/2 km in latitude. To ensure uniform cell dimensions, we converted the coordinates of points |${x}_1$| and |${x}_2$| from the WGS84 geographic reference system to kilometric coordinates in the RDN2008 Italy zone (E-N) CRS EPSG:7794 using the QGIS software (see the Data Availability section).

Following Gasperini et al. (2021), Biondini et al. (2023) and Biondini & Gasperini (2023) we only consider cells in which at least one earthquake of magnitude |${M}_\mathrm{ w} \ge 4.0$|  |$( {z \le 50\ \mathrm{ km}} )$| occurred on the Italian mainland during the period 1600–1959 according to the seismic catalogue CPTI15 (Rovida et al. 2022). Isolated cells not contiguous with the two main tessellations are excluded. The resulting grids are composed of 167 and 170 cells for |${R}_1$| and |${R}_2$|⁠, respectively (Fig. 1). This approach is consistent with the methodology defined by Gasperini et al. (2021) to prevent the overestimation of the forecasting performance by excluding predominantly aseismic areas, such as the island of Sardinia.

To calibrate and pseudo-prospectively test the relative efficiency of various forecasting models, we use the HORUS catalogue (Lolli et al. 2020; see the Data Availability section), which collect earthquakes from 1960 to the present and is homogeneous in terms of magnitudes. To balance statistical resolution and the inherent uncertainties in magnitude measurements, we preliminarily bin the magnitudes to a single decimal digit. Hence, all magnitude thresholds mentioned below must be understood as reduced by half of the bin size (0.05).

During the learning period ranging from 1990 to 2004, the catalogue reports 18 target shocks with |${M}_\mathrm{w} \ge 5.0$| and 8 with |${M}_\mathrm{w} \ge 5.5$| within the double tessellation. In the testing period ranging from 2005 to 2022, it reports 34 target shocks with |${M}_\mathrm{w} \ge 5.0$| and 12 with |${M}_\mathrm{w} \ge 5.5$|⁠.

According to Lolli et al. (2020), the HORUS catalogue can be considered complete for the Italian mainland for |${M}_\mathrm{w} \ge 4.0\ $| since 1960, for |${M}_\mathrm{w} \ge 3.0$| since 1981, for |${M}_\mathrm{w} \ge \ 2.5$| since 1990, for |${M}_\mathrm{w} \ge 2.1$| since 2003 and for |${M}_\mathrm{w} \ge 1.8$| since 2005.

However, for the purpose of determining the b-value, using the method of magnitude differences (van der Elst 2021; Tinti & Gasperini 2024), which is notably resilient to incompleteness, we re-estimated the completeness magnitude (⁠|${M}_\mathrm{ c}$|⁠) using the maximum curvature method (Wiemer & Wyss 2000; Woessner & Wiemer 2005), which is less accurate but preserves more data. The seismicity distribution used for these calculations is illustrated in Fig. 3, which shows the epicentres of earthquakes within the analysis region (R) and the CPTI15 polygon, highlighting the spatial coverage and density of seismic events in the study area. The re-estimated |${M}_\mathrm{ c}$| are |${M}_\mathrm{w} \ge 2.7\ $| since 1960, |${M}_\mathrm{w} \ge 2.0\ $| since 1981, |${M}_\mathrm{w} \ge \ 2.1$| since 1997, |${M}_\mathrm{w} \ge 1.9$| since 2003 and |${M}_\mathrm{w} \ge 1.2$| since 2005 (Fig. 4).

COMPUTATION OF b-VALUE

Up to a few years ago, the determination of the b-value was usually done through a formula, derived many decades ago by Utsu (1965) and Aki (1965), which assumes a continuous exponential distribution of magnitudes. If the magnitudes are binned to only one decimal digit, Utsu (1966) suggested an approximate correction to the formula consisting of subtracting one-half of the binning size (0.05) from the minimum magnitude of completeness of the data set. As the distribution of binned magnitudes is not continuous but discrete, Guttorp & Hopkins (1986) derived the exact formula, using the geometric distribution theory. Such formula, however, was almost totally neglected by the following literature (see Tinti & Gasperini 2024 for a thorough overview).

More recently, van der Elst (2021) showed that the b-value can be determined using the distribution of positive differences between magnitudes of successive shocks and, in this case, the estimation is more resilient to the incompleteness of the data set, particularly if the differences lower than a given trimming threshold |${\rm{\Delta }}{M}_\mathrm{ c}^{\prime}$| are discarded. Tinti & Gasperini (2024) confirmed that the effect of incompleteness is mitigated by using trimmed magnitude differences but not that positive differences are significantly better than negative ones. However, as the so-called b-positive (b⁺) estimator proposed by van der Elst (2021) is largely used in the most recent literature, we adopt it even in this work. For computing the b-value, we use the expression proposed by Tinti & Gasperini (2024):

where |$\delta $| is equal to one half of the magnitude bin size (typically 0.05) and |$\overline {{\rm{\Delta }}M} $| is the average of the positive magnitude differences. Such expression is algebraically equivalent to that proposed by van der Elst (2021) in its eq. (A3) but can be preferred because it does not use the hyperbolic arc-cotangent function, which is not provided in some software environments. Although the b⁺ estimator is more resilient than methods based on the distribution of magnitudes, it is not completely insensitive to incompleteness. Hence, a completeness analysis (Woessner & Wiemer 2005) is needed anyhow before using it.

van der Elst (2021), in analysing real sequences, considered two options: in the first one, he used all catalogued earthquakes and then applied the b⁺ estimator above a minimum positive magnitude difference |${\rm{\Delta }}{M}_\mathrm{ c}^{\prime}$| ranging from 0.3 to 0.4 units, while in the second option, he used only the earthquakes with magnitude above the completeness magnitude |${M}_\mathrm{ c}$| computed as that of maximum curvature of the FMD + 0.2 according to Woessner & Wiemer (2005) and applied the b⁺ estimator to all positive magnitude differences. He obtained similar results in the two cases. However, as the true b-value of real sequences is not known, it is not possible to decide which strategy is the best. Both selection strategies imply a reduction of the number of remaining data available for b-value determination but, according to the computation made by van der Elst (2021) for the 2019 Ridgecrest Sequence, the first strategy (⁠|${\rm{\Delta }}{M}_\mathrm{ c}^{\prime}$| cut and not |${M}_\mathrm{ c}$| cut) would seem to preserve more data.

Tinti & Gasperini (2024), by statistical analyses of simulated data sets, showed that the theoretical b-value is correctly determined when the trimming threshold is equal to the bin size (only |${\rm{\Delta }}M = 0$| discarded), if the cutting magnitude |${M}_\mathrm{ c}$| is assumed to be that of maximum curvature of the FMD (without + 0.2 correction), and when the data set is strongly incomplete data sets (without |${M}_\mathrm{ c}$| cut), if using trimming thresholds four to five times the binning size. Such analyses also indicate that the number of remaining data is 20–30 per cent larger in the former case. Hence in the following, we will adopt the strategy of cutting magnitudes below the maximum curvature of the FMD and considering all positive differences (⁠|${\rm{\Delta }}M > 0$|⁠).

SETTING UP THE BVAL FORECASTING METHOD

We first remove all magnitudes below the point of maximum curvature of the FMD determined as reported in section analysis region and earthquake catalogue. The b-value is then calculated in each tessellation cell starting from the beginning of the catalogue, using a sampling window with a fixed number of events, |${N}_m = 150$| according to Gulia et al. (2024). For each sampling window, we compute the difference vector D obtained by subtracting the magnitude of each earthquake from the magnitude of the preceding one, that is, |${D}_i\ = \ ( {{m}_i\ - \ {m}_{i - 1}} )$| for |$i\ = \ 2,\ 3,\ \ldots,\ {N}_m$|⁠. We then filter out any negative or zero difference and compute the average of the remaining ones denoted as |$\overline {{\rm{\Delta }}M} $|⁠. The latter is used in eq. (1) to calculate the b-value. We create a temporal series of the b-values by moving forward the sampling windows one event at time and recalculating the parameter for each new window until the end of the data set is reached. According to this forecasting model, alarms are triggered each time the b-value falls below a critical threshold value |${b}_{\mathrm{ th}}$|⁠, coming from a larger value. The alarm lasts for a duration |$\Delta t$|⁠, starting from the time of the last earthquake in the sample of |${N}_m$| events. If the b-value in the previous window was already below the threshold, no alarm is issued to avoid repeated warnings for the same low b-value condition. The choice of the optimal threshold |${b}_{\mathrm{ th}}$| is determined in a specific calibration step where various thresholds are tested using the learning data set to find the value that best balances the sensitivity (detecting potential seismic events) and specificity (minimizing false alarms). Using the calibrated threshold, we evaluate pseudo-prospectively the efficiency of the forecasting method by applying it to the testing data set.

TESTING ALARM-BASED FORECASTING METHODS

A target earthquake is considered successfully predicted if it occurs within a cell where one or more alarm windows of duration |${\rm{\Delta }}t$| are active. Conversely, it is considered a failure to predict if it occurs outside any active alarm window. Thus, the prediction of a target earthquake is evaluated according to a dichotomy: predicted or not predicted.

According to Molchan (1990, 1991) and Molchan & Kagan (1992), the miss rate |$\nu $|⁠, that is the fraction of unpredicted earthquakes, is computed as:

where N is the total number of target events and h is the number of target events that were successfully predicted. Additionally, the overall fraction of space-time occupied by alarms is computed as:

where K is the total number of cells, T is the total duration of the forecasting experiment and |$d{c}_i$| is the total time coverage of alarms within each i-th spatial cell. |$d{c}_i$| can be calculated by multiplying the alarm time |${\rm{\Delta }}t$| by the number |$n\ $| of defined alarms and then subtracting the sum of the intersections between the different alarm windows |$\cap {t}_\mathrm{ s}$| as follows:

Following Shebalin et al. (2011), we also calculate a weighted fraction of space-time occupied by alarms, which accounts for the long-term earthquake rate within each cell. Such rates are computed using the Italian historical earthquake data in the period 1620–1959, according to the CPTI15 V4.0 catalogue (Rovida et al. 2020, 2022; see the Data Availability section), following the procedure described in Gasperini et al. (2021). The weighted fraction of space-time occupied by alarms is then computed as:

where |${\bar{\lambda }}_i$| is the long-term average rate of earthquakes |${M}_\mathrm{w} \ge 4.0$| for each i-th cell. Details of weight computations are reported in Table S1 (Supporting Information). This weighting procedure penalizes forecasts in highly seismic areas, where predictions are inherently easier, by increasing the space-time fraction occupied by alarms proportionally to the historical earthquake rate.

The miss rate |$\nu $| and the fractions of space-time occupied by alarms |${\tau }_\mathrm{ u}$| (unweighted) or |${\tau }_\mathrm{ w}$| (weighted) are used to construct the so-called Molchan diagram (Molchan 1990, 1991, 1997; Molchan & Kagan 1992; Molchan & Keilis-Borok 2008). The latter is a graphical tool commonly used to assess the forecasting performance of alarm-based models (see for example the one reported in Fig. 5).

In a paradoxical forecasting method that does not issue any alarm, we have |$\tau = 0$| and it is impossible to forecast any target event. This corresponds to the point |$( {\tau ,\ \nu } ) = ( {0\ \mathrm{ per}\ \textrm{cent},\ 100\ \mathrm{ per}\ \textrm{cent}} )$| in the upper left corner of the diagram of Fig. 5. On the other hand, if the method issues a set of alarms that cover the entire space-time volume, all target earthquakes are certainly predicted. This corresponds to the point |$( {\tau ,\ \nu } ) = ( {100\ \mathrm{ per}\ \textrm{cent},\ \ 0\ \mathrm{ per}\ \textrm{cent}} )$|⁠. The diagonal line connecting these two points of equation |$\nu = 1 - \tau $| represents the prediction performance of a purely random method that simply forecast earthquakes proportionally to the space-time fraction occupied by alarms.

By analysing different alarm durations |$\Delta t$|⁠, we calculated the corresponding values of |$\tau $| and |$\nu $| for each duration, allowing us to plot a set of points on the Molchan diagram forming a path known as the `Molchan trajectory`. Models with trajectories that lie below the diagonal line and are close to the lower left-hand corner of coordinates |$( {\tau ,\nu } ) = ( {0\ \mathrm{ per}\ \textrm{cent},0\ \mathrm{ per}\ \textrm{cent}} )$| show higher predictive capabilities than the random prediction model. In contrast, models with trajectories extending above the diagonal line show lower predictive performance than the random model (see Fig. 5). The ratio between the success rate and the space-time fraction occupied by alarms,

indicates the probability gain compared to the random chance (Kagan 2009). A value of |$G = 1$| aligns with the random chance diagonal, while |$G > 1$| characterises skilled forecasting methods|$.$|

The Molchan diagram provides insights into how the miss rate fluctuates with changes in the fraction of space-time covered by alarms, highlighting the essential trade-off between alarm coverage and minimizing failure rates. However, it remains primarily a visual tool. To quantitatively assess the forecasting performance of models, Zechar & Jordan (2008, 2010) introduced the area skill (AS) score |${a}_\mathrm{f}( \tau )$|⁠. The latter quantifies how much better or worse a model performs compared to a random prediction baseline and is defined as the integral of the success rate |$1 - \nu ( \tau )$| normalized to the alarm space-time coverage |$\tau $|⁠:

The overall AS score, calculated across the entire Molchan trajectory ranging from |$\tau = 0$| per cent to |$\tau = 100$| per cent, quantifies the area above the Molchan trajectory. The AS score is normalized such that its value ranges from 0 to 1, with values closer to 1 indicating better forecasting performance than the reference random model, represented by the diagonal line (Fig. 5), which is characterized by an overall AS score of 0.5.

Zechar & Jordan (2008, 2010) also derived a Gaussian asymptotical estimate of the variance of |${a}_\mathrm{f}( \tau )$| as |${\sigma }^2 = \frac{1}{{12}}N$|⁠, where N indicates the number of target shocks. While the AS score can be computed for any |$\tau $|⁠, Zechar & Jordan (2008, 2010) suggested that the statistical power of the test generally improves with larger values of |$\tau $|⁠. Therefore, they recommended to use |${a}_\mathrm{f}( {\tau = 100\ \mathrm{ per}\ \textrm{cent}} )$| for hypothesis testing to maximize the sensitivity of the analysis.

CALIBRATION AND OPTIMIZATION OF FORECASTING MODELS

BVAL, FORE and ETAS models are applied to forecast target earthquakes with |$M_{ \mathrm{ w}} \ge 5.0$|⁠, which occurred inland in Italy during the learning period 1990–2004. As possible precursory seismicity, also the data in the preceding 15 yr interval from 1975 to 1989 is considered.

Following the approach of Biondini & Gasperini (2023), we used the maximum AS score as the criterion for selecting the optimal forecasting thresholds for each analysed model. Additionally, we also provide information on the fraction of space-time (⁠|${\tau }_{1\mathrm{ yr}}$|⁠) occupied by alarms of one year duration (⁠|${\rm{\Delta }}t = \ $|1 yr) for each analysed threshold. This metric indicates the proportion of time and geographic area where alarms were issued. Even if not used for the optimal threshold selection, this additional information provides a comprehensive view of the operational impact of the chosen thresholds. For example, a high fraction of space-time occupied by alarms may indicate a high number of false alarms, potentially reducing the model practical utility. Conversely, a too-low fraction may suggest that the chosen precursor is overly prudential and may miss target events.

For the BVAL model, in Fig. 6 (with numerical values reported in Table S2 of Supporting Information) we show the behaviour of the AS scores (lines) and of |${\tau }_{1\mathrm{ yr}}$| (grey bars) with the b-value threshold |${b}_{\mathrm{ th}}$| varying from 0.50 to 1.30. The red and blue lines indicate the AS scores as a function of the unweighted and weighted fraction of space-time occupied of alarms, respectively. For both curves the highest overall scores (⁠|${a}_\mathrm{f}( {{\tau }_\mathrm{u} = 100{\rm{\ per\ cent}}} )$|=0.81 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100{\rm{\ per\ cent}}} )=$|0.79, respectively) are obtained for |${b}_{\mathrm{ th}} = 0.90$|⁠. For such threshold the fraction of space-time covered by alarms of duration |${\rm{\Delta }}t = \ $|1 yr is |${\tau }_{1\mathrm{ yr}} = 0.017$| that is 1.7 per cent of the total duration of the learning time interval.

For the FORE model, an alarm is issued in a spatial cell every time a potential foreshock with |$M_{ \mathrm{ w}} = M \pm {\rm{\Delta }}M$| occurs within that cell (Gasperini et al. 2021). The optimal range of potential foreshocks was determined by comparing the AS scores obtained in the learning interval by using magnitude ranges with central values |$4.1 \le M \le 4.8$|⁠, and deviations from the central value |$0.1 \le {\rm{\Delta }}M \le 0.4$| but ensuring that the upper bound is not larger than 4.9 (⁠|$M + {\rm{\Delta }}M \le 4.9$|⁠), so that the magnitude of the foreshocks is lower than that of any target earthquakes (Fig. 7; with values in Table S3 of Supporting Information). The maximum overall AS scores (⁠|${a}_\mathrm{f}( {{\tau }_\mathrm{u} = 100{\rm{\ per\ cent}}} )=$|0.90 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100{\rm{\ per\ cent}}} )$|0.88) are obtained for |$M_{ \mathrm{ w}} = 4.5 \pm 0.3$|⁠. For such ranges, the fraction of space-time covered by alarms of duration |${\rm{\Delta }}t = \ $|1 yr is |${\tau }_{1\mathrm{ yr}} = 0.045$| that is 4.5 per cent of the total duration of the learning time interval.

For the ETAS model, an alarm is issued, in a spatial cell, whenever the expected daily rate (⁠|${\lambda }_{\textrm{daily}}$|⁠) exceeds a given threshold. ETAS parameters were fitted using the maximum likelihood approach (see details in Appendix  A) for the learning period (1990–2004), separately for the two tessellations, R₁ and R₂ (see Table A1). Following Biondini & Gasperini (2023), the expected rate of earthquakes with |${M}_\mathrm{w} \ge 5.0$| is recalculated for each cell whenever an earthquake with |${M}_\mathrm{w} \ge 2.5$| occurs within the entire observation region R = (⁠|${R}_1 + {R}_2)$|⁠. Each cell is thus characterised by its own time-series of expected |${\lambda }_{\textrm{daily}}$|⁠.

The optimal alarm threshold |${\lambda }_{\mathrm{ th}}$| for the expected daily earthquake rate is determined by varying it logarithmically from |$5 \times {10}^{ - 7}$| to |$9 \times {10}^{ - 3}$| (Fig. 8; with values in Table S4 in the Supporting Information). The maximum overall AS scores (⁠|${a}_\mathrm{f}( {{\tau }_\mathrm{u} = 100{\rm{\ per\ cent}}} )=$| 0.88 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100{\rm{\ per\ cent}}} )=$| 0.83) are obtained for |${\lambda }_{\mathrm{ th}} = 1 \times {10}^{ - 5}$|⁠. For such threshold the fraction of space-time covered by alarms of duration |${\rm{\Delta }}t = \ $| 1 yr is |${\tau }_{1\mathrm{ yr}} = 0.337$| that is 33.7 per cent of the total duration of the learning time interval. Such |${\tau }_{1\mathrm{ yr}}$|⁠, definitely larger than the one obtained by optimizing other two methods, highlights, right from this optimization stage, a different behaviour of the ETAS model compared to the other two forecasting models, which might condition its practical applicability.

RESULTS OF THE PSEUDO-PROSPECTIVE TESTING

We used the optimal thresholds and ranges determined in the previous step for applying the simple model BVAL, FORE and ETAS to pseudo-prospectively predict 34 earthquakes with magnitudes |${M}_\mathrm{w} \ge 5.0$| and 12 earthquakes with magnitudes |${M}_\mathrm{w} \ge $| 5.5 occurred during the 18 yr of the pseudo-prospective test period 2005–2022. As possible precursor seismicity of targets events in the testing period 2005–2022, we used also the data in the preceding 18 yr interval from 1986 to 2004.

We use the same optimal thresholds even for applying the two ensemble models EADD and EMUL, obtained by summing and intersecting, respectively, the alarm windows of the two simple BVAL and FORE models.

In Fig 5, the Molchan trajectories of the BVAL model applied to the prediction of target earthquakes |$M \ge 5.0$| are shown. They are obtained by varying the alarm time |$( {{\rm{\Delta }}t} )$| from fractions of seconds to the total duration of the pseudo-prospective prediction experiment (T = 18 yr). The red and blue lines correspond to unweighted (⁠|${\tau }_\mathrm{u}$|⁠) and weighted (⁠|${\tau }_\mathrm{w}$|⁠) fractions of space-time occupied by alarms, respectively (see in Table S5 of the Supporting Information, the numerical values of the plotted trajectories). Both the red and blue Molchan trajectories lie well below the diagonal line representing the predictive performance of a purely random method and even below the 1 per cent confidence level threshold, highlighting their higher predictive ability. The BVAL model fails to predict the total number of target earthquakes with |${M}_\mathrm{w} \ge 5.0$| even with |${\rm{\Delta }}t$|=18 yr, for which only 30/34 (88.2 per cent) are predicted with time-space coverages of |${\tau }_\mathrm{u} \approx $| 33 per cent and |${\tau }_\mathrm{w} \approx $| 51 per cent. With |${\rm{\Delta }}t$| = 1 yr, approximately 68 per cent of target earthquakes are successfully forecasted, achieving space-time coverages of |${\tau }_\mathrm{u}$|= 8.6 per cent and |${\tau }_\mathrm{w}$|= 12.7 per cent. For a |${\rm{\Delta }}t$| = 3 months (0.25 yr), 64.7 per cent of targets earthquakes are successfully forecasted with |${\tau }_\mathrm{u}$| = 2.96 per cent and |${\tau }_\mathrm{w}$| = 4.45 per cent. Using |${\rm{\Delta }}t$| = 1 d, 50 per cent of target earthquakes are successfully forecasted, with |${\tau }_\mathrm{u}$| = 0.05 per cent and |${\tau }_\mathrm{w}$| = 0.09 per cent. The overall AS scores |${a}_\mathrm{f}( {{\tau }_{\mathrm{ u}} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.88 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.85 are significantly larger of the expected score of the random model (0.5), based on the Student's t-test. The Molchan diagram for the BVAL model applied to predict target earthquakes with |${M}_\mathrm{w} \ge 5.5$| is shown in Fig. S1 of Supporting Information. It shows high prediction skills, with overall AS scores of |${a}_\mathrm{f}( {{\tau }_\mathrm{u} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.92 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.90 for the unweighted and weighted fraction of space-time occupied by alarms, respectively. Even with |${\rm{\Delta }}t$| =18 yr, the model forecasts only 11 of the 12 target events (92 per cent), with |${\tau }_\mathrm{u}$| and |${\tau }_\mathrm{w}$| remaining unchanged from the previous case due to consistent alarm issuance. For |${\rm{\Delta }}t$| = 1 yr and |${\rm{\Delta }}t$| = 1 d, 83 per cent (10/12) and 67 per cent (8/12) of target shocks are predicted, respectively. The numerical values of the plotted miss rates |$\nu $| and of the AS scores as a function |${\tau }_\mathrm{u}$| and |${\tau }_\mathrm{w}$| for targets with |${M}_\mathrm{w} \ge 5.0$| And |${M}_\mathrm{w} \ge 5.5$| are reported in Tables S5 and S6 (Supporting Information), respectively.

For the FORE model, Fig. 9 shows the Molchan trajectories for target earthquakes |${M}_\mathrm{w} \ge 5.0$|⁠. As for the BVAL model, also the FORE model overperform the prediction performance of the random method since the Molchan trajectories for both unweighted and weighted fraction of space-time occupied by alarms lies well below the diagonal line and the 1 per cent confidence level threshold. Unlike the BVAL model, the FORE model predicts all the target earthquakes with |${M}_\mathrm{w} \ge 5.0$| for |${\rm{\Delta }}t$|=18 yr, for which |${\tau }_\mathrm{u}$| = 41.5 per cent and |${\tau }_\mathrm{w}$| = 58.2 per cent. For |${\rm{\Delta }}t$|= 1 yr, approximately 79 per cent (27/34) of targets are predicted with |${\tau }_\mathrm{u}$| = 3.4 per cent and |${\tau }_\mathrm{w}$| = 5.5 per cent. With |${\rm{\Delta }}t$| = 3 months, 74 per cent (25/34) of target earthquakes are predicted with |${\tau }_\mathrm{u}$| = 0.95 per cent and |${\tau }_\mathrm{w}$|=1.58 per cent. For |${\rm{\Delta }}t$| = 1 d, 47 per cent (16/34) of target earthquakes are predicted with |${\tau }_{\mathrm{ u}}$| = 0.01 per cent and |${\tau }_\mathrm{w}$| = 0.03 per cent. The overall AS scores |${a}_\mathrm{f}( {{\tau }_\mathrm{u} = 100\ \mathrm{ per}\ \textrm{cent}} )$|=0.95 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.92 are significantly higher than the random model's expected AS score (0.5). The Molchan diagram of the FORE model applied to |${M}_\mathrm{w} \ge $| 5.5 earthquakes (Fig. S2, Supporting Information) show very high prediction skills, with |${a}_\mathrm{f}( {{\tau }_\mathrm{u} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.98 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100{\rm{\ per\ cent}}} )$| = 0.97. The total number of target events (12) are predicted starting from |${\rm{\Delta }}t$| = 9 yr with |${\tau }_\mathrm{u}$| = 24.8 per cent and |${\tau }_\mathrm{w}$| = 36.3 per cent. For |${\rm{\Delta }}t$| = 1 yr and |${\rm{\Delta }}t$| = 1 d the 92 per cent (11/12) and the 58 per cent (7/12) of target shocks are predicted. Detailed numerical values and AS scores are provided in Tables S7 and S8 (Supporting Information).

For the ETAS model, Fig. 10 displays the Molchan trajectories for target earthquakes with |${M}_\mathrm{w} \ge 5.0$|⁠. Like the two previous models, ETAS outperforms the random prediction method, as both the unweighted and weighted Molchan trajectories lie well below the diagonal line and even below the 1 per cent confidence level threshold. It predicts all target events with |${\rm{\Delta }}t$| = 3 d or larger, with |${\tau }_\mathrm{u}$| = 30.1 per cent and |${\tau }_\mathrm{w}$| = 48.1 per cent. With |${\rm{\Delta }}t$| = 1 yr, the fraction of space-time occupied by alarm increase to |${\tau }_{\mathrm{ u}}$| = 33.4 per cent and |${\tau }_\mathrm{w}$| = 53.2 per cent. For |${\rm{\Delta }}t$| = 1 d, 97 per cent (33/34) of target earthquakes are predicted with |${\tau }_\mathrm{u}$| = 20.8 per cent and |${\tau }_\mathrm{w}$| = 33.9 per cent. The overall AS scores are |${a}_\mathrm{f}( {{\tau }_\mathrm{u} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.97 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.95. The Molchan diagram for the ETAS model applied to |${M}_\mathrm{w} \ge $| 5.5 earthquakes (Fig. S3, Supporting Information) show very high prediction skills, with |${a}_\mathrm{f}( {{\tau }_\mathrm{u} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.99 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.98. All 12 target events are predicted starting from |${\rm{\Delta }}t$| = 0.125 d (3 hr), with corresponding |${\tau }_\mathrm{u}$| = 5.0 per cent and |${\tau }_\mathrm{w}$| = 8.0 per cent. Detailed numerical values and AS scores are provided in Tables S9 and S10 (Supporting Information).

For the ensemble EADD model, Fig. 11 shows the Molchan trajectories for target earthquakes |${M}_\mathrm{w} \ge 5.0$|⁠. As for other models, it overperform the prediction performance of the random method. It predicts all the 34 target earthquakes with |${M}_\mathrm{w} \ge 5.0$| from |${\rm{\Delta }}t$| = 16 yr for which |${\tau }_\mathrm{u}$|  |$= 48.4$| per cent and |${\tau }_\mathrm{w} = 67.3$| per cent. For |${\rm{\Delta }}t$| = 1 yr, approximately 82 per cent (28/34) are predicted with |${\tau }_\mathrm{u}$| = 10.2 per cent and |${\tau }_\mathrm{w}$| = 15.1 per cent. For |${\rm{\Delta }}t$| = 3 months, it predicts 76 per cent (26/34) of target earthquakes with |${\tau }_\mathrm{u}$| = 3.5 per cent and |${\tau }_\mathrm{w}$| = 5.3 per cent. For |${\rm{\Delta }}t$| = 1 d, 68 per cent (23/34) of target earthquakes are predicted with |${\tau }_\mathrm{u}$| = 0.06 per cent and |${\tau }_\mathrm{w}$| = 0.11 per cent. The overall AS scores |${a}_\mathrm{f}( {{\tau }_{\mathrm{ u}} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.96 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.94 are significantly higher than the random model expected AS score (0.5). The Molchan diagram of the EADD model applied to |${M}_\mathrm{w} \ge $| 5.5 earthquakes (Fig. S4, Supporting Information) show very high prediction skills, with |${a}_\mathrm{f}( {{\tau }_\mathrm{u} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.97 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.96. The total number of target events (12) are predicted starting from |${\rm{\Delta }}t$| = 6 yr with |${\tau }_\mathrm{u}$| = 31.0 per cent and |${\tau }_\mathrm{w}$| = 45.0 per cent. For |${\rm{\Delta }}t$| = 1 yr and |${\rm{\Delta }}t$| = 1 d, the 92 per cent (11/12) and the 83 per cent (10/12) of target shocks are predicted, respectively. Detailed numerical values and AS scores are provided in Tables S11 and S12 (Supporting Information).

For the ensemble EMUL model, Fig. 12 shows the Molchan trajectories for target earthquakes |${M}_\mathrm{w} \ge 5.0$|⁠. As for other models, it overperforms the prediction performance of the random method since both the Molchan trajectories for unweighted and weighted fraction of space-time occupied by alarms lies well below the diagonal line. Even with |${\rm{\Delta }}t$| = 18 yr, it predicts only 30 of 34 target earthquakes (88 per cent) with |${M}_\mathrm{w} \ge 5.0$|⁠. For |${\rm{\Delta }}t$| = 1 yr, approximately 65 per cent (22/34) are predicted with |${\tau }_\mathrm{u}$| = 1.7 per cent and |${\tau }_\mathrm{w}$| = 3.0 per cent. For |${\rm{\Delta }}t$| = 3 months predicts 62 per cent (21/34) of target earthquakes with |${\tau }_\mathrm{u}$|=0.4 per cent and |${\tau }_\mathrm{w}$|=0.8 per cent. For |${\rm{\Delta }}t$| = 1 d, 29 per cent (10/34) of target earthquakes are predicted with |${\tau }_\mathrm{u}$|=0.003 per cent and |${\tau }_\mathrm{w}$| = 0.01 per cent. The overall AS scores |${a}_\mathrm{f}( {{\tau }_{\mathrm{ u}} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.89 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.86 are significantly higher than the random model's expected AS score (0.5), The Molchan diagram of the EMUL model applied to |${M}_\mathrm{w} \ge $| 5.5 earthquakes (Fig. S5, Supporting Information) show high prediction skills, with |${a}_\mathrm{f}( {{\tau }_\mathrm{u} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.94 and |${a}_\mathrm{f}( {{\tau }_\mathrm{w} = 100\ \mathrm{ per}\ \textrm{cent}} )$| = 0.92. Even with |${\rm{\Delta }}t$| = 18 yr, it predicts only 11 of 12 target earthquakes (92 per cent) with |${M}_\mathrm{w} \ge 5.5$|⁠. For |${\rm{\Delta }}t$| = 1 yr and |${\rm{\Delta }}t$| = 1 d the 83 per cent (10/12) and the 42 per cent (5/12) of target shocks are predicted, respectively. Detailed numerical values and AS scores are provided in Tables S13 and S14 (Supporting Information).

ANALYSIS OF PERFORMANCE OF MODELS BASED ON ALARM DISTRIBUTION

Based on the overall AS score (computed for |$\tau = 100\ \mathrm{ per}\ \textrm{cent}$|⁠), the ETAS model appears to be the best performing model for predicting target earthquake with both |${M}_\mathrm{w} \ge 5.0$| and |${M}_\mathrm{w} \ge 5.5$|⁠. However, a closer examination reveals significant issues associated with its implementation using the alarm-based approach. During the 18-yr testing period, the ETAS model casts more than 1600 000 alarms. These alarms were concentrated in only 134 out of 337 cells, leading to an average of approximately 1.8 alarms per day in each cell. While this high frequency of alarms explain why ETAS can predict all target earthquakes even with very short alarm times (a few days for |${M}_\mathrm{w} \ge $| 5.0 and a few hours for |${M}_\mathrm{w} \ge $| 5.5) it also results in an overwhelming number of alerts. This could be impractical for a real application, where the large number of alarms may complicate the decision-making process and strain the public resources.

Conversely, the BVAL and FORE models issues significantly fewer alarms during the same period, with the BVAL model generating 1796 alarms (less than one each year on average in each cell) and the FORE 782 alarms (about one every three years on average per cell).

The differences between these models are illustrated by their behaviour in Fig. 13 for cell 179, which includes three targets (blue vertical bars): the 2016 August 24 Amatrice earthquake (⁠|${M}_\mathrm{w} = 6.2$|⁠) and two significant aftershocks occurred on 2017 January 18 (⁠|${M}_\mathrm{w} = 5.3\ {\rm{and}}\ 5.7)$|⁠. All the three methods are able to forecast the three target shocks with a minimum |${\rm{\Delta }}t$| of 6 hr for BVAL and FORE, and of about 30 s for ETAS. ETAS alarms (pins) are almost continuous throughout the 18-yr testing period, even when considering short alarm times |${\rm{\Delta }}t$| of just a few hours or minutes. This nearly constant state of alert reflects ETAS's high sensitivity but also highlights its potential for generating ‘alarm fatigue’, where the frequency of alarms diminishes their credibility (Bliss et al. 1995; Graham & Cvach 2010).

By comparison, the alarms issued by the BVAL and FORE models in the same cell are much sparser. While these models may not capture every possible precursor signal, they offer a more balanced approach, issuing alarms selectively and potentially providing more reliable warnings. This contrast highlights the fundamental trade-off in earthquake forecasting between the comprehensive coverage of ETAS and the more targeted and practical utility of BVAL and FORE. This comparison offers useful insights into the practical effects of each forecasting method, helping to determine which model is the most suitable, based on the desired balance between sensitivity and reliability.

As noted by Zechar & Jordan (2008, 2010), the overall AS Score (computed for |$\tau = 100\ \mathrm{ per}\ \textrm{cent}$|⁠) is a valuable criterion to compare forecasting models because its power generally improves with increasing |$\tau $|⁠. Although the ETAS model achieves high overall AS scores, the frequency and duration of the alarms it generates may limit its overall effectiveness. A forecasting method that can achieve a good predictive efficiency with alarms concentrated in relatively short-time intervals (like BVAL and FORE) could be probably more effective and may be better received by the users.

In Fig. 14, we compare the performance of various methods by reporting the AS scores, for targets with |${M}_\mathrm{w} \ge 5.0$|⁠, for |${\tau }_\mathrm{w}$| varying from 0.01 per cent to 100 per cent. For |${\tau }_\mathrm{w}$| < 20 per cent, the higher scores are obtained by the FORE and EADD methods. Even though the ETAS model has the highest score for |${\tau }_\mathrm{w} \ge 20$| per cent, for |${\tau }_\mathrm{w}$|<4 per cent it results even worse than BVAL and EMUL. Note that |${\tau }_\mathrm{w} = 0.01$| per cent means that in each cell there is on average less than one hour of alarm each year, |${\tau }_\mathrm{w} = 0.10$| per cent, less than 9 hr each year, |${\tau }_\mathrm{w} = 1$| per cent, less than 4 d each year, and |${\tau }_\mathrm{w} = 10$| per cent about 36 d.

Similar inferences can be drawn from Fig. 15 where the rate of successful predictions (hit rate, |$h = 1 - \nu $|⁠) is plotted versus |${\tau }_\mathrm{w}$|⁠. In this case, ETAS is the worst model for |${\tau }_\mathrm{w}$| < 2 per cent and becomes the best one for |${\tau }_\mathrm{w}$| > 5 per cent, which mean about 8 and 20 d of alarm each year, respectively, on average in each cell. While ETAS excels in scenarios with extended alarm durations, methods like FORE and EADD may be more suitable for situations where the goal is to maintain high predictive accuracy with minimal disruption to daily life, especially where frequent alarms could lead to public desensitization.

MOLCHAN-SHEBALIN DIAGRAM

Following the methodology introduced by Shebalin et al. (2011), we can compare the miss rates (⁠|$\nu $|⁠) of various analysed models to a skilled reference model (for this study, the ETAS model) playing the role of the purely random prediction on the Molchan diagonal. The skilled reference model is characterized by its space-time fractions occupied by alarms (⁠|${\tau }_\mathrm{ref}$|⁠) and miss rates (⁠|${\nu }_\mathrm{ref}$|⁠). The miss rates (⁠|${\nu }_\mathrm{ref}$|⁠) are plotted along the diagonal line of the Molchan diagram, where |${\tau }_\mathrm{ref} = 1 - {\nu }_\mathrm{ref}$|⁠. For the other models being compared, their expected miss rates (⁠|$\nu $|⁠) must also be plotted against |$\tau = 1 - {\nu }_\mathrm{ref}$|⁠. Since these models may have been computed at different |$\tau $| values compared to the |${\tau }_\mathrm{ref}$| values of the ETAS model, linear interpolation is required to align the |$\tau $| values appropriately. Specifically, for each |${\tau }_\mathrm{ref}$| of the ETAS reference model, we compute the miss rates (⁠|${\nu }_{\mathrm{ int}}$|⁠) for the other models (BVAL, FORE, EADD and EMUL) using the following linear interpolation formula (Biondini & Gasperini 2023):

where |${\tau }_\mathrm{ a}$| and |${\tau }_\mathrm{ b}$| are the fractions of space-time occupied by alarms for the compared models, immediately larger and smaller than |${\tau }_\mathrm{ref}$|⁠, respectively, and |${\nu }_\mathrm{ a}$| and |${\nu }_\mathrm{ b}$| the corresponding miss rates, respectively. We make such interpolation using as |${\tau }_\mathrm{ref}\ $| both |${\tau }_\mathrm{u}$| and |${\tau }_\mathrm{u}$| of the reference ETAS model. The rationale for this interpolation is to provide a consistent basis for comparing the miss rates of different models at the same values of |$\tau $|⁠. Similarly to the Molchan diagrams showed above (Figs 5, and 9 –12), if the miss rates of a compared model |$({\nu }_{\mathrm{ int}})$| for the corresponding |$\tau = 1 - {\nu }_\mathrm{ref}$|⁠, are lower than |${\nu }_\mathrm{ref}$|⁠, such model has a better predictive ability than the reference model. This method provides a rigorous framework for evaluating the predictive capabilities of seismic forecasting models, allowing for a fair comparison based on equivalent conditions. For each model, we also compute a sort of a Molchan–Shebalin area skill (now on MSAS) score, according to eq. (7) where |$\nu $| is substituted by |${\nu }_{\mathrm{ int}}$| and |$\tau $| by |${\tau }_\mathrm{ref}$|⁠. For the reference ETAS model |${\nu }_{\mathrm{ int}} = {\nu }_\mathrm{ref} = 1 - {\tau }_\mathrm{ref}$|⁠, then its AS score is exactly 0.5.

Using the weighted |${\tau }_\mathrm{w}$| of the ETAS model for |${M}_\mathrm{w} \ge 5.0$| as |${\tau }_\mathrm{ref}$|⁠, in Fig. 16 (with values in Table 2) we show that all other models lie below the diagonal representing the ETAS model for |${\tau }_\mathrm{ref}$| < 60 per cent and the FORE and EADD models for |${\tau }_\mathrm{ref}$| < 80 per cent. All models have an overall MSAS score higher than the reference model ETAS. For FORE and EADD models, the score is significantly higher than for the ETAS according to the Student's t-test. The same plot for the |${M}_\mathrm{w} \ge 5.5$| is reported in Fig. S6 of the Supporting Information. The results are like those of Fig. 16 with a slightly better performance of BVAL, FORE, EADD and EMUL with respect to ETAS.

DISCUSSION AND CONCLUSIONS

We developed the BVAL method to forecast potentially damaging target earthquakes (⁠|${M}_\mathrm{w} \ge 5.0$|⁠) in Italy, based on the temporal variations of the slope (b-value) of the frequency–magnitude distribution (Gutenberg & Richter 1944), and we compared it with other forecasting methods based on the statistical properties of seismicity: the FORE method using the occurrence of strong foreshocks (Gasperini et al. 2021), the ETAS method using the epidemic principle that any shock occurred acts as a trigger of the following ones (Ogata 1988), and two ensemble methods obtained by the union (EADD) and the intersection (EMUL) of space time alarm windows of BVAL and FORE.

We were inspired by the hypothesis formulated by Scholz (1968, 2015) and Amitrano (2003) that b-value decreases with increasing differential stress and with the increasing probability of occurrence of strong earthquakes. The same hypothesis brought to the formulation of the Foreshock Traffic Light System (Gulia & Wiemer 2010, 2019, 2021; Gulia et al. 2018) predicting whenever a shock with |${M}_\mathrm{w} \ge 6.0$| is the strongest in a seismic sequence (green light, high b-value) or an even stronger shock is likely to occur close in space in the near future (red light, low b-value).

Our algorithm declares an alarm of duration |${\rm{\Delta }}t$|⁠, which varies from a fraction of second to the entire duration of the experiment, in a square cell with sides of about 42.4 km, every time the b-value in that cell decreases below the threshold |${b}_{\mathrm{ th}} = \ $| 0.9 from a higher value. Such threshold was determined by optimizing the ability of the algorithm to forecast earthquakes with |${M}_\mathrm{w} \ge 5.$|0 of the HORUS catalogue (Lolli et al. 2020) from 1990 to 2004.

If a target earthquake occurs in a cell where an alarm window is active, we count a hit, if a target earthquake occurs outside any alarm window we count a miss. To not miss the prediction of target earthquakes located only slightly outside an alarmed cell, we used a double tessellation with two grids staggered by half the grid side.

We pseudo-prospectively evaluated the predictive efficiency of the forecasting method using the AS score (Zechar & Jordan 2008, 2010) measuring the area above the Molchan trajectory (Molchan 1990, 1991), based on the data of the HORUS catalogue from 2005 to 2022. We obtained the Molchan trajectory by plotting the miss rate |$\nu $| as a function of the fraction of space-time occupied by alarms |$\tau $|⁠, varying the alarm duration |${\rm{\Delta }}t$| from a fraction of second to the total duration of the testing period (18 yr). The lower the miss rates of the Molchan trajectory the larger the AS score and then the more accurate the forecasting method.

We also penalized the forecasting performance of models in highly seismic areas, where random forecasts are easier, by weighting the fraction of space-time occupied by alarms proportionally to the long-term earthquake occurrence rate deduced from the catalogue CPTI15 V4.0 (Rovida et al. 2020, 2022) from 1600 to 1959.

Based on the overall AS score (computed for |$\tau = 100$| per cent) the best performing model would be the ETAS, which achieves a near-perfect AS score (97–99 per cent). While this might indicate ETAS as the most effective model, it might also evidence a critical limitation of the Molchan diagram for testing models, like ETAS, issuing extremely frequent and dense alarms which results in a rapid reduction of the miss rate to zero even for notably short alarms of a few seconds or minutes. This creates an illusion of perfect predictive capability, as all target earthquakes are technically predicted, but at the cost of an overwhelming and nearly constant state of alarm in some cells.

This characteristic of ETAS, while beneficial in a purely statistical sense, undermines its broader applicability. The constant stream of alarms not only risks desensitising the population and creating ‘alarm fatigue’, effective management of emergency responses to natural disasters. In highly alarmed cells, the ETAS model alarms are so continuous that it becomes difficult to distinguish between periods of genuine increased risk and the model default state of alert. This diminishes its overall effectiveness for the operational utility where the goal is to balance predictive accuracy and minimizing unnecessary alerts.

More in general, we can note that high fractions (e.g. >10 per cent) of space-time occupied by alarms may lead to a diminished trust in the model, as the constant presence of alarms could make it harder to discern meaningful signals from the background noise. Hence, we also computed the AS scores of BVAL and other forecasting models even at smaller |$\tau $| down to 0.01 per cent.

If compared with other methods, the BVAL is worse than FORE at any |$\tau $| but better than ETAS for |$\tau < 2 - 4$| per cent. The ensemble model EADD obtained by the sum of alarms issued by both BVAL and FORE is the best performing for space time occupations between 0.04 per cent to 2 per cent, which is the range of time-space occupation for which the forecasting algorithms might be more useful for civil protection purposes.

This preference for other methods with respect to ETAS is also confirmed by the Molchan–Shebalin diagram (Shebalin et al. 2011) in which the role of the random forecast along the Molchan diagonal is played by a reference skilled method (in our case the ETAS method itself). In fact, all other methods describe Molchan trajectories below the reference ETAS diagonal for |$\tau $| < 60 per cent and have the area above the trajectory, larger than 0.5 that corresponds to the ETAS reference model score. In summary, the BVAL forecasting method works rather well so confirming the validity of the basic hypothesis that low b-values might indicate the future occurrence of potentially damaging earthquakes. It is anyhow slightly worse than the FORE method based on the occurrence of strong foreshocks (Gasperini et al. 2021). However, BVAL in additive combination with the FORE improves the efficiency of the latter, particularly at relatively low fractions (<5 per cent) of space-time occupied by alarm that are those that are most likely to be useful for seismic hazard mitigations.

ACKNOWLEDGEMENTS

We thank Rodolfo Console and Aron Mirwald for thoughtful comments and suggestions that helped to improve the paper. This paper benefitted from funding provided by the European Union within the ambit of the H2020 project RISE (no. 821115) and the ETHZ, which in particular fully financed the Postdoc grant of one of the authors (BE). We thank the Istituto Nazionale di Geofisica e Vulcanologia, Italy, grant "Progetto INGV Pianeta Dinamico (NEMESIS)’—code CUP D53J19000170001–funded by Italian Ministry MIUR.

CONFLICT OF INTEREST

The authors declare that there is no conflict of interest regarding the publication of this manuscript.

DATA AVAILABILITY

The data of the HORUS and CPTI15 earthquake catalogues are available from public providers: http://horus.bo.ingv.it and https://emidius.mi.ingv.it/CPTI15-DBMI15/, respectively.

The QGIS geographical software was downloaded from www.qgis.org.

REFERENCES

APPENDIX A: ETAS MODEL FORMULATION AND PARAMETERS FITTING

The ETAS formulation and the background model SVP

The ETAS model, introduced by Ogata (1988), is a stochastic point process used to model aftershock sequences and provide probabilistic estimates of earthquake occurrences. The model addresses the limitation of using a single Omori–Utsu law to describe the temporal decay of aftershock sequences. Guo & Ogata (1997) and Ogata et al. (2003) demonstrated that significant aftershocks often generate their own secondary sequences of aftershocks. As a result, the ETAS model treats each aftershock as a potential trigger for additional aftershocks. The occurrence rate at a given time t is derived from the superposition of time-shifted Omori–Utsu decay functions, each weighted by the magnitude of the parent event (Ogata & Zhuang 2006).

The ETAS model is widely used for earthquake forecasting, including various applications in Italy (e.g. Console et al. 2006; Lombardi & Marzocchi 2010; Falcone et al. 2021; Biondini et al. 2023). The temporal dependence in these models is typically represented by a sum of Omori's law decays, initiated at the times of each earthquake occurrence:

where K, c, p and |$\alpha $| are free parameters, and |$H( {t - {t}_i} )$| is the Heaviside step function, which equals 1 if |$t - {t}_i > 0$| and 0 otherwise. n in the sum represents the total number of earthquakes considered, with each |$i - $|th term in the sum corresponding to an individual earthquake in the sequence. The earthquake productivity is assumed to be an exponential function of the magnitude of the |$i - $|th parent earthquake |${m}_i$| and of the completeness magnitude |${m}_\mathrm{ c}$|⁠.

The spatial decay of the productivity can be described using various probability density function that represent the nucleation kernel (Ogata 1988; Console et al. 2003; Zhuang et al. 2004, 2011; Lombardi & Marzocchi 2010). In this work, we used the spatial probability density function proposed by Ogata & Zhuang (2006), where the smoothing term is an exponential function of the magnitude:

where q, D and |$\gamma $| are free parameters.

The frequency magnitude distribution of earthquakes adheres to the Gutenberg–Richter law (1944):

where |$\beta = b\ {\rm{ln}}10$| and b is a free parameter.

By incorporating a time-invariant background seismicity term |${\lambda }_0( {x,y,m} )$|⁠, the rate density function for the ETAS model is given by:

which explicitly becomes:

The parameter μ represents the ratio between the expected number of independent events of the background seismicity |${\lambda }_0( {x,y,m} )$| and the total number of events.

To describe the background seismicity, we employ the spatially variable poisson (SVP) model. The SVP rate density is given by:

where |${\mu }_0( {x,y} )$| is the spatial rate density of earthquake with |$m \ge {m}_\mathrm{ c}$|⁠, which is a continuous smoothed function of the geographic location. To assess the spatial rate density, we divide each cell of the tessellation into 16 subcells of side |$l = \frac{{30\sqrt 2 }}{4}$| km. The number of earthquakes |${N}_k$| with magnitude |$m \ge {m}_\mathrm{ c}$| is calculated for each subcell. Each |${N}_k$| value, representative of a single subcell, is then smoothed by a Gaussian filter with correlation distance |${{\it{d}}}_\mathrm{ c}$| and normalized to preserve the total number of events as described in Frankel (1995). For each cell, the smoothed |${N}_k$| is given by

where |${\Delta }_{kl}$| is the distance between the centre of the |${k}_{{\rm{th}}}$| and the |${l}_{{\rm{th}}}$| cells. To obtain |${\tilde{N}}_k$| in terms of number of events per unit of time and area, it must be divided by the total duration of the earthquake catalogue and by the area of the cell. The spatial density |${\mu }_0( {x,y} )$| at any given point is calculated taking the weighted average of the four nearest surrounding subcells, where the weights are inversely proportional to the distance from each subcell to the point. To determine the optimal |${d}_\mathrm{ c}$| of the SVP background model, we used the method proposed by Console & Murru (2001). This involves dividing the earthquake catalogue into two subcatalogues of approximately equal temporal length. The value of |${d}_\mathrm{ c}$| is then chosen by maximizing the log-likelihood of one subcatalogue using the smoothed seismicity derived from the other subcatalogue. This analysis is performed separately for both subcatalogues, yielding |${d}_{\mathrm{ c}1}$| = 16.0 and |${d}_{\mathrm{ c}2}$| = 13.0. The optimal correlation distance |${d}_\mathrm{ c}$| is given by the average of these two values, resulting in |${d}_\mathrm{ c}$| = 14.5. Given the parameter |${d}_\mathrm{ c}$|⁠, the spatial density of earthquakes |${\mu }_0( {x,y} )$| for the SVP background model can be evaluated for each cell and point in space.

ETAS parameters fitting

The best fit for the free parameters of the ETAS model is achieved by maximizing the log-likelihood of an inhomogeneous Poisson point process, using the data set of earthquakes occurred in the learning period (1990–2004). The log-likelihood function is given by:

where |$\lambda ( {{t}_i,{m}_i,{x}_i,{y}_i} )$| is the rate density function of the ETAS model (eq. A5); |$( {{t}_\mathrm{ a},{t}_\mathrm{ b}} )$| represent the fitting period; |${m}_\mathrm{ T}$| and |${m}_\mathrm{ u}$| are the lower and upper bounds of the magnitude range, set to 2.5 and 7.5, respectively; and R is the spatial region used for parameter fitting.

The analysis region R consists of two different tessellations, which form a grid of partially overlapping cells. Consequently, the fitting of the free parameters was performed separately for each tessellation, |${R}_1$| and |${R}_2$|⁠, (see Fig. A1). The b-value is estimated independently for each tessellation using the b⁺ estimator introduced by van der Elst (2021) and recently explored by Tinti & Gasperini (2024). To assess the b-value, we use the seismicity data from the learning period (1990–2004) within each tessellation, resulting in b-values of 1.035 for |${R}_1$| and 1.033 for |${R}_2$|⁠. The parameter q is set to 1.5 according to Lombardi & Marzocchi (2010) showing that the static stress change decreases with epicentral distance as |${r}^{ - 3}$|⁠. Similarly, the parameters k, p, c, |$\alpha $|⁠, d and |$\nu $| have been fitted by maximizing the log-likelihood of earthquakes with magnitude |$m \ge 2.5$| occurred within each individual tessellation, |${R}_1$| and |${R}_2,\ $| respectively, during the period 1990–2004 using the interior-point method (Wright 2004). The standard error (SE) for each parameter is computed as the square root of the corresponding diagonal element of the variance–covariance matrix computed as the inverse of the Hessian matrix for log-likelihood function at the maximum (Ogata 1988; Lolli et al. 2009). The optimal fitted parameters with their standard deviations are reported in Table A1. Additionally, the expected number of events (E) and the log-likelihood (⁠|$\mathrm{ ln}L$|⁠), shown in Table A1 are calculated based on the fit to earthquakes with magnitudes |$m \ge 5.0$|⁠. This approach ensures a more robust fit by initially using a broader range of magnitudes (⁠|$m \ge 2.5$|⁠), as described above, to obtain a more reliable parameter estimation with a larger data set. Subsequently, the log-likelihood and the expected number of events for earthquakes with |$m \ge 5.0$| are recalculated with the previously fitted parameters,. This two step-process leverages the extensive data set for parameter estimation and then focuses on higher-magnitude events for assessing model performance.

Using the fitted parameters, the daily rate of target earthquakes (⁠|$m \ge 5.0$|⁠) is calculated for each cell within both the tessellations |${R}_1$| and |${R}_2$| for both the learning period (1990–2004) and testing period (2005–2023). The daily rate for each cell, which refers to the expected number of such events per day within the cell, is obtained by integrating the rate density function |$\lambda ( {t,\ m,\ x,\ y} )$| over the magnitude range |$[ {{m}_\mathrm{ T} = 5.0,\ {m}_\mathrm{ u} = 7.5} ]$| and the spatial cell region C, while excluding the integration over time. The daily rate |${\lambda }_{\textrm{daily}}( {m,x,y} )$| is thus given by: