A fault-based application to model seismicity rates for seismic hazard assessment in the southern Apennines (Italy)

Although fault-based approaches to seismic hazard assessment have been increasingly adopted worldwide, the official Italian hazard model, on which the national building code is based, still relies on a catalogue-based framework, with well-known limitations in capturing the long-term recurrences of large-magnitude events. In this study, we present a fault-based application to model seismicity rates for the southern Apennines (Italy) that incorporates a multi-fault rupture assumption. This area is of particular interest due to its active seismicity and the presence of large dams, for which robust long-term hazard estimates are essential. We use the SHERIFS code to model seismicity rates at the fault system-level, which allow us to explore epistemic uncertainties of fault and seismicity parameters (rupture scenarios, scaling laws, b-values and background seismicity). Our results highlight the key role of rupture models: scenarios allowing multi-fault ruptures outperform single-fault rupture models in terms of agreement with the regional seismicity and paleoseismic rates. Our findings support the inclusion of multi-fault rupture models in PSHA logic trees for the region and emphasize the need for improved fault behaviour characterization in southern Italy.

1 Introduction

In recent decades, fault-based approaches to seismic hazard assessment have gained increasing attention worldwide as an alternative to traditional methods relying exclusively on earthquake catalogues. These latter approaches, based on the original formulation of Probabilistic Seismic Hazard Analyses (PSHA) by Cornell (1968), have completeness limitations related to the combined effect of short earthquake catalogues (usually no more than a thousand years) and the long recurrence intervals of moderate-to large-magnitude earthquakes. These effects hamper the earthquake forecasting capabilities of PSHA, especially in regions with lower deformation rates. Moreover, catalogue-based approaches classically employ distributed seismicity models (area sources, grid sources, etc.) which, although suitable for representing background seismicity, i.e. earthquakes occurring on unmapped or poorly constrained faults, tend to unify or smooth out the seismic hazard over larger regions rather than localizing it along main active fault structures. Conversely, fault-based PSHA allows also the incorporation of fault-specific geological data (e.g., paleoseismology, geomorphology) into the source modelling, enabling a time-extended and more reliable characterization and localization of the seismic activity of a region. This approach provides improved spatial resolution of seismic hazard associated with moderate-to-large magnitude earthquakes (M_w≥6.5–7.0), whose occurrence is strongly controlled by the geometry, slip rates and segmentation of major active faults (e.g., Field et al., 2009; Gerstenberger et al., 2024).

Quantitative comparisons between the two approaches have shown that incorporating geological fault information may significantly modify both the spatial variability and the amplitude of ground-motion estimates. For example, fault-based models have been shown to produce increases in peak ground acceleration (PGA, 10 % probability of exceedance in 50 years) of up to 0.3 g compared to catalogue-based models in regions characterized by low deformation rates and limited instrumental seismicity (e.g., Williams et al., 2023), while differences exceeding 100 % have been observed for sites located close to major active faults (Gómez-Novell et al., 2020b). In California, the transition from area-source models to the fault-based UCERF framework reduced discrepancies between observed and modelled rates of intermediate-magnitude earthquakes ($\mathrm{6.5}\le M\le \mathrm{7.0}$), improving consistency with historical catalogues (Field et al., 2009).

Despite the evolution and implementation of fault-based PSHA approaches in many countries worldwide (e.g., New Zealand, Gerstenberger et al., 2024; California, Petersen et al., 2024) and the advances achieved by large-scale hazard models such as ESHM20 (Danciu et al., 2024), whose results have been incorporated into the second generation of Eurocode 8 (Labbé and Paolucci, 2022), catalogue-based PSHA remains the standard practice, particularly for official seismic hazard assessment and regulatory applications. For example, in Italy the current building code (NTC, 2008, 2018) still relies on an area-source earthquake rupture forecast (Stucchi et al., 2011), as does the current Spanish seismic code (NCSE, 2002).

However, the reliance of catalogue-based hazard models on historical earthquakes can lead to the exclusion of faults with paleoseismic evidence of rupture – and thus capable of generating significant earthquakes – that have not produced events within the historical record. A relevant example, in Italy, is the Mount Vettore Fault, which ruptured during the 2016 Central Apennines earthquake even though it was previously classified as silent (Galadini and Galli, 2003; Galli et al., 2019).

Beyond these limitations, catalogue-based approaches do not account for the complexity of earthquake ruptures observed in nature. Seismological, geological and paleoseismological data show that earthquake ruptures can be very complex, involving the simultaneous activation of faults with different characteristics. Recent complex coseismic ruptures highlighted the need to include the multi-fault earthquakes in PSHA, going thus beyond strict fault segmentation assumptions (e.g., M_w 6.9 Irpinia in 1980, Bernard and Zollo, 1989; M_w 7.1 El Mayor-Cucapah in 2010, Wei et al., 2011; M_w 8.6 Sumatra in 2012, Zhang et al., 2012; M_w 6.5 central Italy in 2016, Peruzza et al., 2016; M_w 7.8 Kaikoura in 2016, Hamling et al., 2017; M_w 7.8 and M_w 7.7 Türkiye–Syria in 2023, Liu et al., 2023). Capturing such multi-fault ruptures is a key challenge in fault-based PSHA: implementations in California with UCERF3 (Field et al., 2009) or New Zealand (Stirling et al., 2012) pioneered the inclusion of such fault rupture complexities in PSHA.

However, incorporating faults into seismic hazard assessment in low-strain regions, such as Europe, is not yet a widespread practice and remains a topic of debate in many countries. A regulated implementation of such fault-based models in Europe remains a challenge mainly due to the larger recurrence intervals of the active faults and therefore to the scarcity of data to characterize such activity. In recent years, some initiatives such as the Fault2SHA Working Group of the European Seismological Commission have been founded to discuss the implementation of fault-based PSHA in Europe (central Italy, e.g., Peruzza et al., 2011; Pace et al., 2018; Valentini et al., 2018, 2019; Valentini, 2020; France, Scotti et al., 2014; Greece, Deligiannakis et al., 2018; and Spain, Gómez-Novell et al., 2020b). Along this line, in the recent MPS19 Seismic Hazard Model of Italy, two fault-based models were included in the earthquake rupture forecast, marking a significant step forward in Europe (Meletti et al., 2021; Visini et al., 2021).

