Paleoearthquakes modelling and effects of uncertainties on probability assessment of next fault ruptures: the case of Central Italy surface faulting earthquakes

SUMMARY

Using the published paleoseismological trenching data for 16 faults in Central Italy, we compile a new data base of surface faulting earthquakes, having a quite stationary temporal distribution since 6000 BCE. By applying a probabilistic aggregation method, we correlate the event ages from distinct trenches on each fault, to construct all possible individual fault rupture scenarios, consistent with geological constraints. These inferred fault time histories are the basis for both individual fault and regional seismic hazard evaluation. We found that the mean recurrence time of each fault goes from about 1 to 4 thousand years for individual faults, whereas the value at regional scale is close to 120 yr. The small size of individual fault data samples does not allow us to infer straightforward information on the fault temporal behaviour, but only to evaluate the reliability of a chosen occurrence model for each fault. Therefore, hazard assessment is carried out by including the uncertainties related to both ages and probability distribution of the interevent times. We find that both these have a large impact on the probabilities of next rupture for individual faults: these depend on basic features of the temporal model and on the relation between the elapsed time and the mean interevent time. At a regional scale, we cannot exclude the simplest possible model, that is, the Poissonian behaviour, that provides quite stable probabilities of future events, close to 27 per cent in the next 50 yr.

1 INTRODUCTION

Paleoseismological studies provide essential information for understanding long-term rupture patterns of faults and for imposing important constraints on the distribution of strong earthquakes over scales of time and magnitude that are useful for seismic hazard assessment. This information is particularly crucial both for regions of low-medium seismic rate, such as Italy, where the average recurrence period of major earthquakes on a fault reaches thousands of years, and in areas with no or poor historical seismic data.

The processing and the analysis of paleoseismological data require a delicate balance between geological and statistical inferences, especially if finalized to seismic forecast and hazard assessment. In particular, robust statistical methods must provide, besides other things, an adequate treatment of all epistemic (and aleatory) uncertainties, such as those associated with ages of paleoearthquakes or related to the unknown temporal behaviour of faults. The identification and treatment of all sources of uncertainty is a very tricky question, especially because, in some cases, the paleoseismological data may depend on expert opinion, which is itself a source of uncertainty and of possible biased information. Moreover, whereas the uncertainties associated with ages of paleoearthquakes are generally addressed by means of Monte Carlo simulations (Ellsworth et al. 1999; Biasi et al. 2002; Parsons 2008; Parsons 2012), the identification and the treatment of epistemic uncertainty about the seismic behaviour of faults open practical questions, since the data for individual faults are, in general, not sufficient to identify a unique occurrence model (Ellsworth et al. 1999).

The simplest interpretation of paleoseismological data, commonly adopted in the literature, usually involves considering the average recurrence time as representative of fault ruptures behaviour, implicitly assuming that faults generate earthquakes with a quasi-periodic frequency and with similar magnitudes over multiple seismic cycles (Nishenko & Buland 1987; Ellsworth et al., 1999; Matthews et al., 2002). If one hypothesizes a broader range of possible recurrence and magnitude distributions, such as time clustering (Kagan & Jackson 1999), fault interaction (Stein et al. 1992; King et al. 1994) or variable magnitudes (Weldon et al. 2004; Jackson & Kagan 2006; Page et al. 2009; Parsons & Geist 2009), then interpreting paleoseismic information becomes more complicated, but probably closer to the behaviour of the natural systems. Such differing models should be reconciled by several factors, such as missing events in geological evidence, indistinguishability of two or more clustered earthquakes in the paleoseismic records or the different resolution between paleoseismological and seismological investigations.

A deeper knowledge of fault rupture mechanisms may be provided by multidisciplinary approaches. An example of these is the combination of expert opinion and statistical methods on geological data given by Cinti et al. (2021). They apply a new quantitative multistep method, based on the statistical analysis of trenching data, yielding rupture scenarios for individual faults, for fault systems, and for the whole region under investigation. By analysing two data sets of paleoearthquakes for 10 faults in the Central Apennines, built by clearing (‘Cleared’ data set) or by maintaining the author's expert opinions about the age and recognition of events (‘Further Constraints’ data set), Cinti et al. (2021) quantify the impact that age uncertainty in paleoseismic records has on the statistical reconstruction of the paleoearthquake histories and on the measurement of the seismic rates. They find that (a) the scenarios from the alternative data sets for individual faults appear quite stable with the rupture histories of the last 7 kyr and show a high degree of completeness for the reconstructed chronology of activity; (b) paleoseismological average regional recurrence agrees with the M6+seismic rate from the historical records; (c) the scenarios derived from both ‘Cleared’ and ‘Further Constrained’ data sets show peaks of potential multiple fault activation, alternating with more ‘quiet’ periods, suggesting cycles of fault ruptures having a time length of 2000–2500 yr, comparable with the timescale of typical individual fault loading.

This study is a follow up of Cinti et al. (2021) and has been carried out with the main aim of characterizing the earthquake occurrences and of constructing earthquake rupture forecasts, both for individual faults and at regional scale. To achieve this goal, we compile a new data set, named Central Apennines Paleoearthquakes (CAP), built directly on the basis of the paleoearthquakes recognition and ages expressed by the authors in their published scientific papers, without removing existing expert opinion but always disaggregating data from different trenches.

