**1. Introduction**

A typical aspect of the observed seismicity in the northern and central Apennines, and in the whole Italian region more generally, is the occurrence of earthquake sequences characterized by multiple, similarly large mainshocks. An example of this behavior is the quantitative model "Every Earthquake Precursory According to Scale" (EEPAS), applied by Rhoades and Evison [1,2,3]. According to their quantitative definition, introduced by Evison and Rhoades [4], swarms are seismic sequences constituted by at least three earthquakes whose magnitudes are linked to each other by empirical rules.

In this study we define as a multiplet a set of two or more earthquakes, with the following conditions: (a) the first event has a magnitude equal to or larger than a given threshold; (b) the others occur within a time difference and distance defined by the Gardner and Knopoff [5] criterion from each other; and (c) within a given magnitude range. This definition is different from that usually applied for common seismicity patterns such as foreshock–aftershock sequences and clusters (e.g., Gentili and Di Giovambattista [6]).

Building upon a previous paper (Console et al. [7]), in which we examined the aspect of multiple mainshocks in central Italy, in this study we aim at verifying if a synthetic catalog reproduces this kind of earthquake clustering. For this purpose, we apply a new version of the simulator algorithm, in which the role of stress transfer among elements

**Citation:** Console, R.; Vannoli, P.; Carluccio, R. Physics-Based Simulation of Sequences with Foreshocks, Aftershocks and Multiple Main Shocks in Italy. *Appl. Sci.* **2022**, *12*, 2062. https://doi.org/ 10.3390/app12042062

Academic Editor: Stefania Gentili

Received: 22 December 2021 Accepted: 8 February 2022 Published: 16 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> Istituto Nazionale di Geofisica e Vulcanologia, 00143 Rome, Italy; rodolfo.console@ingv.it (R.C.); roberto.carluccio@ingv.it (R.C.)

of an expanding rupture is enhanced. Moreover, we give also examples of other spacetime seismic features exhibited by synthetic catalogs, both in short- (days–months) and long-term (years–centuries), some of which were observed in real earthquake catalogs.

In Section 2 we present a brief description of the algorithm used for detecting multiple events in an earthquake catalog, based on the previously cited (a), (b), and (c) criteria. This algorithm is applied for providing a possible metric for comparing real observations with simulations.

Section 3 gives an outline of the seismotectonic model of our study area of northern and central Apennines, along with examples of recent and historical sequences of multiple mainshocks observed in this region.

In Section 4, after a short introduction of the new version of the simulator employed in this study, we show the results obtained applying this simulation code to the above mentioned seismotectonic model of the study area. Having tried three choices for the two main free parameters present in the algorithm, for a total of nine different combinations, we chose one of them by a criterion based on the analysis of the multiplets in the synthetic catalog of 100,000 years. Some features of this preferred simulated catalog are then compared in several ways with a real set of observations lasting only 370 years in the same seismogenic area.

Section 5 reports other results of spatio-temporal analysis of the same 100,000 years simulated catalog that appear to be consistent in reasonable way with real seismicity patterns not strictly related to our study area. In particular, we show that the use of simulators allows testing hypotheses of seismogenic models in a way that is not possible on the basis of real observations, due to lack of completeness and homogeneity of these observations in the long-term.

#### **2. The Algorithm for Identification of Multiple Events**

A special algorithm for the search of multiple events in a seismic catalog was created to use it as metrics in the comparison between the simulator results and the observed seismicity of the studied region. The computer code is "customer-built" and it was already introduced by Console et al. [7]. There is no specific definition of "sequence with multiple main shocks", nor any fixed magnitude values. We give here a brief description for a better understanding of its use. At its first level, the algorithm systematically analyzes time-ordered couples of events to check if they meet some constitutive conditions. Once matching couples are found, they are then used as elements for ordered noncyclic graph construction. These graphs can be 'traversed' to find in them the searched multiple events groups. In accordance with the above definitions, we developed a method based on four criteria for our comparisons among couples of events (Table 1):


