Earthquake Rupture Forecasts for the MPS19 Seismic Hazard Model of Italy

In recent years, new approaches for developing earthquake rupture forecasts (ERFs) have been proposed to be used as an input for probabilistic seismic hazard assessment (PSHA). Zone- based approaches with seismicity rates derived from earthquake catalogs are commonly used in many countries as the standard for national seismic hazard models. In Italy, a single zone- based ERF is currently the basis for the official seismic hazard model. In this contribution, we present eleven new ERFs, including five zone-based, two smoothed seismicity-based, two fault- based, and two geodetic-based, used for a new PSH model in Italy. The ERFs were tested against observed seismicity and were subject to an elicitation procedure by a panel of PSHA experts to verify the scientific robustness and consistency of the forecasts with respect to the observations. Tests and elicitation were finalized to weight the ERFs. The results show a good response to the new inputs to observed seismicity in the last few centuries. The entire approach was a first attempt to build a community-based set of ERFs for an Italian PSHA model. The project involved a large number of seismic hazard practitioners, with their knowledge and experience, and the development of different models to capture and explore a large range of epistemic uncertainties in building ERFs, and represents an important step forward for the new national seismic hazard model.


Introduction
Probabilistic seismic hazard analysis (PSHA) is the standard approach for quantifying the likelihood that ground shaking will exceed various levels during a given time window. The most important consequences of national-scale applications of PSHA include building code updates and planning actions for seismic risk mitigation and emergency response. While the common PSHA methodology remains based on classical approaches [e.g. Cornell 1968, McGuire 2005, significant progress has been made in refining its implementation through the improved knowledge of its two main inputs: the earthquake rupture forecast (ERF) and ground motion model (GMM). Advances in active fault studies and improved earthquake catalogs may help to refine ERF characterization, while the parameterization of new GMMs is constantly enhancing the estimations of the expected ground motion and its variability due to source and ground characteristics.
In the early 2000s, the Istituto Nazionale di Geofisica e Vulcanologia (INGV) produced the current official ERF for Italy using a time-independent (long-term) area source input [Stucchi et al., 2011].
In 2015, the Seismic Hazard Center (Centro di Pericolosità Sismica, CPS) of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) was commissioned by the Italian Civil Protection Department to engage and coordinate

Structure of an ERF input: rates and parameters
An ERF returns the long-term rate of all earthquakes throughout the region above a specified threshold; in MPS19, this value was fixed as Mw ≥ 4.5, based on the observation that earthquakes with a magnitude less than Mw 4.5 are generally unlikely to cause significant damage in Italy. In general, there are many possible approaches to build an ERF, and the overall framework in which all the ERFs are constructed is designed to accommodate many alternatives and, possibly, their future updates or modifications. An ERF is composed of two levels of information: a) long-term rate of earthquake occurrences and b) geometrical parameterizations of the seismogenic sources responsible for the earthquake occurrence.
The long-term rate of earthquake occurrences of the ERF is given as an incremental magnitude-frequency distribution. The ERFs for the MPS19 use a discrete incremental magnitude-frequency distribution defined by the minimum value of magnitude, bin width, and the annual number of events for different bins. More details on the evaluation of the long-term rates of occurrence for the ERFs are provided in the next sections. All the working groups involved in the parametrization of the ERFs evaluated the rate of occurrences in the area defined in Figure   1b, as the area covered by the Italian earthquake catalog, the Catalogo Parametrico dei Terremoti Italiani, CPTI15 [Rovida et al., 2016;2020].
To describe the seismogenic source geometries, we adopted the set of possible typologies provided by the OpenQuake Engine [Pagani et al., 2014]. We thus identified three main source types into which the ERFs for the Italian PSHA model could be grouped. 1) Area sources: this is a very common source type, used in many national and regional PSHA models, in which the source geometry is defined by a closed polygon. 2) Point sources: this is the source type commonly adopted to model distributed seismicity. Note that an area source is actually a container of point sources. 3) Fault sources: this is a less often used source geometry since it implies that a sufficient number of faults have been mapped and characterized in the model area. Prospective earthquake ruptures are constrained on the fault planes. Note that ERFs of this type use point sources to describe off-fault earthquake occurrences.
Along with the source typology, ERF parameters for each seismogenic source include information on the top and bottom depth of the seismogenic layer, a magnitude independent hypocenter distribution, a magnitude-depth independent distribution of nodal planes, a magnitude scaling relation -all ERFs adopted the area-Mw relationship of Wells and Coppersmith [1994] and a rupture aspect ratio equals to 1. All the data related to point sources (ERFs and other earthquake parameters) use a common grid spaced at 0.1° × 0.1°.
Hereinafter we describe some basic assumptions for the definition of these source typologies and for the modeling [detailed information on model settings can be found in Pagani et al., 2014]. In the case of area and fault sources, seismicity is homogeneously distributed over the area zone. Every rupture has a rectangular shape and is centered on the hypocentral position. Hypocenters are described by a magnitude independent distribution. Ruptures are generated by conserving the area computed using the specified magnitude-area scaling relationship and the corresponding value of magnitude and by the attributed aspect ratio. For point and area sources, the extent of the rupture is limited at the top and bottom by two parallel planes representing the upper and lower seismogenic boundaries. Nodal planes are used to define the three-dimensional position of the ruptures while rake is utilized to attribute the correct coefficients of the GMMs. Meletti et al. [2021] compiled a collection of new datasets to be used as inputs for ERF development in MPS19.

Common inputs
They include, among others: (1) an updated historical-instrumental earthquake catalog [CPTI15, Rovida et al., 2016; 2020], (2) an instrumental earthquake catalog later  ERFs for the MPS19 SHM of Italy historical and a statistical approach. This information is supplied along with the catalog. Moreover, to support the working groups involved in the parametrization of the ERFs, we developed a set of common data regarding the maximum magnitude (Mmax), seismogenic depths, and nodal planes for the study area. Working groups could use none, one, or more of this set of common data, depending on their needs (see a detailed description of the ERFs in Section 3.3). The inputs are described as follows.

