Sequence-based hazard analysis for Italy considering a grid seismic source model

Earthquakes are usually clustered in both time and space and, within each cluster, the event of highest magnitude is conventionally identified as the mainshock, while the foreshocks and the aftershocks are the events that occur before and after it, respectively. Mainshocks are the earthquakes considered in the classical formulation of the probabilistic seismic hazard analysis (PSHA), where the contribution of foreshocks and aftershocks is usually neglected. In fact, it has been shown that it is possible to rigorously, within the hypotheses of the model, account for the effect of mainshock-aftershocks sequences by means of the sequence-based PSHA (i.e., SPSHA). SPSHA extends the usability of the homogeneous Poisson process, adopted for mainshocks within PSHA, to also describe the occurrence of clusters maintaining the same input data of PSHA; i.e., the seismic rates derived by a declustered catalog. The aftershocks’ occurrences are accounted for by means of conditional non-homogeneous Poisson processes based on the modified Omori law. The seismic source model for Italy has been recently investigated, and the objective of the study herein presented is to include and evaluate the effect of aftershocks, by means of SPSHA, based on a new grid model. In the paper, the results of PSHA and SPSHA are compared, considering the spectral and return periods that are of typical interest for earthquake engineering. Finally, a comparison with the SPSHA map based on a well-established source model for Italy is also provided.


Introduction
Probabilistic seismic hazard analysis [PSHA, Cornell, 1968;Mc Guire, 2004], in its classical formulation, describes the earthquakes' occurrence through a homogeneous Poisson process (HPP). The HPP is completely characterized by the (time-invariant) rate, that is, the mean number of events in the time unit (usually one year). The main result of PSHA is the rate of earthquakes causing exceedance of an arbitrary threshold of one ground motion intensity measure ( ) at a site of interest. This allows to compute the probability of exceedance of any ground motion intensity measure value in any time interval. ANNALS OF GEOPHYSICS, 64, 2, SE214, 2021; doi:10.4401/ag-8586 OPEN ACCESS In fact, while earthquakes are typically clustered in time and space, PSHA only models the occurrence of mainshocks, which are usually defined as the largest magnitude event in a cluster, and the effect of foreshocks and aftershocks, that is, the earthquakes of the cluster occurring before and after the mainshock, respectively, is neglected. This is required to meet the high desirable assumption of HPP occurrence of earthquakes and also because models involved in PSHA are based on seismic catalogs that usually lack enough data to completely characterize the occurrence of earthquakes during seismic sequences. Thus, for PSHA purposes, seismic catalogs go through a preliminary procedure to remove foreshocks and aftershocks, that is, declustering [e.g., Gardner and Knopoff, 1974].
However, recognizing that clusters, as a whole, occur in the same number of mainshocks alone, Iervolino et al. [2014] developed the so-called sequence-based PSHA or SPSHA that allows including the effect of both mainshocks and aftershocks (neglecting foreshocks) in hazard analysis. The main advantage of this approach is that it retains the HPP hypothesis for the occurrence of mainshock-aftershocks sequences, while aftershocks are accounted for via conditional non-homogeneous Poisson processes (NHPP). SPSHA avoids the issues of non-declustered catalogues such as completeness with respect to aftershocks. In fact, models to account for aftershocks in hazard analysis based on non-declustered catalogs have also been proposed [e.g., Marzocchi and Taroni, 2014]; however, considering the rates from non-declustered catalogs as those of a HPP presents some issues [see Iervolino, 2019, for a discussion].
The current official [i.e., adopted by the Italian building code; C.S.LL.PP., 2018] seismic hazard study of Italy [Stucchi et al., 2011] has been computed via PSHA referring to a seismic source model developed some years ago [Meletti et al., 2008]. The effect of aftershocks, evaluated by means of SPSHA, on that source model has been discussed by Iervolino et al. [2018] showing that, for the studied spectral periods and return periods (the same analyzed in the following), hazard increase due to aftershocks is up to 30% and it is, generally, more significant within or around areas exposed to comparatively higher hazard according to classical PSHA. Recently, a set of new source models (eleven working groups developed ninety-four source models) has been developed for Italy [Meletti et al., 2021]. From these models, a weighted average grid-seismicity model was derived. The objective of the study herein presented is to develop SPSHA using the cited grid-seismicity model, so as to evaluate the effect of aftershocks on the hazard in terms of the ground motion corresponding to exceedance return periods of typical interest for structural engineering. To this aim, in the remainder paper, after presenting the average seismicity model, the basics of PSHA and SPSHA are briefly recalled. Then, selecting two IMs corresponding to two spectral pseudo-accelerations (hereafter identified as spectral accelerations), PSHA and SPSHA are carried out for Italy considering a recent ground motion prediction equation (GMPE). Results are presented in terms of countrywide maps of acceleration values with fixed exceedance return period at each site. The comparison of the corresponding maps, derived by PSHA and SPSHA, allows the quantification of the hazard variations (increments) of SPSHA with respect to PSHA. Then, two specific sites, representative of low and high seismic hazard in the country, are considered. For these sites, a range of spectral periods is investigated comparing uniform hazard spectra (UHS) obtained by PSHA and SPSHA. Moreover, hazard variations are compared for a wide range of return periods. A discussion about the comparison between the results of this work and those based on the source model considered by the official hazard map precedes some final remarks that close the paper.