The values associated to the criteria #1 and #2 are selected by the user and based on needs, expert judgment and/or knowledge of the instrumental–historical seismicity of the area. The relations for time (#3) and distance (#4) thresholds as a function of event magnitudes were derived from Gardner and Knopoff [5] empiric tables and as epicentre distance threshold. Notice that the criterion #3 does not provide for a choice, while the criterion #4 provides for the choice among three different ways of applying the formulas of Gardner and Knopoff [5]. In this study, we used the sum of the radii. Once launched, the algorithm parses the catalog through a cyclic, the three-step analysis procedure is repeated until the end of the file:


Table 1 summarizes the formulas, explaining the rationale by which they were used in this study, and lists our choices for the threshold values.

**Table 1.** Constitutive criteria for our clustering analysis algorithm. For each criterion, the second column shows rationale on which it is based, the third column contains formulas, and the last column contains choices used in this paper.


Even if the algorithm cannot be called "optimal" in principle, since it is based on an arbitrary choice among possible criteria, it is, however, quite effective, and its importance lays in the metrics it represents for comparison among seismic catalogs.

## **3. Seismotectonic Model**

The seismogenic model of the study area straddles northern and central Italy from the large flat area of the Po Plain (to the north) toward the northern flank of the Gran Sasso mountain range, the highest sector of the Apennines (to the south, Figure 1). The study area is wider than that previously studied (Console et al. [8]) through an old version of the simulator and the seismogenic sources now come from the latest 3.3.0 version of the DISS (DISS Working Group [9]).

Historical and instrumental seismicity in the study area is mainly distributed along the axis of the northern and central Apennines chain and, secondarily, in correspondence with its foothills, plains, and coastal areas (Rovida et al. [10]). The causative sources of the earthquakes of these two regions have different parameters and kinematics, as shown by focal mechanisms (Pondrelli et al. [11]), active stress indicators (Mariucci and Montone [12]), geological data (see DISS Working Group [9] and references therein), and active strain data (Devoti et al. [13]). As matter of fact, the GPS data show the crustal extension at a rate of about 3 mm/yr across the Apennines belt and the compression towards the Adriatic foreland (Devoti et al. [13]).

The active extension along the backbone of the Apennines is accommodated by normal faulting, which dominates along the hinge of the chain at shallow crustal seismogenic depth (blue polygons in Figure 1; e.g., Vannoli et al. [14]). The strongest extensional recent earthquakes in the study area occurred during the 2016-2017 central Italy seismic sequence