The use of faults data in earthquake forecasting or seismic hazard assessment is a well-established practice worldwide (Frankel 1995; WGCEP 2003; Stirling et al. 2012; Field et al. 2014; Woessner et al., 2015). In most of these cases, faults are modelled as independent sources following a quasi-periodic model (Youngs & Coppersmith 1985; Nishenko & Buland 1987) or, more generally, a time-dependent renewal process, for which hazard rate varies with time since last event (Lindh 1983; Sykes & Nishenko 1984; Cornell & Winterstein 1988; Ellsworth 1995; Ogata 1999). In cases where one suspects that the time-dependent models are too poorly constrained, the Poisson model, having a constant hazard rate and the memoryless property, provides a ‘reference’ probability calculation that reflects only the long-term rate of earthquakes (WGCEP 2003). In this study we consider two temporal models: a reference Poisson process, predicting a memoryless occurrence of fault ruptures, and a Brownian Passage Time model (Kagan & Knopoff 1987; Ellsworth et al. 1999; Matthews et al., 2002), capable of representing a wide variety of time-dependent behaviours, according to parameters values. Specifically, we compute the probabilities of occurrence for the next event, both for individual faults and at regional scale, and we show the effects of time-dependent and time-independent occurrence hypothesis on earthquake forecast assessment.

The paper is organized as follows. In the next section, we describe the set of CAP data, relative to 16 faults of Central Apennines. In Sections 3 and 4 we apply the procedure by Cinti et al. (2021) to produce rupture scenarios for each fault. Section 5 is devoted to earthquake forecast assessment, for individual faults and at regional scale, made by including all uncertainties related to both ages and rupture temporal models.

2 BUILDING THE CAP DATA SET

Active tectonics in Central Apennines is mostly accommodated by a set of normal faults, active since the Early Quaternary, having a slip rate up to 1 mm yr⁻¹ (Galadini & Galli 2000). This is one of the regions with the highest seismic hazard in Italy (Meletti et al., 2021) for which a wide documentation on the largest earthquakes of last centuries is preserved in the historical catalogue (Rovida et al. 2020).

The region considered in this study (Fig. 1) has been recently struck by two seismic sequences, comprising several destructive earthquakes: the 2009 L'Aquila (main event M_w 6.3; Chiaraluce 2012) and the 2016–2017 Central Italy (Amatrice event M_w 6.1; Norcia event M_w 6.5; Improta et al. 2019) sequences. The largest, documented, earthquakes in the last four centuries occurred at 1703 (Valnerina, M_w 6.9 and Aquilano, M_w 6.7) and at 1915 (Fucino, M_w 7.1), for a total of about 20 events with M_w > 5.5 in the past ∼400 yr (Rovida et al. 2020). Both historical and recent data agree in showing temporally and spatially clustered ruptures along individual faults and/or fault systems in the region (Wedmore et al. 2017).

A most ancient surface-faulting earthquakes, since the Late Pleistocene-Holocene, are documented in the large paleoseismology literature, which allowed us to extend our knowledge on the seismicity of this area, to identify the major active structures and to define their kinematic characteristics (Galadini & Galli, 2000; Cinti et al. 2021 and references therein). In particular, Cinti et al. (2021) built the M6 + rupture histories of 10 individual faults in Central Apennines, on the basis of two data sets of paleoearthquakes ages, differentiated for degree of uncertainties, with the aim to explore the sensitivity of ruptures pattern in time to data ambiguity. The records from published trench data were systematically analysed and reviewed, in order to minimize biases given by the author's expert opinion. The main steps of their review were:

1. identification of paleoearthquakes recognized and dated in individual trenches; application of a disaggregation procedure when data from individual trenches along the same fault were presented as merged;

2. compilation of the ‘Cleared’ data set, obtained by only considering ages of paleoearthquakes derived by laboratory measured ages (i.e. radiocarbon, TL, etc.) or robust stratigraphic correlations;

3. compilation of the ‘Further Constraints’ data set, by using further elements that provide clues to set paleoearthquake age ranges, such as relative dating of the stratigraphic units, archaeological evaluation on pottery findings or other remains, historical considerations, stratigraphic correlation with other nearby trenches, etc.

Cinti et al. (2021) show that the differences of possible scenarios of past earthquakes occurrence, obtained by the statistical analysis of these two data sets, appear negligible, given the uncertainties affecting the paleoearthquake ages, particularly for the older events.

In the light of these results, this work intends to investigate the differences that may result from the use of the native data, published by the authors, by applying only the first of the three steps listed above, that is, the disaggregation of data from different trenches. We consider this step mandatory for two reasons: first, to avoid biases related to the different weights that the authors may have unwittingly assigned to the results of each trench, and, secondly, to allow the automatic aggregation of paleoearthquakes on the same fault, by means of the sweep-line algorithm proposed by Cinti et al. (2021). This automatic aggregation, being independent from any author's personal interpretation, may, in some cases, produce rupture scenarios for the fault, different from that proposed in the published paper, although derived from the same trench data.

In this work, we extrapolate the native data for the same 10 faults investigated by Cinti et al. (2021) and for 6 new faults in the same area (Figs 1 and 2), including the records of paleoearthquakes published until the year 2021. The resulting CAP data set comprises 16 faults, 67 trenches and 198 individual trench paleoearthquakes. For each trench we have from 1 to 7 paleoearthquakes, with a median of three records. Full information on the trench sites and relative publications are provided in Table S1 in the Supporting Material.

In a preliminary step of our study, we test and verify that the results of statistical analysis, reached by Cinti et al. (2021) for 10 of the 16 faults, remain unchanged considering the large uncertainties related to the age of paleoearthquakes (see Appendix  A). Specifically, to evaluate the impact of the expert opinion on the results, we compare the rupture scenarios obtained from the CAP data set and those derived from the ‘Cleared’ data set (intended as the most unbiased data set) by Cinti et al. (2021), for the 10 faults (F1–F10). This analysis shows that the differences in rupture scenarios obtained from the two data sets, have an order of magnitude negligible compared to the events age uncertainties and allows us to use the CAP data set, with no need to deeply clear the published trench data.

3 DEFINITION OF INDIVIDUAL FAULT RUPTURE SCENARIOS