Grid (average) Italian seismic source model
Starting in 2015, a large community of researchers, led by the Italian National Institute of Geophysics and Volcanology (INGV), developed a set of ninety-four seismic source models for Italy. These models have been weighted according to statistical performances in describing the earthquake occurrences and taking into account an experts' elicitation process [Meletti et al., 2021]. Moreover, the models were combined via a logic tree approach, with a set of selected GMPEs [Lanzano et al., 2020], to obtain about six-hundreds branches.
To facilitate the reproducibility of the results without implementing the whole logic tree, a weighted average grid-seismicity model was also developed (see Data and resources section). Such a model, adopted in this study, is represented by a grid of more than eleven thousand point-sources covering and surrounding the country. The model describes the mainshocks occurrence via the so-called annual activity rates (i.e., average number of events per year per magnitude interval); according to them, first, the minimum and maximum moment magnitude ( ) of the generated mainshocks, on each source, are assigned, that is, , and , , respectively (the subscript denotes terms referring to mainshocks for reasons that will be clear in the next section). Then, assuming bins of magnitude Eugenio Chioccarelli et al.
of width equal to 0.1 and values between , and , , the rate of mainshocks with magnitude belonging to the generic i-th bin centered on = is the activity rate, , , ∀ = 1,2, ... . The first bin of all the sources of the grid is centered at 4.5 so that the minimum magnitude of earthquakes is equal to 4.45, while , varies over the grid. Figure 1(a) shows the mean annual number of earthquakes of magnitude larger than 4.45 that, for each pointlike source, is the sum of all the activity rates. In the same picture, the black dotted lines identify the two portions of the territory (norther Italy and Calabria) in which the maximum magnitude of earthquakes is 8.35.
In all the remaining areas, maximum magnitude is 9.05. The grid model is completed by a distribution of rupture mechanisms associated to each point-source.
With the aim of providing a general overview of the activity rates, they are given in Figure 1(b) for all the point sources. The figure allows showing the general trend of the rates: among all the sources, , ₌₄ . ₅ ranges between 4.95E-07 and 8.87E-04 events per year and , ₌₆ . ₅ ranges between 2.60E-09 and 1.55E-05 events per year. Although the discussion of the source model is out of the scope of this paper, it should be noted that maximum values of magnitude are increased with respect to previous seismological source models for Italy; this is also an effect of the tapered magnitude models adopted in some branches of the logic tree. However, the rates associated to the largest magnitudes (and their contribution to the hazard) are small, if not negligible. To give an example, , ₌₈ . ₃, that corresponds to the last magnitude bin for some of the sources, ranges between 7.87E-22 and 1.36E-10 whereas , ₌₉ (when larger than zero) varies between 1.03E-48 and 3.82E-44. Moreover, as discussed in the following, in the analyses carried out in this paper, the earthquake rates with magnitude larger than 6.9 are associated to the 6.9 magnitude bin in order to not extrapolate the ground motion prediction equation.