that struck central Apennines with multiple mainshocks (Table 2). The sequence initiated on 24 August 2016 with the *Mw* 6.2 Amatrice earthquake and was followed on 26 October 2016 by the *Mw* 6.1 Visso earthquake, about 25 km to the north. The largest event, the *Mw* 6.6 Norcia earthquake, occurred on 30 October 2016 and nucleated between the source regions of the two previous mainshocks (e.g., Michele et al. [15]; Rovida et al. [10]). Low-magnitude earthquakes of this sequence still occur today (http://terremoti.ingv.it/ accessed on 22 November 2021). This seismic sequence activated a circa 80 km long, NNW-SSE trending, low-angle multiple fault systems (IDs 127 and 128 in Figure 1). These fault systems exhibit complex ruptures and are the easternmost normal faults of the central Apennines, just west of where compressional activity prevails (e.g., Basili et al. [16]; Bonini et al. [17]; Di Bucci et al. [18]; DISS Working Group [9]).

**Figure 1.** Fault systems and earthquakes. Forty-three DISS (version 3.3.0) seismogenic fault systems are divided into 198 quadrilaterals that best approximate the DISS composite sources, and they are labeled with last three numbers of their DISS-IDs (DISS Working Group [9]). They are shown in accordance with their kinematics (extensional in blue, compressional in red, strike-slip in green), and have colored circles associated with their upper-left corners. Epicentres of the earthquakes from 1650 to 2020 A.D., with *Mw* ≥ 5.5, within a 5 km buffer from faults (dotted line) are shown by black circles. Main shocks of Table 2 are labeled in black (Rovida et al. [10]).

The active compression in the Adriatic foreland is mainly accommodated by thrust faulting (e.g., Vannoli et al. [19], Vannoli et al. [20]). Thrust faulting is widespread along the external fronts (red polygons in Figure 1) and propagates from the inner and coastal areas towards the offshore (to the east) and the Po Plain (to the north). Strongest recent compressional earthquakes of the study area occurred in Emilia during the 2012 sequence. This sequence began with the 20 May *Mw* 6.1 earthquake and was followed on 29 May 2012 by the *Mw* 5.9 earthquake; therefore, it is characterized by two similarly large mainshocks (see also Figure 1 in Console et al. [7]). The causative faults systems of the 2012 sequence are the external arcs of the most advanced and buried portions of the northern Apennines (IDs 103 and 51 in Figure 1; e.g., Vannoli et al. [20]).

Therefore, the seismogenic model of northern and central Apennines includes onshore and offshore seismogenic sources characterized by both extensional and compressive kinematics (DISS Working Group [9]). In addition, dextral strike-slip faulting is present in the southernmost study area, at the northern border of the Gran Sasso ridge (green polygons in Figure 1). Generally, the transverse structures are faults inherited from older

tectonic phases that cut the Adriatic foreland areas accommodating the segmentation of the thrust fronts and the outward propagation of the fold and thrust belts (e.g., Zampieri et al. [21]). Specifically, these strike-slip sources are high-angle, ENE–WSW-trending faults bounding the central Apennines thrust fronts and the southern part of the Apennines basal decollement. They are relatively deep (having 15–20 km of maximum depth), with shear zones that affect the Adriatic foreland (IDs 135 and 134). The western source (ID 135) is believed to be responsible for the seismic sequence that includes two relatively similar large mainshocks that occurred on 5 September 1950 ( *Mw* 5.7) and 8 August 1951 ( *Mw* 5.3; see Table 2).

In summary, the earthquake sequences characterized by at least two similarly large mainshocks are rather common in the study area, affect compressional, extensional, and strike-slip environments, and are very different from the sequences made up of a single large earthquake followed by aftershocks of decreasing magnitude. Figure 2 shows the epicentres of the CPTI15 catalog from 1650 to 2020, with *Mw* ≥ 5.0, and the colored lines connect the multiple main shocks events recognized by the algorithm described in the text and reported in column "Csum" of the Table 2. In the same Table the column "C1st" shows the results of the algorithm applying the criterion 4a.

**Table 2.** Largest sequences in real catalog with at least two main shocks of past 370 years (1650– 2020; magnitude and locality from CPTI15). Results of algorithm for detecting sequences with multiple main shocks in study area are shown in last two columns. Column "Csum" shows results of algorithm applying criterion 4c, while column "C1st" criterion 4a (Y: simulated; N: not simulated). Kin: Kinematics; N: normal; S: strike-slip; T: thrust; n.a.: not applicable; \* inferred (faults and kinematics responsible for historical earthquakes are inferred)



**Table 2.** *Cont.*

The seismogenic model upon which we applied the simulator code was derived from the Composite Seismogenic Sources (CSS) of DISS, version 3.3.0 (DISS Working Group [9]). The CSSs are parameterized crustal faults based on regional surface and subsurface geological data, and they are believed to be capable of producing *Mw* ≥ 5.5 earthquakes. We converted the 43 CSSs identified in the study area into 198 quadrilaterals specifically developed for this study, and this is consistent with all the geometrical and kinematics parameters supplied for the CSSs (Figure 1). The Table S1 in the Supplementary Material reports the list and the parameters of the 198 quadrilaterals recognized in the study area. Figure S4 shows a sketch of a quadrilateral fault segmen<sup>t</sup> and the description of its geometrical parameters.

**Figure 2.** Epicenters of CPTI15 catalog from 1650 to 2020, with *Mw* ≥ 5.0. Colored lines join mainshocks of the same sequence recognized by algorithm and described in text and in Table 1. The colors indicate the number of mainshocks for each sequence (the three sequences consisting of two mainshocks are shown in light green, and so on; see histogram in inset).

#### **4. Simulation of the Seismicity**

By means of a newly developed version of our physically based simulation code Console et al. [22] and references therein), we compiled synthetic earthquake catalogs lasting 100,000 years for events of magnitude ≥ 4.2 within the polygonal area depicted in Figure 1. In this version, no application is performed on the State and Rate formulation, but we adopted an enhanced role of the static Coulomb stress transfer between every ruptured element of the fault model and all the other elements in the surrounding faults. In this new version of the code, the magnitude distribution of the simulated catalog is controlled by two free parameters to be selected by the user (Console et al. [7,23]):


