An updated area-source seismogenic model (MA4) for seismic hazard of Italy

A new probabilistic seismic hazard model, called Modello di Pericolosità Sismica 2019 (MPS19), has been recently proposed for the Italian territory, as a result of the efforts of a large national scientific community. This model is based on 11 groups of earthquake rupture forecast inputs and, particularly, on 5 area-source seismogenic models, including the so-called MA4 model. Data-driven procedures were followed in MA4 to evaluate seismogenic parameters of each area source, such as upper and lower seismogenic depths, hypocentral-depth distributions, and nodal planes. In a few cases, expert judgement or ad hoc assumptions were necessary due to the scarcity of data. MA4 consists of 20 seismicity models that consider epistemic uncertainty in the estimations of the completeness periods of the earthquake catalogue, of maximum magnitude values and of seismicity rates. In particular, five approaches were adopted to calculate the rates, in the form of the truncated Gutenberg–Richter frequency–magnitude distribution. The first approach estimated seismicity rates using earthquakes located in each area source, while the other approaches firstly calculated the seismicity rates for groups of areas considered tectonically homogeneous and successively partitioned in different ways the values to the area forming each group.

The results obtained in terms of seismic hazard estimates highlight that the uncertainty explored by the 20 seismicity models of MA4 is at least of the same order of magnitude as the uncertainty due to alternative ground motion models.

1 Introduction

A recent project led by the Seismic Hazard Centre (Centro di Pericolosità Sismica, CPS) of the Italian Istituto Nazionale di Geofisica e Vulcanologia (INGV) and funded by the Italian Civil Protection Department produced a new time-independent probabilistic seismic hazard assessment (PSHA) model for Italy, called Modello di Pericolosità Sismica 2019 (MPS19; Meletti et al., 2021). The model consists of a suite of earthquake rupture forecasts (ERFs) and ground motion models (GMMs), described in Visini et al. (2021) and Lanzano et al. (2020), respectively, that are based on updated and new data acquired in the last decade after the release of the current reference Italian seismic hazard model in 2004–2006 (MPS04, Modello di Pericolosità Sismica 2004; Stucchi et al., 2011), which is currently adopted as seismic input in the Italian building code (NTC, 2018).

In particular, MPS19 consists of 564 alternative seismic hazard models (i.e. logic-tree branches) obtained by combining 11 groups of ERFs, each made by a different number of sub-models (for a total of 94 ERFs) to consider the epistemic uncertainty inside each group, with a set of six GMMs (three for active shallow crustal regions, two for subduction zones and one for volcanic areas). In terms of seismic source typologies, 5 groups of ERFs out of 11 consider area sources: 2 are based on smoothed seismicity calculated on a grid of points; 2 combine faults sources with background seismicity; and 2 derive earthquake rates from geodetic data over a grid of points. The ERFs are based on updated and new historical, geological, geodetic and palaeoseismological datasets collected over the last 15 years, since the realization of MPS04, for the Italian territory and its conterminous areas.

In this paper we describe one of the five area-source ERFs that is the so-called MA4 model (i.e. area-source model no. 4, modello ad aree no. 4 in Italian), based on seismogenic zoning ZS16 (zonazione sismogenetica 2016). ZS16 represents the evolution of previous area-source models proposed in the last 30 years as unique area-source input in seismic hazard assessment performed in Italy; these models are all based on the same seismotectonic approach to seismogenic zoning described by Meletti et al. (2000).

The ZS16 of the MA4 model incorporates a number of different parameters for each defined area source: (a) geographical boundaries, (b) top and bottom depth of the seismogenic layer, (c) hypocentre distribution, and (d) style of faulting. For each area source of ZS16, MA4 uses five alternative frequency–magnitude distributions, providing the annual rates of all earthquakes with M_w≥4.5, which is the threshold magnitude adopted in MPS19 (Meletti et al., 2021; Visini et al., 2021).

In the following, we first briefly introduce the input data considered for developing the ZS16 and the MA4 model and then describe the methods used to define the geometry of area sources and to estimate, for each of them, the top and bottom depth of the seismogenic layer, hypocentre distribution, style of faulting, and annual rates of earthquake occurrence. Finally, seismic hazard estimates computed using the MA4 model are presented and discussed.

2 Chronology of seismogenic zonings and PSHA models for Italy

The MA4 ERF is based on the area-source approach for defining seismicity parameters. The choice operated in this work was to update previous zonings designated for previous seismic hazard assessment in Italy. Since 1990, several seismogenic zonings have been released, adopting the seismotectonic approach described in Meletti et al. (2000). In the first half of the 1990s, the first zoning adopted in a seismic hazard computation for the whole national territory was ZS4 (zonazione sismotettonica vers. 04, Meletti et al., 2000, shown in Fig. 1), used for the PS4 (pericolosità sismica vers. 04) hazard model (Slejko et al., 1998). ZS4 was delineated interpreting the seismicity in terms of tectonic regimes and, as a second-order criterion, the spatial variation in seismicity; the model was constituted by 80 area sources.

Figure 1 (a) Comparison of ZS4 by Meletti et al. (2000) and ZS9 by Meletti et al. (2008). ZS9 is the seismogenic zoning of the current Italian seismic hazard model (MPS04) by Stucchi et al. (2011). The sketch in (b) shows the relationships between the seismogenic zoning and the seismic hazard (PSH, probabilistic seismic hazard) models developed in Italy. Note that only for the last model, ZS16, did we distinguish the seismogenic zoning from the earthquake rupture forecast model.

At the beginning of the 21st century, an updated version of the area-source zoning (ZS9; Meletti et al., 2008, shown in Fig. 1) was defined for the elaboration of the MPS04 seismic hazard model (Stucchi et al., 2011). ZS9 resulted from modifications, merges and eliminations of the numerous areas delineated in the previous zoning of ZS4, as well as from the introduction of new ones. The goal of ZS9 was to build a model consistent with new data collected at the time of its development.

Since most of that knowledge was considered still reliable during the development of ZS9, this later was built without introducing substantial novelties as regards the general kinematic framework on which ZS4 was based. In some cases, groups of area sources of ZS4 were merged on the basis of the characteristics of the kinematic domain to which each of the area sources was attributed. In the meantime, the geometry of the area sources was modified according to the changed seismotectonic knowledge. Most importantly, in ZS9 area sources were designed strictly enveloping the seismogenic sources that were at that time known and defined in the DISS database (Database of Individual Seismogenic Sources; Basili et al., 2008). In ZS4, on the contrary, the areas extended over the known seismogenic sources, including regions where faults were not mapped, according to what was thought to be a cautionary criterion. We should consider that, for PSHAs using an area-source model, the seismicity rates computed using faults and earthquakes located inside the area are equally spaced in a grid of point sources where each point has the same seismicity occurrence properties (i.e. rate of events generated). The arbitrary increase in the surface of some of the area sources of ZS4, therefore, led to a reduction in the hazard estimate of PS4 in those areas. ZS9 was then developed by constraining the geometry of the area sources to the location of seismogenic faults and historical and instrumental earthquakes, avoiding arbitrary extensions of the dimensions of the area sources. Figure 1 shows the ZS9 model, consisting of 36 area sources, together with ZS4.