Sequence-based probabilistic seismic hazard analysis
For a chosen intensity measure threshold , PSHA allows computing the rate of mainshocks producing the exceedance of at the site, , . Such a rate, which characterizes the HPP describing the occurrence of the earthquakes causing > at the site, with reference to the grid source model just described, can be computed according to Equation (1): 3 SPSHA for Italy using a recent source model

Eugenio Chioccarelli et al.
In the equation, is the number of point-like seismic sources affecting the hazard of the site (e.g., within a certain distance from the site), , , is the rate of mainshocks from source = 1,2,... within the magnitude bin centered on , [ > | , ] is the conditional probability that the mainshock, with magnitude and generated by the point-like source with distance from the site equal to , causes the exceedance of ; i.e., provided by a GMPE.
As demonstrated in Iervolino et al. [2014], the average number of seismic sequences (mainshocks and following aftershocks) that in one year cause at least one exceedance of at the site, that is , can be computed via SPSHA, modelling aftershocks' occurrence via a conditional NHPP, following Yeo and Cornell [2009]. In fact, can be computed via Equation (2): In the equation, is the probability that the maximum intensity among the aftershocks, ∪ , following the mainshock of features = , = , does not exceed . Thus, the product of the probabilities in the curly brackets is the probability that is exceeded neither in mainshocks nor in subsequent aftershocks. [ ∪ ≤ | , ] is provided by Equation (3), in which the subscript ( ) denotes terms referring to aftershocks: ( 3) where [ > | , ] is the conditional probability, provided by a GMPE, that is exceeded given an aftershock of magnitude = and source-to-site distance = . The term , | , is the joint probability density function (PDF) of magnitude and distance of aftershocks, conditional on the magnitude and distance of the mainshock.
Indeed, assuming the conditional independence of the and random variables, the latter distribution can be written as , | , = | · | , , where | is the conditional PDF of aftershocks magnitude (i.e., derived from the model of Gutenberg and Richter, 1944), and | , is the conditional PDF of the distance of the site to the aftershocks. The aftershocks magnitude is bounded by a minimum magnitude, ,, and the mainshock magnitude; i.e., ∈ ,, . Given the location of the site for which hazard is being computed, the aftershocks' location depends on mainshock's magnitude and location (although distance is indicated in the equations; see the next section for details) and thus the aftershocks source-to-site distance, ∈ ,, ,, depends on magnitude and location of the generating mainshock. Finally, , , , are the parameters of the modified Omori law which probabilistically defines the process of occurrence of aftershocks after the mainshock [see Yeo and Cornell, 2009], and is the (assumed) duration of the aftershocks sequence.