In this study we model fault seismicity rates in the southern Apennines using a relaxed segmentation framework, i.e., considering multi-fault rupture scenarios. To our knowledge, this is the first study in the area to adopt such an approach for seismic hazard assessment. The Irpinia region, severely affected by the 1980 M_w 6.9 earthquake, the third deadliest in Italy (Peruzza, 2018), hosts several large earth dams located near active but poorly characterized faults, highlighting the need for accurate estimates of fault activity rates.

The study is designed as a comprehensive sensitivity test, aimed at observing how various model configurations of fault ruptures, fault and seismicity parameters impact fault system-level seismicity rates for seismic hazard assessment.

Geological and seismological setting

The study area, covering about 50 000 km², is centred around the Irpinia district, one of the most seismically active regions of the southern Apennines (southern Italy; Fig. 1). The southern Apennines are part of an east-verging thrust belt related to the west-dipping subduction of the Apulian lithosphere beneath the Tyrrhenian domain (Doglioni et al., 1996). Since the Late Quaternary, slab rollback and trench retreat have promoted extensional tectonics along the Apennine chain, resulting in a system of large, segmented NW–SE-striking normal faults accommodating crustal stretching (Doglioni et al., 1996). At the same time, strike-slip to oblique-slip faulting is observed in the outer sectors of the belt and within the Adriatic foreland. These structures are commonly interpreted as the reactivation of inherited crustal discontinuities within the Apulian foreland, occurring at depth (typically ∼ 15–30 km; Latorre et al., 2023) and partially decoupled from the shallow extensional regime affecting the chain (generally confined within the upper ∼ 10–15 km of the crust; e.g., Boncio et al., 2007; Vannoli and Burrato, 2018). Recent studies (e.g., Adinolfi et al., 2015) have proposed that the coexistence of extensional and strike-slip faulting may reflect a vertically partitioned seismotectonic system, with shallow normal faulting and deeper strike-slip deformation. However, the relationship between these deformation regimes, and whether they are mechanically coupled or spatially overlapping, remains debated.

Figure 1 Geographical distribution of the M_w≥4.0 earthquakes occurred in the southern Apennines between 1019 and 2017 ($\mathrm{4.0}\le M_{\mathrm{w}}\le \mathrm{7.2}$; reference catalogue: CPTI15 v4.0, Rovida et al., 2022) and the fault traces of the 35 seismogenic sources used for this study (references: Valentini et al., 2017; DISS Working Group, 2025). Earthquakes are colour-coded by magnitude class (M_w bin = 1.0), with larger symbols representing stronger events; major historical earthquakes are labelled by year. Faults are colour-coded by tectonic domain. Latitude and longitude grid shown in WGS84 coordinates.

This tectonic framework explains the subdivision of the study area into two main domains or fault systems. The strongest seismic events of the region (Fig. 1) mostly occurred in one of the domains, a segmented belt of large normal NW–SE striking faults, running along the chain axis, hereafter “SubArea 1” (DISS Working Group, 2025). In the other system, i.e., the foreland, faulting develops along E–W right-lateral strike-slip to oblique-slip faults, hereafter “SubArea 2” (DISS Working Group, 2025).

The strongest events (Fig. 1) reported in the Parametric Catalogue of Italian Earthquakes – CPTI15 v4.0 (Rovida et al., 2022) for this region are the 1456 Sannio-Irpinia sequence (M_w 7.2), which is one of the most destructive events that took place in the Italian peninsula, the 1688 (M_w 7.1) and the 1857 (M_w 7.1) earthquakes. The probably most famous destructive earthquake, the M_w 6.9 that hit the Irpinia district in 1980, had devastating effects in a large area of the southern Apennines, causing thousands of casualties and high damage to buildings. Latorre et al. (2023) observe that the highest earthquake concentration in the southern Apennines consistently occurs along the faults responsible for the 1980, M_w 6.9 Irpinia earthquake. Close to this epicentre, it is also worth to mention the 1694 (M_w 6.7) Irpinia-Basilicata earthquake, considered the ancestor of the 1980 (DISS Working Group, 2025), the 1851 (M_w 6.5), and the 1930 (M_W 6.7) Irpinia earthquake.

2 Methods, data and model parameters

2.1 The SHERIFS method

In this study, we employ SHERIFS V1.1 (Seismic Hazard and Earthquake Rate In Fault Systems), an approach introduced by Chartier et al. (2017) and further developed as an accompanying code by Chartier et al. (2019), designed to model seismicity rates in active fault systems using geological fault data as inputs. The approach treats fault systems as a whole, explicitly accounting for multi-fault rupture scenarios in the estimation of seismicity rates. So far, SHERIFS has been tested in several regions worldwide including SE Spain (Gómez-Novell et al., 2020a, b), the Marmara Region in Türkiye (Chartier et al., 2021), the Tibetan Plateau (Cheng et al., 2021) and the Levant Fault (El Kadri et al., 2025).

To derive seismicity rates from geological slip-rate data, SHERIFS requires several inputs: (i) the 3D geometry and slip rates of the fault system, (ii) a set of multi-fault rupture scenarios based on user-defined criteria (e.g., inter-fault distance), and (iii) a target Magnitude-Frequency Distribution (MFD) specified at the fault-system level (e.g., a Gutenberg-Richter law). The epistemic uncertainties of the input parameters in SHERIFS are explored through an iterative process that performs the computation n times exploring n different values of the input parameters within their uncertainty ranges. The result is the average and uncertainty dispersion cloud of all samples explored.

