Articles | Volume 11, issue 11
Geosci. Model Dev., 11, 4383–4397, 2018
Geosci. Model Dev., 11, 4383–4397, 2018

Development and technical paper 01 Nov 2018

Development and technical paper | 01 Nov 2018

Bayesian earthquake dating and seismic hazard assessment using chlorine-36 measurements (BED v1)

Bayesian earthquake dating and seismic hazard assessment using chlorine-36 measurements (BED v1)
Joakim Beck1, Sören Wolfers1, and Gerald P. Roberts2 Joakim Beck et al.
  • 1Computer, Electrical and Mathematical Sciences & Engineering (CEMSE), King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia
  • 2Department of Earth and Planetary Sciences, Birkbeck College, University of London, WC1E 7HX, UK

Correspondence: Joakim Beck (


Over the past 20 years, analyzing the abundance of the isotope chlorine-36 (36Cl) has emerged as a popular tool for geologic dating. In particular, it has been observed that 36Cl measurements along a fault plane can be used to study the timings of past ground displacements during earthquakes, which in turn can be used to improve existing seismic hazard assessment. This approach requires accurate simulations of 36Cl accumulation for a set of fault-scarp rock samples, which are progressively exhumed during earthquakes, in order to infer displacement histories from 36Cl measurements. While the physical models underlying such simulations have continuously been improved, the inverse problem of recovering displacement histories from 36Cl measurements is still mostly solved on an ad hoc basis. The current work resolves this situation by providing a MATLAB implementation of a fast, automatic, and flexible Bayesian Markov-chain Monte Carlo algorithm for the inverse problem, and provides a validation of the 36Cl approach to inference of earthquakes from the demise of the Last Glacial Maximum until present. To demonstrate its performance, we apply our algorithm to a synthetic case to verify identifiability, and to the Fiamignano and Frattura faults in the Italian Apennines in order to infer their earthquake displacement histories and to provide seismic hazard assessments. The results suggest high variability in slip rates for both faults, and large displacements on the Fiamignano fault at times when the Colosseum and other ancient buildings in Rome were damaged.

1 Introduction

A fundamental problem in earthquake science is the paucity of reliable earthquake records including multiple large-magnitude earthquakes on individual faults. This hinders a more advanced understanding of earthquake recurrence, which is a prerequisite to forecasting future earthquakes. A promising approach to address this problem is in situ chlorine-36 (36Cl) cosmogenic exposure dating of active normal faults (Zreda and Noller1998; Mitchell et al.2001; Schlagenhauf et al.2010). This approach is based on the fact that earthquakes progressively exhume bedrock fault planes and thereby expose the bedrock to an increasing amount of cosmic radiation, which is the dominant source of 36Cl production in rock. The resulting characteristic 36Cl concentration profiles along fault planes therefore provide information about the timing and intensity of past earthquakes.

A comprehensive mathematical model of 36Cl production was provided in Gosse and Phillips (2001) and later formed the basis of a MATLAB code that computes 36Cl concentration profiles from temporal sequences of ground displacements (Schlagenhauf et al.2010). Manual attempts to find best fits have subsequently been used for various faults (Benedetti et al.2002; Palumbo et al.2004; Schlagenhauf et al.2010, 2011; Yildirim et al.2016). Manual fits, however, can be deceived by local minima and do not impart information on possible alternative solutions and the resulting uncertainties in the inference of earthquake histories. Indeed, the complexity of the model and the abundance of uncertain parameters results in a highly nonlinear and non-convex problem, so that statistically reliable claims cannot be made without thorough modeling of prior beliefs, and proper accounting for parameter and model uncertainties.

To address these issues, Cowie et al. (2017) used a Bayesian Markov-chain Monte Carlo (MCMC) sampler (Robert and Casella2004) to generate ensembles of plausible solutions. However, the candidate space considered was restricted by the assumption of equally spaced and sized displacements in active slip time periods. Also, priors on the placement and number of active time periods were not discussed. Furthermore, their approach did not incorporate uncertainties in important parameters such as the 36Cl spallation and muonic production rates, the colluvial wedge mean density, which is difficult to measure in the field, and the timing of the demise of the Last Glacial Maximum (LGM) after which slip was preserved on the fault plane, which is only known imprecisely. Finally, their implementation of the MCMC algorithm was not fully described, and no evidence of convergence was shown.

The new Bayesian MCMC method proposed in the current work adopts the Brownian motion model for earthquake recurrence (Matthews et al.2002) to form a candidate space for earthquakes that arise and the associated prior probabilities. Besides improving the recurrence of prior earthquakes, this enables forecasts of future earthquakes, which is crucial for the subsequent task of regional seismic hazard assessment (Pace et al.2006, 2014). Furthermore, we forgo the unrealistic assumption in Cowie et al. (2017) that all displacements are of equal size and instead allow sizes to lie in a fault-dependent range. We employ parallel tempering (Woodard et al.2009) to avoid premature convergence to local rather than global optima due to a lack of explorative capabilities. We verified the correct implementation of our MCMC algorithm by the test of Talts et al. (2018) and measure global convergence with the diagnostic of Gelman and Rubin (1992). Since no model of 36Cl production can be completely accurate, we include an estimation of model discrepancy (Rougier2007; Rougier et al.2013; Brynjarsdottir and O'Hagan2014), which results in more realistic confidence bands and may help guide various research efforts to improve the modeling of 36Cl accumulation. We also account for uncertainties in important parameters of the model that were previously held fixed. More specifically, we assume uncertain production rates, attenuation lengths, and colluvium density, and we infer the time during the demise of the LGM at which the effect of erosion became outpaced by the ground displacements. We apply our algorithm to a synthetic case to verify the identifiability of past earthquakes in our model.

Finally, we study earthquake displacement histories for the Fiamignano and Frattura faults in the Italian Apennines. The results provide new evidence of slip-rate variability in normal faults in the Italian Apennines. The timing of slip-rate episodes can be reconciled with the historical record of earthquakes and damage to the Colosseum and other ancient buildings in Rome. Beyond proposing a new approach to 36Cl earthquake dating, we extend our model to allow for Bayesian regional probabilistic seismic hazard assessment.

We supplement this paper with an easy-to-use MATLAB code of the proposed Bayesian MCMC method for earthquake dating and regional probabilistic seismic hazard assessment.

2 Simulation of 36Cl concentrations along fault scarps

Typical faults scarps are shown in Fig. 1a–c. These are characterized by two slopes that were originally joined, but are now offset across a geological fault due to surface slip events during earthquakes. The original planar slope (see Fig. 1d) was formed during the LGM, which for southern Europe ended around 20 000 to 12 000 years before present (20–12 ka), through intense erosion on the upthrown side of the fault that was subject to freeze-thaw action (frost shattering) and sedimentation of the liberated slope debris (colluvium) on the downthrown side of the fault. Figure 1d–f also show how the morphology of the scarp changes through time across the LGM to post-glacial transition.

Before the demise of the LGM (see Fig. 1d), erosion via freeze-thaw action rapidly removed the surface uplifted by earthquakes (Allen et al.1999; Peyron et al.1998). During the demise of the LGM, decreasing erosion rates may have allowed fault slip to be preserved as a scarp, depending on the relative rates of erosion and fault slip. However, after the demise the LGM (Fig. 1e), the slip was preserved due to low erosion rates relative to the fault slip rates (Roberts and Michetti2004; Cowie et al.2013, 2017). The type of slip is illustrated in Fig. 1h, which shows decimeter-scale slips produced during two earthquakes at the Mt. Vettore fault on 24 August 2016 (Mw 6.1) and on 30 October 2016 (Mw 6.6). The Mw 6.6 surface rupture formed in 2–4 s during coseismic slip, recorded by Global Navigation Satellite System receivers placed either side of the fault before the earthquake (Wilkinson et al.2017).

The preserved fault plane can be sampled above ground and in excavated trenches. The 36Cl concentration in a fixed sample of limestone bedrock evolves according to

(1) d g ( t ) d t = ϕ ( t ) - λ g ( t ) ,

where g(t) is the 36Cl concentration (atoms g−1), ϕ(t) is the production rate (atoms g−1 yr−1), and λ is the decay rate of 36Cl (yr−1). The main 36Cl production pathways in pure limestones are spallation of 40Ca, muon capture by 40Ca, and thermal neutron capture by 35Cl (Marrero et al.2016). These processes depend on cosmic radiation flux, which is attenuated by the surrounding environment composed of air, colluvium, and rock (see Fig. 1g). As a consequence, the production rate ϕ(t) in a rock sample through time is strongly influenced by earthquake-induced changes to the surrounding environment. Roughly speaking, samples taken from the trench have experienced 36Cl production at low rates due to shielding by overlying colluvium and neighboring limestone bedrock (see Fig. 1g), whereas aboveground samples have experienced early subsurface production before exhumation as well as subsequent surface production. This results in a characteristic 36Cl concentration profile along the fault plane, which captures its history of ground displacements. Note that samples must be taken from a fault plane away from post-glacial alluvial fans and eroded gullies where exposure is also influenced by post-glacial erosion and sedimentation (see Fig. 1f)

Figure 1Images and illustrations of scarp evolution and radiation environment.


Formulae for surface production through the aforementioned processes as well as the associated attenuation, or shielding, factors can be found in Gosse and Phillips (2001) and Schlagenhauf et al. (2010). More than a decade of contributions to the modeling of 36Cl production (Lal1988; Phillips et al.1996; Mitchell et al.2001) have been assembled into a MATLAB code that computes 36Cl concentrations in a set of bedrock samples for given sequences of scarp displacements and event times; the code is provided in the supplement of Schlagenhauf et al. (2010). The code solves Eq. (1) by a first-order finite-difference scheme, and uses various simplifications for calculating the shielding factors; for example, it assumes that the cosmic ray flux decays exponentially with the depth of a sample beneath the colluvial wedge.

As part of this work, we provide a MATLAB code that circumvents some of the approximations and simplifications of Schlagenhauf et al. (2010), in particular in the computation of shielding factors. We use exact solutions of Eq. (1), piecewise exponentials, which is possible under the assumption of a constant-in-time production rate. These changes improved the predicted 36Cl concentration by around 5 % in our numerical experiments. Furthermore, we introduce an offline phase during which we pre-compute a database of shielding factors for a sparse grid (Barthelmann et al.2000) of possible values of the model parameters. By interpolating between these factors, we were able to accelerate the computations during the inverse problem by 2 orders of magnitude. Finally, an improvement of particular relevance to our case studies is that we also consider events before the end of the LGM. We approximate the effects of the LGM on erosion in a binary manner, by assuming a single point in time, Tinit, before which erosion immediately eroded scarps formed by earthquakes, and after which erosion stopped completely.

The provided MATLAB code calculates 36Cl concentrations in a set of bedrock samples for given sequences d=(d1,,dN) and s=(s1,,sN) of scarp displacements and event times and the following fault site properties.

  • Geometric description (see Fig. 1g):

    • dip of the lower slope, α ();

    • fault plane dip, β ();

    • dip of the upper slope of the footwall, γ ();

  • Geological properties:

    • rock mean density, ρrock (g cm−3);

    • colluvial wedge mean density, ρcoll (g cm−3);

    • spallation production rate of 36Cl by fast secondary neutrons at the surface from 40Ca, Ψsp (atoms g−1 yr−1);

    • slow muon capture production rate of 36Cl at the surface from 40Ca, Ψμ (atoms g−1 yr−1);

    • neutron apparent attenuation length, Λsp, and muon apparent attenuation length, Λμ, for a horizontal unshielded surface (g cm−2);

    • chemical compositions (ppm) of the rock in each sample and of the soil and pebbles in the colluvium.

  • Further properties:

    • the influence of the geomagnetic field is specified through scaling factors for fast neutrons and slow muons, see the supplement of Schlagenhauf et al. (2010);

    • the time during the demise of the LGM at which the effect of erosion became outpaced by the ground displacements, Tinit (yr−1).

3 Bayesian inference of displacement histories using MCMC

In this section, we present a Bayesian MCMC algorithm that solves the inverse problem of the previous section, i.e., the inference of past earthquakes from 36Cl measurements y=(y1,,yM) at M sample sites along a single fault scarp. More specifically, we compute posterior distributions for the vectors of fault displacements and event times, d=(d1,,dN) and s=(s1,,sN), respectively (and in particular for the number of events, N). In addition, we treat the values of z:=(Tinit,Ψsp,Ψμ,Λsp,Λμ,ρcoll) as uncertain parameters and include them in the inference problem.

Denoting the vector of true 36Cl concentrations by y*, and the output of the computer model of Eq. (2) for given values of x:=(d,s,z) by g(x), we first observe the equation

(2) y = y * + ϵ = g ( x ) + δ + ϵ ,

where ϵ is the measurement error, i.e., the discrepancy between the true and the measured concentrations, and δ:=y*-g(x) is the model error, i.e., the discrepancy between the model output and the true concentrations (Rougier et al.2013; Brynjarsdottir and O'Hagan2014). We assume that the measurement errors are independent and normally distributed, ϵN(0,σ2) for a given vector σ=(σ1,,σM) of positive real numbers. The model error is assumed to be proportional to the measurements, δN(0,(ρy)2), and we include the value of ρ>0 in the inference problem.

To describe the Bayesian inference method, let us denote the vector of all unknowns by θ:=(x,ρ,ϕ), where ϕ contains auxiliary variables that are described in Sect. 3.1 below. If we have a prior belief on the value of θ, described by a probability distribution Pθ, then the posterior distribution for θ given the measurements y can be found using Bayes' rule,

(3) P θ | y P y | θ P θ ,

with the distribution of y conditioned on a fixed value of θ, or likelihood, given by


To gain information about the posterior distribution Pθ|y, we employ a MCMC method (Robert and Casella2004), which generates samples that, roughly speaking, behave as if they were drawn from the posterior and thus can be used to approximate statistical properties thereof. More specifically, we use a Metropolis–Hastings MCMC method, which generates a sequence (chain) (θk)k=1 of random samples, where an initial sample is taken from the prior distribution and each successive sample is generated from its predecessor by means of random proposal functions and an acceptance step that guarantees that the distribution of θk converges to the desired distribution Pθ|y as k→∞. To accelerate this convergence, we employ a parallel tempering approach that simultaneously generates L>1 chains (θk(l))k=1, 1lL with progressively flattened likelihood distributions,


and randomly swaps states between neighboring chains such that the resulting samples of the first chain (θk(1))k=1, which uses κ1=1, are still distributed according to Pθy (in the limit). This has been shown to accelerate the exploration of the state space in cases where the posterior distribution has multiple local maxima (Woodard et al.2009).

To fully specify our inference method, it remains to describe the prior distribution Pθ and the proposal functions that are used for sample generation.

Figure 2Sample path of stochastic process used for earthquake event time generation. Sample path generated with X0=0.16, Tmin=-30 kyr, t1=-10 kyr, ν=(15 kyr, 3 kyr), τ=0.5, and Δt≈0.007 kyr. The corresponding earthquake times are represented by dots on the time axis, the switch point t1 is represented by a red vertical line.


Figure 3Measured and sampled 36Cl concentrations for the synthetic case with 146 36Cl measurements


3.1 The prior distribution Pθ

We describe the prior distributions of different components of θ with the understanding that separately described components are assumed to be stochastically independent.

For the parameter ρ, which controls the relative model error, we assume a uniform prior distribution, ρU[0,ρmax=0.1].

To describe a prior for the earthquake times s (and in particular the number of events N), we extend our state space by a stochastic process in [-Tmin,0] that is used to model earthquake occurrence (see Fig. 2). In doing so, we follow Matthews et al. (2002), where the time between successive earthquakes is modeled by an inverse Gaussian distribution, which is also called the Brownian passage-time (BPT) distribution since it describes the time required for a Brownian motion to reach a certain threshold. To account for large-scale time periods with differing average earthquake frequencies, we extend the model in Matthews et al. (2002) by including a Poisson (Λ Tmin)-distributed number J of switch points (tj)j=1J, where tjU([-Tmin,0]) and ΛU([10-4,10-3]). This is equivalent to the assumption that the times between successive switch points are exponential random variables whose mean 1∕Λ is inferred in the interval [103,104]. We then define the drift a(t) and volatility b(t) of the Brownian motion in the intervals of the resulting partition of [-Tmin,0] by


where the average interarrival times ν=(νj)j=0J are independent and identically distributed according to an inverse gamma distribution with mean mU([200,2000]) and shape parameter αU([1,10]), and where τU([0,1]) controls the short-scale recurrence variability (with the above definitions, the standard deviation of interarrival times in the jth subinterval is τνj). Numerically, given values of ν and τ, we use forward Euler–Maruyama time stepping (Higham2001) with time discretization Δt≈15 (yr−1) to simulate the resulting stochastic process on [-Tmin,0], with the modification that we reset the process to 0, and generate an earthquake by extending d and s, each time it reaches the threshold 1. In formulae, we let

(5) X ^ i + 1 := X i + a ( i Δ t ) Δ t + b ( j Δ t ) W i + 1 , X i + 1 := X ^ i + 1 if  X ^ i + 1 < 1 , 0 otherwise.

The values WiN(0,Δt) together with ν and τ and the initial value X0U([0,1]) form the vector of auxiliary variables ϕ referred to above, which fully determines the process (Xi)i=0Tmin/Δt, which in turn determines the vector of earthquake times s.

For the prior distribution of the displacement vector d conditioned on the number N of earthquake events, we use a uniform prior on the hypercube [dmin,dmax]N conditioned on the requirement that n=1Ndn=Hsc, where Hsc (see Fig. 1) is the present height of the fault scarp, and dmin and dmax are fault-dependent bounds on displacement sizes.

Figure 4The synthetic case with 146 36Cl measurements: the accumulated displacement (a), the mean earthquake intensity (b), and an event scatter plot (c) showing true events (blue) and events from posterior samples (with gray scale to indicate frequencies).


Figure 5Gelman–Rubin diagnostic for the synthetic case with 146 36Cl measurements


Finally, we assign prior distributions to the components of z: TinitU([12000,20000]), ΨspN(48.8,1.7), ΨμN(190,19), ΛspU([180,220]), and ΛμU([1300,1700]), where the first three are based on results of Allen et al. (1999), Stone et al. (1996), and Heisinger et al. (2002), respectively, and the remaining are chosen to be uniform around values taken from Schlagenhauf et al. (2010). It is difficult to accurately measure the value of ρcoll in the field due to compaction of the sediment during excavation of the sample trench at the base of the fault scarp, so we adopt the relatively wide prior ρcollU([1.2,1.8]) for all case studies considered in this work.

Figure 6Measured and sampled 36Cl concentrations for the synthetic case with 16 36Cl measurements


Figure 7The synthetic case with 16 36Cl measurements: the accumulated displacement (a), the mean earthquake intensity (b), and an event scatter plot (c) showing true events (blue) and events from posterior samples (with gray scale to indicate frequencies).


3.2 MCMC proposal functions

To explore the state space, we design a number of proposal functions and apply a subset of these before each rejection step in the MCMC algorithm. For each component of z as well as for the values of ρ, ν, τ, and X0, we include a global proposal from the corresponding prior distribution as well as a local proposal based on a normal distribution around the current value. To propose the partition of [-Tmin,0] and the corresponding piecewise constant drift and volatility coefficients, we employ reversible jump MCMC (Green1995), which allows for the application of MCMC to variables whose state space contains subspaces of different dimensionality. To propose new values of the variables (Wi)i=1Tmin/Δt, which drive the process (Xi)i=1Tmin/Δt, we again use a global proposal, which redraws all values independently, as well as a Brownian bridge-type local proposal that redraws the values within a random subinterval of [-Tmin,0] from their prior distribution conditioned on maintaining their sum. If a proposal changes the number of earthquakes, the earthquake displacement vector d (which then has a different size) is re-sampled too, and if a proposal changes the number of switch points, the vectors of drift and volatility coefficients are resampled as well.

Figure 8Location map of the central Apennines and the Fiamignano and Frattura sample sites. Holocene active faults and historical ruptures adapted from Roberts and Michetti (2004), Pace et al. (2006), Cowie et al. (2017), and Mildon et al. (2017).


4 Application to synthetic 36Cl data

In this section we apply our Bayesian MCMC method to synthetic 36Cl data. To generate these data, we drew d and s from the prior distributions described in Sect. 3.1 with dmin=10 (cm) and dmax=110 (cm) and applied the computer model of Sect. 2 using the chemical compositions of 146 rock samples from the Fiamignano fault (Cowie et al.2017) and the values z=(Tinit=-19000, Ψsp=49.5,Ψμ=200, Λsp=195, Λμ=1700,ρcoll=1.6). The remaining parameters, which are not part of the inference problem, were chosen as Hsc=2705, Htr=115, ρrock=2.7, α=23, β=42, and γ=33, based on the true values of the Fiamignano fault that were measured in the field. Finally, we perturbed the 36Cl concentration values according to Eq. (2), with standard deviation σ=0.025g(x) for the measurement error and ρ=0.03 in the model error term. The realizations of d and s are given in the Supplement.

We ran our MCMC method using the priors specified in Sect. 3.1 for the uncertain parameters, i.e., TinitU([12000,20000]), ΨspN(48.8,1.7), ΨμN(190,19), ΛspU([180,220]), and ΛμU([1300,1700]), ρcollU([1.2,1.8]).

The results presented in this section are based on two independent MCMC chains, each consisting of L=20 parallel tempering levels. We performed 1 374 462 MCMC iterations of each chain, resulting in a combined number of 2 474 032 samples in the first levels of the two independent chains after a 10 % burn-in period. Among these samples were 843 735 distinct scenarios, whose repetitions correspond to rejected proposals and reflect their statistical weight. Finally, to accelerate post-processing and to decrease memory consumption, we only saved each fifth scenario together with its number of repetitions (thinning).

Figure 9Damage to the Colosseum in Rome.


Figure 3 shows that the medians of the 36Cl concentrations of the posterior samples provide a good fit to the synthetic measured concentrations. More importantly, Fig. 4a shows that the posterior median of the accumulated displacement is close to the true value throughout the entire time span. This indicates that past earthquake activity can be recovered despite measurement errors (here 2.5 %), model errors (here 3 %), and parameter uncertainties. Figure 4 (b, c) shows the posterior mean earthquake intensity, i.e., the average displacement per time (computed for bins of width ∼500 years), and a scatter plot of all earthquake events of the posterior samples.

In both plots, the erosion during the LGM has removed information and so leads to almost uniform posterior intensities across time and through event sizes during the LGM, which means that individual events in that period cannot be recovered. Nevertheless, including such events in our approach can be considered a way to obtain realistic prior 36Cl concentrations at the end of the LGM, Tinit, as opposed to the alternative to start with zero concentrations or imposing ad hoc pre-exposure times (Schlagenhauf et al.2010).

Figure 5 shows the convergence diagnostic of Gelman and Rubin (1992). The potential scale reduction factor (PSRF) associated with the accumulated displacement values of Fig. 4 is less than 1.08.

Since the collection and chemical analysis of 36Cl samples involve time-consuming and costly fieldwork and lab work, in particular the accelerator mass spectrometry (AMS) analysis, it is worthwhile knowing whether fewer than 146 samples can provide similar insights. To answer this question, we repeated our computations with a subset of 16 rock samples (see Fig. 6).

Figure 10Fiamignano fault: measured and sampled 36Cl concentrations.


Figure 11Fiamignano fault: the accumulated displacement (a), the mean earthquake intensity (b) and an event scatter plot (c) showing events from posterior samples (with gray scale to indicate frequencies).


Figure 7 shows that this reduced dataset leads to similar results as the complete dataset (cf. Fig. 4). These results are based on 421 022 MCMC iterations of each chain and using the same number of chains, burn-in, and thinning as before, which led to a maximal PSRF value of 1.02.

5 Application to normal faults in the Italian Apennines

The Italian Apennines contain many examples of bedrock scarps on active normal faults and it has been suggested that rates of slip produced by repeated earthquake rupture since the LGM can be used to investigate seismic hazards and the mechanics of continental deformation (Piccardi et al.1999; Roberts and Michetti2004; Cowie et al.2013). The active normal faults work to extend the continental crust at the present day in formerly thickened crust of the Alpine–Apennines collision zone. The extension started at 2–3 Ma, as evidenced by dated sediments in extensional sedimentary basins formed by fault activity (Cavinato and Celles1999; Roberts et al.2002). The basins occur on the downthrown side of the faults, and contain accumulations of sedimentary layers that are usually hundreds to a few thousand meters thick. When added to the offset implied by the uplifted mountains on the upthrown side of the faults, total fault offsets are up to 2–2.5 km, which means long-term rates of vertical motion across the faults are in the order of 0.1 cm yr−1.

In this section we apply our method to the Fiamignano and Frattura faults of the Italian Apennines, and show its applicability to regional probabilistic seismic hazard assessments.

Figure 12Fiamignano fault: detailed view of posterior displacement history and earthquake intensity. The vertical blue lines show the timings of two events that damaged the Colosseum.


5.1 Fiamignano fault

The Fiamignano fault is located on the southwest flanks of the Apennines, approximately 60 km northeast of Rome. Based on a study of 34 damaging earthquakes recorded in Guidoboni et al. (2007), Galli and Molin (2014), and Vittori (2015), intensity VIII damage is expected at distances of less than ∼40–60 km of an epicenter from a Mw∼6.0–6.5 event, and the Fiamignano fault is one of the few Apennine normal faults with that epicentral distance to Rome (see Fig. 8) thus making it a plausible source of the historic earthquake damage.

Figure 13Frattura fault: measured and sampled 36Cl concentrations.


Figure 14Frattura fault: the accumulated displacement (a), the mean earthquake intensity (b) and an event scatter plot (c) showing events from posterior samples (with gray scale to indicate frequencies).


Figure 15Frattura fault: detailed view of posterior displacement history and earthquake intensity.


Rome was damaged in earthquakes of intensity VIII during at least three earthquakes in the 5th, 9th, and 14th centuries. Other earthquakes damaged the city in AD 801, 1091, 1231, 1279, 1298, 1328 (Guidoboni et al.2007). The Colosseum was damaged in AD 484 or 508 based on a stone inscription where a person known as “Decius Marius Venantius Basilus,” a prefect of the city, declares that he directly paid for restoration works after an earthquake; the uncertainty in age is because two Decius' who were consul, one in AD 484 and another in AD 508. Damage consisted of collapse of the colonnades in the summa cavea (upper seating section for plebian spectators) with major damage to the arena and podium. Collapse of the outer rings of the Colosseum is sometimes attributed to an AD 847 earthquake, when the nearby church of Santa Maria Antiqua was abandoned after earthquake damage (Vittori2015). An earthquake also damaged Rome on 9 September AD 1349, one of three large earthquakes in the Apennines on that day. This earthquake is thought to be linked to the Fiamignano fault based on the observation that taxes were reduced in the vicinity of the Fiamignano fault due to misery and depopulation after the earthquake (Guerrieri et al.2002). Damage and collapse reported in 14th century accounts are summarized in Galli and Molin (2014), where the earthquake on that day is described as “the strongest seismic shaking ever felt in Rome” (Galli and Molin2014). It is postulated that there was an abrupt collapse of the southern external ring of the Colosseum during this earthquake, as observed in Fig. 9, because a bull fight was hosted there in AD 1332, suggesting the Colosseum was intact, and an announcement that the collapsed stones were for sale was made in the second half of the 14th century. Furthermore, an intact Colosseum can be seen on a coin from AD 1328, whereas a damaged Colosseum is displayed on a 15th century image (Galli and Molin2014). Minor damage to the Colosseum also occurred in an AD 1703 earthquake (Vittori2015).

One of the goals of this work is to investigate whether the Fiamignano fault is a candidate for the historical damage accounts in Rome based solely on 36Cl analysis. 36Cl data from the Fiamignano fault have previously been collected and analyzed in Cowie et al. (2017), where geomorphic and structural field mapping as well as laser and radar datasets was presented as evidence that their site was exhumed by tectonic slip as opposed to exhumation by erosion. Cowie et al. (2017) tried to infer an earthquake history for the fault, and the results suggest rapid displacement between approximately 2000 years ago and AD 1349. However, the results were biased by the modeling constraint that the most recent earthquake was enforced to be at AD 1349 and by the use of constant inter-event times between slip-rate change points and constant displacement sizes.

Figure 16Posteriors of next earthquake time for Fiamignano (a) and Frattura (b) according to BED with dmin=10 and dmin=50, and according to a hybrid BED-BPT approach (with dmin=10 in the BED part).


In our analysis, we use the same site parameter values and priors as in Sect. 4 except that we now use dmax=300, motivated by the large displacements for similar faults presented in Wells and Coppersmith (1994). We performed 1 311 016 MCMC iterations of each chain and used the same number of chains, the same burn-in, and thinning as in the synthetic case. A maximal PSRF value of 1.03 using two independent MCMC chains indicates convergence. The minimal weighted root mean squared error (wRMS) among all scenarios was 1.42. The time required for the offline calculations was ∼15 min and the total time was ∼12 days. In this work, we performed all numerical experiments on Intel Xeon E5-2680 v2 at 2.80 GHz (20 cores) with MATLAB (2016a).

The sampled posterior 36Cl concentrations fit well to the measurements as observed in Fig. 10.

Our results agree with Cowie et al. (2017) in that we find clear evidence of slip-rate variability and the slip rate increasing through time before cessation of slip in the last 600–700 years. Although the present 28 m offset in the plane of the fault occurred at an average slip rate of ∼0.2 cm yr−1, we find rapid slip of 1–1.4 cm yr−1 centered around 1 ka (see Fig. 11). A detailed view for the period from 5 ka to present is shown in Fig. 12, where we observe that both of the two major earthquake events to have been observed in AD 847 and 1349 fall within this region of high intensity, see Fig. 12, and potentially large displacements indicated by dark gray regions in Fig. 11c. In fact, even the AD 801, 1091, 1231, 1279, 1298, and 1328 earthquakes fall within the high-intensity region at times when our findings suggest relatively low displacement sizes. This deserves further investigation, and we hope it leads to paleoseismic studies of offset late Holocene sediments on the Fiamignano fault that may be able to verify or refute these possibilities.

In conclusion, our findings suggest that slip is highly clustered in time and that the Fiamignano fault is a plausible source of the AD 847 and 1349 earthquakes associated with damage to the Colosseum and other ancient buildings known from the historical record.

5.2 Frattura fault

36Cl data for the Frattura fault have been collected and analyzed in Cowie et al. (2017), where geomorphic and structural field mapping as well as laser and radar datasets was presented as evidence that their site was exhumed by tectonic slip as opposed to exhumation by erosion. Unlike for the Fiamignano fault, the 36Cl data were collected sparsely (15 samples), similar to the situation in our synthetic case with a reduced amount of data shown. The site-specific parameters for the Frattura fault are α=25, β=53, and γ=28; fault scarp height Hsc=1570, trench depth Htr=130; and the rock mean density ρrock=2.7. Here the maximum displacement size is again chosen conservatively as dmax=300; though we mainly expect displacement sizes less than 100 cm based on the data of Wells and Coppersmith (1994). As in the previous cases, we let ρcollU([1.2,1.8]).

There are no known historical earthquakes for this fault consistent with records from towns and cities nearby, such as Sulmona and Pescasseroli, for which earthquake records extend back to Roman times. Indeed, the record is thought to be complete since AD 1349 for magnitudes larger than Mw5.8; see Guerrieri et al. (2002), Guidoboni et al. (2007), Pace et al. (2006), and Roberts and Michetti (2004) for discussions of the completeness period.

We performed 522 560 MCMC iterations of each chain and used the same number of chains, burn-in and thinning as in the synthetic case. A maximal PSRF value of 1.003 indicates convergence of the MCMC algorithm. The minimal wRMS among all scenarios was 1.78. The same settings are used for the offline calculations (∼15 min) and the total time was ∼3 days.

Figure 14 suggests a scarp age at approximately 19 ka consistent with expected slope stabilization ages, which may well be the marker for the time of the demise of the LGM shown in Fig. 1. After ∼19 ka, Figs. 14 and 15 indicate highly variable slip rates throughout the investigated time domain: initially a constant moderately high slip rate until ∼15 ka, a decreasing slip rate until ∼8 ka, a sudden peak in slip from ∼5 to 3 ka, and a low occurrence probability of earthquake events in the past ∼2500 years. The lack of earthquakes in the past ∼2500 years is consistent with historical earthquake records, which show no major earthquakes in the vicinity.

The displacement histories of the Fiamignano and Frattura faults provide insights into slip-rate variability. The relatively short-lived bursts of high activity observed in our results occur at different times on faults in the same tectonic setting. This suggests that this is not a regional pulse of synchronous high slip, but probably related to the dynamics of interaction between each fault and its neighbors (Cowie et al.2012), and that slip is highly clustered in time, which has important implications for forecasting seismic hazards as discussed below.

5.3 Application to regional probabilistic seismic hazard assessment

Our results show that calculating probabilistic seismic hazard is considerably more challenging than previously thought. Typical fault-based and time-dependent seismic hazard models are based on the BPT distribution and require specification of the most recent earthquake time, the mean inter-event time, and the coefficient of variation (CV) value (standard deviation of inter-event times divided by the mean inter-event time) to factor in variability in the inter-event time (Pace et al.2016). However, our results show that in addition to the variability in inter-event times around a constant slip rate, faults show heightened activity and quiescence over time periods lasting a few millennia relative to the longer-term deformation rate. The differences in slip rate between time of heightened activity (>1 cm yr−1) and quiescence (<0.1 cm yr−1) are dramatic. These two timescales of slip-rate variability are not considered by current methods for calculating probabilistic seismic hazard (Pace et al.2006, 2016; Tesson et al.2016). To address this omission, we extend our more complex earthquake recurrence process – a Brownian motion with time-varying drift and noise – into the future, which enables us to sample more realistic next earthquake event times. More specifically, for each of our posterior samples, we continue the simulation of the corresponding Brownian motion sample path (such as that shown in Fig. 2) until another earthquake is generated. For this purpose, we use the sample-dependent values of the parameters describing the behavior of the Brownian motion: the average time 1∕Λ between large-scale changes in slip activity, the small-scale recurrence variability τ, and the hyperparameters α and m that describe the distribution of average interarrival times in new slip activity periods. Thus, our approach is a natural extension of the state-of-the-art methodology that accounts for large-scale slip-rate variability and informs the resulting more complex model by the full displacement history at a given fault site.

Posteriors of the next earthquake time for the Fiamignano and Frattura faults are shown in Fig. 16. For comparison, we show results based on dmin=10 and dmin=50. As expected, the assumption that only earthquakes with a slip larger than 50 cm can occur at a given fault reduces the total number of earthquakes and thus the predicted hazard of an impending future earthquake. Furthermore, we include the results of a hybrid BED-BPT approach that uses the results of our MCMC algorithm but neglects different timescales of slip-rate variability in future event time simulations. (Since typical historical earthquake records are too sparse to provide realistic estimates of inter-event times and CV values, a pure BPT approach based on historical earthquake records would vastly underestimate seismic hazard.) For this purpose, we randomly generate a future event time from the BPT distribution, conditioned on the event occurring in the future, for each sample of our MCMC algorithm, using the mean inter-event time, CV value, and the most recent event time of that sample. Figure 16 shows that the effect of predicting future events with the BPT distribution is similar to that of using dmin=50. While the latter underestimates the number of past earthquakes, the former underestimates the chance of a fault becoming active despite recent inactivity.

An important finding is that the values differ between the two faults, and with measurement data from more faults it would be of great interest to map such posteriors across entire regions such as that shown in Fig. 8. We expect site-specific values for these posteriors to change on a length scale of 20–30 km, the length of individual faults. Thus, our approach facilitates high spatial resolution seismic hazard mapping by implicitly including individual seismic sources (active/capable faults), the long-term slip rates, and importantly the two timescales of slip-rate variability (millennial-scale heightened activity or quiescence, and inter-earthquake time variability).

However, a note of caution is that these results are probably only meaningful for near-future prediction. A physical basis that explains the cause of the slip-rate variability resulting from fault interaction would further improve probabilistic hazard analysis, as mentioned in Tesson et al. (2016). During and after rupture, stress is transferred onto neighboring faults, and this is thought to produce a temporal variation in slip rate on faults that manifests itself in terms of temporal earthquake clustering (Scholz2010; Cowie et al.2012; Mildon et al.2017). Interaction allows the fault systems to share out the work associated with deforming a region between different faults so that only a few of the total number of active faults take part in regional deformation on millennial timescales (Cowie et al.2012). Further work, including 36Cl results from many faults in the same region, is needed to elucidate such interaction.

6 Conclusions

This work provides a validation of the 36Cl approach to inference of earthquake occurrence from the demise of the LGM until present. We propose a Bayesian MCMC method for the study of earthquake displacement histories, with applications to regional probabilistic hazard assessment. The method improves on the 36Cl modeling in Schlagenhauf et al. (2010), the Bayesian inference for 36Cl earthquake recovery in Cowie et al. (2017), and the regional probabilistic hazard assessment in Pace et al. (2006) and Tesson et al. (2016). After demonstrating identifiability in the inverse problem through a synthetic case study, we present probabilistic earthquake displacement histories for the Fiamignano and Frattura faults in the Italian Apennines. We obtain highly variable slip rates at both faults, in agreement with earlier studies of the region. At the Fiamignano fault, our findings suggest slip in earthquakes at times when the Colosseum and other ancient buildings in Rome were damaged. Conversely, at the Frattura fault, our result is consistent with the fact that no large earthquakes were reported since Roman times.

Code availability

The MATLAB code Bayesian Earthquake Dating (BED) v1.0.1 of our Bayesian MCMC method for 36Cl earthquake dating and regional probabilistic seismic hazard assessment and the data required to reproduce our results are available as Supplement and at (Beck et al., 2018). Future releases of BED can be found at (last access: 22 August 2018; Beck et al., 2018).


The supplement related to this article is available online at:

Author contributions

JB and SW developed and implemented the MCMC algorithm. All authors contributed to the underlying statistical modeling and the 36Cl simulation itself. JB was the main author of the manuscript with SW authoring Sect. 3. GPR contributed with an in-depth description of the physical problem with the corresponding figures, co-wrote the introduction, collated the information on the historical accounts, and provided geological and seismological interpretation of the results.

Competing interests

The authors declare that they have no conflict of interest.


This work was supported by NERC directed grant NE/J017434/1 “Probability, Uncertainty and Risk in the Environment”, NERC standard grant NE/I024127/1 “Earthquake hazard from 36-Cl exposure dating of elapsed time and Coulomb stress transfer”, NERC standard grant NE/E016545/1, “Testing Theoretical Models for Earthquake Clustering using 36Cl Cosmogenic Exposure Dating of Active Normal Faults in Central Italy”, and the KAUST Office of Sponsored Research (OSR) under award no. URF/1/2584-01-01. We thank the many participants in these grants for discussions on earthquakes and 36Cl, although the views expressed in this paper are our own and any misconceptions are our sole responsibility. Francesco Iezzi kindly provided the field photograph in Fig. 1h. Eutizio Vittori kindly advised on interpretations of earthquake damage to buildings in Rome.

Edited by: Thomas Poulet
Reviewed by: Bruno Pace and one anonymous referee


Allen, J. R. M., Brandt, U., Brauer, A., Huntley, B., Keller, J., Kraml, M., Mackensen, A., Mingram, J., Negendank, J. F. W., Nowaczyk, N. R., Watts, W. A., Wulf, S., Zolitschka, B., Hubberten, H.-W., and Oberhänsli, H.: Rapid environmental changes in southern Europe during the last glacial period, Nature, 400, 740–743, 1999. a, b

Barthelmann, V., Novak, E., and Ritter, K.: High dimensional polynomial interpolation on sparse grids, Adv. Comput. Mathe., 12, 273–288, 2000. a

Beck, J., Wolfers, S., and Roberts, G. P.: Bayesian Earthquake Dating (BED) (Version 1.0.1), Zenodo,, 2018. 

Benedetti, L., Finkel, R., Papanastassiou, D., King, G., Armijo, R., Ryerson, F., Farber, D., and Flerit, F.: Post-glacial slip history of the Sparta fault (Greece) determined by 36Cl cosmogenic dating: evidence for non-periodic earthquakes, Geophys. Res. Lett., 29, 1–4,, 2002. a

Brynjarsdottir, J. and O'Hagan, A.: Learning about physical parameters: the importance of model discrepancy, Inverse Prob., 30, 114007,, 2014. a, b

Cavinato, G. and Celles, P. D.: Extensional basins in the tectonically bimodal central Apennines fold-thrust belt, Italy: response to corner flow above a subducting slab in retrograde motion, Geology, 27, 955–958, 1999. a

Cowie, P., Scholz, C., Roberts, G. P., Walker, J. F., and Steer, P.: Viscous roots of active seismogenic faults revealed by geologic slip rate variations, Nat. Geosci., 6, 1036,, 2013. a, b

Cowie, P. A., Phillips, R. J., Roberts, G. P., McCaffrey, K., Zijerveld, L. J. J., Gregory, L. C., Faure Walker, J., Wedmore, L. N. J., Dunai, T. J., Binnie, S. A., Freeman, S. P. H. T., Wilcken, K., Shanks, R. P., Huismans, R. S., Papanikolaou, I., Michetti, A. M., and Wilkinson, M.: Orogen-scale uplift in the central Italian Apennines drives episodic behaviour of earthquake faults, Tech. Rep. 44858, 2017. a, b, c, d, e, f, g, h, i, j

Cowie, P. A., Roberts, G. P., Bull, J. M., and Visini, F.: Relationships between fault geometry, slip rate variability and earthquake recurrence in extensional settings, Geophys. J. Int., 189, 143–160, 2012. a, b, c

Galli, P. A. and Molin, D.: Beyond the damage threshold: the historic earthquakes of Rome, B. Earthq. Eng., 12, 1277–1306, 2014. a, b, c, d

Gelman, A. and Rubin, D. B.: Inference from iterative simulation using multiple sequences, Stat. Sci., 7, 457–472, 1992. a, b

Gosse, J. C. and Phillips, F. M.: Terrestrial in situ cosmogenic nuclides: theory and application, Quaternary Sci. Rev., 20, 1475–1560, 2001. a, b

Green, P. J.: Reversible jump Markov chain Monte Carlo computation and Bayesian model determination, Biometrika, 82, 711–732, 1995. a

Guerrieri, L., Pascarella, F., Silvestri, S., and Serva, L.: Evoluzione recente dea paesaggio e dissesto geologico-idraulico: primo risultati in un'area campione dell'Appennino Centrale (valle del Salto–Rieti), Boll. Soc. Geol. Ital, 57, 453–461, 2002. a, b

Guidoboni, E., Ferrari, G., Mariotti, D., Comastri, A., Tarabusi, G., and Valensise, G.: Catalogue of Strong Earthquakes in Italy (461 BC-1997) and Mediterranean Area (760 BC-1500), available at: (last access: 22 October 2018), 2007. a, b, c

Heisinger, B., Lal, D., Jull, A., Kubik, P., Ivy-Ochs, S., Knie, K., and Nolte, E.: Production of selected cosmogenic radionuclides by muons: 2. Capture of negative muons, Earth Planet. Sci. Lett., 200, 357–369, 2002. a

Higham, D. J.: An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev., 43, 525–546, 2001. a

Lal, D.: In situ-produced cosmogenic isotopes in terrestrial rocks, Annu. Rev. Earth Planet. Sci., 16, 355–388, 1988. a

Marrero, S. M., Phillips, F. M., Caffee, M. W., and Gosse, J. C.: CRONUS-Earth cosmogenic 36 Cl calibration, Quaternary Geochronol., 31, 199–219, 2016. a

Matthews, M. V., Ellsworth, W. L., and Reasenberg, P. A.: A Brownian model for recurrent earthquakes, B. Seismol. Soc. Am., 92, 2233–2250, 2002. a, b, c

Mildon, Z. K., Roberts, G. P., Faure Walker, J. P., and Iezzi, F.: Coulomb stress transfer and fault interaction over millennia on non-planar active normal faults: the Mw 6.5–5.0 seismic sequence of 2016–2017, central Italy, Geophys. J. Int., 210, 1206–1218, 2017. a, b

Mitchell, S. G., Matmon, A., Bierman, P. R., Enzel, Y., Caffee, M., and Rizzo, D.: Displacement history of a limestone normal fault scarp, northern Israel, from cosmogenic 36Cl, J. Geophys. Res.-Solid Earth, 106, 4247–4264, 2001. a, b

Pace, B., Peruzza, L., Lavecchia, G., and Boncio, P.: Layered seismogenic source model and probabilistic seismic-hazard analyses in central Italy, B. Seismol. Soc. Am., 96, 107–132, 2006. a, b, c, d, e

Pace, B., Bocchini, G. M., and Boncio, P.: Do static stress changes of a moderate-magnitude earthquake significantly modify the regional seismic hazard? Hints from the L'Aquila 2009 normal-faulting earthquake (Mw 6.3, central Italy), Terra Nova, 26, 430–439, 2014. a

Pace, B., Visini, F., and Peruzza, L.: FiSH: MATLAB Tools to Turn Fault Data into Seismic-Hazard Models, Seismol. Res. Lett., 87, 374–386, 2016. a, b

Palumbo, L., Benedetti, L., Bourles, D., Cinque, A., and Finkel, R.: Slip history of the Magnola fault (Apennines, Central Italy) from 36 Cl surface exposure dating: evidence for strong earthquakes over the Holocene, Earth Planet. Sci. Lett., 225, 163–176, 2004. a

Peyron, O., Guiot, J., Cheddadi, R., Tarasov, P., Reille, M., de Beaulieu, J.-L., Bottema, S., and Andrieu, V.: Climatic reconstruction in Europe for 18,000 yr BP from pollen data, Quaternary Res., 49, 183–196, 1998. a

Phillips, F. M., Zreda, M. G., Flinsch, M. R., Elmore, D., and Sharma, P.: A reevaluation of cosmogenic 36Cl production rates in terrestrial rocks, Geophys. Res. Lett., 23, 949–952, 1996. a

Piccardi, L., Gaudemer, Y., Tapponnier, P., and Boccaletti, M.: Active oblique extension in the central Apennines (Italy): evidence from the Fucino region, Geophys. J. Int., 139, 499–530, 1999. a

Robert, C. P. and Casella, G.: Monte Carlo Statistical Methods, Springer Texts in Statistics, 2004. a, b

Roberts, G. P. and Michetti, A. M.: Spatial and temporal variations in growth rates along active normal fault systems: an example from the Lazio–Abruzzo Apennines, central Italy, J. Struct. Geol., 26, 339–376, 2004. a, b, c, d

Roberts, G. P., Michetti, A. M., Cowie, P., Morewood, N. C., and Papanikolaou, I.: Fault slip-rate variations during crustal-scale strain localisation, central Italy, Geophys. Res. Lett., 29, 1–4,, 2002.  a

Rougier, J.: Probabilistic inference for future climate using an ensemble of climate model evaluations, Clim. Change, 81, 247–264, 2007. a

Rougier, J., Sparks, S., and Hill, L. J.: Risk and Uncertainty Assessment for Natural Hazards, Cambridge University Press, 2013. a, b

Schlagenhauf, A., Gaudemer, Y., Benedetti, L., Manighetti, I., Palumbo, L., Schimmelpfennig, I., Finkel, R., and Pou, K.: Using in situ Chlorine-36 cosmonuclide to recover past earthquake histories on limestone normal fault scarps: a reappraisal of methodology and interpretations, Geophys. J. Int., 182, 36–72, 2010. a, b, c, d, e, f, g, h, i, j

Schlagenhauf, A., Manighetti, I., Benedetti, L., Gaudemer, Y., Finkel, R., Malavieille, J., and Pou, K.: Earthquake supercycles in central Italy, inferred from 36Cl exposure dating, Earth Planet. Sci. Lett., 307, 487–500, 2011. a

Scholz, C. H.: Large earthquake triggering, clustering, and the synchronization of faults, B. Seismol. Soc. Am., 100, 901–909, 2010. a

Stone, J., Allan, G., Fifield, L., and Cresswell, R.: Cosmogenic chlorine-36 from calcium spallation, Geochim. Cosmochim. Acta, 60, 679–692, 1996. a

Talts, S., Betancourt, M., Simpson, D., Vehtari, A., and Gelman, A.: Validating Bayesian Inference Algorithms with Simulation-Based Calibration, arXiv preprint arXiv:1804.06788, 2018. a

Tesson, J., Pace, B., Benedetti, L., Visini, F., Delli Rocioli, M., Arnold, M., Aumaître, G., Bourlès, D., and Keddadouche, K.: Seismic slip history of the Pizzalto fault (central Apennines, Italy) using in situ-produced 36Cl cosmic ray exposure dating and rare earth element concentrations, J. Geophys. Res.-Solid Earth, 121, 1983–2003, 2016. a, b, c

Vittori, E.: Archaeoseismic evidence in Roma, Italy (pre-congress field trip), in: 6th INQUA International Workshop on Active Techonics Paleaoseismology and Archaeoseismology, 19–24 April 2015, 1–25, 2015. a, b, c

Wells, D. L. and Coppersmith, K. J.: New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement, B. Seismol. Soc. Am., 84, 974–1002, 1994. a, b

Wilkinson, M. W., McCaffrey, K. J., Jones, R. R., Roberts, G. P., Holdsworth, R. E., Gregory, L. C., Walters, R. J., Wedmore, L., Goodall, H., and Iezzi, F.: Near-field fault slip of the 2016 Vettore M w 6.6 earthquake (Central Italy) measured using low-cost GNSS, Scientific Reports, 7, 4612, 2017. a

Woodard, D. B., Schmidler, S. C., and Huber, M.: Conditions for rapid mixing of parallel and simulated tempering on multimodal distributions, The Ann. Appl. Prob., 19, 617–640, 2009. a, b

Yildirim, C., Ersen Aksoy, M., Akif Sarikaya, M., Tuysuz, O., Genc, S., Ertekin Doksanalti, M., Sahin, S., Benedetti, L., Tesson, J., and Team, A.: Deriving earthquake history of the Knidos Fault Zone, SW Turkey, using cosmogenic 36Cl surface exposure dating of the fault scarp, in: EGU General Assembly Conference Abstracts, Vol. 18, p. 14251, 2016. a

Zreda, M. and Noller, J. S.: Ages of prehistoric earthquakes revealed by cosmogenic chlorine-36 in a bedrock fault scarp at Hebgen Lake, Science, 282, 1097–1099, 1998. a

Short summary
Seismic hazard assessment requires records of earthquake recurrence with many slip events. Current data from paleoseismology on individual faults are sparse and do not provide stable estimates of earthquake recurrence. We propose a statistical model-based method to study timings of earthquakes over the past few millennia. The results agree with historical earthquakes for faults in the Italian Apennines, and can aid future studies of fault interactions over multiple earthquake cycles.