Mmax
We based the assessment of Mmax on the estimate of the maximum observed earthquake in the historical seismicity record from CPTI15. First, we divided the study area into different tectonics domains (#1-#19 in Figure   1d), based on (i) seismotectonic data by Meletti et al., [2021], (ii) a regionalization proposed by Delavaud et al.
[2012], for GMMs applicability and, (iii) the superzones (as in the European seismic hazard model, ESHM2013) for the evaluation of the maximum magnitude [Woessner et al., 2015]. Secondly, we assigned the earthquakes in CPTI15 to the respective tectonic domains where they are located. Based on the average error of magnitude estimates for earthquake occurrences pre-and post-1980, we introduced a minimum threshold of uncertainty for the Mw evaluation of 0.2 and 0.3, for the instrumental and historical part of the catalog, respectively. Therefore, the maximum observed magnitude in a tectonic domain is the largest magnitude observed including the uncertainty ("Mw obs + uncertainty"). In Figure 2, we show "Mw obs + uncertainty" and the corresponding value without uncertainty, "Mw obs ", for the 19 tectonic domains.
Then, following Woessner et al. [2015], we introduced a minimum value of Mmax to each tectonic domain (Mw tect ): 6.5 for the active crustal areas in peninsular Italy, 6.0 for the Tyrrhenian tectonic domain, and 5.6 for the volcanic area of Etna. The first value of Mmax, attributed to a tectonic domain, Mw max1 , thus is the largest between Mw tect and Mw obs + uncertainty.
In lack of better qualitative and quantitative data coverage, we defined a second set of Mmax, Mw max2 , with uniform magnitude increments to account for epistemic uncertainties to be of 0.3, except for the Etna volcanic domain where we defined it to be 0.  Figure 2. It must be considered that the Mmax of the CSS (Mw css ) is estimated in different ways: by empirical regression, from the literature, or by associating the largest historical earthquake of the given area to the CSS. In the latter case, however, given the different development of DISS 3.2.1 and CPTI15, there is not necessarily a correspondence between the magnitude associated with a CSS in DISS 3.2.1 and the one provided by CPTI15. It is important to note that this is the case for tectonic domain 2 in Figure 2 (the location of tectonic domain 2 is shown in Figure   1d), where Mw css exceeds the Mw max2 . While Mw css was evaluated using the Mw of the 1690 earthquake as derived from Guidoboni et al. [2007;, the Mw of the same earthquake in CPTI15, including uncertainty, is ~6.5 and it is located in the tectonic domain 3. In this case, we decided to anchor Mmax to the Mw tect which is also in line with the most recent catalog. In the tectonic domain 12 (location in Figure 1d), the highest Mw css derives from the coseismic uplift modeling from Meghraoui et al. [2004] of the Mw 6.8 2003 earthquake. In this case, Mw max2 agrees with the original estimate of the magnitude of the 2003 earthquake; therefore, we assumed Mw max2 as a reasonable value.

Seismogenic depths and hypocentral distribution
To estimate the seismogenic depth parameters for the earthquakes above the magnitude threshold relevant for the PSHA, we used the instrumental catalog by Gasperini et al. [2016]. We assumed the upper and lower limits of the seismogenic layer to correspond with the 5th and 95th percentiles of the distribution of earthquake hypocenter depths within the tectonic domains shown in Figure 1d [according to Boncio et al., 2009], and the prevalent hypocentral depths to correspond with the modal values of the depth distribution. To estimate these values from the instrumental catalog, we removed earthquakes with a fixed depth (i.e., 0 km, 5 km, and 10 km) that represent ~10% of the total. We used the earthquakes that can be associated with the crustal seismicity only; therefore, we removed earthquakes with a focal depth > 40 km. We excluded the Etna area (#17 in Figure   1d) from this computation, as we used an ad-hoc ERF input (see Section 3.4). We removed earthquakes with hypocentral depth < 10 km located within polygon #17, to avoid counting volcanic earthquakes in the tectonic domain #13. For the volcanic Ischia area (#19 in Figure 1d), we adopted "expert judgment" values of 0 km, 4 km, and 1 km respectively for the upper seismogenic depths, lower seismogenic depths, and hypocentral depth, mainly based on the parameters from the 2017 Mw 3.9 seismic sequence [De Novellis et al., 2018].
For each of the 17 tectonic domains (all except the volcanic domain of the Etna #17 and Ischia areas #19), we calculated the 5 th and 95 th percentile of the hypocentral depth distributions, and we considered these values as raw estimates for the upper and lower boundaries of the seismogenic layer. Using the hypocenters within the upper and lower boundaries, we calculated the unimodal and bimodal distributions that best approximate the observed distribution. In the case of a bimodal distribution, weights to the two peaks are given proportionally to their density. We compared the Akaike Information Criterion [AIC, Akaike, 1974] index of both distributions to select which of the two models best approximates the observed depth distributions of the 17 domains. We used different magnitude thresholds, Mw 2.5, 3.5, and 4.5, to calculate percentiles and modal values.  Data are for the tectonic domain #6 (location in Figure 1d).

Nodal Planes
We used data from Roselli et al. [2018] to assign, within the common grid of points, the observed percentage of normal, reverse, and strike-slip earthquakes, and the average distribution of the P (pressure), T (tension), and N (neutral) axes for each of these types of earthquakes. The method is a modification of the Cumulative Moment Tensor (CMT) introduced by Kostrov [1974], where all data are weighted according to their spatial distance from the grid cell. The dataset used in Roselli et al. [2018] is composed of focal mechanisms from the European -Mediterranean RCMT catalog [Pondrelli, 2002;Pondrelli and Salimbeni, 2015] and the Italian CMT dataset from 1976 to the present  combined with present-day stress indicators from the latest stress map release for Italy [Montone et al., 2012;Montone and Mariucci, 2016]. The After a careful classification based on the interval-range of plunges of P, T, and N stress axes [following Zoback, 1992], the dataset consists of 392 focal mechanisms subdivided into 179 normal, 128 reverse, and 85 strike-slip faulting events.
The grid dataset is represented by strike, dip, and rake of the two nodal planes and by a weight proportional to the observed seismic moment. A solution obtained by two strike-slip vertical faults (north-south and east-west oriented) and four normal 60° dipping faults and reverse 30° dipping faults (north-south oriented faults dipping to east and west, and east-west oriented faults dipping to north and south) was assigned to the points without an estimate in Roselli et al. [2018]. The gridded dataset is given in Appendix 3.

National ERFs: Eleven inputs and their uncertainties covering the shallow crustal sources
In 2015, a panel of experts appointed by the Civil Protection Department and INGV engaged the scientific community through a public call for proposals of original ERFs to be used in the MPS19 project. Researchers responded by organizing themselves into working groups and eventually providing eleven independent ERFs, covering the whole Italian territory, and an additional ERF specifically designed for the Etna volcano. The description of the eleven ERFs is organized in groups according to the seismogenic source typologies (see Section 3.1). The abbreviations of the ERFs names follow a nomenclature in Italian that we reuse here for consistency with the description of the MPS19 model reported in Meletti et al. [2021]. Five area-source ERFs (MA for "Modelli ad Area") for which seismicity rates are based on the statistics of the earthquake catalog; two point-source ERFs based on the smoothed seismicity approach (MS for "Modelli basati su Smoothed seismicity"), that consider the spatial distribution of moment rates depending on the density distributions of earthquakes; two fault-source ERFs (MF for "Modelli basati su Faglie") based on three-dimensional seismogenic sources identified through tectonic and geophysical evidence, which include background seismicity represented by point sources; two point-source ERFs based on geodetically-derived seismicity rates (MG for "Modelli derivati dalla Geodesia") in which that consider the spatial distribution of the geodetic strain rate field.
Worthy of note, most of these ERFs (six out of eleven) are zonation-free, implying an effort to overcome the limitations of the classical area sources, and two of them, the two MG for instance, include innovative approaches because they are essentially based on instrumentally measured ground-surface deformation. Table 1 lists the proposed ERFs and a summary of their main features.
Epistemic uncertainty, i. e., how the outcome varies as a function of the variability of the input parameters of the ERF, is explored within each ERF (Table 1). This step is crucial to estimate the overall epistemic uncertainty in the final model, which includes the uncertainty of each ERF and the uncertainty among the quasi-independent ERFs (even though they have been developed independently from each other, many of them use the same amount of information). As stated above, this paper is not aimed at describing in detail every single ERF. Notwithstanding, a brief description of the methodologies of each ERF, along with the indication of the common inputs used to build them, is given in the following and aims to provide a sufficient degree of knowledge of the characteristics of the ERFs. Figure 5 shows the mean rates of occurrence of earthquakes with Mw greater than or equal to 4.5 for the eleven national ERFs described in the next section.  EQ rates replaced in macro-areas = macro-areas (shown in Figure 7) into which earthquake rates of the ERF were replaced by the average earthquake rates.

MA1
MA1 is a seismotectonic zonation built by integrating new and updated seismotectonic data into the so-called SZ9 zonation scheme by Meletti et al. [2008]. The common inputs used by MA1 are: CPTI15, instrumental catalog, completeness periods, Mmax, and dataset of focal mechanisms. The area sources in MA1 are generally more detailed than those in SZ9, and a few new zones have been introduced, including areas previously not considered seismogenic; only where the differences between the new area sources and those of the SZ9 were found negligible, in terms of geographical boundaries and seismotectonic characteristics, the original geometries of the SZ9 zones were used. MA1 evaluated hypocentral depths, seismogenic layers, and nodal planes for each zone, starting from the common data available. After a first test, where each zone had its independent b-value [Santulin et al., 2017], individual seismicity rates have been computed for each zone, using the declustered CPTI15 catalog with the provided completeness, through two different approaches: (i) zones were grouped into macro-areas to evaluate the b-value and then this is used to calculate the seismicity rates within the zones according to a Gutenberg-Richter (GR) distribution (the zones are kept independent for the a-value computation); (ii) observed seismicity rates for each zone.
MA1 explores uncertainties on the maximum magnitude, completeness approaches (historical and statistical), and approaches for seismicity rates (GR and observed).

MA2
MA2 is a two-level area source model. Each level is constructed independently from the other. The first level, named Level 0, includes 24 large zones whose definition is inspired by the so-called superzones of the Woessner et al. [2015] model and aims at identifying areas with common geodynamic and tectonic characteristics reflected in the current stress regime, also considering the presence of volcanoes. Because of their large size, these zones often include different tectonic regimes. Therefore, seismotectonic data are used to derive the prevalence of the associated fault rupture parameters, such as the relative importance of different faulting styles and orientation of nodal planes, the depth of the seismogenic layer, and the prevalent hypocentral depth. The second, named Level 1, is a finer subdivision of Level 0 that includes 54 zones. It is primarily based on the presence or absence of CSSs in DISS 3.2.1. Where present, the CSS guided the geometric description of the source zone and the refinement of the associated fault rupture parameters mentioned above. Occasionally, like in the northern Apennines and the Calabrian Arc, some zones were included to incorporate tectonic features only partially addressed by the DISS 3.2.1. In areas of moderate seismic activity, where the characterization of the CSS is scarce, the Level 1 zones are based on geologic, tectonic, and seismicity knowledge at the regional scale. The common inputs used by MA2 are CPTI15 and DISS 3.2.1. For the seismicity annual rate calculations, MA2 adopts the same approach described below for MA3.

MA3
MA3 is a two-level area source model. The first, Level 0, is a 2D subdivision of the crustal space in large regions (initially set at 25.000 km 2 ), each of them having consistent geotectonic properties inspired by concepts and definitions of plate tectonics (Bates and Jackson, 1980), and already used in earthquake hazards studies [Delavaud et al., 2012;Selva et al., 2016;Basili et al., 2021]. The second, Level 1, is a subdivision of the regions in Level 0, considering geotectonic properties at a local scale, including the presence of volcanoes.
For the seismicity annual rate calculations, Level 0 is used to estimate β (as 2/3 of the GR b-value) and the corner moment of a tapered-Pareto frequency-moment distribution by Kagan [2002a], based on the maximum-likelihood algorithm by Kagan and Jackson [2000]. Several Level 0 areas were merged according to their tectonic category to maximize the number of earthquakes in case of analogous tectonic processes. Level 1 is then used to estimate the annual rate (a-value) by anchoring the theoretical magnitude-frequency distribution to the CPTI15. The relations by Kanamori and Brodsky [2001] were used for the conversion between seismic moment and magnitude. Regarding the earthquake catalog, MA3 used the declustered version, considering the completeness, with a threshold magnitude set to Mw ≥ 4.8 and a bin size of 0.2 magnitude units. To limit the model to tectonic crustal processes, the events classified as volcano-tectonic ones and the deep earthquakes associated with subduction processes were removed.
To explore the epistemic uncertainty, MA3 developed eight different frequency-moment distribution models, two adopting either type of completeness (historical and statistical) and four weighted alternatives for the anchoring.
The first uses a fixed magnitude value of Mw = 5.5; the second uses the magnitude bin with the relatively highest rate; the third uses the class that yields the lowest residuals; the fourth uses the magnitude value that yields the lowest RMS. The assigned weights were defined as 0.3, 0.1, 0.3, and 0.3, respectively.

MA4
The area sources of MA4 represent an update of the so-called SZ9 zonation scheme (Meletti et al., 2008). To estimate the expected seismicity rates, MA4 were grouped the zones into macro-areas, according to their tectonic regime. Assuming that the distribution of the earthquake sizes follows a truncated GR model, MA4 applies the Weichert [1980] approach to estimate the parameters of the model for each macro-area.
Subsequently, four approaches were used to distribute the activity rates from a macro-area to the zones included. The first approach was based on the average ratio of the observed rates of occurrence in the zones in respect to the macro-area they belong to. The second one used the slope of truncated GR distribution of the macro-area and the annual rate of exceedance of the minimum magnitude of completeness of each zone, ensuring that the resulting distributions add to the one of the macro-area. The third approach maximized the log-likelihood of the observed and forecasted number of earthquakes in each zone. The fourth approach minimized the RMS of the observed rates of occurrence of each area. In addition, Weichert [1980] was also applied to the single zones, with upper and lower bounds of the slope value of the model.
MA4 explored uncertainties of completeness (historical and statistical), maximum magnitude, and activity rates.

MA5
MA5 uses the same SZ as MA4 but a different approach to obtain the seismicity rates. The common inputs used by MA5 are CPTI15 and Mw max2 . MA5 jointly estimates the completeness time and the seismicity rate of the complete part of the datasets (mainshocks in CPTI15) of each seismogenic zone, which are constituted by events of magnitude exceeding a varying threshold [Rotondi and Garavaglia, 2002]. Assuming that the main Francesco Visini et al.
earthquakes follow a homogeneous Poisson process with constant rates and that the beginning of the complete part of the catalog corresponds to a change in the Poisson parameter, the issue of completeness time and seismicity rate estimation was translated into a change-point problem which combines stochastic simulation methods with a Bayesian hierarchical model [Pievatolo and Rotondi, 2008]. In detail, MA5 assumes that the process has two rates, one in the incomplete part and another in the complete one; following the Bayesian approach, both rates and the time when the process changes were considered as random variables, and the method computes, in addition to their estimates, the entire probability distribution of these variables as well.
Hence, MA5 provides a more informative result than those of the most used methods. For example, it was possible to get more changepoints among which to choose the best one based on cultural and historical information or the associated probability [Rotondi and Garavaglia, 2002]. The events of the parametric catalog associated with each seismogenic zone were partitioned into subsets taking 4.45 as the initial magnitude threshold and then increasing this value gradually by 0.1 bin size, up to the maximum magnitude observed in that zone. The analysis was repeated in each subset, and the estimated rates were smoothed by a power function to give the rate of bins without observed events as well and to extrapolate the rate value up to the maximum magnitude. Moreover, following the Bayesian approach, MA5 is able to exploit prior information contained in other databases, allowing it to analyze even scant sets of data. MA5 used the best estimate of the completeness times and the seismicity rates given by this approach.

MF1
MF1 is a fault-based ERF entirely relying on the geological data from the DISS 3.2.1. Other datasets, such as the common inputs CPTI15 and focal mechanisms, stress orientations [Carafa and Barba, 2013], and crustal models [Laske et al., 2013], were used to constrain various parameters of the model but not the seismicity rates.
Although MF1 was based on fault geologic information only, an empirical proximity function was specifically developed to predict off-fault seismicity, i.e., the occurrence of earthquakes of any size away from the faults.
These characteristics ensured the independence of the seismicity rates predicted by the ERF from the earthquake catalog, thereby providing a valid alternative to explore the epistemic uncertainty in PSHA. The basic idea was that potentially seismogenic faults could provide a representation of the seismicity over a time span much longer than that of earthquake catalogs. The seismic moment rate of a seismogenic fault is derived from the geologic moment rate by multiplying the shear modulus, by fault length and width, by slip rate, and by a coefficientoften referred to as seismic efficiency -that determines how much of the geologic rate is converted into the seismic rate. Using the moment conservation principle, the seismic moment rate of any seismogenic fault is related to the annual earthquake rate of a tapered GR magnitude-frequency distribution, using the formulations by Kagan [2002a;2002b] and assuming a single β equal to 0.667 (corresponding to GR b-value = 1) for all the faults. To evaluate the off-fault seismicity rates, MF1 followed an empirical approach aimed to capture the natural distribution of the observed earthquakes around the faults and the respective location uncertainties. To this aim, an empirical proximity function representing the frequency at which an earthquake has occurred at any distance from a mapped fault was developed. Then, the earthquakes predicted by the calculated magnitude-frequency distributions, as previously described, are symmetrically spread around proportionally to the magnitude and distance pairs dictated by the empirical proximity function. MF1 considered uncertainties in the slip rates assuming they also incorporate different levels of the seismic efficiency.

MF2
MF2 [Murru et al., 2020] combines the seismic rates obtained by different datasets from the common inputs to cover all the Italian territory: DISS 3.2.1, CPTI15, the instrumental catalog, seismogenic depths, hypocentres, and nodal planes solutions. For the historical catalog, both statistical and historical criteria of completeness [Meletti et al., 2019] were used. MF2 adopted, for the whole of Italy, a tapered GR relation with a uniform b-value (0.99) and a uniform corner magnitude (Mw 7.3), both estimated by the seismic catalogs. MF2 is composed of two different source types merged into a single ensemble ERF: the first realized using the individual seismogenic sources in the DISS 3.2.1, the second using a smoothed seismicity approach over a grid of points. For each individual seismogenic source, MF2 considered the following information: location, dimension, and slip rate.
The use of geological information overcomes the problem of the short time length of the available catalogs.
The seismic moment rate of each individual source was computed using the seismic-moment conservation principle [Field et al., 1999], which allows deriving seismic moment rate from long-term slip rate and geometry of the seismogenic source. Then, balancing the seismic moment rate over a tapered GR magnitude-frequency distribution from magnitude Mw ≥ 4.5, seismicity rates were evaluated for each individual source. In addition, using the smoothed-seismicity approach proposed by Frankel [1995], with the correction for the time-varying completeness magnitude [Hiemer et al., 2014], and considering as input the combined historical and instrumental catalogs, the rates of occurrence for Mw

MS1
MS1 is developed using smoothed seismicity methods for long-term earthquake rate space-time forecasts of future Italian seismicity for MPS19. The common inputs used by MS1 are: the CPTI15, seismogenic depths and hypocentres, and nodal planes solutions. A forecast of the probabilities of (Mw ≥ 4.5) earthquakes is calculated from ensemble smoothed seismicity models derived from CPTI15 (1000-2014) and the 1981-2016 instrumental earthquake catalog by Gasperini et al. [2016]. Two different smoothed seismicity methods are adopted for the MS1 construction following: 1) the well-known and widely applied fixed smoothing [Frankel, 1995], such that the kernels used to smooth catalog-derived seismicity rates are invariant to spatial variations in seismicity rate and 2) the adaptive smoothing method [Helmstetter et al., 2007], which determines the smoothing distance for each earthquake from its distance to neighbouring earthquake epicenters. The optimized fixed smoothing distances and adaptive neighbour numbers are identified by ranking the ability of all the trial smoothed seismicity models to forecast the location of earthquakes in the recent part of the earthquake catalog with a likelihood parameter.
Finally, MS1 is created as an ensemble model that combines two different but equally weighted smoothing models (adaptive and fixed) from the two earthquake catalogs through a logic tree to improve the forecast capability. MS1 uses the common input data for seismogenic layers, hypocentral depths, and nodal planes.
MS1 considers the uncertainties in the model constructions using different percentage combinations of the models obtained from two different smoothing approaches [for details see, Akinci et al., 2018].

MS2
MS2 followed the Woo [1996] approach and used a zone-free method solely based on the use of the earthquake catalog CPTI15. The common inputs used by MS2 are: the CPTI15, seismogenic depths and hypocentres, and nodal planes solutions.
Proxy seismogenic sources were created from the epicentral locations of the events that were smoothed according to their fractal distribution in space. A square sub-grid of point sources (approximate dimensions 10 x 10 km, with step 0.2 km) was defined around each node of a calculation grid (0.1° × 0.1°), where the activity rates of each point source were calculated from the density and proximity of events lying within the considered Francesco Visini et al.
ERFs for the MPS19 SHM of Italy magnitude range. The contribution of each earthquake to the seismicity of the region was smeared over all subgrid points falling within an epicentral distance that depends on the magnitude of the event itself. In doing so, Woo's (1996) method uses concepts from fractal geometry [e.g., Kagan and Jackson 2000], and in contrast to the rigid zonation of Cornell's approach (which assumes a uniformly distributed seismicity within each seismogenic zone), it depicts seismicity that is spatially non-uniformly distributed. The magnitude-dependent relationship is defined by a kernel function. Finally, the activity rates associated with each node of the calculation grid were computed by summing the activity rates of the point sources multiplied by the area of the sub-cells.
To consider the uncertainties associated with some assumptions, MS2 explored eight combinations resulting from kernel shapes, kernel directionalities, and completeness periods.

MG1
MG1 is a geodetically-derived seismicity ERF based on the analysis of the common data of 919 GPS horizontal interseismic velocities and assuming that the future seismic moment release will scale linearly with the rate of present-day strain accumulation [Ward, 1998;Ward, 2007]. The strain rate tensor field has been calculated on a regular 0.1° × 0.1° grid using the VISR software [Shen et al., 2015], considering the variable station spacing for the optimal smoothing parameters, and then converted into rates of seismic moment accumulation following Kostrov [1974]. The rate of seismic moment accumulation density was finally converted to earthquake rate density (or earthquake potential) under the assumption that seismic moment distributes into earthquake sizes that follow a tapered G-R distribution of given b-value and corner magnitude [Kagan, 2002b]. MG1 uses two alternative approaches, Ward [1998] and Savage and Simpson [1997], to treat the non-uniqueness of the relationship between scalar strain rates and seismic moment. The effect of uncertainties and spatial variations in the seismogenic thickness has been introduced by considering alternative realizations of the seismogenic thickness, ranging between 7 and 13 km [from Chiarabba and De Gori, 2016]. The forecasted seismicity rate has been scaled to the seismic moment release of the declustered CPTI15 catalog in the interval 1787-2015 to account for the effects of declustering and aseismic components. Following the approach described above, the seismicity rates forecasted by MG1 satisfy the overall 1787-2015 CPTI15 moment release while respecting the spatial variations shown by the geodetic strain rate field. Finally, MG1 also used the common input data of nodal planes.

MG2
MG2 uses common data of GPS interseismic horizontal velocities and stress azimuth data (Carafa and Barba, 2013) to determine both the interseismic and long-term strain-rates and velocities on a finite element grid by means of the NeoKinema code version 5.2 [Bird and Liu 2007;Bird 2009;Bird and Carafa 2016]. Short term geodetic and non-tectonic (e.g., landslide movements) transients are modeled as principal factors in the covariance matrix  before running NeoKinema, and the strain rate in Calabria is calculated only for the upper plate, under the assumption of locked subduction [Carafa et al., 2018]. The methodology to compute models of the long-term earthquake rates density for shallow events from the long-term strain rate follows Bird and Liu [2007], Bird et al. [2010], and Bird and Kreemer [2015]. The parameters of the Tapered GR frequency-magnitude relations and seismic coupling were computed at a regional scale. In detail, the upper crustal earthquakes are assigned from the CPTI15 catalog to their presumed causative faults from the DISS 3.2.1 database. Then, the events are grouped into three sub-catalogs, corresponding to the compressional, extensional, and strike-slip fault classes, and the seismicity parameters, including the seismic coupling factor, are computed [see Carafa et al., 2017 for details]. Hence, the seismic moment rate for each cell of the grid is calculated in a two-step procedure: first, the interseismic strain-rates of the portion of the finite elements falling in the individual grid cells are converted to long-term tectonic moment rates; second, the long-term tectonic momentrates are converted to seismic moment rates. MG2 adopted three alternatives for smoothing parameter of SHmax and GPS and four alternatives of coupled thickness.

Common ERFs: Volcanic, Subduction, and External sources
We developed three ERFs for (a) the Etna area (tectonic domain 17 in Fig. 1d); (b) the subduction shallow interface seismicity and deep intra-slab seismicity of the Calabrian Arc spanning from the Ionian Sea to the southern Tyrrhenian Sea across the Calabria region; and (c) the seismogenic sources located outside the area of the CPTI15 catalog. The three ERFs are described in the following.
The ERF for Etna was adapted from Azzaro et al. [2017] and Peruzza et al. [2017]. This ERF includes a rather high degree of detail and complexity, motivated by the huge amount of geological, seismological, and geodetic data for the area. Originally, this ERF consists of three alternatives that explore the impact of different seismic source typologies (areas, faults, and background) and earthquake occurrences (time-dependent and timeindependent). It also considered the topographic surface, which influences the site-source distances in this high relief area. For the purposes of application to the MPS19 project, the time-dependent model was discarded, and the focus was directed exclusively on the Poisson model. Also, the topography was not incorporated for the sake of maintaining homogeneity with the other models on a national scale, and therefore the sources were all modeled with respect to the sea level. From the three available alternatives, we selected the one that combines individual sources (seismogenic faults) and gridded seismicity. The individual sources model the strongest seismicity (Mw≥ 4.5) using a historical approach, and therefore with the rates estimated by the maximum magnitudes observed and relative average times of recurrence, whereas the gridded seismicity models the smaller earthquakes (3.5 ≤ Mw < 4.5) with point sources with variable depth.
The ERF input for the subduction of the Calabrian Arc uses the slab geometry from Maesano et al. [2017] that has already been used for tsunami hazard from seismic sources [Basili et al., 2021;Scala et al., 2020]. We Moho [Grad et al., 2009] of the upper plate, respectively. The seismic moment rate for the interface was calculated using the convergence rates derived from Davies et al. [2018] and Carafa et al. [2018], multiplied by the shear modulus, the interface area, and the seismic efficiency coefficient. The seismic moment rate of the slab interface was related to the annual earthquake rate by applying the moment conservation principle to a tapered GR magnitude-frequency distribution using the formulations by Kagan [2002a;2002b] and assuming a single β equal to 0.667 (corresponding to GR b-value = 1). For the intraslab seismicity, we resampled the three-dimensional geometry of the slab with a 10-km-spaced grid of points. The slab thickness was taken as equal to the crustal thickness (seafloor to Moho) measured in the undeformed part of the subduction (Ionian Sea), which then constrained other seismogenic parameters. In particular, the expected nodal planes that were assigned to the grid points have earthquake ruptures at 45° angles with respect to the local slab dip, and the rupture size is confined within the subducting slab. The expected seismicity rates were calculated according to the Weichert [1980] approach using, as input, the deep earthquakes in the CPTI15 catalog, the statistical approach for estimating completeness, and a magnitude-frequency distribution of the type Tapered, with corner magnitude equal to Mw 7.5. The total magnitude-frequency distribution of the seismicity rate was uniformly distributed on the grid points. The ERF input for the interface was added to the ERF MF1 and MG2, which, in their development, had been able to separate seismicity due to the interface from the crustal tectonic seismicity.
To parametrize sources located outside the area of the CPTI15 catalog, in order to minimize the edge effect, we started from the ESHM2013 European hazard model [Woessner et al., 2015], which for extent, methodology and typologies of ERFs, was deemed useful for this purpose. The ESHM2013 model consists of three types of ERFs: (i) areas; (ii) faults and background (FSBG); (iii) smoothed seismicity and faults [SEIFA, Hiemer et al., 2014].
To avoid as much as possible the overlapping of MPS19 sources with the external ones, we resampled, where necessary, the three ERFs and used only the sources located outside the national ERFs. First, the area source and the FSBG ERFs were resampled according to a regular grid spacing 0.2° × 0.2°. Then, to each point of this grid, we assigned a weighted average seismicity rate given by the mean of the seismicity rates of the three ERFs. This was made in accordance with ESHM2013 for the return periods between 475 and 2475 years, whose weights are 0.5, 0.2, and 0.3 for areas, FSBG, and SEIFA, respectively. Finally, for each of the eleven national ERFs, we included only the points located outside the area covered by the ERF and within 100 km from the outer edge.

Francesco Visini et al.
In Figure 6, we show the areas covered by the three common ERFs in respect to the area of interest of the eleven national ERFs.

The final ERFs, their weighing, and the ensemble modeling
We adopted the following procedure to weigh the eleven national ERFs in the framework of MPS19. We subjected each of the eleven ERFs to an elicitation procedure and consistency tests. The elicitation procedure aims at investigating if an ERF may be considered plausible and sound from a group of independent experts through a formal experts' elicitation. Consistency tests were performed to evaluate whether the average rate and uncertainty obtained for each ERF (by a weighted combination of the alternatives of every single ERF) describes satisfactorily the spatial distribution of seismicity and the total number of events in the last centuries [see, for instance, Werner et al., 2010]. In particular: (i) we integrated the CPTI15 with the earthquakes that occurred in 2015 and 2016, whose magnitude values were extracted from the instrumental catalog by Gasperini et al. [2016], and (ii) we adopted two pairs of time and magnitude of completeness: Mw ≥ 6 from 1660 and Mw ≥ 5.6 from 1787. The consistency tests were the N-test and the S-test. The N-test verifies the consistency of models with the total number of observations, while the S-test verifies the consistency of models with the spatial distribution of observations. These tests have already been used to check the performances of medium-term ERFs in the Italian region . A model that passes both tests explains satisfactorily the overall number of target earthquakes and their location.
We finally applied two scoring methods to the N-and S-tests, the pari-mutuel gambling score, and the Loglikelihood score [Scherbaum et al., 2009;Zechar and Zhuang, 2014] to finally obtain a weight for each ERF.
The weight (W) of each ERF is a combination of the elicitation (W1) and the results of the consistency tests (W2): where C is a normalization factor that sets W in the range [0,1]. The weight W1 is based on expert opinion elicited from a pool of experts that analyzed the details of the implementations of all ERFs. Similar to Basili et al. [2021], the ERF weights have been evaluated by applying an Analytic Hierarchy Process [AHP, Saaty 1980Saaty , 1998Goepel 2013]. The weight W2 instead favours ERFs with strong compatibility with a subset of the input data, specifically selected for testing purposes. More details of the elicitation procedure, consistency tests, and weighting are described in Basili et al. [2021].
All the eleven national ERFs passed the elicitation and consistency tests. The adopted weights W, as well as the W1 and W2, are given in Table 1.
As described in Meletti et al. [2021], we carried out consistency tests also on macro-areas to investigate potentially different performances of the ERFs depending on the area, given the different reliability of the input data used to build each ERF. To define the regions to be tested, we based on the following criteria: avoid taking more than one completeness zone, or at least limit the overlap with a second completeness zone as much as possible; avoid the proliferation of the number of regions; define regions to be as large as possible; avoid maximizing the result of the ERF closely related to the seismic catalog; define regions to be as homogeneous as possible from the tectonic viewpoint; avoid or minimize the offshore part included in the regions; separate the regions with a relatively higher amount of geological data, earthquakes, and geodetic stations from the relatively poorer ones. Following these criteria, we defined the seven macro-areas shown in Figure 7. To handle the case when sources of an ERF span multiple macro-areas (e.g.., ERFs MA1 to MA5 and MF1), we first resampled ERFs using a grid of points 0.1° × 0.1° spaced, and then applied consistency tests to the ERFs for the seven macro-areas.
The results of these regional tests are listed in Table 1 and indicated if an ERF input failed to represent the observations. In detail, MA2, MA3, and MA5 failed the S-test for macro-area 7; MF1 failed the S-test for macroareas 1 and 3; MS2, MG1, and MG2 failed the N-test for macro-area 4; and MG2 failed the N-test for macro-area 3 and the S-test for macro-area 6. In macro-areas where one or more ERFs failed a test, we replaced the seismicity rates of those ERFs by the weighted average seismicity rates obtained from the combination of the seismicity rates of the ERFs that passed the consistency tests in that region. This implied, for example, that the subduction interface ERF was weighed out together with the portion of MF1 and MG2 contained in the macro-area 6, whereas the intraslab ERF remained with all the other ERF.
We named the ERFs after this regional weighing procedure as "final" ERFs. Finally, we extended back to Mw 3.5 the rates of occurrence of the ERF in the Ischia island, as already done for Etna. Ischia is characterized, in fact, by volcanotectonic seismicity, with earthquakes with magnitude in the range ~3.5 -5, resulting in severe ground shaking, as testified by the effects associated with the Md 4.0 2017 Casamicciola earthquake [Nappi et al. 2018;Cubellis et al., 2020].

ERFs for the MPS19 SHM of Italy
We used the left tail of the magnitude frequency distribution of each final ERF to extrapolate the earthquake occurrences down to Mw 3.5. Figure 8 shows the mean rates of occurrence of earthquakes with Mw larger or equal than 4.5 of the eleven national ERFs, after replacing earthquake rates with the average values for the ERFs indicated in Table 1.
After having defined the final ERFs and their weights (Table 1), we estimated the epistemic uncertainty through a quantitative evaluation of the dispersion of the forecasts given by the different final ERFs. In particular, following the probabilistic framework by Marzocchi and Jordan [2017], the forecasts provided by the eleven final ERFs and their different parameterizations are meant to sample the epistemic uncertainty for details, see also Meletti et al., [2021]. To build the inputs used for MPS19, we added Etna, the subduction, and the external region inputs described in Section 3.4 to the national-scale ERF. To show the epistemic uncertainty resulting in the evaluation of the seismicity rates, in Figure 9, we plot the seismicity rates computed by the final eleven ERFs and resulting in 94 alternatives. In Figure 9, it is also possible to note the effect of the replacement of the original earthquake rates with the average solution.

Discussion
The ERFs, along with the ground motion models, are the key points for the input in PSHA. One of the most innovative aspects of MPS19 concerns the quantification of the epistemic uncertainty, following the hierarchy of uncertainties defined in Marzocchi and Jordan [2017]. Specifically, we sampled the epistemic uncertainty through a set of forecasts of the Italian seismicity, which is given by the ensemble of the alternative interpretations of each ERF [Marzocchi et al., 2015].
The ensemble accounts for eleven alternative ERFs, five of which are based on area sources, two on active faults plus background (or off-fault) seismicity, two on smoothed seismicity, and two on geodetic-derived strain rates, for a total of 94 forecasts. Although to a lesser extent, the epistemic uncertainty is also explored through the differences within the groups of ERFs. For example, the different geometric realizations of the ERFs based on area sources result from the application of significantly different principles.
With respect to the ERF currently used for the official seismic hazard model of Italy [ZS9 in Meletti et al., 2008], the described ERFs represent an important step forward. ZS9 was based on a single seismotectonic scheme, and its epistemic uncertainty only dealt with the completeness of the earthquake catalog and the earthquake rate estimations treated by four logic-tree branches. Moreover, we are conscious that the epistemic uncertainty of the ensemble does not represent the full distribution of all source parameters, but we certainly enlarged our view with the inclusion of alternative area sources and ERFs not based on earthquake catalogs.
The number of ERFs and the difference in their approaches made use of updated and harmonized earthquake parameters [e.g., Rovida et al., 2016] et al., 2017], and required the revision or development of tools for analyzing input data and computing expected seismicity rates. The ERFs for MPS19 can be considered innovative under several aspects; moreover, the main point we highlight here is that the ERFs were developed within a community-based effort that promoted advanced and original methods for earthquake modeling. We submitted ERF methodologies to an experts' opinion procedures to verify their scientific robustness and tested the consistency of each ERF with respect to the observed earthquake occurrences to avoid inputs that under-or overestimate observations. We replaced earthquake rates for those ERFs that do not pass the tests in some regions. Figure 9 displays, for the seven macro-areas tested, the earthquake rates generated by the ERFs, for three thresholds of magnitudes (Mw 4.5, 5.5, and 6.5), compared with the observed rates calculated from the CPTI15 using both statistical and historical magnitude-dependent completeness intervals by Meletti and Marzocchi [2019]. The ERFs for which we replaced earthquake rates are shown by a yellow square. Figure   9 qualitatively shows the relationship between the observed and forecasted seismicity rates by each ERF. It is interesting to note how, for the ERFs extensively based on the seismic catalog for rate estimation (e.g., ERF from MA1 to MA5, MS1, and MS2), Figure 9 represents, in a certain way, an index of the approach's maturity. The areabased ERFs, for example, return forecast ranges that include the observed rates. It is also important to note the consistency and, therefore, in some ways the maturity, of the fault-based ERFs, whose forecasted seismicity rates tend to replicate the observations, apart from a few areas where data are particularly scarce (e.g., the west coast of peninsular Italy). The geodetic-based ERFs were used here for the first time to produce nationwide models. These   Figure 7. Circles indicate the seismicity rates for each alternative of each ERF, for magnitude greater than or equal to 4.5, 5.5, and 6.5. Yellow squares indicate the seismicity rates of the average model that replace the ERFs. Blue and red lines indicate the observed seismicity rates according to the historical and statistical completeness intervals.
ERFs certainly differ more than the others from the observed ones. Indeed, this can be considered an added value, as such inputs allow us to see the ongoing deformation in a different and integrative way from the seismic catalog.
Furthermore, the problem of estimating seismicity rates from a single point of view, the seismic catalog, for example, can lead to overfitting situations that do not allow us to consider anything more than what has been observed. Surely, both geodetic and fault ERFs need more time to be weighed with more incisiveness, but it should be emphasized that, along with the development of geodetic and fault ERFs, the tests for weighing need to be improved, too.
Regarding the ERF development, which is the core of this work,  has shown, for example, that a model of elastically unloading seismogenic faults within a viscously deforming lithosphere can use GPS and stress direction data to determine how the extension is partitioned between fault slip and bulk lithosphere permanent strain. Geodetic velocities are collected in the interseismic tectonic stage with steady secular deformation; thus, the main improvement of geodetic models rests in the estimation of long-term average slip rates for active faults of Italy using GPS measurements.
Depending on the scale of the application, fault-based ERFs may suffer from the inhomogeneous distribution of the available data, which tends to increase the overall epistemic uncertainty of the model [e.g., Basili et al., 2013] as the size of the study area increases. To cope with these limitations, so far most analyses have been developed at a regional scale, for example in central Italy [e.g., Pace et al., 2010;Peruzza et al., 2011;Visini and Pace, 2014]. The availability of detailed databases of fault geometries and slip rates [e.g., Valensise and Pantosti, 2001;Boncio et al., 2004;Basili et al., 2008;Burrato et al., 2008;Kastelic and Carafa, 2012;Maesano et al., 2015;Valentini et al., 2017;Valentini, 2020;Faure Walker et al., 2021] is paramount to build alternative PSHAs. The definition of the 3D geometries of faults and the discrimination of the tectonic components from other aseismic processes remain however the inherent challenges of such databases and underline the need to explore the uncertainties that arise from them. For some faults, for example, slip rates may bear the effects of sediment compaction [e.g., Scrocca et al., 2007;Maesano et al., 2015], gravitational movements [Kastelic et al., 2017;Di Naccio et al., 2019], or aseismic behaviour . Moreover, fault-based ERFs could potentially benefit from the recent development of multi-fault rupture modeling [Field et al., 2013;Chartier et al., 2019;Verdecchia et al., 2018;Visini et al., 2020;Valentini et al., 2020;Scotti et al., 2021]. Although certainly more mature than the others, areabased ERFs can also benefit from other improvements; the choice of criteria for drawing the areas, and in particular, the balance between expert judgment and data-driven approaches can still be further explored.
Finally, we underline that all ERFs include a systematic analysis of the most important uncertainties (epistemic uncertainty), according to experts. The testing phase did not highlight the existence of blatant unknown unknowns (ontological error); however, as for most hazard analyses in all fields, the discovery and reduction of unknown unknowns remain an important challenge.

Conclusions
The MPS19 project by Meletti et al. [2021] aimed at building a new seismic hazard model for Italy, responding to the requirements of the Italian Civil Protection Department to build a time-independent model based on newly available data acquired in the last decades and on reliable, consolidated, and testable components. In this paper, we presented the ERFs used in MPS19. In particular, we defined eleven ERFs for the entire Italian territory, with the addition of (a) an ERF for the volcanic region of Mt. Etna; (b) an ERF for the deep seismicity related to the subduction in the southern Tyrrhenian Sea; (c) an ERF for the seismogenic sources located outside the area covered by the CPTI15 catalog.
For each of the eleven ERFs, we explored the epistemic uncertainties, so in the end, we collected a total of 94 ERFs for MPS19.
We used all the recently published data, and we implemented the most advanced practices by introducing several novelties for defining the area source, smoothed seismicity, faults, and geodetic ERFs.
We tested the ERFs at both a national scale as well as in macro-areas to check the variation of the forecasting capability in space and as a function of the number of input data. We remark that if an ERF does not pass both consistency tests in a macro-area, its overall earthquake rate is replaced by the weighted average of the other consistent ERFs, preserving the geometries, depth, and focal mechanisms of the tested ERF.
The entire approach is a first attempt to build an inclusive model of ERFs, involving a large number of seismic Francesco Visini et al.
hazard practitioners with their models, to capture and explore a broad range of epistemic uncertainty. While some phases can be further improved, the overall MPS19 project represents an important step forward for a national seismic hazard model assessment and such an approach can be applied worldwide.