The first step of our study consists in the generation of a list of possible rupture histories for each of the 16 individual faults, by applying the multistep method, proposed by Cinti et al. (2021), on the ‘Original’ data set. This method searches for all possible intersections among age series of paleoearthquakes, found in different trenches along the same fault, by means of an ad-hoc developed sweep line algorithm (Bentley & Ottmann 1979; Cinti et al. 2021). The idea behind this algorithm is to imagine that a line moves along the time axis, stopping at and reporting intersections of paleoearthquake ages, as they are encountered (fig. 3 of Cinti et al 2021). The constraints of the algorithm are that (a) the event probability density function (PDF) is uniform over age ranges, (b) each identified fault rupture scenario satisfies all trench observations, in terms of number and temporal sequence, (c) each paleoearthquake contributes for only one specific age interval in each scenario, (d) two paleoearthquakes of the same trench cannot define the same event of a fault scenario and (e) for each scenario the intersections between observations of different trenches are maximized (i.e. the age ranges are defined by all paleoearthquakes belonging to an intersection encountered by the line; see Cinti et al. 2021 for details). The so-defined fault scenarios are, then, validated and selected by checking if they are consistent with all constraints provided by paleoseismic investigation (assuming that the whole fault ruptures at once).

Generally, the probabilistic aggregation processes provide a preferred correlation of trench event ages by (a) quantitatively measuring the amount of overlap between different site PDFs of paleoearthquakes and (b) by multiplying the probabilities of the site PDFs, to generate a final PDF for each earthquake (DuRoss 2011; Gómez-Novell et al. 2023). The aggregation approach used here has different characteristics. In fact, we prefer to adopt an alternative strategy, with the aim of creating as few restrictions as possible. First, we assign a uniform probability distribution to the maximum age range of the event, identified by the C14 calibrated data with two standard deviations of uncertainty. The inclusion of PDFs for limit dates (when they are known) undoubtedly would improve the statistical analysis of paleoearthquakes. Here, we consider the more general case in which the PDFs of calibrated dates are not published and, then, we do not include the ‘calibration’ probabilities in our analysis. The proper probabilistic way to represent a state of ignorance is a uniform PDF and so we assume it for the time of fault rupture, giving the same chance to the true numerical date of limiting deposit ages inside the calibration range. It is a reliable assumption in cases of historical or well constrained and not overlapping limiting ages, such as in most of our paleoseismological data (Table S1). Moreover, previous studies showed that the choice of PDF inside the 2σ uncertainty ranges of paleoearthquakes is relatively insignificant (Fraser et al. 2010). However, the adoption of a uniform distribution allows us to rule out any bias on results coming from more sophisticated, but subjective, PDFs, without giving up a probabilistic treatment of data. Indeed, when you impose subjective PDFs for the trench limit ages (such as, e.g. the normal distribution assumed by Benedetti et. al 2013 or by the PEACH software, Gómez-Novell et al. 2023), you do not ensure a more robust analysis, since you associate arbitrary probabilities to possible limit dates, that may significantly affect overall integration. Moreover, if you truncate the calibration ranges to 1 sigma level (Gómez-Novell et al. 2023), you are heavily reducing (of more than 30 per cent) the real uncertainty on the dates at the limits, acting on the whole scenario planning process. However, the event detection based on an overall PDF is not devoid of subjective choices, such the ‘preferred correlation’ of DuRoss et al. (2011) or the ‘prominence threshold of peaks’ of Gómez-Novell et al. (2023). Finally, even if the integration of event PDFs is correctly obtained and analysed, the inferred scenario represents the most probable one, not the only one possible. A proper account of uncertainties must consider all chances, not only the most likely ones.

By applying our aggregation method, we produce fault rupture scenarios obtained by all possible intersections between event ages. Then, we check that they respect initial chronologic and stratigraphic constraints. Finally, we apply a selection on this basis. Therefore, we do not build a scenario based on the highest overlap between event ages (DuRoss 2011) or on the picks identification of the mean PDF (Gómez-Novell et al. 2023), since this could lead to a misconstruction of fault histories, rather than to a reduction in event time uncertainties. Instead, we consider also apparently less likely configurations, which are, however, consistent with the geological constraints. As an example, we consider the chance that also highly overlapping trench ages refer to different earthquakes, differently from methods based on PDFs overlap. In this way, we fully consider the uncertainty on fault histories, without preferring scenarios that we subjectively think are more likely, but only ensuring consistency between selected scenarios and chronological and stratigraphic constraints.

By applying our procedure, we select a total of 33 regional scenarios, obtained by considering all possible combinations of scenarios of 16 individual faults (ranging from 1 to a maximum of 6, Table 1 and Fig. 3). These fault and regional scenarios allow us to obtain information on the past surface faulting earthquakes (M6+) in the studied region, including all the uncertainties affecting the paleoearthquakes ages.

4 CHANGE POINT ANALYSIS OF INTEREVENT TIMES FOR THE WHOLE CAP DATA SET

The combination of fault data collected in the CAP data set allows us to infer information about the temporal distribution of the interevent times (IETs) of the whole region. All statistical analyses shown and discussed in this and in following sections are based on simulations of individual fault rupture histories. These last are obtained by randomly choosing a scenario and by randomly selecting a time from each age interval of that scenario for each fault (Table 1).

To obtain and analyse the statistical measures on the temporal distribution of IETs we follow the method described in Cinti et al. (2021) and summarized below. First, we repeat the following steps 10⁶ times:

• for each fault, we simulate a rupture history for each fault;

• we obtain a regional rupture scenario by combining the 16 fault rupture histories;

• we compute the regional IETs for each couple of simulated occurrence times and combine the i-th IET with the midpoint T_i of two simulated occurrence times, to which it refers.

Then, we sort the times T_i obtained in all 10⁶ simulations together and we compute the median value (M_IET) and the width of 90 per cent confidence interval (Δ90 per cent_IET) of associated IETs, in windows of 10³ T_i values, moved forward by one, through the whole data base of simulations. These temporal series represent the average temporal behaviour and the related uncertainty of regional IETs, considering the uncertainty about the times of individual-fault ruptures.