Hazard analysis input models
The source model described in Section 2 is adopted for both PSHA and SPSHA. The adopted GMPE for modelling the propagation of both mainshocks and aftershocks is the one proposed by Bindi et al. [2011], which is fitted on Italian records and is one of the best performing models for Italy [Lanzano et al., 2020]. The GMPE refers to Joyner-Boore distance [ , Joyner and Boore, 1981] as metric of source-to-site distance. Assuming that each point-like source is representative of the epicenter of the earthquakes generated by that source, the R jb distance is obtained converting the epicentral distance according to Montaldo et al. [2005]. (Such a conversion is applied for magnitudes up to 6.9 in accordance to what discussed in the following).
The chosen GMPE is defined within the 4.0 -6.9 magnitude and 0 -200 km distance ranges. Therefore, to avoid extrapolations, earthquakes characterized by magnitude and source-to-site distance outside these ranges are neglected. However, as pertaining to magnitude, to still consider the contribution of high-magnitude earthquakes, in each point source, the activity rates associated to magnitudes higher than 6.9 are added to the rate of magnitude 6.9 (i.e., concentrated at 6.9).
To provide a quantitative measure of the modification introduced on the magnitude distribution of generated mainshocks, the values of the activity rate centered at magnitude 6.9, , ₌₆ . ₉, were computed. In the original model, considering the whole grid of point sources of Figure 1 concentrating rates of higher magnitudes in this bin, , ₌₆ . ₉ ranges between 2.80E-09 and 2.85E-05. Therefore, this correction is considered viable. In all the analyses discussed in the following, class A soil (i.e., rock), according to the classification of the adopted GMPE, was assumed.
SPSHA also requires models for describing the aftershocks' occurrence given a mainshock of known magnitude and location. More specifically, the parameters describing the aftershocks occurrence in Equation (3) are those of Lolli and Gasperini [2003]: = -1.66, = 0.96, = -0.03 (in days), = 0.93, which refer to generic (i.e., typical) Italian sequences. All the mainshocks may generate aftershocks for which minimum magnitude is assumed equal to 4.0 (that is, , < ,).
As pertaining to the location of aftershocks, the model of Utsu [1970] is adopted: it assumes that aftershocks may occur with equal probability in a circular area centered on the mainshock location (it was verified in Iervolino et al., [2014], that a particular choice of the shape of the aftershock source does not significanlty affects the results). The size of such an area, , depends on the magnitude of the mainshock ( = ) via Equation (4), in squared kilometers: = 10 ⁻⁴ . ¹.
The duration of the aftershocks sequence, that is in Equation (3), was considered equal to 90 days according to Iervolino et al. [2018], where it has been proven that extending longer this period does not have a significant impact on the results.

Analysis and results
All the analyses were performed with the REASSESS V2.0 software [Chioccarelli et al., 2019]. Two types of results are presented in this section. First, Section 5.1 discusses the hazard maps for four return periods ( ) of particular structural engineering interest (50, 475, 975 and 2475 years) and two spectral intensity measures that are the peak ground acceleration, , and the spectral acceleration at 1 second vibration period, (1s). The maps are derived through both PSHA and SPSHA and are compared to evaluate, at a large scale, the hazard increment due to aftershocks. Then, two sites are considered in Section 5.2, for a more detailed analysis. In fact, UHS from PSHA and SPSHA are evaluated for spectral accelerations in a range of vibration periods ( ) up to four seconds. Finally, the hazard increments are also quantified for and (1s), considering return periods up to ten thousand years.

Hazard maps
PSHA and SPSHA were performed on a regular grid of about ten thousand inland sites. For each site, and (1s) values corresponding to the four selected return periods were computed. Results are shown as maps in Figure   2, that refers to , and Figure 3, that refers to (1s). Each figure is made of eight panels divided into two rows: the first one shows the results of PSHA, while the second row shows the results accounting for the aftershocks' effect on the exceedance rate, via SPSHA. Moreover, each column identifies a return period, i.e., 50, 475, 975 or 2475 years (yr), from left to right. In both figures, the sites of L'Aquila and Milan, considered in the following section, are also identified (black triangles). Figure 3 show that, for each return period, hazard increment due to aftershocks is comparatively larger in sites for which the mainshock hazard is large: this is particularly clear comparing map (b) with map (f), that is = 475yr, and map (c) with map (g), that is = 975yr. This can be intuitively motivated because sites subjected to stronger or more frequent mainshocks are those in which significant seismic sequences are more likely to be observed.