One of the main characteristics of SHERIFS is that the approach treats the prescribed slip rate of the fault system as a budget. Through an iterative procedure, this slip rate budget is distributed among faults according to the magnitudes they can accommodate, following analytical relationships that link seismic moment rate to seismicity rates. The slip rate budget is spent until the pre-imposed MFD target is reached and the slip rate budgets are consumed. In some cases, the target is reached before the budget of all faults is exhausted and therefore there is remaining slip rate budget that is not converted into seismicity rates. This remaining budget, namely Non-Main-Shock slip (hereafter NMS) in SHERIFS, is indeed an artifact of the modelling. However, according to the authors, small percentages of NMS (generally below 30 %–40 % of the total budget of the fault), can be interpreted as a natural phenomenon (e.g., post-seismic slip or creep) as geological slip rates might not be exclusively resulting from the coseismic phase, but also from other post-seismic relaxation processes. Conversely, NMS percentages above 30 %–40 % are less easily attributable to the natural phenomena described. Instead, they most likely indicate that the combination of input hypotheses used (fault parameters, target MFD and rupture hypotheses) is not optimal in the SHERIFS framework and that they should be reconsidered (Chartier et al., 2019).

In our model, the epistemic uncertainty is explored in two ways:

1. Model branches: we define 13 modelling branches (Fig. 2), covering variables such as rupture scenarios, scaling relationships, b-values, and background seismicity and buffer areas.

2. Random samples: we run 10 random samples per branch that explore different values of slip rate on faults within their uncertainty ranges.

Figure 2 Structure of the main model branches (blue) used to explore epistemic uncertainty in the fault-based model, and of the consistency check criteria (violet) adopted to assess their performance.

The inputs used in our modelling, which are discussed in detail in the following sections, are archived and openly accessible in a Zenodo repository (Alessandrini et al., 2026).

2.2 Input data and model parameters

2.2.1 Fault sources

We use 35 faults of the southern Apennines for our study. All fault identification names (thereafter IDs), and their main characteristics are listed in Table 1, while their geometries are shown in Fig. 3. Fault traces, 3D geometry (dip and seismogenic depth), slip rates and kinematics are mainly extracted from the Individual Seismogenic Sources (ISS) within the Database of Individual Seismogenic Sources – DISS 3.3.1 (DISS Working Group, 2025), a geo-referenced archive of information on seismogenic sources and of their seismogenic potential (Basili et al., 2008). This database is a reference in Italy and many recent seismic hazard studies have employed this database for seismic source modelling (e.g., Stucchi et al., 2011; Woessner et al., 2015). We include some additional active faults reported by Valentini et al. (2017), that are not part of the DISS (respectively faults 11, 15, 16, 26, 27, 28, 29, 30, 32, 33, 34, 35).

Table 1 List of faults employed in the models. Faults highlighted by a star (^*) are the integration of Valentini et al. (2017) while the other faults come from the DISS 3.3.1 database (DISS Working Group, 2025). For each of them, we list here the name, the fault ID used in the model, the dip angle, the kinematics (S = strike-slip, N = normal), the seismogenic depth range (km), length and the slip rate range (mm yr⁻¹). Paleoseismic data are available for faults reported in the DISS database and here marked by a double star (${}^{**}$).

Figure 3 Representation of the 35 fault sections, coloured based on their mean slip rates (references: Valentini et al., 2017; DISS Working Group, 2025).

Fault geometries and segmentation are directly adopted from these datasets, which define individual seismogenic sources based on geological, geophysical and seismotectonic criteria, including geological mapping, seismicity distribution, geomorphic evidence and geodetic observations. Potential discontinuities in fault traces, particularly within the strike-slip systems, may reflect the presence of unmapped or blind structures. For example, additional faults are identified in the ITHACA catalogue (ITaly HAzards from CApable faults; ITHACA Working Group, 2019). However, these structures generally lack the quantitative parameters required for earthquake rate modelling, such as slip rates and seismogenic depths. For this reason, only faults with sufficiently constrained geometric and kinematic parameters are included in the present modelling framework. The adopted segmentation therefore reflects the current state of knowledge on active faults for seismic hazard applications in the region.

The slip rates we assigned to each fault come from the DISS 3.3.1 database and by Valentini et al. (2017) and are constrained by geological and geodynamic observations. Variability in slip rates across adjacent faults is consistent with the heterogeneous accommodation of strain within fault systems, which depends on factors such as fault geometry, orientation and mechanical properties. This behaviour has been documented in several tectonic settings, including the Apennines, and is reflected in regional fault databases (e.g., Faure Walker et al., 2021).

2.2.2 Segmentation rules and rupture scenarios

We explore fault-rupture scenarios to compute system-level seismicity rates under different hypotheses of fault connectivity within the fault system. We define three fault rupture scenarios for the whole fault system (Table 2), each representing an incremental increase in fault connectivity and, consequently, in multi-fault rupture length. These scenarios are exploratory hypotheses, and as such, other configurations could be tested.

Table 2 List of the three fault rupture scenarios, characterized by incremental multi-fault ruptures. For each scenario the combinations of the maximum ruptures are listed, and the expected maximum magnitude range is provided. Fault sections (from 1 to 35) are shown in Fig. 3.

• Set_0. In this hypothesis, only single fault section ruptures are allowed. The segment lengths reported in the literature for each fault define the maximum rupture length, and multi-segment ruptures involving neighbouring faults are not allowed.

• Set_1. In this scenario a subset of multi-fault ruptures is allowed based on geometric compatibility rules (i.e., faults with similar strike, dip direction, seismogenic depth), kinematics compatibility, as well as a fault distance threshold. In terms of kinematic compatibility, we follow the cumulative rake change criterion suggested by Milner et al. (2013), which limits the total rake variation within a rupture to 180° to prevent kinematic incompatibility. This means that only faults with the same kinematics are allowed to rupture together. In terms of distance for rupture propagation, we select 5 km as the threshold for rupture propagation, consistent with previous studies on multi-fault rupture connectivity, including the UCERF-3 framework (e.g., Milner et al., 2013; Field et al., 2014; Scotti et al., 2020), where similar inter-fault jump distances are commonly used. Based on the above criteria, Set_1 permits ruptures between neighbouring faults up to a maximum of 4 faults rupturing together (Fig. 4) with rupture lengths up to 75 km, which falls well within the typical rupture lengths observed in Gerstenberger et al. (2024), usually below a few hundred kilometres.