The development of ZS16 was driven by the choice to update the area-source model of ZS9 only where new data suggest different interpretations. To summarize, as shown in Fig. 1, ZS16 updates ZS9 and constitutes the base for the MA4 ERF. MA4, along with other ERF models, was used as one of the inputs of MPS19 (Meletti et al., 2021; Visini et al., 2021).

3 Input data

Area sources for PSHA represent regions with seismicity spatially uniform in terms of earthquake occurrence rates, maximum magnitude, expected rupture mechanism and so on. In our model, mapped active faults played a major role in defining the boundaries of the area sources; however we integrated geological data with historical and instrumental seismicity, as well as with geophysical data, including geodetic strain field, maximum horizontal stress (Sh_max) orientation and focal mechanisms, to define the geographical borders of the areas, the prevailing style of faulting, the upper and lower bounds of the characteristic seismogenic depth, and the distribution of hypocentral depths.

To determine the boundaries and the seismic parameters of the area sources, we collected and analysed several seismotectonic datasets (Fig. 2), some of which were compiled in the framework of MPS19 (Meletti et al., 2021) to be used as common inputs for the development of all the ERFs. Among these datasets, we used a historical earthquake catalogue (Catalogo Parametrico dei Terremoti Italiani, version 1.5, hereinafter CPTI15; Rovida et al., 2016, 2020); an instrumental earthquake catalogue (Gasperini et al., 2016; Lolli et al., 2020); version 3.2.1 of the Database of Individual Seismogenic Sources (DISS 3.2.1; Basili et al., 2008; DISS Working Group, 2018); a harmonized GPS velocity model for the Mediterranean area (Devoti et al., 2017); and other geological and geophysical data, available for specific regions and for the whole territory, as described in the following. Data used to draw boundaries of area sources of ZS16, earthquake catalogues and a brief description of the area sources are available on request from the corresponding author.

Figure 2 Main datasets (a, b, c) used to build the ZS16 (d) seismogenic zoning. (a) Earthquake epicentres from the CPTI15 catalogue (the green polygon represents the area covered by the CPTI15, as described by Rovida et al., 2016, 2020); (b) focal mechanisms of earthquakes with M_w≥4 (Pondrelli et al., 2020); (c) composite seismogenic sources from the DISS3.2.1 database (Basili et al., 2008; DISS Working Group, 2018); (d) seismogenic zoning ZS16 proposed in this study.

3.1 Earthquake catalogues

CPTI15 v1.5 lists 4389 earthquakes with moment magnitude M_w≥4.0 or macroseismic intensity ≥5 that occurred in the Italian and neighbouring areas (Fig. 2a) in the period 1000–2014. The catalogue provides epicentral locations and homogeneous M_w estimates derived from both macroseismic and instrumental data. The catalogue takes advantage of the wealth of macroseismic intensity data related to both historical and recent earthquakes collected in the Italian Macroseismic Database DBMI15 (version 1.5; Locati et al., 2016). The parameters of 43 % of the earthquakes are calculated from such intensity data with the Boxer algorithm (Gasperini et al., 1999, 2010) with an updated calibration (Rovida et al., 2020). Following Gasperini et al. (2012), instrumental magnitudes are moment tensor solutions complemented with proxy M_w obtained from magnitude estimates in other scales according to Gasperini et al. (2013) and Lolli et al. (2014, 2015, 2018). For the sake of homogeneity, the M_w of modern earthquakes with both macroseismic and instrumental magnitudes is the mean of the two estimates.

Although different methods for identifying mainshocks are available in the literature, within the MPS19 project the widely used and tested procedure by Gardner and Knopoff (1974) with the space and time windows defined therein was selected. The procedure resulted in a catalogue of 3353 mainshocks, corresponding to 76 % of the whole CPTI15, which was used in all the ERFs of MPS19.

To define the seismogenic layers and the depth distributions of the earthquakes in the area sources, we also considered an instrumental catalogue with homogeneous M_w determination for the period 1981–2015 that contains about 256 000 events without any lower-magnitude threshold (Gasperini et al., 2016). An updated and slightly different version of this catalogue was later published by Lolli et al. (2020), with a thorough description of the input data and selection criteria. The catalogue contains the instrumental part of CPTI15, complementing it for magnitude <4.0.

3.2 Focal mechanisms