Figure 2 and
Hazard increment over the country due to aftershocks is quantitively represented in Figure 4: in each map, the percentage increments computed as the difference between the results of SPSHA and PSHA, divided by results of PSHA, are presented. Maps from (a) to (d) are the increases in terms of for the four considered return periods; maps from (e) to (h) are the increases in terms of (1s).
Referring to the maps for , the following considerations can be made: (i) all the increases are larger than 7% and in the most of the sites the increases are between 10% and 15%; (ii) the number of sites in which increases are lower than 10% decreases with the increasing return period; (iii) the increases between 15% and 20% are mostly observed along the southern Apennines for = 475yr and = 975yr and in the south-eastern Sicily (i.e., an area also interested by volcanic activity) for ≥ 475yr. Hazard increments for (1s) are more homogeneous (in the adopted color scale): they are between 4% and 10% for all sites and considered return periods. However, as shown in the next section, increments for (1s) maintain the dependency on .

Site-specific analysis
The comparison between PSHA and SPSHA results for specific sites is reported in this section. The considered sites are L'Aquila (central Italy, 13.42°E and 42.34°N) and Milan (northern Italy, 9.12°E and 45.46°N). The sites are selected because they are representative of high and low hazard in the country, respectively (see the maps shown in the previous section).
First, the site of L'Aquila is considered. Figure 5(a) identifies its geographical location. In Figure 5(b) the UHS for the four considered return periods are shown: black lines are results of SPSHA while grey lines are for PSHA. According to the GMPE, the range of vibration periods is 0-4s. Figure 5(c) reports the hazard increase, as already defined for Figure 4, for the four return periods considered before and varying the vibration period. As shown, for vibration periods lower than 0.4s, the increases are larger than 8% for any return period. More specifically the largest increase for each spectrum is about 19%, 20%, 18%, and 16% for equal to 50, 475, 975 and 2475yr, respectively.
For vibration periods larger than 1.0s, the increases are relatively constant, ranging between 5% and 8% for the investigated return periods.
It is interesting to extend the analysis to other return periods. To this aim, Figure 5(d) represents the hazard increases as a function of the return period ranging from zero to 10000 years. The figure shows that the trend for is non-monotonic: increments are between 9% and 15% with a maximum corresponding to about 220 years.
On the other hand, increments for (1s) monotonically decrease with the increasing return period, ranging between 6% and 8% for the return periods of main earthquake engineering interest.
The site of Milan is considered in Figure 6. Similarly to the previous case, Figure 6(a) shows the geographical location of the site; UHS are reported in Figure 6(b) for the four return periods of interest while corresponding hazard increases are reported in Figure 6(c). Hazard increments are nonmonotonic for spectral periods lower than 0.5s. For = 50yr, the maximum increase is about 12% and occurs at vibration period equal to 0.2s. For = 475yr, the maximum increase is about 13% at the vibration period equal to 0.1s. Finally, for return periods equal to 975 and 2475 years, maximum values correspond to 0.07s spectral period and are about 14%. Similar to L'Aquila, for spectral periods larger than 1.0s, the hazard increases are relatively constant and within 4% and 7% for all the considered return periods.

7
SPSHA for Italy using a recent source model This is because spectral ordinates associated to large spectral periods are generally affected by larger magnitudes than those affecting low spectral periods [e.g., Iervolino et al., 2011]. Thus, the contribution of aftershocks, which have magnitude lower than the corresponding mainshock in the framework of SPSHA, decreases as the spectral period increases. Figure 6(d) shows the hazard increases as a function of the return periods. For larger than 50 years, the maximum hazard increase for is about 10% for equal to about 1100 years, whereas for (1s) the maximum is about 7% for equal to 100 years. Eugenio Chioccarelli et al.