The temporal trend of both M_IET and Δ90 per cent_IET shows a decrease from about 6000 BCE to nowadays (Fig. 4). The temporal fluctuations of largest IETs amplitude, on a timescale of centuries-thousands of years, may represent a period of lack of data or of lower seismic activity. Indeed, the more difficult definition of the stratigraphy and the unavailability of good datable material going back in time increase both the chance of missing events and the uncertainty of paleoearthquakes, due to the lower resolving power of events capturing. Therefore, the differences in the distribution of the oldest and newest paleoearthquakes, common to all faults (Fig. 2), cannot be neglected.

To rigorously identify turning points in the temporal behaviour of regional IETs (i.e. the decrease of data quality and of completeness, going back in time), a change point analysis (Lavielle 2005; Killick et al. 2012) is performed on the inferred temporal distribution of regional IETs. This is done by applying the (base 10) logarithmic scale to M_IET and Δ90 per cent_IET, given their large-scale fluctuations. A change point at about 6000 BCE is recognized for both M_IET (5830–5790 BCE) and Δ90 per cent_IET(5960–5920 BCE). It divides the history of IETs in two periods in which the median recurrence time (computed as the median of M_IET time histories) is close to 120 and 480 yr, for the more and the less recent periods, respectively (Table 2, Fig. 4). The analysis of Δ90 per cent_IET provides a second change point close to 200 BCE (240–210 BCE), after which the uncertainties about both geologic data and interpretative issues notably decrease. The median values of Δ90 per cent_IET for the three periods (from the most recent to the oldest), identified by two change points, are 200, 370 and 1380 yr.

From the change point analysis, we infer that the CAP data have a substantially lower degree of completeness and a smaller uncertainty since 6000 BCE. So, we decided to use for our statistical study only this younger part of the data set that can be considered the most robust for detecting rupture patterns in time. For faults F7, F12 and F16, we also include one paleoearthquake identified by an interval time over 6000 BCE (Table 1), but most likely occurred within this time. Note that, by excluding data older than 6000 BCE, the differences between scenarios may significantly decrease, up to vanish, as for the Montereale fault (F5; see Table 1). In Fig. 4(c) we show a zoom of the chronogram of scenarios identified by the Sweep Line method (Fig. 3) for the selected period.

5 RECURRENCE MODELS OF INDIVIDUAL FAULTS AND FOR REGIONAL SEISMICITY

One model for probabilistic seismic hazard analyses is the time-independent Poisson model (Petersen et al. 1996; Frankel et al. 1997), predicting a constant (in time) occurrence rate. This is considered appropriate when no information, other than the average seismic rate, is known. In case of studies of individual faults, the Poisson model provides a ‘baseline’ conservative estimate of probabilities of fault ruptures, particularly useful when the time-dependent features of fault behaviour are poorly constrained. From a mathematical point of view the adoption of this model only requires the computation of the mean of interevent times (μ). This can be estimated from slip rate and fault dimension, if the earthquake chronology is unknown, or directly from the recurrence intervals, if the timing of past faulting events is known (Mc Calpin 1996). In this last case, the uncertainty on events ages must be adequately addressed.

More complex probability models account for different theories, put forward by the scientific community, such as the periodic behaviour predicted by the elastic rebound theory (Reid 1910), the quasi-periodic recurrence of earthquakes (Bakun & Lindh 1985; Nishenko 1985), or the clustering model of earthquake supercycles (Sieh et al. 2008). Generally, all these models are statistically represented by renewal processes, for which the probability of occurrence depends on the time elapsed since the last event, such as lognormal (Nishenko & Buland 1987), Weibull (Hagiwara 1974) or Brownian Passage Time (BPT; Kagan & Knopoff 1987; Ellsworth et al. 1999; Matthews et al., 2002). If a time-dependent behaviour is supposed, the coefficient of variation, besides μ, plays a key role in characterizing earthquake recurrence in major tectonic faults, since it measures the degree of clustering or periodicity.

Faults with fewer than 10 IETs do not provide enough information to reliably determine the shape of the probability density function for fault earthquake recurrence (e.g. Ellsworth et al. 1999). However, it is possible that enough trenching-based IETs can be obtained on faults characterized by very high rates of activity. In the study area, the longest rupture history on an individual fault was obtained for the Paganica-S.Demetrio fault (F8) that has a maximum number of seven events (six IETs), since 6000 BCE. To overcome this limitation, we arbitrarily choose to employ two models as representative of the whole class of renewal models and of different possible time-dependent behaviours: the Poisson model, which yields time-independent rupture probabilities with a long-term rate μ, and the Brownian Passage Time (BPT) distribution, which is specified by an additional parameter, the aperiodicity α, coincident with the coefficient of variation and measuring the degree of clustering (α > 1) or periodicity (α < 1; α = 1 indicates an aperiodic process) of fault ruptures. The application of these models in paleoseismology is complicated by two factors, the small sample of paleoearthquakes for individual fault (N) and the event's age uncertainty, both leading to a significant bias in estimated parameters, particularly in periodicity (Kempf & Moernaut 2021). Therefore, both these points must be adequately dealt, to correctly infer information on fault behaviour.

The analysis of faults in the studied area of the Central Apennines, by Poisson and BPT models, is structured as follows. First, we attempt to extract useful knowledge on the earthquake recurrence mechanism, by applying well-known statistical measures on simulated histories of each fault, based on the available scenarios (Table 1 and Fig. 3). Then, we inspect both the overall variability and the sensitivity to data and model uncertainties of probabilities of future fault ruptures.

5.1 Statistics of recurrence intervals

The mean and the aperiodicity of IETs are the basic ingredients to compute earthquake probabilities by renewal models. The estimation of these parameters must be carried out with extreme caution for those faults with IETs affected by large uncertainty and not enough to characterize the interevent time distribution (Nomura et al 2011; Kempf & Moernaut 2021).

A common metric to measure the overall regularity of occurrences in a renewal process is the coefficient of variation (CoV), also known as aperiodicity (α), of IETs, defined as the ratio between the standard deviation σ divided the mean μ.