• Set_2. The distinct feature of Set_2 is the adoption of a larger distance threshold of 10 km for rupture propagation, while maintaining the requirement of compatible faulting mechanism. This configuration is intentionally more speculative and is introduced to explore an upper-bound range of rupture behaviours and their impact on seismicity rates. The choice of extending the distance threshold beyond 5 km is motivated by observations of complex multi-fault ruptures in the study area, such as the 1980 M_w 6.9 Irpinia earthquake, which involved multiple adjacent fault segments (faults 21–22–23) separated by distances exceeding 5 km. Under this configuration, up to 8 faults can rupture together (Fig. 4) resulting in rupture lengths up to 184 km.

Figure 4 Maximum extent (number of faults) of the multi-fault rupture scenarios in Set_1 and Set_2 explored in this study. Faults with the same colour are allowed to rupture together in a same rupture within that particular scenario set.

In order to study the relationship between magnitude and rupture area, we employ the empirical scaling relationships of Wells and Coppersmith (1994), Leonard (2010), and Thingbaijam et al. (2017). These scaling laws, widely used in PSHA applications, have broad applicability across tectonic settings and are implemented in the current version of the SHERIFS code. Leonard (2010) has been employed in regional studies in Italy and is used for the parametrization of national databases (e.g., DISS Working Group, 2025), supporting its suitability for the study area. Wells and Coppersmith (1994) provide a robust and widely adopted model that does not require specific adjustments between continental extensional and compressional regimes. In addition, Thingbaijam et al. (2017) provides coefficients specifically derived for normal faults, which represent the dominant faulting style in the Apennines, allowing a better representation of this tectonic regime. The maximum magnitude (M_max) range is set by the scaling relationship, and it is shown for each scenario in Table 2. A comparative evaluation of their performance within the study area is presented in Sect. 4.3, where their implications for the modelling results are discussed.

2.2.3 Model MFD shape

To model the magnitude-frequency distribution (MFD) of seismicity on individual faults, we adopt a Gutenberg-Richter (G-R; Gutenberg and Richter, 1944) relationship, which is widely used in PSHA and is consistent with the Italian national seismic catalogue (e.g., CPTI; Rovida et al., 2022). To explore the influence of epistemic uncertainty, we test two different b-value ranges: a broader interval (0.9–1.1) based on the national CPTI catalogue (Rovida et al., 2022), and a narrower one (0.93–0.96) supported by previous studies focused on this region (Gulia and Meletti, 2007; Meletti et al., 2008).

2.2.4 Background seismicity

One common challenge in fault-based PSHA is selecting a magnitude threshold for a confident assignation of seismicity to the known active faults, as small-to-moderate earthquakes rarely cause surface faulting and therefore may not be related to the known active faults. Because the classical selection of a cut-off threshold can be somewhat arbitrary, SHERIFS deals with this issue by incorporating magnitude-dependent ratios that determine the proportion of the seismicity that will be assigned to the faults, usually increasing with magnitude. That is, for each magnitude a user-defined proportion – ratio – of the seismicity rates of the faults is subtracted and assigned to the background as a distributed source. This approach allows for a smoother distribution of seismicity between faults and background and is consistent with real observations.

We explore two alternative scenarios for background seismicity:

• BG_1, a fault-only model where all seismicity is attributed to the known active faults.

• BG_2, a hybrid model where a magnitude-dependent fraction of seismicity is assigned to the background buffer (Fig. 5), outside the mapped faults.

Figure 5 Buffer area (10 km around the faults, dashed line) used in the BG_2 scenario to define the background seismicity and to extract the catalogue for the consistency check. Earthquakes from the CPTI15 v4.0 catalogue (Rovida et al., 2022) are shown for reference. Internal buffer areas are also displayed for the two structural domains (SubArea 1, blue shaded buffer and SubArea 2, pink shaded buffer).

To define the BG_2 scenario, we set up a buffer area of 10 km around the faults (Fig. 5), where background seismicity is assigned. This 10 km distance is consistent with background zone sizing proposed by Chartier et al. (2019) (e.g., 10 or 20 km from faults). In order to evaluate the sensitivity of the results to this modelling choice, we also tested alternative buffer distances of 5 and 15 km. The comparison between buffer sizes is further discussed in Sect. 4.4, and the corresponding models and figures are provided in the Zenodo repository (Alessandrini et al., 2026).

The ratios for each magnitude bin (Table 3) are custom set up using a magnitude-dependent criterion, based on the assumption that larger earthquakes are more likely to be associated with known, geomorphologically recognized faults, whereas smaller events may occur in secondary unmapped or buried faults. Consequently, BG_2 assigns higher background ratios to lower magnitudes.

Table 3 Ratio of seismicity occurring on known faults for each background scenario considered in the model (M_w bin = 0.5).

2.2.5 Model setup

Each branch of the model is computed 10 times using different random values (samples) within the uncertainty ranges of the different input parameters (e.g., fault slip rates, b-values) to explore the epistemic variability of the seismicity rates linked to the input parameter uncertainties. This results in the computation of 1080 different models, 360 for each rupture scenario branch.

2.3 Comparison with observations

2.3.1 Regional historical and instrumental seismicity

