What Is the Effect of Seismic Swarms on Short-Term Seismic Hazard and Gutenberg-Richter b-Value Temporal Variation? Examples from Central Italy, October–November 2023

: A seismic hazard can be quantified by using probabilities. Modern seismic forecasting models (e


Introduction
Seismic sequences are usually of two types.The first one is the classical "mainshockaftershock" sequence (MAs), where an initial strong event E1 gives birth to its own progeny of triggered events, generally occurring within the triggering area of the mainshock, from the first generation on.The starting event E1 is usually the largest of the sequence, but it may also happen that a stronger offspring E2 occurs.In such a case, either one simply acknowledges the possibility that an aftershock stronger than the mainshock can occur, or the initial shock E1 is re-labeled as foreshock together with all the events preceding E2, the latter becomes the new mainshock, and all its triggered events represent the new progeny ("foreshock-mainshock-aftershock" sequence type).
The second kind of earthquake sequence is instead named "seismic swarm" (SW), and is characterized by numerous series of small-to-moderate events occurring in a local area within a relatively short period of time.
The definition of both the two types of sequences presents a delicate point, about which there is no clear agreement in the literature.Regarding MAs, the main drawback is related to the fact that there is no general rule to identify an event as a mainshock.The labeling of events as "foreshocks", "mainshocks" and "aftershocks" is completely subjective and not always accepted by the scientific community, since no difference in their physical mechanisms has been proven (e.g., [1]).For example, the most used probabilistic model in statistical seismology, the Epidemic-Type Aftershock Sequence (ETAS) branching process [2,3], distinguishes only between triggered and background events.The latter represents the stationary component of the seismic process, which generates the family of offsprings just because is the first event in time (not the strongest at all).The delicate point in the definition of SWs relies instead on the fact that there is not a precise minimum number of small events for the sequence to be characterized as a swarm type, nor a reference maximum magnitude threshold the events should have.It is also important to stress that, in both cases, it is very difficult to understand into which of the two types an ongoing seismic sequence will evolve; the labeling of MAs or SWs is possible, just a posteriori, after the end of the sequence itself (not to mention the yet non-shared and non-rigid definition of the "end of the sequence" itself).
Another term of comparison between MAs and SWs concerns the seismic hazard in the short term.The need to deliver quasi-real-time earthquake forecasts is becoming more and more important nowadays to help communities prepare for possibly destructive seismic sequences.A leading method in this regard is the probabilistic approach, which allows us to quantify how likely the occurrence rate of events in a given space-time-magnitude domain will be.Operational Earthquake Forecasting (OEF [4]) systems have been developed in several countries worldwide (e.g., Italy [5]; New Zealand [6]; United States [7]; Israel [8]) to produce these short-term probabilities in real time, thus allowing us to highlight (and interpret) their evolution in space and time.More precisely, OEF systems regularly (e.g., daily) deliver the forecast of earthquakes over a specific spatiotemporal window of interest (e.g., one week in 0.1 • × 0.1 • cells) and above a given threshold of magnitude, based on some probabilistic models typically used in statistical seismology, such as ETAS.Evaluating the evolution of the delivered probabilistic forecasts in near-real time may give insights into the difference between MAs and SWs, as well as into how they affect the short-term probabilistic rate.In general, the models used in OEF to date are expected to better catch the spatiotemporal changes in seismicity in the case of MAs, since the occurrence of a strong initial event more easily influences the probabilistic rates by inducing a large increase.However, there is no condition/hypothesis behind the OEF models that would impede the short-term seismic hazard from also changing when SWs occur.In this paper, we contribute to this discussion by investigating how earthquake swarms affect the short-term probabilistic seismic rate in the case of three real seismic sequences in Central Italy.Save for the drawbacks previously discussed, we ascribe all these sequences to the swarm type.
To complement the analysis, we also investigate the effect of earthquake swarms on short-term seismic hazard in terms of the temporal evolution of the b-value parameter of the Gutenberg-Richter frequency magnitude distribution (FMD).The b-value controls the slope of the FMD in a semi-log plot: the larger this parameter, the smaller the proportion of strong shocks with respect to the little ones.Great interest has been shown in the literature about the possibility of interpreting the spatiotemporal variability of the b-value in terms of the physical process, both before the occurrence of a strong event and during an ongoing seismic sequence, thus considering the capability this parameter could have in acting as a proxy for the average stress condition in the crust.A heated debate has been sparked among seismologists on this matter.
On the one hand, anomalies and transient variations in the b-value have been observed both preceding strong shocks (e.g., [9,10]) and during earthquake sequences (e.g., [11,12]), these changes being explained by some evidence of correlation between the b-value and differential stress, material heterogeneity, or fluid diffusion [13][14][15][16].In light of these results, some authors proposed to analyze temporal changes in the b-value in terms of precursors for an upcoming large event: examples are the traffic light classification by [17] and the combination of ionospheric GPS-TEC data by [18].
On the other hand, high uncertainty is related to the estimation of the b-value, which is indeed prone to several potential sources of bias, such as catalog incompleteness, magnitude binning, and the size of the earthquake catalog [19][20][21][22].Obviously, a wrong estimation inevitably leads to apparent variations that do not have any real physical meaning.
In this paper, we investigate the temporal evolution of the b-value by considering the weighted likelihood method [23], which bases the estimation of this parameter on the full history of data (likelihood estimation), with positive weights depending on the time lag between the first event in the catalog and all the other events [24].Uncertainties are also determined by considering the Gaussian confidence interval [25].This method does not require subjective choices in the estimation procedure, thus allowing us to obtain robust b-value time series that can be properly evaluated in terms of "different effect of MAs and SWs on short-term seismic hazard".

The Seismic Catalog
To study the effect of earthquake swarms on short-term seismic hazards, we focused on three catalogs from Central Italy, the latter being one of the most active seismic areas in this country.More precisely, we considered the events inside circles with a 30 km radius, centered at the geographical coordinates of (1) the municipality of Sora (Frosinone, Lazio region), i.e., Longitude 13.617 • and Latitude 41.717 • ; (2) the municipality of Monte Cavallo (Macerata, Marche, Italy), i.e., Longitude 13.001 • and Latitude 42.992 • ; and (3) the municipality of Lucoli (L'Aquila, Abruzzo, Italy), i.e., Longitude 13.339 • and Latitude 42.291 • .The temporal interval considered was 1 January 2018-27 November 2023.We also set a maximum depth of 40 km, to overcome issues related to the likely unreliable hypocentral location of deeper events, and a minimum threshold magnitude ML of 1.0.The seismic map is given in Figure 1, along with time vs. magnitude plots.All these three areas experienced strong earthquakes in recent centuries, according to the CPTI15 historical seismic catalog [26].
In recent months, these three areas have experienced swarm-like seismic activity, characterized by a high number of events, all with relatively small magnitudes.More precisely, since 1 October 2023, the recorded events with ML 1+ are seventy-eight in the Sora catalog (two events with ML 2.5+, maximum magnitude ML 2.8), three hundred and sixty-three in the Monte Cavallo catalog (five events with ML 2.5+, maximum magnitude two events with ML 2.9) and fifty-nine in the Lucoli catalog (three events with ML 2.5+, maximum magnitude ML 3.7).From Figure 1b-d, it is possible to note an increase in the seismicity in the right part of the plots.Figure 1b shows a huge increase in the number of events, while in the two other zones, this increase is partially masked by the higher background seismicity level.
Since the events' magnitudes are expected to follow the Gutenberg-Richter frequency distribution [27], we estimated the completeness threshold for the three considered catalogs by means of the Lilliefors test [28,29] for magnitudes' exponentiality.The determination of the proper completeness magnitude M c is actually a delicate point to address, as several factors may produce artifacts biasing its estimate and consequently the estimate of the b-value.Several methods have been proposed in the literature to estimate M c , such as the Maximum Curvature (MaxC) or the 90% and 95% goodness of fit of the FMD to the exponential distribution (e.g., [30][31][32]).However, none of them is specific for testing the exponential property of the magnitudes, well acknowledged in statistical seismology.The canonical statistical test in this scope is the Lilliefors, which is indeed a modification of the classical nonparametric one-sample Kolmogorov-Smirnov, where data are tested for coming from an exponential distribution instead of the classical Gaussian one [28,33].Despite being less conservative than, for example, the MaxC method, the Lilliefors test allows us to obtain a statistically robust result.Because of the closeness of the three areas, all in Central Italy, and due to the comparable seismic activity, we also stress that we estimated a single completeness magnitude for all the cases, obtaining a value of 1.5.

The Italian Operational Earthquake Forecasting System
As anticipated in the Introduction, in order to study the evolution of these swarmtype sequences and to investigate potential scenarios after the temporal window of analysis, we followed two main methods.The first one is related to the Operational Earthquake Forecasting (OEF) Italian system [5], which produces real-time short-term earthquake forecasts in each 0.1° × 0.1° cell of a spatial grid covering the entire Italian territory, according to the standards of the Collaboratory for the Study of Earthquake Predictability (CSEP, [34]).At midnight every day and after the occurrence of any ML 3.5+ event, the OEF-Italy system delivers the weekly expected rate of events with ML 4+ and 5.5+, with Modified Mercalli Intensity (MMI) VI+, VII+, and VIII+.The forecast is probabilistic and based on the ensemble combination of three models typically used in statistical seismology: the Epidemic-Type Aftershock Sequence (ETAS) proposed by [35], the Epidemic-Type Earthquake Sequence (ETES) by [36], and the Short-Term Earthquake

The Italian Operational Earthquake Forecasting System
As anticipated in the Introduction, in order to study the evolution of these swarm-type sequences and to investigate potential scenarios after the temporal window of analysis, we followed two main methods.The first one is related to the Operational Earthquake Forecasting (OEF) Italian system [5], which produces real-time short-term earthquake forecasts in each 0.1 • × 0.1 • cell of a spatial grid covering the entire Italian territory, according to the standards of the Collaboratory for the Study of Earthquake Predictability (CSEP [34]).At midnight every day and after the occurrence of any ML 3.5+ event, the OEF-Italy system delivers the weekly expected rate of events with ML 4+ and 5.5+, with Modified Mercalli Intensity (MMI) VI+, VII+, and VIII+.The forecast is probabilistic and based on the ensemble combination of three models typically used in statistical seismology: the Epidemic-Type Aftershock Sequence (ETAS) proposed by [35], the Epidemic-Type Earthquake Sequence (ETES) by [36], and the Short-Term Earthquake Probability (STEP) by [37].Common to all these models is the Omori-Utsu law for the temporal decay of aftershocks [38,39] and the decreasing exponential Gutenberg-Richter law for the magnitudes.Both ETAS and ETES belong to the class of self-exciting Hawkes processes of branching type and are based on the rationale that seismic events cluster in space and time and evolve in successive generations independently of the others.Instead, the STEP model first recasts the forecasts in terms of strong ground-shaking probability; then, it combines a time-dependent model of earthquake occurrence, based on faults and historical data, with (i) a Reasenberg and Jones of clustering model [40], (ii) a sequence-specific model, and (iii) a spatially heterogeneous model.The combination of ETAS, ETES, and STEP is weighted according to the past performance of the three models considered.As input, OEF uses events with ML 2.5+ and is continuously updated according to the catalog recorded by the seismic surveillance system of the INGV.The OEF-Italy system is shown to deliver globally reliable forecasts [41,42], and is therefore the first tool we used in this study to investigate the evolution of the considered swarm-type seismic sequences.

The Gutenberg-Richter b-Value Estimation
The second method we followed is instead related to the well-known b-value parameter of the Gutenberg-Richter law, which, as anticipated in the Introduction, controls the proportion of larger shocks with respect to the smaller ones.In this paper, we analyzed the temporal variation of the b-value by means of the weighted likelihood method [23], which allowed us to estimate this parameter based on the full history of available data.The weights for the likelihood are defined as in [24]: each weight W i (i = 1, . .., N events) is positive and expressed as an exponential kernel depending on the time lag between the i-th event and the first event in the catalog: where k is a normalizing constant and α is a "forgetting parameter", which quantifies the amount of relevant past information.The parameter α is found through a maximum likelihood estimate (MLE) approach.The weighted likelihood estimation b of the b-value, properly corrected to account for magnitudes' binning (see also [43]), is then obtained as follows: where M i are the events' magnitudes, M c is the completeness magnitude, and ∆M is the magnitude binning (in our case, 0.1).The uncertainty was finally determined by considering the normal approximation, with the b-value standard error σ equal to b/ √ N [19,25,44].The difference between the weighted likelihood method and the classical [25] MLE approach is straightforward, with the first approach being a generalization of the latter.In the MLE case, each event has the same importance in the estimation.In the weighted likelihood case, the importance is related to the temporal distance of events used in the estimation: the larger the distance between the time of the event and the actual time (∆t i in Equation ( 1)), the smaller the importance in the estimation.The advantage of this method is that we do not need subjective choices for the building of the b-value time series: in the classical rolling-window approach, a fixed number of events is used to build the time series (usually 50, 100, or 200), and this number is a subjective decision of the researcher.Here, the parameter α, which rules the importance of past information in the estimation of the actual b-value, is estimated through an MLE approach, avoiding subjective choices on the number of events used.

OEF Probabilities
The OEF probabilities of MMI VI+, VII+, or VIII+ and ML 4+ or 5.5+ events, for the three areas under analysis, are given in Table 1.Specifically, we report the background rate computed during the last 5 years for each case, together with the maximum peaks of probability since October 2023, as well as the weekly probabilities delivered on the last date considered in our catalogs, i.e., 27 November 2023.The OEF-Italy graphical interfaces relative to the case ML 5.5+ are also given in Figure S1 of the Supplementary Material (the top, middle, and bottom panels are for the Sora, Monte Cavallo, and Lucoli catalogs, respectively).The results show that, for all the catalogs considered, there is a small increase from the background rates to the highest peaks observed during the swarm activity since October 2023.Precisely, from the zero increase obtained for Monte Cavallo in the cases of MMI VI+ and ML 5.5+, we observe a maximum increase of about four times in the MMI VII+ and MMI VIII+ rates for Lucoli.The latter is indeed the catalog for which the OEF rates are shown to be more influenced by the seismic activity, although very slightly.Interestingly, the swarm sequence in Lucoli is not the one containing the highest number of ML 2.5+ events among the catalogs considered here (we recall that only these events are accounted for, in computing the OEF rate).That is to say, one ML 3.7+ event (as in Lucoli) weighs more than a higher number of smaller events, provided they all have an ML of 2.5+, including two ML 2.9 (as in Monte Cavallo).Finally, Sora shows consistency for the five cases (all rates increase by about 1.5 times).
It is important to stress that, after the occurrence of the strongest events in the three considered areas during the entire operating period of OEF (from 2009), the OEF rates actually increased considerably.Regarding Sora, the maximum peaks were obtained during the 2009 L'Aquila sequence, which indeed also affected this area due to its proximity to the Abruzzo region, and during the 2013 sequence, the rates increased by factors from 20 to 35, reaching, for example, an MMI VIII+ rate of 2 × 10 −3 and an ML 4+ rate of 0.05.The same sequence induced an increase in the OEF rates relative to Lucoli from two to three orders of magnitude.On 6 April 2009 at 8:00 a.m., OEF released an 80% probability of this area having an ML 4+ event during the following week.A 90% probability was finally released for MMI VI+ and ML 4+ events after the Norcia Mw 6.5 earthquake in the area around Monte Cavallo; here, OEF rates increased by factors from 100 to 250 during the Central Italy sequence.Overall, comparing the OEF maxima probabilities delivered after the strongest events experienced since 2005 with those obtained during the current swarm activity, we observe an increase factor IF from the background that is about 20 times higher for Sora, 60 times higher for Lucoli, until 250 times higher for Monte Cavallo.The factors of increase for all the cases are given in Table S1 of the Supplementary Materials.
In conclusion, we can deduce that the occurrence of a swarm-type sequence does not significantly increase the OEF-Italy system probabilities, and the observed increases are smaller than one order of magnitude (i.e., less than 10 times) with respect to the background probabilities.

Temporal b-Value Variations
The results of the estimation of the temporal b-value variations from 1 January 2019 to 27 November 2023 are shown in Figure 2.These results show three different situations for the three catalogs.In the case of Sora (Figure 2a), the b-value shows no significant variation for the whole time period.This also reflects the result that this catalog shows, with respect to the other catalogs, the smallest OEF increase factor IF from the background when comparing the maxima probabilities delivered after the strongest events experienced since 2005 with those obtained during the current swarm activity.In conclusion, the b-value variations for all three catalogs do not show significant deviation from the background value during the seismic swarms of October and November 2023, i.e., estimated values considering the uncertainty are similar to the background value.Overall, the results obtained for the b-value are consistent with those Regarding the Monte Cavallo catalog, we have a larger number of events during the analyzed time period, and we can observe a fluctuation in the b-value.In the last month, it is possible to observe a decrease in the b-value; however, this decrease led to the b-value being near the background value (i.e., the value computed using all the events in the complete part of the catalog).Recalling that this parameter is estimated by means of the magnitudes (usually the local ML), the evidence just discussed agrees with the fact that Monte Cavallo shows, for the OEF ML4+ and ML5.5+ rates, an IF∼1 during the current swarm activity and experienced some of the highest IFs after the largest shock (see Table S1 of the Supplementary Materials).
For the Lucoli catalog, we observe moderate fluctuations of the b-value, and generally, these fluctuations are not significantly different from the background value.Also in this case, we can observe a small decrease in the b-value in the last part of the catalog; nevertheless, its decrease is not significant and is constrained within the uncertainty range.Recalling that Lucoli is indeed the catalog for which we observed a slightly higher influence of the swarm activity on the OEF rates, we could have expected to obtain a more evident fluctuation of the b-value in this case.However, we have to recall that Lucoli experienced the lowest number of events during the current swarm activity: 84% and 24% less than in Monte Cavallo and Sora, respectively.Furthermore, it experienced the highest magnitude, i.e., one ML 3.7, which is quite but not significantly higher than the single ML 2.8 and the pair of ML 2.9 that occurred in Sora and Monte Cavallo, respectively.Now, recalling that a sufficient statistic for the b-value estimate is indeed the mean magnitude [19,25,45], it seems reasonable to have obtained variations for this parameter, in the case of Lucoli, which are lower than but comparable to those for Monte Cavallo and slightly higher than for Sora.
In conclusion, the b-value variations for all three catalogs do not show significant deviation from the background value during the seismic swarms of October and November 2023, i.e., estimated values considering the uncertainty are similar to the background value.Overall, the results obtained for the b-value are consistent with those obtained for the OEF rates, in terms of their relative variations.Moreover, since the b-value time series of Monte Cavallo and Lucoli shows general fluctuations in recent years (Figure 2a,b), we underline the difficulty of identifying significant b-value variations in such a type of time series.Statistical b-value fluctuations, not related to any physical phenomena, mask the possible identification of increased/decreased b-value patterns if we have a limited number of events, as in our swarm seismicity case.Therefore, the study of b-value variations of aftershock sequences is clearly different from our situation, since the number of available events is dissimilar, at most hundreds in the swarm case vs. thousands for aftershock sequences.
Finally, we stress that for all the catalogs, the maximum observed magnitude in the considered time period is smaller than ML 4.5; therefore, we do not expect short-term aftershock incompleteness problems related to our 1.5 completeness threshold.We also verified this argument by observing that the sequential number of events shows a homogeneous pattern when plotted versus the magnitudes (heuristic technique [46]; see Figure S2 of the Supplementary Materials).We recall that short-term aftershock incompleteness could lead to a biased estimation of the b-value [22].For the sake of comparison, we finally performed the estimation of the b-values for the three catalogs by considering also the classical method of a rolling window with a fixed number of 100 events.As expected, results are consistent with those obtained above, but the estimations present a higher volatility (see Figure S3 of the Supplementary Materials).

Conclusions
Seismic swarms are characterized by numerous series of small-to-moderate events occurring in a local area within a relatively short period of time.Some of these swarms can lead to strong seismic events, as in the case of the 2009 L'Aquila earthquake in Central Italy (Mw 6.1 [47]).In this work, we analyze the effect of three seismic swarms on the probability released by the Italian Operational Earthquake Forecasting system (OEF) and on the temporal b-value variations.The study of b-value variations is important because its significant decrease could be related to an increase in the probability of a large impending earthquake.Our results highlight how the effect of swarms in OEF is very limited.No matter how many events occur, if they remain "small" (i.e., smaller than 4.0), the increase in the weekly probability of strong events in that area is less than one order of magnitude (i.e., 10 times).
The b-value temporal variations show, for two out of three time series, several fluctuations over recent years.These fluctuations include both increases and decreases in the b-value, making the interpretation of a deviation from the background value during the seismic swarms difficult, and consequently making it impossible to infer something from this parameter.
We can therefore conclude that, although the Central Italy zone is a well-monitored seismic area, with a completeness magnitude of ML 1.5 and with the OEF Italian system operating for more than 10 years, the occurrence of the three studied seismic swarms is not indicative of a significant variation in the probability of large impending events.Since both the OEF system and the studies of b-value variations are similar in other parts of the world, we can generalize our results.Seismic swarms containing about 50-200 events, and with a maximum observed magnitude smaller than 4.0, do not significantly increase the probability of strong impending events for OEF and b-value methodologies.
These findings should push seismologists to find other ways of studying seismic swarms if they want to extract information related to future large events, for example, investigating new spatio-temporal and magnitude statistical relations or using machine or deep learning techniques (e.g., [48,49]).Such techniques adopt optimized methods to reveal the temporal seismic structure, as they are truly trained from data.

Supplementary Materials:
The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/geosciences14020049/s1, Figure S1: OEF probabilities of ML 5.5+ events: top, middle, and bottom panels for Sora, Monte Cavallo, and Lucoli catalogs, respectively; Figure S2: Magnitudes versus sequential number of events: top, middle and bottom panels respectively for Sora, Monte Cavallo and Lucoli catalogs; Figure S3: b-value time series (thick lines), with the associated uncertainty (+/− σ, thin lines), estimated by means of a rolling window with a fixed number of 100 events: panels (a-c) respectively for Sora, Monte Cavallo and Lucoli catalogs; Table S1: Factor of increase in the OEF probabilities for the three considered catalogs, from the background to both the maximum probability during the swarm current activity and the maximum probability after the occurrence of the corresponding past strongest events.

Figure 1 .
Figure 1.Panel (a): epicenters of the events considered in this work; black for the Sora catalog, blue for the Monte Cavallo catalog, and red for the Lucoli catalog.Panels (b-d): time vs. magnitude plot for the events in the Sora, Monte Cavallo, and Lucoli catalogs, respectively.

Figure 1 .
Figure 1.Panel (a): epicenters of the events considered in this work; black for the Sora catalog, blue for the Monte Cavallo catalog, and red for the Lucoli catalog.Panels (b-d): time vs. magnitude plot for the events in the Sora, Monte Cavallo, and Lucoli catalogs, respectively.

Geosciences 2024 , 11 Figure 2
Figure 2. b-value time series (thick lines), with the associated uncertainty (+/− σ, thin lines), for the Sora catalog (a), the Monte Cavallo catalog (b), and the Lucoli catalog (c).Horizontal lines represent the b-value computed for all catalogs (considered the background value).

Figure 2 .
Figure 2. b-value time series (thick lines), with the associated uncertainty (+/− σ, thin lines), for the Sora catalog (a), the Monte Cavallo catalog (b), and the Lucoli catalog (c).Horizontal lines represent the b-value computed for all catalogs (considered the background value).

Table 1 .
OEF probabilities for the three considered catalogs.