Establishing a rigorous nomenclature for the type of BPT models is challenging, but, by following Salditch et al. (2020), we define a BPT model as (a) ‘strongly periodic’ if 0≤α<0.5α (σ<μ/2); b) ‘weakly periodic’ 0.5≤α<1 (μ/2≤σ<μ); (c) ‘aperiodic’ α∼1 (σ∼μ); (d) ‘clustered’ if α>1 (σ >μ). The terms ‘strongly periodic’ and ‘weakly periodic’ indicate the degree of regularity of earthquakes occurrence. For clustered cases, the value of α affects also the intercluster distances, which are larger for higher α.

Parameter α is not sufficient to fully understand the origin of intermittency in the occurrence rate. Specifically, it ignores the order in which events occur and if shorter interevent times tend to cluster or are separated by longer interevent times. In order to limit this deficiency, Goh and Barabási (2008) introduced the memory coefficient, that is, a normalized form of autocorrelation calculated between consecutive interevent times. Unfortunately, the memory coefficient is useless in renewal modelling of small samples (Kempf & Moernaut 2021), such as the paleoseismology records.

A second benchmark, to investigate the temporal features of fault ruptures, is given by the Bayesian Information Criteria (BIC; Akaike 1974; Schwarz 1978)

where N is the sample size, L is the likelihood of the model and k is the number of parameters of the model, that is, 1 for Poisson (μ) and 2 for BPT (μ and α). It is a selection method, under the maximum-likelihood estimation framework, identifying as preferred the model with the minimum BIC value and it is particularly recommended when the number of parameters is high relative to sample size (more than 30 per cent). In our study, we evaluate the ability of BIC to identify the correct model, between Poisson and BPT, for each fault, by computing the quantity ΔBIC = BIC_BPT–BIC_POISSON (therefore ΔBIC > 0 and ΔBIC < 0 if the preferred model is Poisson and BPT, respectively).

Being aware that statistical modelling of a few data is highly uncertain and that statistical measures, such as α and BIC, are poorly informative for small data sets, we quantify the performance of these statistical tools by means of synthetic samples, having 6 data at most (the maximum number of IETs in our faults, since 6000 BCE), generated by Poisson and BPT models with different parameters. Specifically, we simulate samples, having a size N from 3 to 6, by Poisson and BPT models, predicting random, clustering, periodic and aperiodic behaviours, and estimate μ and α (by applying the maximum-likelihood estimator) and ΔBIC for each of them. The results, shown in Fig. 5 and described below, are independent of the value of μ used for simulations and, therefore, we fix μ =1000 yr, whereas α goes from 0.1 to 2.0, with a step of 0.1. By this analysis we infer some guidelines about the chances of correctly identifying the type of recurrence model, listed below.

• In cases of clear periodic BPT models (0.1< α <0.5), the estimated values of μ (marked by |$\hat{\mu }$| in the following) approach the true values also for small data samples. For N = 6 the relative width of the 90 per cent confidence interval goes from 5 per cent (α = 0.1) to 30 per cent (α = 0.5). The corresponding values for N = 3 go from 10 per cent to 50 per cent. On the other side, the values of |$\hat{\mu }$| may be strongly biased for Poisson, aperiodic or clustered models (Fig. 5a).

• The estimated values of α (marked by |$\hat{\alpha }$| in the following) are < 1 for both periodic and clustered models, whereas |$\hat{\alpha }$| > 1 is very unlikely for periodic models (probability < 1 per cent for α<0.5 and < 10 per cent for 0.5<α<0.7; Fig. 5b). Therefore, we infer two points: (1) values of |$\hat{\alpha }$| > 1 indicates the periodic behaviour as unlikely and (2) it is highly probable that small samples do not allow us to correctly identify a Poisson/aperiodic and a clustering behaviour. This is due to the less intuitive meaning of μ and to the larger correlation between |$\hat{\mu }$| and |$\hat{\alpha }$| for clustered models, respect to periodic BPT one.

• The ΔBIC scarcely identifies the right model for 3 ≤ N ≤ 6, when spanning from about −5 to 2 (shaded area in Fig. 5c), whereas smaller and larger values indicate BPT and Poisson, respectively, as the preferred model, at 90 per cent confidence level.

With these guidelines in mind, we explore the possibility of getting information on fault rupture behaviour, in the studied area of Central Apennines, by means of Poisson and BPT models. The temporal behaviour of individual faults is often represented by the mean of interevent times, μ, meant as a measure of the activity of a fault. In our case, this parameter is affected by two sources of uncertainty: the uncertainty on paleoearthquake times, from which it is computed, and the estimation error, which is as larger as smaller is the sample size (Kempf & Moernaut 2021). The first source of uncertainty may be investigated by estimating μ on 10⁶ simulated rupture histories, obtained by randomly choosing the scenario and a time from each age interval of that scenario, using paleoearthquakes since 6000 BCE (Table 3). We obtain values on average going from about 1 to 4 thousand years, but the uncertainty on these values depends on resolution of event ages (Table 1): we identify small confidence intervals of |$\hat{\mu }$| for some faults (F3, F4, F6, F8, F14), but their width may reach 1–2 orders of degree for the other faults (Table 3). Anyway, we must remember that the significance of |$\hat{\mu }$| depends on the (true) distribution of IETs: whereas it well represents the expected earthquake's ‘return period’ for periodic occurrences, it may misconstrue the true variability of earthquake interevent times for different temporal models (Fig. 5a).

To try to correctly interpret the temporal behaviour of individual fault ruptures, we estimate α and apply the statistical measure ΔBIC on faults having at least four earthquake ages (and, therefore, three IETs), for all scenarios and since 6000 BCE; these are: F3, F4, F6, F8, F11, F14, F15. This analysis is structured as follows:

• we simulate 10⁶ rupture histories for each fault, by randomly choosing the scenario and a time from each age interval of that scenario;