We adopt the Parametric Catalogue of Italian Earthquakes – CPTI15 v4.0 (Rovida et al., 2022) as a regional seismic catalogue to compare the synthetic MFDs modelled in this study and to evaluate the models that best fit such seismicity rates. CPTI is nowadays the main tool for seismic hazard models in Italy, as it homogenises the magnitudes of pre-instrumental (historical) and instrumental events (Rovida et al., 2020). The catalogue contains all known Italian earthquakes from year 1000 AD. All the available events for this work have been selected within a radius of about 150 km encompassing the study area in the southern Apennines. The resulting selection (Fig. 1) spans a time interval from year 1019 to 2017, and a moment magnitude range of $\mathrm{2.9}\le M_{\mathrm{w}}\le \mathrm{7.2}$. The minimum magnitude for our analysis is set at M_w 4.0, and the completeness magnitudes to obtain the GR distributions for the CPTI catalogue are reported in Rovida et al. (2020). For the comparison with the models, only the events falling within the buffer area in Fig. 5, which is regarded as the area of influence of the active faults, are considered.

2.3.2 Paleoearthquake rates

In faults with available paleoseismological studies and, as such, estimates of paleoearthquake recurrence intervals, we compare the modelled participation rates of that fault with the estimated paleoearthquake rates from paleoseismological studies (paleoearthquake rate = 1 $/$ paleoearthquake recurrence interval). We considered that the paleoseismological earthquake rates (paleorates) data represent earthquakes of M_w≥6.0–6.5, based on the statistical fact that most surface rupturing earthquakes (i.e., those observed in paleoseismology) have a magnitude equal or higher than the said threshold (Baize et al., 2020; Nurminen et al., 2022).