Aftershock's contribution according to MPS04
In this section a comparison with results of Iervolino et al. [2018] is reported. This latter study discussed SPSHA for the source model adopted in the current seismic hazard model for structural design in Italy, which is referred to as MPS04 [see Stucchi et al., 2011, for details]. The MPS04 source model was characterized by areal seismic source zones (shown in Figure 7) with prevalent rupture mechanism associated to them. Morevoer, MPS04 featured a logic tree with sixteen branches; in Iervolino et al. [2018] only one of them was considered. It was based on activity rates, computed from a version of the seismic catalog older than that at the basis of the rates discussed in 9 SPSHA for Italy using a recent source model Section 2, and the considered GMPE was that of Ambraseys et al. [1996], which refers to surface-waves magnitude 1 .
Such differences suggest that a direct comparison between hazard from the previous and current source model is challenging, as discussed in Meletti et al. [2021]; however, it is possible to compare hazard increments due to aftershocks, because in the two cases, increments are self-consistent.
Differences between the hazard increases computed in Iervolino et al. [2018] and those of Figure 4 are reported in Figure 7 (for a grid of about four thousand sites adopted in the cited paper); in the figure, areal source zones for MPS04 are also reported. Three main conclusions can be drawn from the comparison: (1) in the most of cases, differences are within ∓5%; (2) in some cases, they are between 5% and 10% (e.g., the south-eastern Sicily and the southern Apennines) or lower than -5% (e.g., in the northwestern region and only in the case of ); (3) as expected, seismogenic zones have significant influence on the geographical distribution of the differences (e.g., in very few sites, at the boundary of seismogenic zones and only for , differences are between 10% and 15%) 2 .

Conclusions
PSHA, in its standard version, accounts for mainshocks, neglecting the contribution of the other events of seismic clusters; i.e., foreshocks and aftershocks. The, recently developed, SPSHA allows to account for the effects of also aftershocks on the exceedance rate of a ground motion intensity measure at a site of interest, avoiding the issues of using a non-declustered seismic catalog. Because a new set of seismic source models has been recently developed for Italy, from which a weighted average grid-source model has been derived, this study assessed the hazard increases due to the aftershocks' when SPSHA is employed.
1 Ambraseys et al. [1996] provides the largest horizontal acceleration response ordinates whereas the model adopted in the previous sections provides the geometrical mean of horizontal components. Thus, the hazard maps shown in the two papers cannot be directly compared and a conversion model needs to be applied [e.g., Beyer and Bommer, 2006]. 2 Sardinia was not considered in MPS04. Nationwide hazard maps were computed, on rock site conditions, via PSHA and SPSHA referring to and (1s), and considering four return periods between 50 and 2475 years. Hazard increment was quantified showing that it is dependent on the site, the exceedance return period and the spectral period of interest. The increment in terms of is between 7% and 15%, with the highest value recorded along the southern Apennines for = 475yr and = 975yr, and in the eastern Sicily for ≥ 475yr. The increases in terms of (1s) are generally lower, ranging between 4% and 10%, for all sites and considered return periods.
Two specific sites, L'Aquila and Milan, were analyzed in larger details. The analysis of single sites allowed deepening the trends of hazard increments for other vibration periods and return periods with respect to those discussed at national scale. Comparing the uniform hazard spectra from PSHA and SPSHA, it was shown that the hazard increases are larger for vibration periods between zero and about 0.5s and tend to be constant for vibration periods larger than 1.0s. The effect of return period was also discussed, showing hazard increment up to 10000 years. Results are different for the two considered spectral periods: the hazard increases for are larger than those for (1s) and, especially for the site of L'Aquila, show significant variations with .
Finally, a comparison of SPSHA-related hazard increments in the case of the previous source model of the country, currently used for seismic design and based on areal seismic source zones, was discussed. Comparable effect of aftershocks on hazard, with respect to the average source model considered herein, was found.

Data and resources.
The weighted average grid-seismicity model adopted in this paper is available at http://wpage.unina.it/iuniervo/papers/MPS19_weighted_average_model.xlsx. All the other data used in this study come from the listed references.