• for each simulated record, we estimate Poisson and BPT parameters, μ and α, and compute the measure ΔBIC;

• we try to infer information on the recurrence model of each fault, by comparing the estimated parameters with the range of values obtained from rupture history simulations (Fig. 5).

The empirical distributions of |$\widehat {\alpha \ }$| and ΔBIC for the selected seven faults allow us to infer a few, not conclusive, indications on their temporal behaviour (Fig. 6). For each fault and for each simulated time history, we compare the inferred values of |$\widehat {\alpha \ }$| and ΔBIC with the theoretical distribution, obtained for different BPT and for Poisson models by means of simulations (Fig. 5), in order to decide which of them may be considered as likely. Specifically, we verify if each estimated parameter is inside or outside the 90 per cent confidence interval given for that parameter by a model (i.e. if the model is rejected or not, on the basis of that parameter) and, therefore, we compute the percentage of rejections on all simulated histories. By this way, we build a reliability curve indicating the degree of consistency between a prefixed rupture model (Poisson or BPT) and fault data (including their uncertainty). Clear periodic behaviours (α<0.8) are partially (i.e. only for part of scenarios) rejected (at 90 per cent confidence level) for faults Fucino (F3), Monte Marine (F6), Paganica-San Demetrio (F8), Norcia (F14) and Aremogna-Cinquemiglia (F15), with a decreasing proportion of rejections as α approaches 1 (red lines, Fig. 6). In particular, a strong periodic model (α < 0.5) turns out to be very (i.e. for all simulated histories) unlikely for the Paganica-San Demetrio (F8) fault. The Trasacco-Luco fault (F4) stands out, among all faults, for the narrow range of |$\widehat {\alpha \ }$| values, close to 0.2, suggesting a periodic behaviour. About one half and two-thirds of scenarios provide values of |$\widehat {\alpha \ }$| and ΔBIC not explained by a Poisson and a clear clustered (⁠|$\alpha $|>1.2) model, respectively. The interpretation of results for the Vettore-Bove (F11) is particularly difficult, since values of |$\widehat {\alpha \ }$|and ΔBIC are partially inconsistent with any model, including Poisson.

In summary, the analysis of paleoearthquakes of individual faults shows that the reliability of a model is highly sensitive to uncertainties about the rupture histories. In any case, although results do not allow us to conclusively choose a model over the others, we may rank possible models for all investigated faults, basing on their reliability on simulated histories (Fig. 6).

5.2 Individual fault rupture probabilities: uncertainty and sensitivity analyses

The statistical analysis of paleoearthquakes, described in Section 5.1, do not provide conclusive information on the temporal behaviour of fault ruptures. This is not surprising considering the small number of paleoearthquakes available for individual faults (Ellsworth et al. 1999). Therefore, the impact of both recurrence model and data uncertainties must be explored, to achieve a credible probabilistic fault rupture hazard assessment.

Assuming a renewal behaviour of fault ruptures, represented by a recurrence model M, the most effective fault hazard indicator is the probability P_M(ΔT|T_L) that one or more earthquakes will occur during an interval of interest ΔT, conditional to the time from the last event T_L. The computation of this probability requests multiple input data: the occurrence model and its parameters and the last event occurrence time. Since all of these are affected from uncertainty, an effective analysis of P_M(ΔT|T_L) must include all possible choices of them.

Even if we have identified more likely models for some faults, here we prefer to give the same chance to each model for two reasons: (1) the examined seven faults show different, more likely, types of behaviour and this result deserves more investigation and (2) we have no indication to prefer a specific model for the remaining nine faults (with less of four events), not examined in the previous subsection. In light of this, the fault hazard analysis, applied on all 16 faults, is structured as follows:

• we redact a list of possible models M, for each fault, by considering different versions of Poisson and BPT models, obtained by varying both μ, in the empirical 90 per cent confidence interval (Table 3), and the aperiodicity α, uniformly from 0.1 to 2.0 (with a step of 0.1);

• we compute 10⁶ values of P_M(ΔT|T_L), by varying the model M and the last event T_L (if it is unknown), for ΔT = 50 yr. We randomly vary the scenario for Roccapreturo (F10) fault, the sole case in which T_L is different in two scenarios (Table 1).

The main products of this analysis are an overall uncertainty estimation of P_M(ΔT|T_L), including all sources of uncertainty (last event, model, parameters), and a sensitivity analysis, quantifying the impact of a single source of uncertainty on the probability P_M(ΔT|T_L).

The first result is that, if we include all sources of uncertainty, P_M(ΔT|T_L) mostly varies by several orders of magnitude and is picked around 1 per cent–10 per cent for all faults (Fig. 7). The sensitivity analysis of P_M(ΔT|T_L) to model M shows that the 16 faults may be divided in two classes, depending if they have largest probabilities for periodic (F1, F2, F5, F7, F9, F10, F12, F13, F15, F16) or clustered (F3, F4, F6, F8, F11, F14) models (Fig. 8). Faults of first class have probabilities that might approach 100 per cent, whereas the faults of the second class have probabilities not exceeding 20 per cent. This maximum probability threshold is valid for all faults, under the hypothesis of a Poissonian/aperiodic or clustered behaviour, with the exception of Pettino (F7), for which probabilities given by a Poisson model reach to 30 per cent. The division of faults in two classes is mainly given by the relationship between the mean interevent time and the occurrence of the last event. The higher probabilities for periodic faults are obtained in cases in which the last elapsed time is close to |$\hat{\mu }$|⁠, whereas the higher probabilities in case of a clustering behaviour are given for faults ruptured by recent earthquakes.

To clarify how the occurrence model drives the identification of faults most prone to rupture, we build three scenarios for all 16 faults, by doing following choices:

• we fix μ to the median of distribution obtained by simulated histories (Table 3);

• we choose three models: the Poisson, a weakly periodic BPT (with a = 0.5) and a clustered BPT (with a = 1.5);