The seismogenic model adopted in the simulation algorithm is depicted in Figure 1, and the slip rates assumed for each fault segmen<sup>t</sup> are the highest values of the range reported by the DISS database (DISS Working Group [9]; Table S1 of the Supplementary Material).

We carried out a set of tests to investigate the effect that the two above described free parameters have on the magnitude distribution of the output catalogs, letting the S–R parameter assume the values 0.1, 0.2, and 0.3, and the A–R parameter the values 2, 5, and 10, respectively. The results of these tests are reported in Figure S1 of the Supplementary Material. Each 100,000 years catalog was divided in 270 groups of 370 years (with the purpose of simulating many instances of the real catalog), counting the number of multiplets contained in each of them. Table 3 reports the averages and the standard deviation for the 270 elements population. Then, for each of the nine cases, we carried out the same analysis on 50 randomized catalogs obtained from the 100,000 years original ones by shuffling the origin time of all earthquakes by a random permutation. Finally, the average and the standard deviation of the ratios between the total number of multiplets in the original catalogs and those obtained from the respective randomized catalogs was computed (Table 4).

On the basis of the above-mentioned tests, although the largest number of multiplets is provided by the couple of parameters 0.1 and 2, we chose the simulation obtained with the values 0.2 and 10 for the S–R and A–R parameters, respectively, which gives the largest ratio of multiplets. Figure 3 shows the results of this simulation with the 13,845 earthquakes having *Mw* ≥ 5.0, evidencing the fault segments where the number of simulated earthquakes is higher. Our simulation algorithm does not produce any seismic activity outside the borders of the faults considered in the seismogenic model.

A comparison of seismic features detected in the CPTI15 catalog, in the time interval 1650–2020, and the 100,000 years simulated catalog for the study area is shown in Table 5. For example, in this table we can compare the rate of earthquakes with *Mw* ≥ 5.0 in the simulated catalog (0.138/yr) with the corresponding rate of earthquakes with *Mw* ≥ 5.0 in the real catalog (0.573/yr). This circumstance is justified by the adoption of relatively high values of the S–R and A–R free parameters, which favor the growth of nucleated ruptures and accordingly produce a relatively large quantity of strong earthquakes. Moreover, we should take into account the fact that the source model adopted in our simulation does not include the numerous small sources, capable of producing only *Mw* ≤ 5.5 earthquakes.

Table 5 shows a comparison of seismic features detected in the CPTI15 catalog, in the time interval 1650–2020, and the 100,000 years simulated catalog for the study area. In Figure 4, we show the Magnitude–Frequency Distribution (MFD) of the simulated catalog, compared with that of the 1650–2020 CPTI catalog for events above the completeness threshold magnitude of 5.0.