Paleoseismic data, including annual paleoearthquake rates, are available for three faults activated during the 1980 M_w 6.9 Irpinia earthquake: Pescopagano (21), Colliano (22), and San Gregorio Magno (23), thanks to several studies (e.g., Pantosti and Valensise, 1990; D'Addezio et al., 1991; Pantosti et al., 1993).

3 Results

Our results show that the choice of rupture scenario is the primary factor driving differences in the modelled seismicity rates, whereas the other explored branches have a much more limited impact. For this reason, to enable a direct comparison of the different rupture scenarios, we display our results for each rupture scenario showing a single branch for the scaling relationship, b-value, background seismicity and buffer area (“Leonard (2010)”, “b-value = 0.9–1.1”, “BG_1”, “buffer area = 10 km”, respectively). The results for all the other branches are available in the Zenodo repository of the paper (Alessandrini et al., 2026). The quality of our models is tested through the evaluation of their fit with the regional earthquake catalogue MFD, the overall model performance based on the SHERIFS NMS proportion, and the agreement with paleoseismic slip rates.

3.1 Modelled MFDs: fit with the earthquake catalogue

In the upper three panels of Fig. 6, we show the fit between the synthetic MFDs, i.e., the modelled seismicity rates, and the regional seismic catalogue for the different rupture scenarios. In each scenario, the synthetic MFD is formed by combining 10 different MFDs, generated from the 10 random samples explored in the modelling.

Figure 6 Columns show different rupture scenarios. Above: comparison between the modelled MFDs (green) and the seismicity rates from the historical and instrumental regional catalogues within the study area (red). Solid green line is the mean MFD and green patches represent the uncertainty (16–84 percentiles). Below: NMS (expressed as a percentage) for each random sample explored in the model.

Results show that the choice of the rupture configuration directly influences the shape of the synthetic MFDs. The intermediate multi-fault rupture scenario Set_1 and the large multi-fault rupture scenario Set_2 exhibit the best fit with the catalogue, as the modelled rates with SHERIFS overlap with the catalogue-derived MFD. In particular, Set_1 displays a strong agreement with the catalogue for magnitudes up to M_w 6.0, while Set_2 tends to slightly underestimate the MFD across this magnitude range, but the modelled and observed rates still remain in good agreement when accounting for the respective variability.

Conversely, for Set_0 there is no match between the modelled and catalogue MFDs: the synthetic MFD lies consistently below the catalogue-derived curve across the entire magnitude range, even when considering their respective uncertainty bounds.

3.2 Non-Main-Shock slip

The choice of a rupture configuration also influences how the slip rate budget is consumed in the different SHERIFS iterations. When part of this budget remains unspent at the end of the calculation, it is referred to as non-mainshock slip. The NMS is expressed as a percentage and represented through histograms for each of the 10 random samples in each scenario (Fig. 6, lower panels). Following Chartier et al. (2019), NMS percentages above 30 %–40 % in a model indicate poor model performance. Accordingly, we use the NMS calculated by SHERIFS as an indicator of model quality: lower NMS correspond to better performance, as the models are able to spend the prescribed slip rate in the computation.

The NMS is also useful to differentiate between models with otherwise similar fits with the catalogue, as in the case of Set_1 and Set_2 (Fig. 6, upper panels). In our results, Set_1 exhibits the lowest NMS values among all models (mean value of 15.5 %), with 90 % of its samples falling below the 40 % threshold. Set_2 shows slightly higher values, but its mean value (23.9 %) remains well within the acceptable range. Conversely, Set_0 displays consistently high NMS values, ranging from (85 %–97 %). These results confirm that Set_1 and Set_2 yield the most robust performances, in line with their closer agreement with the catalogue-derived seismicity rates (Fig. 6).

3.3 Paleoearthquake rates on faults

Here we compare the participation rate curves modelled with SHERIFS in each of the three fault sections with the paleorates documented (faults 21, 22, 23; Fig. 3), for each rupture scenario, including their magnitude and rate uncertainty ranges (Fig. 7). The participation rates are derived from all ruptures involving these fault sections.

Figure 7 Comparison of the modelled paleoearthquakes rates on faults and paleoseismic data. Green solid curves are the modelled mean cumulative participation rates rupturing the specific fault section; dashed lines are the dispersions of individual samples; purple points are the paleorates measured on the fault trench; purples boxes are the paleodata uncertainty ranges. The location of each fault (rows) is highlighted in red within the fault maps, for the three different scenarios (columns).

While all three rupture scenarios tend to underestimate the paleoseismic rates for fault 21, Set_1 and Set_2 closely match the paleoearthquake rate mean values and their uncertainty ranges for faults 22 and 23. Set_1 only slightly overestimates the central value for fault 22, but remains within the confidence intervals. In contrast, Set_0 consistently and significantly underestimates the measured rates for all three faults, confirming its weaker performance in line with the previous quality tests.

4 Discussion

4.1 Segmentation models in the southern Apennines

The results of our analysis show that the impact of different rupture scenarios leads to markedly different behaviours in terms of seismicity rates, MFD shape, and expected maximum magnitudes. These results have implications for the segmentation models to be considered for PSHA in the southern Apennines.

Both Set_1 and Set_2, which allow ruptures involving multiple adjacent faults, show a good fit to the catalogue. Set_1, in particular, stands out as the most balanced configuration, showing strong consistency with the observed seismicity and producing a maximum expected magnitude (M_w 7.4) closer to historical records with respect to Set_2 (M_w 7.6). Conversely, individual segment ruptures show lower consistency with the observations for the study region. The good agreement of multi-fault rupture earthquake rates with the seismic and paleoseismic records in the region indicates that such multi-fault rupture scenarios are statistically consistent with observations in the southern Apennines, given the modelling set up and fault input data used. These results are also in line with historical observations, such as the 1980 Irpinia earthquake, which involved the rupture of multiple fault segments (21–22–23; Fig. 3, Table 1). In this sense, this study goes beyond simply generating alternative rupture sets and instead evaluates which classes of rupture complexity are required to reconcile independent observations at the regional scale.

Multi-fault rupture scenarios also have important implications for the maximum expected magnitude in the region, allowing earthquakes up to M_w 7.6 in the most permissive configurations (Set_2). Although such magnitudes are not represented in the historical catalogue (CPTI; Rovida et al., 2022), their occurrence cannot be rejected based on the available geological and seismological constraints adopted in this study. It is important to remark that SHERIFS is not a physics-based model (e.g., it does not simulate stress transfer or rupture propagation dynamics) to evaluate whether the fault system can accommodate such large magnitudes. Instead, it relies on statistical metrics derived from observations, which might not fully capture the specific characteristics of the southern Apennines fault system. Nevertheless, the explored rupture scenarios are not arbitrarily generated. Their definition is constrained by fault geometry, segmentation patterns, rupture jump distances, kinematic compatibility and seismogenic depth constraints, as discussed in Sect. 2.2.2.

Despite this and as observed in other tectonic regions (e.g., California, Field et al., 2014; New Zealand, Quigley et al., 2017; Taiwan, Chang et al., 2023), relaxing segmentation can give enhanced performance of seismicity models for PSHA. While the rupture scenario models explored in this study are a sample of the possible rupture scenario models to be explored, our results give valuable insight for future seismic hazard assessments in the region, suggesting that rupture complexity is a first-order control on long-term seismicity rates in the southern Apennines.

4.2 Tectonic domains impact

The modelled faults in the study area span two tectonically distinct domains: SubArea 1, characterised by extensional normal faulting, and SubArea 2, dominated by strike-slip structures (Fig. 1). These tectonic domains not only differ in fault kinematics but also in their seismic expression: SubArea 1 shows a denser clustering of seismicity, while SubArea 2 appears notably more quiescent. To investigate how these differences may influence model behaviour, we analysed the two sub-areas independently, that is, we computed the models considering only the faults and seismicity catalogue within each sub-area (the background seismicity buffers specific to each sub-area are shown in Fig. 5).

The results show that the two tectonic domains exhibit different statistical performances under the tested rupture hypotheses. In particular, we observe that the modelled rates of SubArea 1 faults (Fig. 8) return considerably better fits with the catalogue rates than SubArea 2 (Fig. 9), which systematically underestimates the catalogue. In SubArea 1, Set_1 is the configuration with the best overall performance for its agreement with the catalogue for M_w>6 and considerably low NMS percentages. The other sets underperform Set_1: Set_0 has higher NMS percentages, although holding mean values within the 40 % threshold; Set_2 slightly underestimates the catalogue rates. These results, especially the agreement of Set_1 and the NMS values across models, are consistent with the full model analysis (both sub-areas combined) and suggest that the extensional domain of the southern Apennines (SubArea 1) contributes substantially to the overall seismicity rates reproduced by the model. Conversely, the results of SubArea 2 indicate that the currently available fault data do not fully capture the seismicity rates observed in the catalogue in this area. While the lower seismic activity in SubArea 2 could suggest a role of tectonic behaviour, the modelling framework adopted in this study is not intended to resolve tectonic controls on seismicity. Rather, it evaluates the statistical consistency of alternative rupture hypotheses with available observations, which are inherently influenced by the quality and completeness of input data.

Figure 8 SubArea 1 models: columns show different rupture scenarios. Above: comparison between the modelled MFDs (green) and the seismicity rates from the historical and instrumental regional catalogues within the study area (red). Below: NMS (expressed as a percentage) for each random sample explored in the model.

Figure 9 SubArea 2 models: columns show different rupture scenarios. Above: comparison between the modelled MFDs (green) and the seismicity rates from the historical and instrumental regional catalogues within the study area (red). Below: NMS (expressed as a percentage) for each random sample explored in the model.

The observed underestimation of catalogue rates in SubArea 2 may be explained by two principal factors. First, the faults in this region generally have lower slip rates (Fig. 3), and in some cases, relatively long lengths. Since SHERIFS treats slip rate as a budget that is spent on larger magnitudes first (due to its moment rate formulation), these faults tend to release most of their seismic moment through fewer, higher-magnitude events. As a consequence, the earthquake rates across the full MFD are reduced, particularly at lower magnitudes.

A second factor may lie in the limited geological constraints available for several faults in SubArea 2. Many slip rates in this region come from regional geodynamic models rather than direct geological constraints (see DISS Working Group, 2025), and paleoseismic data are notably scarce. Although some structures – particularly in the Gargano Promontory – are well studied and associated with significant seismicity, the overall quantity and resolution of geological data in SubArea 2 remain more limited than in SubArea 1, likely due to the lower number of large earthquakes historically recorded in the region.

Additionally, certain faults in SubArea 2 are modelled with long, simplified traces that may not fully represent their internal segmentation. This generalisation could lead to an overestimation of rupture length in the models, ultimately affecting seismicity rate calculations. Similar effects have been observed in previous SHERIFS applications (e.g., Gómez-Novell et al., 2020a), where simplified fault geometries impacted the agreement between modelled and observed seismicity. From this perspective, the poorer performance of SubArea 2 should not be regarded solely as a limitation of the model, but rather as an informative result. Neither the current fault model, as represented in DISS, nor the wide range of rupture configurations explored in this study are able to satisfactorily reproduce the observed catalogue rates in this sector. This mismatch suggests that important geological constraints may still be incomplete or insufficiently resolved, such as slip-rate estimates, fault geometries, segmentation patterns, or potentially unmapped active structures.

Overall, the comparison between SubArea 1 and SubArea 2 highlights the critical role of geological constraints in controlling the robustness of fault-based PSHA models and helps identify priorities for future fault characterization efforts in the southern Apennines.

4.3 Scaling relationship impact

Among the three scaling relationships tested, the differences between Leonard (2010) and Wells and Coppersmith (1994) are minor, with both yielding consistent and reliable results across all rupture scenarios. The Leonard (2010) scaling relationship was ultimately selected as the reference scaling relationship in this study (Fig. 6), as it delivers a slightly better fit to the MFDs, lower NMS percentages, being also a more recent formulation applicable to a wider range of tectonic settings.

In contrast, the scaling law by Thingbaijam et al. (2017), while producing a closer match to the regional MFD in the single-fault scenario (Set_0), results in only minor differences across the three rupture sets (Fig. 10), making it less effective in distinguishing between different rupture behaviours. This, together with a generally broader dispersion and a systematically higher proportion of non-main-shock slip, limits its applicability in our modelling framework.

Figure 10 Rows correspond to different scaling laws. Columns 1–3 show the comparison between modelled MFDs (green) and seismicity rates from the historical and instrumental regional catalogues within the study area (red), for the three rupture scenarios. Column 4 displays the NMS (expressed as a percentage) for each random sample explored under all rupture scenarios.

This behaviour may be attributed to the nature of the Thingbaijam et al. (2017) scaling law, which instead of geological data is based on source modelling inversion datasets. This scaling law tends to estimate smaller magnitudes for a given rupture area compared to Wells and Coppersmith (1994), as it can be seen by the maximum magnitudes of the models depending on the scaling relationship. Because SHERIFS operates under the preservation of the seismic moment assumption (budget), the reduction of M_max in the Thingbaijam et al. (2017)-based models is compensated in the whole MFD, which leads to more similar earthquake rates across all models. This formulation seems less effective for our fault system and therefore has not adopted as the reference model.

4.4 Impact of MFD shape and background seismicity

Concerning the MFD shape assumption, the choice of a Gutenberg–Richter distribution fits the observed seismicity well and proves to be a valid and robust statistical representation. This choice is also consistent with long-standing applications of the GR law to the instrumental seismicity in Italy (e.g., Stucchi et al., 2011; Rovida et al., 2020; Meletti et al., 2021). Exploring a wider or a narrower b-value range has also a limited impact in our results (Fig. 11). Generally, the wider b-value range of 0.9–1.1 describes slightly better the regional seismicity and shows slightly lower NMS percentages. This is because the MFD shape is less restrictive in these models and allows the modelling to better capture the rate variability of the catalogue.

Figure 11 Rows correspond to different b-values. Columns 1–3 show the comparison between modelled MFDs (green) and seismicity rates from the historical and instrumental regional catalogues within the study area (red), for the three rupture scenarios. Column 4 displays the NMS (expressed as a percentage) for each random sample explored under all rupture scenarios.

Consistent with previous studies that acknowledge the relevance of distributed background seismicity models to account for events not clearly linked to known faults (e.g., Stirling et al., 2002; Chartier et al., 2019; Gerstenberger et al., 2024), we configured BG_2 to include a magnitude-dependent fraction of seismicity within a 10 km buffer zone surrounding the mapped faults Set_0 is the one showing the larger differences related to background seismicity, allowing to partially compensate the model underestimation of the earthquake catalogue rates (Fig. 12). Conversely, the effect of BG_2 on multi-fault rupture models (Set_1 and Set_2) is minimal, and in some cases slightly worsens the fit to the observed MFD. Additionally, BG_1 yields slightly lower NMS values (Fig. 12), with only marginal differences. This indicates that introducing background seismicity does not produce any significant improvement in our modelling results and, as such, the relative performance between rupture scenario sets remains invariant.

Figure 12 Rows correspond to different background seismicity scenarios. Columns 1–3 show the comparison between modelled MFDs (green) and seismicity rates from the historical and instrumental regional catalogues within the study area (red), for the three rupture scenarios. Column 4 displays the NMS (expressed as a percentage) for each random sample explored under all rupture scenarios.

To evaluate the sensitivity of these results to the definition of the background area, we additionally tested buffer distances of 5 and 15 km (available in the Zenodo repository; Alessandrini et al., 2026). The results obtained with these alternative configurations are consistent with those derived for the 10 km buffer, confirming that the overall outcomes are not significantly influenced by the selected buffer size. Moreover, a 5 km buffer from fault traces is not optimal because it does not always capture the whole surface projection area of the faults and, as such, might exclude some seismicity related to the faults at study, while the 15 km buffer might include seismicity out of the known fault systems.

Although a distributed seismicity component is widely considered essential in modern PSHA frameworks, especially to account for unmapped faults, the similarity between fault seismicity rates and earthquake catalogue indicates that most of the recorded seismicity within the buffer is likely linked to the active fault structures considered. Nonetheless, the introduction of BG_2 also adds an additional layer of epistemic uncertainty, as the magnitude-dependent ratios adopted are based on expert judgment and remain inherently arbitrary. Given that our tests indicate a limited impact on overall model performance, we prefer to adopt the simpler BG_1 configuration, which avoids unnecessary complexity without compromising the quality of the results.

4.5 Model weighting and implications for probabilistic seismic hazard

The main challenges in fault-based PSHA applications adopting relaxed segmentation frameworks concern: (i) the definition of possible earthquake rupture scenarios and (ii) the estimation of their associated rates of occurrence. While the approach employed in this study does not aim to resolve the physical plausibility of multi-fault ruptures in the southern Apennines – which would require physics-based investigations – it provides a systematic framework to structure and explore epistemic uncertainties related to rupture configurations and earthquake rates within a PSHA context. In this framework, the approach does allow to narrow some of these uncertainties through consistency analyses with available observational data (instrumental seismicity and paleoseismic rates). At the same time, our hypothesis-based setup allows preserving the probabilistic nature of PSHA, as the results are intended to be implemented within a logic-tree scheme.

It is important to note that this study represents only a subset of the many rupture scenarios and seismicity parameters that could be explored for the southern Apennines and is therefore not intended to confirm or rule out other configurations. We also acknowledge the limitations in model quality check, as the available datasets – namely the regional seismic catalogue and the sparse paleoseismic records – may be insufficient to fully assess model performance.

Table 4 Summary table of the main outcomes discussed for each branch explored.

In this sense, the main outcome of this study is to provide a modelling framework that enables the use of geological-based criteria to assign weights to rupture hypotheses within a PSHA logic tree. The modelling framework we propose follows a logic-tree structure that can be directly implemented within PSHA engines such as OpenQuake, as SHERIFS is explicitly designed to explore epistemic uncertainties in fault parameters and rupture configurations and to propagate them into probabilistic hazard calculations. The weights derived for the different model branches remain region-specific, as they are based on comparisons with local observational constraints that are not directly transferable to other tectonic settings. However, the proposed weighting rationale and the general methodological framework can be used in other active regions worldwide.

Based on the performance of the different models across all consistency tests (summarized in Table 4), Set_1 emerges as the most robust configuration and should be assigned the highest weight in a PSHA logic tree, given its strong agreement with the catalogue MFD. Set_2, while slightly less aligned with the catalogue MFD, provides the best results in terms of NMS and paleoseismic fit, and should therefore receive a comparable weight. In contrast, Set_0 shows consistently poorer performance and limited agreement with observational data, and should be assigned a lower weight.

As for the other modelling branches explored, while their influence is less pronounced than that of the rupture scenarios, they can still serve as additional criteria for developing hazard models. Among the scaling relationships tested, Leonard (2010) is our preferred choice due to its consistent performance across scenarios. Wells and Coppersmith (1994) yields very similar results and can therefore be considered a valid alternative. In contrast, the relationship proposed by Thingbaijam et al. (2017) generally results in poorer model performance, with greater dispersion in the MFDs and higher NMS percentages.

Regarding the b-value, we find that adopting a broader range (0.9–1.1) results in better agreement between the modelled and observed seismicity, making it our preferred option. As for background seismicity, BG_1 (where all seismicity is attributed to the faults) is favoured, as it explains the observed seismicity well and yields model performances comparable to BG_2 for logic tree weighting purposes. Additionally, BG_1 avoids the added layer of arbitrariness required by BG_2, thus reducing unnecessary complexity in the modelling framework. Sensitivity tests performed using alternative background buffer distances (5, 10 and 15 km) show only minor variations in model performance, with the 10 km configuration representing the preferred compromise between spatial association of seismicity with mapped faults and the obtained results.

5 Conclusions

This study explores the impact of different variables on fault-based seismicity rate modelling in the southern Apennines, using a relaxed segmentation approach that allows multi-fault ruptures. The region, affected by the 1980 M_w 6.9 Irpinia earthquake, includes several critical infrastructures, highlighting the need for accurate seismic hazard assessment.

We computed seismicity rates from geological data in 35 seismogenic sources employing the SHERIFS algorithm through an exploration tree that combines different rupture scenarios and seismic parameter assumptions (e.g., scaling relationships, b-value, background seismicity, buffer area) to explore intrinsic epistemic uncertainties in fault data. We evaluated and quantified model performance by comparing the modelled seismicity rates with the regional earthquake catalogue data and paleoearthquake rates from paleoseismic studies.

Our results indicate that rupture assumptions strongly control the model outcomes: scenarios restricted to single-fault ruptures are systematically less consistent with observations than more complex rupture configurations. Multi-fault rupture scenarios are required within the explored modelling framework to reconcile geological and seismological observations at the regional scale. Importantly, these findings are consistent with evidence provided by the 1980 Irpinia earthquake, suggesting that similar rupture complexity cannot be neglected in long-term earthquake rupture forecasts for the southern Apennines.

Beyond the comparison of alternative rupture models, this study provides insight into the relative importance of different sources of epistemic uncertainty. In particular, rupture assumptions exert a stronger influence on model performance than background-seismicity assumptions, b-values and scaling relationships, while the contrasting behaviour of SubArea 1 and SubArea 2 highlights the importance of fault-system characteristics and geological constraints in controlling model robustness. In this sense, the study does not aim to reduce epistemic uncertainty, but rather to identify, structure, and evaluate the assumptions that most strongly affect fault-based PSHA results, thereby identifying priorities for future fault characterization and seismotectonic investigations in the Southern Apennines.

Nevertheless, this modelling has some limitations. The limited availability and quality of fault data, together with the necessary simplification of fault geometries, result in poorly constrained input parameters. In addition, the regional seismic catalogue and the limited number of paleoseismic constraints may not fully capture the long-term behaviour of the fault system. Finally, the modelling framework relies on statistical consistency with observations rather than physics-based rupture modelling, which could further refine rupture scenario selection and seismicity-rate estimates.

Addressing these aspects will be essential for the next generation of fault-based seismic hazard models in southern Italy, where the high seismicity and the presence of large populations and critical infrastructures demand accurate estimates of large-magnitude events, including those with long recurrence intervals.