• we choose the central time of the most recent paleoearthquakes as the occurrence time of the last event T_L, with the exception for the historical events. This choice was done only for calculation purposes and does not indicate a preference among the possible values.

Fig. 9 shows the probabilities P_M(ΔT|T_L) for all 16 faults and for ΔT = 50 yr, for the three models. We report a single value for Roccapreturo (F10) fault, since we obtain very close results by separately considering the two different values of T_L provided by two scenarios. The largest probabilities for the Poisson model are obtained for Fucino (F3), Paganica-San Demetrio (F8) and Sulmona (F16) faults, having the smallest estimated average interevent times, |$\hat{\mu }$| (Table 3). The periodic model identifies Monte Ocre (F2), Aremogna Cinque-Miglia (F15) and Sulmona (F16) as the most probable faults, since the chosen values of elapsed times are larger than the average interevent times. A particular attention must be paid for Sulmona (F16), having a well-constrained elapsed time (close to about 1900 yr), higher than 90 per cent of all possible values of μ (Table 3). The periodic model provides a low value (close to 0 per cent) for Trasacco-Luco (left histogram for F4 on Fig. 9), but it is worth considering that this fault was possibly only partially ruptured by the last event (1915 CE), sympathetically with Fucino (F3), and, therefore, forthcoming earthquakes are possible, also for a periodic behaviour. For this reason, we recalculate the probabilities P_M(ΔT|T_L) by considering the penultimate earthquake (508 CE) as the most recent one and finding a value of 5 per cent for a periodic BPT (right histogram for F4 on Fig. 9). Finally, the clustered model identifies Fucino (F3), as the fault most prone to rupture, considering the combination of the low average recurrence and the recent last event. As said above, all values of P_M(ΔT|T_L) reported in Fig. 9 were computed by choosing the values of T_Land μ. Different values of T_L and μ may change both the values and the order of probability values for different models, in agreement with what is shown in Figs 6, 7 and 8.

5.3 Regional rupture probabilities: uncertainty and sensitivity analyses

The information on the paleoearthquakes occurring on the 16 faults, allow us also to estimate the regional probability of future surface rupture earthquakes, given the uncertainties on data and temporal behaviour of individual fault ruptures. To this aim, we generate composite synthetic paleoseismic histories since 6000 BCE, by transferring the synthetic histories of ruptures on individual faults to a single timeline. In this way, we generate 10⁶ regional synthetic histories without considering the possibility of interfault triggering of earthquakes. Table 4 summarizes basic information on these simulated histories, on results of their statistical analysis by renewal models and on forecast computations. The 90 per cent of generated synthetic histories contain from 47 to 52 earthquakes. The 90 per cent confidence interval of median interevent time for the composite histories ranges from 90 to 170 yr. In order to test if the BPT model may be preferred to the Poisson one, we adopt two measures: ΔBIC and the Lilliefors test statistic. Both of them indicate the Poisson model as a reliable distribution for the regional seismicity, including all sources of uncertainty involved in dating determinations, since the percentage of cases in which they are out the 99 per cent confidence interval given by Poisson model is below 1 per cent and, then, compatible with randomness.

Since we do not have any evidence of a time-dependent, regional, earthquake occurrence, we adopt the Poisson model to compute the probability to have one surface faulting earthquake, at least, in the region, within ΔT yr. The resulting probabilities span in quite narrow ranges, close to about 6 per cent, 30 per cent and 50 per cent for ΔT = 10, 50 and 100 yr, respectively, starting from 2023. These probabilities do not take into account the contribution of possible unknown faults and, therefore, they may be a lower bound of the actual values.

6 DISCUSSION AND CONCLUSIONS

The overview of results reached in this study leads to the following issues.

• The method, adopted here and presented by Cinti et al. (2021), allows us to use the native trenches data for production of fault rupture scenarios, regardless of the subjective interpretation, provided by the authors, correlating paleoseismic data to earthquake timing. This is an important step towards a comprehensive knowledge of fault rupture mechanism and an objective regional earthquake probability assessment, based on paleoseismological information and the associated uncertainties. Results inferred without recognizing this uncertainty may potentially be flawed. Furthermore, a comprehensive uncertainty assessment must be as free as possible of any subjective judgment, to increase the transparency of data and to prevent erroneous conclusions.

• The CAP data after 6000 BCE, in the studied area of Central Apennines, show quite stable statistical properties and, then, provide a singular opportunity for statistical investigations on both fault ruptures and regional seismicity. Quantitatively assessing the completeness of paleoseismic data is a challenging issue, since it involves the completeness of trench stratigraphy, the preservation of evidence and the precision of geochronology. However, any assessment on fault rupture temporal features must include a completeness evaluation, to rule out bias in interpretation. Even if we cannot exclude that some events are unrecorded, the proposed data can be used as a baseline to infer reliable statistical properties of fault and regional seismicity behaviour.