Another metric for comparing our simulations with the real observations is given by the numbers of multiplets counted in the same time interval of 370 years (third line of Table 5) and the mean ratio between these numbers and the corresponding numbers calculated on a set of randomisations (last line of Table 5): these randomisations should effectively destroy the presence of clustering relation among the events. The obtained mean number could in our opinion represent the degree of "clustering" of the catalogs. The value of the ratio systematically greater than one supports the hypothesis that the occurrence of sequences containing multiple mainshocks is not just a casual circumstance. Even if this procedure can give different results changing internal criteria, these criteria are not changed while applying the procedure to the two catalogs to be compared. Figure 5 shows the distribution of the ratios between the number of multiplets identified in the 100,000 years synthetic catalog and 500 randomizations of the same catalog. The average ratio is 2.13 ± 0.24, which denotes a good agreemen<sup>t</sup> between the production of multiplets of the simulated catalog with respect to that of the observations (see also Table 5).

In Figure 4, we show the MFD of the simulated catalog, compared with that of the 1650–2020 CPTI catalog for events above the completeness threshold magnitude of 5.0. This figure shows that the MFD of the simulated catalog does not follow a straight line as expected according to the Gutenberg–Richter law, but exhibits a change in its slope in the magnitude range 5.7 ≤ *Mw* ≤ 7.0, where the b-value decreases dramatically. This circumstance is again due to the selection of the S–R and A–R free parameters, and the boostered role of the Coulomb stress transfer adopted in this particular study, which enhances the growth of nucleated ruptures, producing a sort of characteristic earthquake model.

**Figure 3.** Map of 13,845 simulated earthquakes with *Mw* ≥ 5.0, obtained from 100,000 years simulation. Point opacity is proportional to number of epicenters reported in output synthetic catalog in each cell of fault.

We also performed a comparison between the annual seismic moment rate released by the earthquakes of the simulated catalog and the observed ones. Adopting Hanks and Kanamori [24] magnitude–seismic moment conversion formula, we computed the total seismic moment released in the simulated catalog of 100,000 years. The sum is equal to 0.518 · 10<sup>22</sup> Nm, i.e., a seismic moment rate of 0.518 · 10<sup>17</sup> Nm/year. In a similar way, we computed the seismic moment of all earthquakes of M > 5.0 listed in the observational catalog from 1650 to 2020. A value of 1.35 · 10<sup>19</sup> Nm is acquired, implying a seismic moment rate of 0.365 · 10<sup>17</sup> Nm/year (Table 5). In conclusion, the seismic moment rate released by the simulated catalog is about 1.4 times larger than that of the observed seismicity. This could be explained by the uncertainties in the slip-rate values assumed in our seismogenic model.

**Figure 4.** Cumulative (yellow line) and density (blue dots) Magnitude–Frequency Distributions (MFD) of *Mw* ≥ 5.0 earthquakes of 100,000 years simulated (**left panel**) and observed (**right panel**; CPTI15 from 1650 to 2020 AD) catalogs. Straight dotted lines show best-fit of cumulative distributions.

We should also take into account the limited size of the earthquake catalog considered in the comparison of the observed seismic moment rate with that obtained from simulations. In fact, the duration of the 1650–2020 catalog (370 years) is shorter than the recurrence time on any of the fault segments reported in Table S1 and Figure 1. It is reasonable to hypothesize that this time window, upon which 22 events with *Mw* ≥ 6.0 have occurred, was characterized by a moderate seismic activity in our study area, without a significant contribution of large magnitude events. In contrast with that situation, in the 17th and 18th centuries large magnitude earthquakes occurred in central-southern Italy, outside our study area.

In the same way as we prepared Figure 2, we also plot in Figure 6 the epicentres of the 100,000 years simulated catalog with *Mw* ≥ 5.0. The comparison with Figure 2 shows that the simulated catalog is characterized by a scarce presence of sequences with a number of mainshocks larger than 2.