To collect a representative dataset useful to define the styles of faulting of each area source, we started from the Italian CMT (centroid moment tensor) dataset (Pondrelli et al., 2006; Italian CMT dataset available at https://doi.org/10.13127/rcmt/italy). It includes all available moment tensors for the Italian Peninsula and surrounding areas from 1976 to present with a minimum M_w of 4.0, collected from the global CMT catalogue (Ekström et al., 2012, and references therein) and the European–Mediterranean RCMT (regional centroid moment tensor) catalogue (Pondrelli and Salimbeni, 2015; https://doi.org/10.13127/rcmt/euromed). To reach the best homogeneity in terms of spatial distribution, we included in our dataset the moment tensors obtained through seismic data inversion for a few M_w ≥ 4.0 earthquakes that occurred in the Alpine region, not available in the CMT or RCMT Catalogues, selected from the German Research Centre for Geosciences (GFZ) and ETH Zurich datasets (Saul et al., 2011, and Bernardi et al., 2004, respectively). In addition, to get a dataset with a longer time coverage, we considered first-motion polarity focal solutions for relevant events that occurred before the digital era, such as the 1968 Belice (Sicily) earthquakes. They have been selected from the EMMA database (database of Earthquake Mechanisms of the Mediterranean Area; Vannucci and Gasperini, 2004). When multiple focal mechanisms were available for a single event, the choice followed the quality evaluation given in the EMMA database, where a “preferred” solution is defined (see Pondrelli et al., 2020, for details). The final entire dataset, shown in Fig. 2b, consists of 995 focal mechanism solutions for events in the magnitude M_w range between 4 and 7.

3.3 Active faults

Active faults played an important role in defining the boundaries of the area sources; to this aim we consulted databases referring to a different scale of resolution, from the national to local scale, to include both the general seismotectonic picture and the details in the boundaries.

The DISS database (Basili et al., 2008; DISS Working Group, 2018) is a fundamental product for interpreting the relationships between faults and earthquakes in Italy. DISS 3.2.1 contains 127 individual seismogenic sources (defined as a simplified and three-dimensional representation of a fault plane; individual seismogenic sources are assumed to exhibit “characteristic” behaviour with respect to rupture length/width and expected magnitude), 188 composite seismogenic sources (defined as simplified and three-dimensional representations of crustal faults containing an unspecified number of seismogenic sources that cannot be singled out, Fig. 2c), 35 debated seismogenic sources and 3 subduction zones. All sources are based on geological/geophysical data and cover the whole Italian territory and portions of adjacent countries and seas.

At the national scale, we also considered the Structural Model of Italy (CNR, P. F. GEODINAMICA, 1990) and the seismotectonic model by Meletti et al. (2000). The latter was used as a guide for identifying homogeneous domains of active tectonics in Italy.

In some regions, we integrated the above datasets with data from local detailed geological–structural investigations to define the boundaries of the area sources, for example: Delacou et al. (2004) and Sue et al. (2007) for northwestern Italy; Collettini and Barchi (2002), Boncio et al. (2004), Papanikolaou and Roberts (2007), Lavecchia et al. (2007a), Faure Walker et al. (2010, 2012), Visini (2012), Tesson et al. (2016), and Valentini et al. (2017) for central and southern Italy; and Lavecchia et al. (2007b), Catalano et al. (2010), Billi et al. (2010), Visini et al. (2010) and Mastrolembo et al. (2014) for Sicily.

3.4 Other geophysical data

As a proxy for evaluating the thickness of the crust and defining zones with similar seismogenic thickness, we used the Moho maps by Solarino and Cassinis (2007) and Di Stefano et al. (2011) and the heat flow maps by Della Vedova et al. (2001).

We also considered the regional strain rate fields for the Mediterranean area derived from GPS data (Devoti et al., 2017) and the maximum horizontal stress Sh_max orientation (Mariucci and Montone, 2020) to qualitatively check the homogeneity of the strain rate values and of the Sh_max orientations within the area sources.

4 The MA4 seismogenic model

MA4 is based on seismogenic zoning ZS16, which represents the evolution of previous area-source models proposed in the last 30 years as unique area-source input in seismic hazard assessment performed in Italy. The criteria for defining ZS16 are described in Sect. 4.1; methodologies for the calculation of parameters of the area source, useful for PSHAs, are described in Sect. 4.2 (top and bottom depth of the seismogenic layer), Sect. 4.3 (hypocentre distribution) and Sect. 4.4 (style of faulting). MA4 also models seismicity rates for each area source of ZS16 (Sect. 4.4) and includes epistemic uncertainties in the assessment of the completeness intervals of the earthquake catalogue (Sect. 4.4.1), maximum magnitude (Sect. 4.4.2) and alternative frequency–magnitude distributions (Sect. 4.4.3).

4.1 The ZS16 seismotectonic zoning

Although area sources are widely used for national and international PSHA, there are no standard objective approaches for defining their boundaries. We acknowledge the criteria defined in previous studies (e.g. Giardini, 1999; Meletti et al., 2008; Wiemer et al., 2009; Vilanova et al., 2014; Danciu et al., 2018) to set guidelines for the delineation of area sources in order to describe the correlation between active faults, earthquakes and other geophysical inputs. To update the existing reference national zoning scheme ZS9, we applied the following criteria.

• a

  Start from the area sources of the ZS9 model.

• b

  Be consistent with the general background delineated by the geodynamic model proposed by Meletti et al. (2000); i.e. an area source should belong to a unique tectonic zone (active shallow crustal, volcanic or subduction zone in the specific Italian case).

• c

  Incorporate all recent advances in the understanding of the active tectonics of the territory and in the distribution of seismogenic sources modelled in the DISS 3.2.1 database and other active fault compilations at the national and regional scale (see Sect. 3.3). In particular, for defining area-source boundaries that primarily follow the surface projection of mapped active faults, an area source should not interrupt a normal or reverse fault system unless major differences are observed (changes in stress orientation and/or changes in crustal depth); for strike-slip faults, boundaries should be parallel to the strike of the faults, and the area source should contain the faults.

• d

  Incorporate information derived from the investigation of the most recent seismic sequences that struck Italy after the compilation of ZS9, namely the 2009 L'Aquila, 2012 Emilia and 2016 Amatrice–Norcia sequences.

• e

  Be consistent with the spatial pattern of seismicity depicted by the CPTI15 earthquake catalogue. Area sources whose borders were drawn using mapped active faults (point c) should not cross spatial clusters of earthquakes which are attributed to the same faults. As earthquake locations, both macroseismic and instrumental, are affected by uncertainties, a spatial shift between the possible causative faults and the epicentre can occur.

• f

  Consider for the definition of the boundaries the pattern of seismicity, focal mechanisms, geodetic strain field, Sh_max and heat flow data.

• g

  Account for the variation in the style of faulting and tectonic regime with depth; therefore multiple area sources can overlap on the volume domain.

• h

  Cover the entire Italian territory, as required by MPS19.

Applying these criteria to the data described in Sect. 3, we defined the seismogenic zoning shown in Fig. 2d, consisting of 48 active shallow crustal area sources and 2 area sources corresponding to the Campanian and Mt Etna volcanic districts (i.e. area source nos. 31 and 49, respectively). For the deep seismicity related to the Tyrrhenian subduction intraslab, no area sources were defined, because MPS19 adopted a separate ad hoc ERF for modelling such seismicity. Finally, it is worth noting that three area sources (nos. 19, 20 and 25) showed a different kinematics for shallow and deep seismicity, so they have been split, and computations for them have been done apart (Supplement 2; see Pondrelli et al., 2020, for details).

4.2 Top and bottom depth of the seismogenic layer and hypocentral distributions

Seismogenic depths for each area source of ZS16 were estimated using the instrumental catalogue by Gasperini et al. (2016). In particular, we assumed the upper and lower limits of the seismogenic layer as corresponds to the 5th and 95th percentiles of the depth distribution of the earthquakes inside each area (e.g. Boncio et al., 2009; Stucchi et al., 2011), and we modelled the depth distributions of hypocentres with the peaks of unimodal and bimodal distributions that best approximate the observed values.

To estimate these values, we first removed the earthquakes with fixed hypocentral depth (i.e. 0, 5 or 10 km), which represents ∼10 % of the total. We then considered only the earthquakes whose depth is shallower than the Moho, based on the crustal models by Solarino and Cassinis (2007) and Di Stefano et al. (2011). We calculated the (rounded) 5th and 95th percentiles of the hypocentral-depth distributions, assumed as the upper and lower limits of the seismogenic layer, respectively. As regards area source nos. 19, 20 and 25, showing a different kinematics for shallow and deep seismicity (Pondrelli et al., 2020), we divided them into “shallow” and “deep” area sources (i.e. nos. 19s and 19d and 20s and 20d, i.e. nos. 25s and 25d). Based on crustal thickness and rheological properties, we adopted ad hoc values of depth for these area sources: for shallow area sources, we assumed depth values of 0, 15 and 10 km, respectively, for the top and the bottom of the seismogenic layer and for the hypocentral depth; for deep area sources, we assumed depth values of 15, 30 and 23 km, respectively, for the top and the bottom of the seismogenic layer and for the hypocentral depth.

For depths within the top and the bottom of the seismogenic layer, we computed the modal values, standard deviation and log likelihood of the unimodal and bimodal distributions that best fit the observed values (see Fig. 3 for an example). We evaluated and compared the AIC (Akaike information criterion) index of the unimodal and bimodal distributions to select the best model for the hypocentral-depth distribution of each area. In the case of a unimodal distribution, we used the modal value as representative of the hypocentral depth, while for bimodal distributions, we assigned weights to both modal values by using their mixing proportion value in the bimodal distribution. To evaluate the stability of the results with respect to the number and the magnitude of the considered events, we calculated the upper and lower seismogenic depths and the modal values of the distributions for different minimum magnitudes (from M_w 2.5 to 4.5) and compared the resulting depth estimates with the depth of the composite seismogenic sources of DISS 3.2.1 inside the area. In Fig. 3, we show an example of the results obtained for area source no. 24. We retain that the modal values can be representative of the hypocentral depths; however, future research could better detail correlations among depth distributions, magnitude and kinematics.

Figure 3 Hypocentral-depth distributions (grey bars) for different threshold magnitudes (reported on top of each panel along with the number of considered earthquakes) for area source no. 24. Black lines correspond to the 5th and 95th percentiles, assumed to be the upper and lower seismogenic depths (round value); blue curve and line represent the unimodal distribution and its modal value; red curves and lines represent the bimodal distribution and its two modal values. Solid lines indicate the best model between uni- and bimodal distributions; dashed lines indicate the other model. The right panel shows the depth ranges of the composite seismogenic sources (CSSs) of DISS 3.2.1 inside the area. The depth of the Moho, from Solarini and Cassini (2007) and Di Stefano et al. (2011), is approximately 30–35 km.

In the Mt Etna region, we assigned earthquakes with hypocentral depth <10 km to the volcanic domain (area source no. 49 in Fig. 2d) and earthquakes with hypocentral depth ≥10 km to the underlying active crustal area sources (nos. 44, 45 and 46). For the Campanian volcanic area (no. 31 in Fig. 2d), we adopted depth values of 0, 4 and 1 km, respectively, for the top and the bottom of the seismogenic layer and for the hypocentral depth, mainly based on the parameters of the 2017 M_w 3.9 Ischia earthquake (De Novellis et al., 2018). Supplement 1 lists the parameters derived for all the area sources for minimum magnitude M_w of 2.5, which appears to be the most appropriate threshold value to ensure a significant number of earthquakes for all the areas.

4.3 Style of faulting

Pondrelli et al. (2020) defined the criteria to parametrize the styles of faulting of expected earthquake ruptures and to evaluate their representativeness in each area source. Using available seismic moment tensors for relevant events (M_w≥4.5), first-motion focal mechanisms for less recent earthquakes and also geological data on past activated faults, we collected nearly a thousand data points for seismic events that occurred in the last ∼100 years in the Italian Peninsula and surrounding regions, as described in Sect. 2.2. On this dataset we applied in each seismic area source a procedure that starts with the separation of all available focal mechanisms into the three main tectonic styles, following the rake-based criteria given in Akkar et al. (2014): a rake between − 135 and − 45^∘ defines a normal solution, while one between 45 and 135^∘ is a reverse solution, and all the rest are classified as a strike-slip one. We summed all data within each group, i.e. the normal, the reverse and the strike-slip one. Using the results of the summations, in each area source we identified, when possible, a nodal plane, considering the different percentages of styles of faulting and, where necessary, including a total or a partial (even in terms of tectonic style) random source contribution. Following these steps, we obtained the resulting styles of faulting for each area source reported in Fig. 4 and also in Supplement 2 (Pondrelli et al., 2020). As stated above, in a few cases, changes in tectonic style with depth were identified (area source nos. 19s and 19d, 20s and 20d, and 25s and 25d).

Figure 4 Expected style of faulting for each area source (modified from Pondrelli et al., 2020). Plain circles represent random seismic sources. White circles represent 100 % random, while black, purple and yellow circles correspond to reverse, normal and strike-slip random sources, respectively. Colours for cumulative focal mechanisms follow the same criteria. Focal mechanisms with a grey background or plain circles with darker colours are the sources for deeper layers. Black numbers are the percentages of contribution to the final sources when their sum is the expected style of faulting.

4.4 Annual rates of earthquake occurrences

To estimate the expected seismicity rates of each area source, we adopted a time-independent (i.e. Poisson) model for earthquake occurrence. We assumed that the distribution of the earthquake magnitudes follows the truncated Gutenberg–Richter (TruncGR; Ordaz, 2004) model that has three parameters: Λ₀, which is the cumulative number of earthquakes per unit time equal to or larger than the magnitude threshold or minimum magnitude (M_min) and smaller than the upper (or maximum) magnitude (M_u) and the slope (β, $\mathit{\beta }=\mathrm{2}/\mathrm{3}b$). The TruncGR distribution is the Pareto distribution with the probability density function truncated at both ends. Its cumulative density function related to moment magnitude is

$$\begin{array}{}\text{(1)}& {\mathrm{\Lambda }}_{\left(M\right)}={\mathrm{\Lambda }}_{\mathrm{0}}{\displaystyle \frac{e^{-\mathit{\beta }{M}_{\min }}-e^{-\mathit{\beta }M}}{e^{-\mathit{\beta }{M}_{\min }}-e^{-\mathit{\beta }{M}_{\mathrm{u}}}}},\end{array}$$

where β and Λ₀ were derived from the declustered CPTI15 catalogue (see Sect. 3.1) by adopting the completeness time intervals described in the following (Sect. 4.4.1) and applying a maximum-likelihood fit based on Weichert (1980). M_u, the upper magnitude, is described in the Sect. 4.4.2. The minimum magnitude for the application of Weichert (1980), M_min, is estimated following the maximum curvature approach.

4.4.1 Completeness time intervals

Two independent sets of completeness time intervals for the CPTI15 catalogue were defined according to (i) the historical approach of Stucchi et al. (2004, 2011) and (ii) the statistical method proposed by Albarello et al. (2001).

The historical approach determines the complete intervals analysing the local history of a set of sample localities. Based on this knowledge, the years from which it is unlikely that earthquake effects of a given intensity are not recorded in the local historical sources were determined. The catalogue can be considered complete for earthquakes of the same epicentral intensity (I₀) located at or near the analysed locality. Stucchi et al. (2004, 2011) assessed the starting year of completeness for intensity ≥6 at 18 sample localities and then extrapolated them to the area sources they belong to and to others with similar history and seismotectonic features. As this approach is independent from the catalogue, we used the same completeness intervals of Stucchi et al. (2011) for I₀≥6, adapting the five macroregions defined therein to ZS16. In these macroregions we also evaluated the completeness intervals for I₀ 4–5, 5 and 5–6 with the same criteria, and we assessed all those of a newly introduced offshore region (“Sea” in Fig. 5a). The completeness intervals determined in this way were then applied to M_w bins defined according to the new empirical conversion relation between the epicentral intensity and magnitude of CPTI15 (Rovida et al., 2020). The estimated M_w bins are of 0.23 M_w units, which corresponds to the difference between the M_w values obtained from the discrete I₀ values, including also uncertain intensity assessments as “half degrees”, e.g. 6–7 equal to 6.5. This choice avoids the uneven concentrations in the same bin of M_w values derived from the conversion of different discrete epicentral intensity values of historical earthquakes.

Figure 5 (a) Macroregions adopted for evaluating the completeness time intervals for the earthquake catalogue; (b) tectonic domains (blue polygons) used to calculate β overlapped to the ZS16 area sources (dashed black polygons); (c) map of the values of $M_{{\mathrm{w}}_{\max \mathrm{1}}}$, $M_{{\mathrm{w}}_{\max \mathrm{2}}}=M_{{\mathrm{w}}_{\max \mathrm{1}}}+\mathrm{0.3}$.

The statistical completeness intervals were assessed using the procedure of Albarello et al. (2001) for the same macroregions defined for the historical approach. The method is very sensitive to the number of earthquakes considered in each area and magnitude bin, and to ensure the stability of the results, we selected magnitude bins of 0.46 M_w units, which means grouping in the same class integer and intermediate intensity values (e.g. intensity 6–7 together with 6). The results were then applied to the same bins of 0.23 M_w units defined for evaluating the historical completeness. The completeness starting years for each I₀ $/$ M_w bin and macroregion defined according to the two approaches are shown in Supplement 4. Figure 6 shows a comparison, in each macroregion, of the completeness periods and magnitudes assessed with the historical and statistical approaches. The complete catalogues resulting from the historical and the statistical assessments are made of 2496 and 2603 earthquakes, i.e. 63 % and 66 % of CPTI15 (excluding earthquakes of the Mt Etna volcanic area and subduction events of the Calabrian Arc), respectively.

Figure 6 Plot of time vs. completeness magnitude, defined according to both the historical and statistical approach, in the six macroareas shown in Fig. 5a. Grey circles represent the earthquakes in CPTI15. Hist: historical approach, Stat: statistical approach.

The two methods, although with differences due to their assumptions, provide comparable results within each macroarea. The difference in the number of events in the complete periods arises from the intervals assessed with the statistical approach for low-magnitude bins, which contain more events than the highest ones, which are longer than those obtained with the historical approach. On the contrary, the historical approach determines longer complete periods for the highest-M_w bins, which are less rich in events (Fig. 6).

In conclusion, also taking into account the declustering procedure mentioned above (Sect. 3.1), the catalogues used for calculating annual rates of earthquake occurrences within the MA4 model contain 1800 and 1888 events, obtained, respectively, with the historical and the statistical completeness assessment.

4.4.2 Maximum magnitude

For the definition of the maximum magnitude we used the estimates provided for MPS19, described in Visini et al. (2021), which were based on the estimate of the maximum observed earthquake in the earthquake record from CPTI15. The Italian area was divided into 18 tectonic domains (Fig. 5b), and the earthquakes listed in CPTI15 were assigned to them, according to their location. Based on the average error in magnitude estimates for earthquakes occurring before and after 1980, a minimum value of uncertainty for the M_w evaluation of 0.3 and 0.2 was introduced for the historical and instrumental portion of the catalogue, respectively. The maximum observed magnitude inside a tectonic domain is the largest magnitude observed, including the uncertainty ($M_{{\mathrm{w}}_{\mathrm{obs}}}+$ uncertainty). In low-seismicity regions (for example, those located along the Tyrrhenian coast and in the north of Italy), albeit with a relative abundance of low-magnitude earthquakes, the earthquake record only spans over a few hundred years, and this time interval may be much smaller than the recurrence time of moderate earthquakes. As a result, the maximum observed magnitude in these regions relies on a poor historical record not useful to constrain the upper magnitude. In these situations, Visini et al. (2021) adopted analogies to similar tectonic features (Wheeler, 2009) or cautionary approaches based on expert judgement of minimum values of maximum M_w. Then, following Woessner et al. (2015), the minimum value of maximum magnitude $M_{{\mathrm{w}}_{\max }}$ was assigned to each tectonic domain ($M_{{\mathrm{w}}_{\mathrm{tect}}}$), namely 6.5 for all active crustal areas, 6.0 for the Tyrrhenian tectonic domain and 5.6 for the volcanic area of Mt Etna. Two values of maximum magnitude were then assigned to each tectonic domain: (i) $M_{{\mathrm{w}}_{\max \mathrm{1}}}$ (Fig. 5c), which is the largest value between $M_{{\mathrm{w}}_{\mathrm{tect}}}$ and $M_{{\mathrm{w}}_{\mathrm{obs}}}+$ uncertainty, and (ii) $M_{{\mathrm{w}}_{\max \mathrm{2}}}$, which results by uniformly incrementing $M_{{\mathrm{w}}_{\max \mathrm{1}}}$ by a cautionary value of 0.3 to account for epistemic uncertainties, except for the Mt Etna volcanic domain, where $M_{{\mathrm{w}}_{\max \mathrm{2}}}$ is equal to $M_{{\mathrm{w}}_{\max \mathrm{1}}}$. Furthermore, where active faults are known, Visini et al. (2021) compared the magnitude estimates of the composite seismogenic source (CSS) by DISS 3.2.1 with the $M_{{\mathrm{w}}_{\max \mathrm{1}}}$ and $M_{{\mathrm{w}}_{\max \mathrm{2}}}$ from the earthquake catalogue in each area source. As already pointed out in Visini et al. (2021), caution should be taken when comparing the earthquake-catalogue-based and the DISS-based magnitudes. The $M_{{\mathrm{w}}_{\max }}$ of the CSS (namely $M_{{\mathrm{w}}_{\mathrm{css}}}$) is computed based on empirical regressions, from the literature or by associating the largest historical earthquake of the area to the CSS. In the latter case, however, there is not necessarily a correspondence between the magnitude associated with a CSS in DISS 3.2.1 and the one provided by CPTI15, as a result of the different and independent development of the two databases. For two tectonic domains $M_{{\mathrm{w}}_{\mathrm{css}}}$ exceeds the $M_{{\mathrm{w}}_{\max \mathrm{2}}}$. For tectonic domain 2 (Fig. 5b), $M_{{\mathrm{w}}_{\mathrm{css}}}$ was evaluated using the M_w of the 1690 earthquake as derived from Guidoboni et al. (2007, 2019), whereas the M_w of the same earthquake in CPTI15, including uncertainty, is ∼6.5, and it is located in tectonic domain 3. In this case, Visini et al. (2021) decided to anchor $M_{{\mathrm{w}}_{\max \mathrm{1}}}$ to $M_{{\mathrm{w}}_{\mathrm{tect}}}$, considering that M_w of 6.5 corresponds to the magnitude estimated in the most recent CPTI15 catalogue. In tectonic domain 12, $M_{{\mathrm{w}}_{\mathrm{css}}}$ is the M_w 6.8 2003 earthquake, based on the modelling from Meghraoui et al. (2004) of the coseismic uplift. Because $M_{{\mathrm{w}}_{\max \mathrm{2}}}$ agrees with the original estimate of the magnitude of the 2003 earthquake, Visini et al. (2021) kept the $M_{{\mathrm{w}}_{\max \mathrm{2}}}$ (see Visini et al., 2021, for details).

$M_{{\mathrm{w}}_{\max \mathrm{1}}}$ and $M_{{\mathrm{w}}_{\max \mathrm{2}}}$ are used as upper bounds (M_u) of the earthquake recurrence of various magnitudes specific to each seismogenic source.

4.4.3 Seismic-rate determination

To calculate annual rates of earthquake occurrences for the active shallow crustal areas, we first excluded from the CPTI15 the events belonging to the southern Tyrrhenian subduction (i.e. those located in that area with hypocentral depth larger than 40 km), and we imposed a minimum of 10 earthquakes and at least 2 non-empty classes of magnitude (0.1 bin size) in each area source to derive stable β values; otherwise, we assumed β to be equal to 2.3 (b=1).

We used five different approaches to calculate the seismicity rates for the area sources, which are based on two different assumptions.

The first assumption is that β varies across the areas; then a first approach (named approach i) consists in the classical estimation of the parameters Λ₀ and β (see Eq. 1) directly for the area sources, from the earthquakes belonging to each source with magnitudes between M_min and the maximum magnitude ($M_{{\mathrm{w}}_{\max \mathrm{01}}}$ or $M_{{\mathrm{w}}_{\max \mathrm{02}}}$; “statistical” or “historical” completeness according to the branch of the logic tree).

The second assumption is that β is stable over groups of area sources with similar seismotectonic features. The rationale behind this assumption is to use a sample of earthquakes as rich as possible to estimate a robust β value while maintaining a relationship between the slope of TruncGR and the seismotectonic characteristics of a region (for example area sources characterized by extension). Indeed, for area sources characterized by a “low” rate of seismic activity, the β (or the b) value can be biased and results as an artefact of the low number of data available (e.g. Marzocchi et al., 2020). To maintain β variations from area to area, however, we used the approach i (described above), in which the β value is actually estimated for each area source. Under this second assumption, we first calculated β and Λ₀ for groups of area sources, hereinafter defined as a “macroarea”, using earthquakes belonging to each macroarea with magnitudes between M_min and the maximum magnitude (as for the approach i, $M_{{\mathrm{w}}_{\max \mathrm{01}}}$ or $M_{{\mathrm{w}}_{\max \mathrm{02}}}$; statistical or historical completeness are taken according to the branch of the logic tree). To this purpose, we used the tectonic regions shown in Fig. 5b. As the next step, we assessed the recurrence parameters in each area source within a macroarea. Keeping with our objective to only change Λ₀ according to the level of seismic activity of an area, we defined four different approaches for partitioning Λ₀ to the sources belonging to the same macroarea. In particular, defining Λ_0as and Λ_0ma, respectively, as the values of Λ₀ of the area source and of the macroarea, the four approaches are the following.

• ii

  In each macroarea, we compared the observed seismicity rates of occurrence of the macroarea with the ones of the area sources. We compared seismicity rates of magnitudes ranging between M_min (estimated for the macroarea) and M_u. For each area source and for each magnitude, we calculated the ratio between the observed seismicity rates of occurrences of the macroarea and of the area source. For each area source, we calculated the average of these ratios. Finally, we used the average ratios of all the area sources of a macroarea to proportionally scale Λ_0ma to Λ_0as.

• iii

  To build the TruncGR of each area source in a macroarea, we used β and M_u of the macroarea, whereas Λ_0as is firstly assumed to be the observed rate of exceedances for magnitudes equal to or larger than M_min proper of the area source. Then, we scaled the expected rates of the TruncGR of each area source in a macroarea to ensure they sum up to the TruncGR of the macroarea they belong to.

• iv

  For each area source in a macroarea, the TruncGR uses β and M_u of the macroarea. The Λ_0as is computed at the M_min of the macroarea, and it is the value that maximizes the log likelihood of the observed and expected (forecasted) number of earthquakes in each area source. In particular, for each magnitude of an area source, we computed the expected number of earthquakes by multiplying the seismicity rate of occurrence (not in the cumulative form) with the corresponding interval of completeness. Then, for each area, we sum the log likelihood obtained for each magnitude and find for the Λ_0as that allows for obtaining the best log likelihood.

• v

  For each area source in a macroarea, the TruncGR was computed using β and M_u of the macroarea. As regards Λ_0as, for each area source, we used the observed rates of exceedance for magnitudes equal to or larger than the M_min proper of each area source and computed Λ_0as at the M_min as the value that minimizes the root mean square of the observed vs. expected seismicity rates (given by the TruncGR).

For area nos. 19, 20 and 25, we successively partitioned the five Λ_0as values between the shallow and deep sources, according to their relative percentage of the number of earthquakes with M_w≥2.5 (values are in Supplement 3).

Figure 7 shows an example of the frequency–magnitude distributions calculated directly for area source nos. 24 and 36 and for macroarea nos. 6 and 11, respectively. Whereas in areas of a relatively high rate of seismicity (as for the extensional areas in the central Apennines, area source no. 24), evaluations of β values are approximately stable with both assumptions (between approach i and approaches ii to v), in an area characterized by a low rate, the differences in the estimates of β values can lead to variations in the expected rates modelled by the TruncGR, as for area source no. 36 in the strike-slip regions in the Apulian foreland. As we cannot know a priori what is the best solution between the β-value variable at the level of the single areas and at the level of large tectonic domains, we embedded the five approaches in a logic-tree approach (described in the next section and shown in Fig. 8). In Supplement 3, we listed the parameters of the TruncGR distributions associated with each area source.

Figure 7 Example of the frequency–magnitude distributions calculated for area source (AS) nos. 24 and 36. The frequency–magnitude distributions of the macroarea (MA) they belong to are also shown. The area source and the macroarea are shown in Fig. 5b.

5 Seismic hazard estimation using MA4

Seismic hazard was calculated over the whole Italian territory (including sites located within 5 km outside the borders), adopting the MA4 seismogenic model. For this purpose, 52 area sources were used, considering area source nos. 19, 20 and 25 in the form shallow and deep and discarding area source no. 49 (Mt Etna). In fact, MPS19 defined additional ad hoc ERFs for three specific regions: (a) the Mt Etna volcanic area, which replaces the ERF for area source no. 49; (b) the subduction shallow interface seismicity and deep intraslab 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 catalogue (see Meletti et al., 2021; Visini et al., 2021).

Alternative choices and interpretations about the key elements were embedded in a logic-tree structure (Kulkarni et al., 1984; Coppersmith and Youngs, 1986; Senior Seismic Hazard Analysis Committee, 1997), which is the conventional tool to capture the epistemic uncertainty associated with the input elements of a PSHA model. The adopted logic tree, shown in Fig. 8, consists of a first branching level accounting for the two alternative evaluations of the catalogue completeness time intervals described in Sect. 4.4.1, i.e. one based on historical information (Stucchi et al., 2004, 2011) and one on the statistical approach by Albarello et al. (2001). Then, a second branching level considers the two alternative sets of maximum magnitude M_u (Eq. 1) values (i.e. $M_{{\mathrm{w}}_{\max \mathrm{1}}}$, $M_{{\mathrm{w}}_{\max \mathrm{2}}}$), described in Sect. 4.4.2, and a third level accounts for the five approaches adopted for calculating the frequency–magnitude distributions of each area source (see Sect. 4.4.3). Since we have no specific reason to prefer or to differently weight the 20 resulting ERF branches, the same weight (i.e. $\mathrm{1}/\mathrm{20}$) was assigned to each of them. Following Woessner et al. (2015) and Danciu et al. (2018), we considered the uncertainties on hypocentral depths and focal mechanisms as aleatory, and we used distributions and weights described in the previous sections to model them. Two further branching levels account for alternative choices of GMMs, as selected and applied in MPS19, which are three GMMs for active shallow crustal regions (Bindi et al., 2011, 2014; Cauzzi et al., 2015, with associated weights equal to 0.45, 0.32 and 0.23, respectively) and one model for volcanic areas (Lanzano and Luzi, 2020). The selected GMMs provide estimates of ground shaking in terms of the geometric mean of the horizontal components. In particular, the 3 models adopted for shallow crustal regions derive from a pre-selection of 16 candidate GMMs performed over nearly 1000 models published in the literature. Then the performance of the pre-selected GMMs was evaluated through the comparison with accelerometric records available for Italy and Europe using different scoring methods; a final rank of the 16 candidate GMMs was done, and the 3 best performing models were selected; finally, the weights of the 3 GMMs were assigned, combining the results of the scoring and the weights obtained through an experts' elicitation (see Lanzano et al., 2020, for details). Besides showing the best relative hindcasting skill, the three selected GMMs use different metrics for the distance (i.e. Joyner–Boore distance R_JB for Bindi et al., 2011; hypocentral distance R_hyp for Bindi et al., 2014; and distance from the rupture plane R_rup for Cauzzi et al., 2015) and are calibrated on different datasets (i.e. Italian, European and global, respectively). As a result, we obtained a final logic tree made of 60 branches (Fig. 7).

Figure 8 The logic-tree scheme adopted in this study. ZS16 is the seismogenic zoning. The completeness time intervals for the CPTI15 catalogue were defined according to the historical approach of Stucchi et al. (2004, 2011) and the statistical method of Albarello et al. (2001). M_u $M_{{\mathrm{w}}_{\max \mathrm{1}}}$ and M_u $M_{{\mathrm{w}}_{\max \mathrm{2}}}$ are the two sets of values adopted for the maximum magnitude, described in Visini et al. (2021). Lowercase Roman numerals i to v identify the five approaches used to calculate the annual seismic rates (see Sect. 4.3). The last two nodes concern the GMMs used for active shallow crustal regions and volcanic areas.

Seismic hazard was calculated for rock-site conditions (i.e. EC8 site category A or Vs₃₀≥800 m s⁻¹) using the OpenQuake engine platform (Pagani et al., 2014). We recall here that, for the computation of seismic hazard, OpenQuake discretizes every area source into a regular grid of points, each representing the longitude and latitude of the centre of a rupture. Every rupture has a rectangular shape and is centred on the hypocentral distribution parameterized in ZS16. Ruptures are created using the magnitude–area scaling relationship of Wells and Coppersmith (1994) and by conserving an aspect ratio of 1. The aspect ratio is adjusted if the rupture plane would grow beyond the upper and lower seismogenic depths, in order to keep the area.

Figure 9 illustrates the spatial distribution of mean peak ground acceleration (PGA) values obtained by applying the weighting scheme in Fig. 8, for 10 % and 2 % probabilities of exceedance in 50 years, also referred to as 475- and 2475-year return periods (RPs), respectively. The map for PGA at 10 % probability of exceedance in 50 years (Fig. 9a) shows the highest hazard estimates (PGA ≥0.25 g) in three areas, i.e. in the northeast of Italy, along the central Apennines and in the southern Apennines. Hazard levels ranging between 0.1 and 0.25 g characterize the majority of coastal areas of the Adriatic and Tyrrhenian seas, the southern part of the Po Plain, Sicily, and large parts of southern Italy. The northwest of Italy, Sardinia and part of Apulia show hazard values generally lower than 0.075 g.

Figure 9 Maps of mean values of PGA at 10 % (a) and 2 % (b) probability of exceedance in 50 years. Locations of the three cities selected for detailed analyses are shown.

Figure 10 Hazard curves for PGA for the three selected cities. The curves represent the mean hazard level (black line), the hazard resulting from each of the 60 branches (realizations, grey lines), and the uncertainties expressed through the 16th and 84th percentiles (pct, red lines). The legend in the centre panel refers to all panels.

Figure 11 UHS for 10 % (lower spectra, bright colours) and 2 % (upper spectra, pale colours) probability of exceedance in 50 years for the three selected cities. The mean spectra and the 16th and 84th percentiles are reported. The legend in the first panel refers to all panels.

The map for PGA at 2 % probability of exceedance in 50 years (Fig. 9b) shows that PGA ≥0.5 g characterizes almost entirely the Apennines, part of the southern Po Plain and the northeast of Italy. The southeast of Sicily shows hazard levels ranging between 0.3 and 0.5 g, whereas PGA ranging between 0.2 and 0.3 g characterizes the rest of Sicily and the majority of the Adriatic and Tyrrhenian coastal areas. Hazard levels lower than 0.1 g are obtained in the northwest of Italy, in the southern part of Apulia and in Sardinia.

Figures 10 and 11, respectively, show hazard curves for PGA and uniform hazard spectra (UHS) for spectral periods from 0.1 to 4 s for the cities of Milan, L'Aquila and Syracuse (see locations in Fig. 9), chosen for exemplificative purposes. Figure 10 illustrates the variability in the expected ground motions in PGA, showing the mean hazard level (black line), the hazard curves resulting from each of the 60 realizations (grey lines) and the uncertainties expressed through the 16th and 84th percentiles (red lines). Figure 11 show the mean and the 16th and 84th percentiles of UHS for 10 % and 2 % probability of exceedance in 50 years.

Figure 12 shows the spatial distribution of the coefficient of variation (CoV) of PGA values for 10 % and 2 % probabilities of exceedance in 50 years. Recalling that GMMs were weighted and the 60 branches do not have the same weights, we calculated the CoV as the weighted standard deviation divided by the weighted mean. CoV <0.2 covers a large part of the Italian territory, indicating a low uncertainty in the expected acceleration. CoV >0.2 characterizes the areas with the lowest PGA and the southeastern sector of Sicily. The latter corresponds to the area sources where the β values show the largest differences between the macroarea and the single area-source approaches.

Figure 12 Spatial distribution of the coefficient of variation (CoV) of PGA values for 10 % (a) and 2 % (b) probabilities of exceedance in 50 years.

The CoV values in Fig. 12 contain uncertainty related to both the ERF and the GMMs. We then investigated the relative contribution of both components of epistemic uncertainty at three selected localities. For each site, hazard curves for PGA were plotted by distinguishing the 3 groups composed of 20 realizations that use the same GMM for active shallow crustal regions. In Fig. 13, the three groups are shown by light-red lines for the realizations that adopt Bindi et al. (2011), and light-green lines and light-blue lines are for those using Bindi et al. (2014) and Cauzzi et al. (2015), respectively. Figure 13 also shows the mean probability of exceedance (POE) in 50 years for the three groups. It can be seen that most of the seismic hazard curves from the different GMMs are overlapping, and the uncertainty due to the GMMs is of a similar order of magnitude as the uncertainty related to the ERF. For Milan and L'Aquila, the mean curves of the three groups overlap in a wide range of POE, especially for POE <10 % in 50 years. For Syracuse the Bindi et al. (2011) group intersects the Cauzzi et al. (2015) one at POE of ∼10 %–11 % in 50 years; for PGA <0.02 g, the realizations returning the highest POE are those using the latter GMM, and the realizations returning the lowest POE are those using Bindi et al. (2011). This trend is inverted for PGA >0.05 g or POE <∼12 % in 50 years.

Figure 13 Hazard curves for PGA for the three selected cities. Panels (a), (b) and (c) show the mean hazard curves (black lines) and the hazard curve of each of the 60 branches (realizations, rlzs, coloured per GMM with pale colours). The legend in panel (b) refers to all panels. Panels (d), (e) and (f) show the mean hazard curves obtained as a weighted average of the realizations per each GMM (coloured per GMM with bright colours). The legend in panel (e) refers to all panels. Key for all panels is as follows. BeA11: Bindi et al. (2011), BeA14: Bindi et al. (2014), CeA15: Cauzzi et al. (2015).

To focus on the relative contributions of ERF and GMM uncertainties, in Fig. 14 we illustrate the POE for a series of levels of PGA, by distinguishing the branches that used a particular GMM. In particular, in Fig. 14, the POE for each level of PGA shown is normalized to the maximum, in order to have all values ranging from values >0 and ≤1, for ease of visualization. Normalization does not affect the next considerations, as we are not interested in describing absolute POEs but relative contributions. Figure 14 allows, in a qualitative way, for describing the relative importance of the uncertainty due to ERF modelling in respect to the GMM contribution. We plotted branches that use the same GMM with the same colour: red for those branches that used Bindi et al. (2011), green for Bindi et al. (2014) and blue for Cauzzi et al. (2015). If the uncertainty due to GMMs impacted more than the ERF modelling, we would expect a cluster of values grouped by the same colour (GMM). This is not the case shown in Fig. 13; on the contrary, the symbols of different GMMs are distributed to cover, approximately, the same ranges. The overlaps of values occur at all the PGA values shown and for the three localities; we interpret this as an indication of the importance of ERF modelling that can contribute in an equal, if not even larger, manner than the GMMs.

Figure 14 Relative contributions of epistemic uncertainty in ERF and GMM hazard estimates. Branches that use the same GMM are shown with three colours, red, green and blue. The x axis shows the POEs normalized to their maximum values for each PGA level. The legend in the centre panel refers to all panels. Key for all panels is as follows. BeA11: Bindi et al. (2011), BeA14: Bindi et al. (2014), CeA15: Cauzzi et al. (2015).

6 Discussion

The MA4 seismogenic model in Italy represents a challenge for the construction of an area-source-based model that takes into account the variety of seismotectonic environments, which include spatial and depth variations in the main style of faulting. The MA4 model is based on a seismotectonic zoning (ZS16), defined according to a list of criteria specifically defined to this purpose. However, we cannot exclude the fact that the proposed zoning still contains some a priori bias or, simply, some area sources could not reflect the actual tectonics. Objective criteria to delineate area sources with a quality ranking of the basic data would be an additional step (e.g. Wiemer et al., 2009; Vilanova et al., 2014).

As an example, the debate on “large” vs. “small” areas concerns subjective choices. Small areas are designed to capture changes in seismicity at the local scale (e.g. <20 km), but, in our opinion, these changes are better highlighted with fault-based models or a smoothed seismicity approach. Our idea of zoning aims at individuating the areas that have a homogeneous behaviour from a seismotectonic point of view, focusing on the similarities rather than on the differences. The experience in seismogenic zoning in Italy since the early 1990s (e.g. the ZS4 model by Meletti et al., 2000) suggests that small area sources, such as those representing single faults, can be appreciated by structural geologists but make defining their earthquake recurrence parameters very difficult because of the lack of available data that can even produce apparent differences in seismicity distribution at the local scale. Nevertheless, the size of an area source delineating a pattern of high seismicity should be sufficiently small; otherwise the high rate of seismicity is distributed over a larger area source and thus is dangerously reduced due to the effect known as “spatial smearing” (National Research Council, 1988). Therefore, seismic hazard results are different if the same quantity of seismicity is assigned to sources of different size: the smaller the source area, the higher the resulting hazard estimates and vice versa.

The ZS16 seismotectonic zoning was recently also adopted in the new European Seismic Hazard Model (ESHM20) as the reference area-source model for Italy (Danciu et al., 2021).

To determine the activity rates of area sources from earthquake data, we initially used the concept of macroarea to evaluate some parameters, such as the maximum magnitude and the β value of the TruncGR relation, which were then assigned to all the area sources belonging to the macroarea. This approach serves to prevent biased spatial variations in the b value due to poor datasets. The uncertainty in magnitude and epicentral location was assumed to be minimal with respect to the other sources of uncertainty in the estimation of seismic hazard, such as maximum magnitude, completeness time intervals or the choice of the GMMs. However, uncertainty related to the magnitude–location–depth estimation of the earthquakes propagates throughout the evaluation of the completeness periods and thus impacts the estimation of seismic rates.

MA4 explored a range of sources of epistemic uncertainty and is constituted by 20 alternative ERFs that consider different options in terms of completeness time intervals, maximum magnitude and earthquake rate estimation. We observed that the uncertainty related to this set of ERFs is comparable to or even higher than the uncertainty related to alternative GMMs. Although this observation is not generalizable to other ERFs, it contributes to the discussion on the relative importance of ERFs and GMMs in the overall uncertainty affecting seismic hazard estimates. The GMMs adopted in MPS19 were selected according to statistical criteria and elicitation procedures (Lanzano et al., 2020). We did not test the single ERFs (or branches) of MA4 with statistical procedures, but the mean seismic rates and the associated uncertainty were positively checked against observations (i.e. number of earthquakes occurred in the last centuries) before entering MPS19 (see Meletti et al., 2021, for details). In our opinion, however, trimming ERFs according to their performance against observations contributes to reducing the epistemic uncertainty but could result in a selection of models that is only guided by the earthquake occurrence realization observed in the past.

7 Conclusions

In this study, we presented a new seismogenic model for Italy, called MA4, that is part of a community-based effort that led to the development of a new seismic hazard model for Italy (MPS19; Meletti et al., 2021). MPS19, in fact, involved more than 150 Italian researchers at various stages of the project, and many of them were involved in the building of their earthquake rupture forecast (ERF). The finally adopted 11 groups of ERFs were composed of a number of alternatives that explore the epistemic uncertainty in seismicity modelling (see Visini et al., 2021). In this framework, MA4 represents 1 of the 5 area-source seismicity models included in the set of ERFs.

The MA4 model is based on an update (ZS16) of the previous seismogenic zoning ZS9 (Meletti et al., 2008) adopted by the current Italian reference seismic hazard model MPS04 (Stucchi et al., 2011). The new seismogenic zoning consists of 48 active shallow crustal area sources and 2 area sources corresponding to the Campanian and Mt Etna volcanic districts (the last of which was not used for computation here shown).

We used five different approaches to calculate the expected seismicity rates for the area sources, based on two different assumptions. The first assumption is that the truncated Gutenberg–Richter relation parameter β varies across the areas, and the second assumption is that β is stable over groups of area sources with similar seismotectonic features. For this reason, we determined the seismicity rates by first calculating the truncated Gutenberg–Richter parameters for groups of area sources (macroareas) considered tectonically homogeneous and successively partitioned in different ways the values to the area sources, forming each macroarea.

In conclusion, MA4 consists of 20 seismicity models that consider epistemic uncertainty in the estimations of the completeness periods of the earthquake catalogue, of maximum magnitude values and of seismicity rates. The results show that the uncertainties explored by the 20 different seismicity models is at least of the same order of magnitude of the uncertainty due to the different GMMs. We encourage future studies to account for such uncertainty propagation.