• By means of simulated histories, we find that the mean interevent time of a single fault goes from about 1 to 4 thousand years, whereas it is about 120 yr at regional scale. This last estimation accords with the geodetic measurements (D'Agostino 2014), giving an average recurrence rate of 30–75 yr, for events with magnitude M_W 6.5+, in a 400 km-long section of the central-southern Apennines, comprising our study area (about 120 km-long section). A comparable seismic rate is obtained by complete historical data, giving 4 earthquakes with M_W 6.5 + in the last 500 yr (Rovida et al. 2020).

• The statistical measures applied on seven faults having four or more events (i.e. at least three IET), with the aim to detect a preferred temporal behaviour of individual faults, allow us to consider as unlikely a periodic behaviour for five of them, among which the Pagamica-San Demetrio (F8) having the largest number of known paleoearthquakes. This is not a meaningless result, since it shows that concepts such as ‘recurrence period’ might misinterpret the actual temporal behaviour of faults. On the other hand, this is not a novelty: despite geological observations supporting quasi-periodic recurrence on some faults (e.g. Berryman et al. 2012; Scharer et al., 2017), some studies, by analysing long‐term geological records, suggest more complex patterns of earthquake recurrence, which cannot be explained by the standard earthquake cycle model (e.g. Crone et al. 1997, 2003; Clark et al. 2015; Calais et al. 2016; Taylor‐Silva et al. 2019; Chen et al. 2020; Griffin et al. 2020). Many of these, in agreement with our results, conclude that the regular occurrence is observed for high activity‐rate faults, whereas non-periodic activity is typical for faults in slow deforming and intraplate regions.

• The building of a reliability curve between fault data and rupture models (Fig. 6) is a novelty for the statistical study of short (and uncertain) paleoearthquake records, that do not provide a sufficient basis on which to confidently prefer a model over the others. The assignment of a probability to each model opens new perspectives for a more focused knowledge of fault behaviour and a more constrained calculation of probability of future events that will be realized in future works. Differences in outcomes among faults probably result from different degrees of information, available for each of them, but might also highlight the difficulty or the impossibility to describe the recurrence of large earthquakes by a single model (Wang et al. 2024).

• The reliability of a specific recurrence model is not related to the uncertainty about the ages of events. For example, both F3 and F4 have quite narrow event age intervals (Table 1), but the reliability of each recurrence model is rather different (Fig. 6, solid red lines). On the other hand, the event ages of F15 are very uncertain, but the reliability of each recurrence model is like F3. The probabilities of next fault rupture mostly range from 1 per cent to 10 per cent for all faults (Fig. 7) and are below 30 per cent, in case of a non-periodic model (Fig. 8). Larger probabilities, which can approach 100 per cent, are obtained for periodic behaviour and for elapsed times close to mean interevent time. This proves that probabilities of next events are highly sensitive to model and parameter uncertainties, besides to the lack of information on fault histories. This should be considered when we build future fault rupture scenarios, since neglecting the epistemic uncertainties potentially leads to misleading inference.

These probabilities do not consider the possible effect of stress changes, as a consequence of the coseismic dislocation of nearby faults (Stein et al. 1992; King et al. 1994) or of viscoelastic response of the Earth's layers (Pollitz 1992; Piersanti et al. 1997; Verdecchia et al. 2018), quantitatively assessed through the evaluation of the so-called Coulomb Failure Function [Reasenberg and Simpson, 1992]. Therefore, a further variability of the earthquake occurrence for individual faults should be considered. However, the prospective use of such a kind of model showed that the increase in earthquake predictability may be weakened by a limited knowledge of the fault geometry and of the other source parameters (Woessner et al. 2011; Steacy et al. 2014).

• At a regional scale, we cannot exclude a Poissonian behaviour, that provides rather stable probabilities of future surface faulting events, close to 27 per cent for the next 50 yr. But, again, we must be aware that all findings obtained in this study results from the hypothesis of no interaction between faults. Addressing this chance could reveal new details on regional seismicity, but requires a careful evaluation of fault-interaction mechanism and of statistical techniques to highlight them and, therefore, will be the object of future work.

• The probabilities presented here are not intended to represent absolute estimates; they strongly depend on choices about the occurrence model and data treatment, but help us to understand the sensitivity of regional and fault hazard assessment to all sources of uncertainty. This may drive future research and set novel trends in this field.

ACKNOWLEDGMENTS

The authors are grateful to reviewers (G. Roberts and O. Gómez-Novell) and the associate editor (M. Segou) for their constructive criticisms and valuable comments, which were of great help in revising the manuscript and strongly improved the quality of the paper.

DATA AVAILABILITY

All paleoseismic data, used in this study, are provided in Table S1 of Supplementary Information.

REFERENCES

APPENDIX A: COMPARISON OF THE ‘CLEARED’ AND CAP DATA SETS BY THE CHANGE POINT ANALYSIS

To evaluate how much authors expert opinion affects paleoseismic data, we compare the ‘Cleared’ and the CAP data sets, by means of the change point analysis of interevent times (IETs), limited to the 10 faults (F1–F10) considered by Cinti et al. (2021). The change point analysis allows us to find times at which some statistical properties of a time-series, such as mean and variance, change abruptly. The method consists in the following steps (Lavielle 2005; Killick et al. 2012):

• choose a point that divides the series into two sections and compute the mean of the signal for each section;

• compute the sum of residual (squared) error from its local mean for both sections, to measure how much the property varies for each section;

• vary the location of the division point inside the series and choose the time that minimizes the total residual error as change point.

When the number of change points is unknown, you must add a penalty term to residual error to avoid overfitting. We set this term to 0 to return all possible change points.

For data with several fluctuations and varying in a broad range, it is unlikely to reliably discriminate between noise and signal and the method is prone to raise false change points. This may be resolved by applying a logarithmic scale to data, making the data easier to comprehend.

The temporal trends of both the median value (M_IET) and 90 per cent confidence interval (Δ90 per cent_IET) of IETs are very similar for two data sets and, as a result, the related change points are very close (Fig. A1 and Table A1). For both data sets we find a change point of both M_IET and Δ90 per cent_IET between about 6000 BCE and 4000 BCE. The average values before and after the change point of M_IET are 180 and 790 yr, for the cleared data, and 270 and 770 yr, for the original data, indicating a similar behaviour of two time-series. A second change point is detected for Δ90 per cent_IET, close to 300 BCE and 600 BCE, for the original and the cleared data, respectively, indicating a reduction of data uncertainties for the most recent event ages (Table A1).

By comparing all these results, some interesting features may be noted for the ‘Cleared’ data, such as the largest uncertainty about the change point of M_IET and the older age of the second change point of Δ90 per cent_IET (Table A1). All these features are explained by the larger age uncertainties of ‘Cleared’ data, with respect to the CAP ones, due to fewer constraints adopted during its compilation. Anyway, no significant difference is obtained for the median recurrence times. Therefore, given the high degree of uncertainty of paleoseismic data and the similarity of IET distributions, we consider the two data sets not significantly different for hazard assessment purposes.