**Figure 5.** Histogram of ratios between number of multiplets identified in 100,000 years simulated catalog and 500 randomizations of same catalog.

**Figure 6.** Representation of sequences with multiple mainshocks in 100,000 years simulated catalog, with *Mw* ≥ 5.0. Colored lines show multiple mainshocks recognized by algorithm described in text, with colors indicating respective number of earthquakes for each of them (see histogram in log scale in the inset).

**Table 3.** Average number of multiplets in 100,000 years simulated catalogs in groups of 370 years.


**Table 4.** Ratio between total number of multiplets in original 100,000 years simulated catalogs and average of respective randomized catalogs.


**Table 5.** A comparison of seismic features detected in CPTI15 1650–2020 catalog and 100,000 years simulated catalog (S–R = 0.2 and A–R = 10) for study area.


#### **5. Long- and Short-Term Features of the Simulated Seismicity**

A detailed analysis of the simulated 100,000 years catalog allows the detection of interesting spatiotemporal features showing similarities with analog features existing in the observations.

The following stacking procedure was adopted to highlight if systematic and coherent time features occur before or after "strong" earthquakes: (1) We take into account earthquakes of the simulated catalog with a magnitude greater than 5.2; (2) for each of those events, occurring at time *ti*, a time interval around it (*ti* − Δ*t*, *ti* + Δ*t*) is considered and subgroups of events falling inside that interval and with an epicentral distance less than Δ*r* are added to a stacking list. Their occurrence times are stored as counted relative to *ti*, i.e, with times ranging from − Δ*t* to + Δ*t*; (3) Once the stacking list was filled with all subgroups times, the resulting ( − Δ*t*, Δ*t*) interval is divided into a proper number of bins and events occurrences for any of the bins are counted and reported in the scatter plots.

Long-term seismicity patterns before and after a mainshock are shown in Figure 7. This figure shows the stacked number of *Mw* ≥ 4.2 earthquakes that preceded and followed an *Mw* ≥ 5.2 earthquake within an epicentral distance of 20 km. Here, we may note an acceleration of seismic activity some centuries before a mainshock, a modest quiescence starting 50 years before the mainshock and a strong aftershock occurrence in the following five years. After this aftershock phase, a trend of long-term quiescence recovering in some centuries is noted. In Figure S2 of the Supplementary Material we report the same kind of plots for all the nine combinations of free parameters, showing the same trends, with minor variations.

As far as the short-term features are concerned, a clear foreshock and aftershock pattern of the duration of some weeks before and after a magnitude *Mw* ≥ 5.2 event is visible in the stacking plot of Figure 8a. With the same time scale, Figure 8b shows a clear trend of b-value decreasing before a mainshock of *Mw* ≥ 5.2 and recovering to the average value just after it. Note that the large scattering of the b-values has no real physical meaning, but it is simply due to the limited number of earthquakes on which the b-value is calculated. However, this scattering is much smaller just before and after the mainshocks, when the earthquakes rate is much larger. In Figure S3 of the Supplementary Material, we report the same kind of plots for all the nine combinations of free parameters, showing the same trends, with minor variations. This feature was observed in natural sequences as, for instance, Montuori et al. [25], Papadopoulos et al. [26], Gulia and Wiemer [27].

**Figure 7.** Stacked number of *Mw* ≥ 4.2 earthquakes that preceded and followed an *Mw* ≥ 5.2 earthquake within an epicentral distance of 20 km in 100,000 years simulated catalog.

**Figure 8.** (**a**) Stacked number of *Mw* ≥ 4.2 earthquakes that preceded and followed an *Mw* ≥ 5.2 earthquake within an epicentral distance of 50 km in 100,000 years simulated catalog, zooming on a time scale spanning only 0.1 years (36.5 days); (**b**) average b-value in time bins of 0.365 days before and after an earthquake of *Mw* ≥ 5.2 containing at least 10 events.
