Development and technical paper 01 Nov 2018
Development and technical paper  01 Nov 2018
Bayesian earthquake dating and seismic hazard assessment using chlorine36 measurements (BED v1)
 ^{1}Computer, Electrical and Mathematical Sciences & Engineering (CEMSE), King Abdullah University of Science and Technology (KAUST), Thuwal 239556900, Saudi Arabia
 ^{2}Department of Earth and Planetary Sciences, Birkbeck College, University of London, WC1E 7HX, UK
 ^{1}Computer, Electrical and Mathematical Sciences & Engineering (CEMSE), King Abdullah University of Science and Technology (KAUST), Thuwal 239556900, Saudi Arabia
 ^{2}Department of Earth and Planetary Sciences, Birkbeck College, University of London, WC1E 7HX, UK
Correspondence: Joakim Beck (joakim.beck@kaust.edu.sa)
Hide author detailsCorrespondence: Joakim Beck (joakim.beck@kaust.edu.sa)
Over the past 20 years, analyzing the abundance of the isotope chlorine36 (^{36}Cl) has emerged as a popular tool for geologic dating. In particular, it has been observed that ^{36}Cl 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 ^{36}Cl accumulation for a set of faultscarp rock samples, which are progressively exhumed during earthquakes, in order to infer displacement histories from ^{36}Cl measurements. While the physical models underlying such simulations have continuously been improved, the inverse problem of recovering displacement histories from ^{36}Cl 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 Markovchain Monte Carlo algorithm for the inverse problem, and provides a validation of the ^{36}Cl 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.
 Article
(2505 KB) 
Supplement
(279 KB)  BibTeX
 EndNote
A fundamental problem in earthquake science is the paucity of reliable earthquake records including multiple largemagnitude 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 chlorine36 (^{36}Cl) cosmogenic exposure dating of active normal faults (Zreda and Noller, 1998; 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 ^{36}Cl production in rock. The resulting characteristic ^{36}Cl concentration profiles along fault planes therefore provide information about the timing and intensity of past earthquakes.
A comprehensive mathematical model of ^{36}Cl production was provided in Gosse and Phillips (2001) and later formed the basis of a MATLAB code that computes ^{36}Cl 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 nonconvex 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 Markovchain Monte Carlo (MCMC) sampler (Robert and Casella, 2004) 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 ^{36}Cl 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 faultdependent 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 ^{36}Cl production can be completely accurate, we include an estimation of model discrepancy (Rougier, 2007; Rougier et al., 2013; Brynjarsdottir and O'Hagan, 2014), which results in more realistic confidence bands and may help guide various research efforts to improve the modeling of ^{36}Cl 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 sliprate variability in normal faults in the Italian Apennines. The timing of sliprate 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 ^{36}Cl earthquake dating, we extend our model to allow for Bayesian regional probabilistic seismic hazard assessment.
We supplement this paper with an easytouse MATLAB code of the proposed Bayesian MCMC method for earthquake dating and regional probabilistic seismic hazard assessment.
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 freezethaw 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 postglacial transition.
Before the demise of the LGM (see Fig. 1d), erosion via freezethaw 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 Michetti, 2004; Cowie et al., 2013, 2017). The type of slip is illustrated in Fig. 1h, which shows decimeterscale slips produced during two earthquakes at the Mt. Vettore fault on 24 August 2016 (M_{w} 6.1) and on 30 October 2016 (M_{w} 6.6). The M_{w} 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 ^{36}Cl concentration in a fixed sample of limestone bedrock evolves according to
where g(t) is the ^{36}Cl concentration (atoms g^{−1}), ϕ(t) is the production rate (atoms g^{−1} yr^{−1}), and λ is the decay rate of ^{36}Cl (yr^{−1}). The main ^{36}Cl production pathways in pure limestones are spallation of ^{40}Ca, muon capture by ^{40}Ca, and thermal neutron capture by ^{35}Cl (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 earthquakeinduced changes to the surrounding environment. Roughly speaking, samples taken from the trench have experienced ^{36}Cl 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 ^{36}Cl 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 postglacial alluvial fans and eroded gullies where exposure is also influenced by postglacial erosion and sedimentation (see Fig. 1f)
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 ^{36}Cl production (Lal, 1988; Phillips et al., 1996; Mitchell et al., 2001) have been assembled into a MATLAB code that computes ^{36}Cl 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 firstorder finitedifference 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 constantintime production rate. These changes improved the predicted ^{36}Cl concentration by around 5 % in our numerical experiments. Furthermore, we introduce an offline phase during which we precompute 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, T_{init}, before which erosion immediately eroded scarps formed by earthquakes, and after which erosion stopped completely.
The provided MATLAB code calculates ^{36}Cl concentrations in a set of bedrock samples for given sequences $\mathbf{d}=({d}_{\mathrm{1}},\mathrm{\dots},{d}_{N})$ and $\mathbf{s}=({s}_{\mathrm{1}},\mathrm{\dots},{s}_{N})$ 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 ^{36}Cl by fast secondary neutrons at the surface from ^{40}Ca, Ψ_{sp} (atoms g^{−1} yr^{−1});

slow muon capture production rate of ^{36}Cl at the surface from ^{40}Ca, Ψ_{μ} (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, T_{init} (yr^{−1}).

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 ^{36}Cl measurements $\mathbf{y}=({y}_{\mathrm{1}},\mathrm{\dots},{y}_{M})$ at M sample sites along a single fault scarp. More specifically, we compute posterior distributions for the vectors of fault displacements and event times, $\mathbf{d}=({d}_{\mathrm{1}},\mathrm{\dots},{d}_{N})$ and $\mathbf{s}=({s}_{\mathrm{1}},\mathrm{\dots},{s}_{N})$, respectively (and in particular for the number of events, N). In addition, we treat the values of $\mathbf{z}:=({T}_{\text{init}},{\mathrm{\Psi}}_{\text{sp}},{\mathrm{\Psi}}_{\mathit{\mu}},{\mathrm{\Lambda}}_{\text{sp}},{\mathrm{\Lambda}}_{\mathit{\mu}},{\mathit{\rho}}_{\text{coll}})$ as uncertain parameters and include them in the inference problem.
Denoting the vector of true ^{36}Cl concentrations by y^{*}, and the output of the computer model of Eq. (2) for given values of $\mathbf{x}:=(\mathbf{d},\mathbf{s},\mathbf{z})$ by g(x), we first observe the equation
where ϵ is the measurement error, i.e., the discrepancy between the true and the measured concentrations, and $\mathit{\delta}:={\mathbf{y}}^{*}\mathbf{g}\left(\mathbf{x}\right)$ is the model error, i.e., the discrepancy between the model output and the true concentrations (Rougier et al., 2013; Brynjarsdottir and O'Hagan, 2014). We assume that the measurement errors are independent and normally distributed, $\mathit{\u03f5}\sim \mathcal{N}(\mathrm{0},{\mathit{\sigma}}^{\mathrm{2}})$ for a given vector $\mathit{\sigma}=({\mathit{\sigma}}_{\mathrm{1}},\mathrm{\dots},{\mathit{\sigma}}_{M})$ of positive real numbers. The model error is assumed to be proportional to the measurements, $\mathit{\delta}\sim \mathcal{N}(\mathrm{0},(\mathit{\rho}\mathbf{y}{)}^{\mathrm{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 $\mathit{\theta}:=(\mathbf{x},\mathit{\rho},\mathit{\varphi})$, 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,
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 Casella, 2004), 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) $({\mathit{\theta}}_{k}{)}_{k=\mathrm{1}}^{\mathrm{\infty}}$ 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 $({\mathit{\theta}}_{k}^{\left(l\right)}{)}_{k=\mathrm{1}}^{\mathrm{\infty}}$, $\mathrm{1}\le l\le L$ with progressively flattened likelihood distributions,
and randomly swaps states between neighboring chains such that the resulting samples of the first chain $({\mathit{\theta}}_{k}^{\left(\mathrm{1}\right)}{)}_{k=\mathrm{1}}^{\mathrm{\infty}}$, 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.
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, $\mathit{\rho}\sim \mathcal{U}[\mathrm{0},{\mathit{\rho}}_{max}=\mathrm{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 $[{T}_{min},\mathrm{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 passagetime (BPT) distribution since it describes the time required for a Brownian motion to reach a certain threshold. To account for largescale time periods with differing average earthquake frequencies, we extend the model in Matthews et al. (2002) by including a Poisson (Λ T_{min})distributed number J of switch points $({t}_{j}{)}_{j=\mathrm{1}}^{J}$, where ${t}_{j}\sim \mathcal{U}\left(\right[{T}_{\text{min}},\mathrm{0}\left]\right)$ and $\mathrm{\Lambda}\sim \mathcal{U}\left(\right[{\mathrm{10}}^{\mathrm{4}},{\mathrm{10}}^{\mathrm{3}}\left]\right)$. 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 [10^{3},10^{4}]. We then define the drift a(t) and volatility b(t) of the Brownian motion in the intervals of the resulting partition of $[{T}_{\text{min}},\mathrm{0}]$ by
where the average interarrival times $\mathit{\nu}=({\mathit{\nu}}_{j}{)}_{j=\mathrm{0}}^{J}$ are independent and identically distributed according to an inverse gamma distribution with mean $m\sim \mathcal{U}\left(\right[\mathrm{200},\mathrm{2000}\left]\right)$ and shape parameter $\mathit{\alpha}\sim \mathcal{U}\left(\right[\mathrm{1},\mathrm{10}\left]\right)$, and where $\mathit{\tau}\sim \mathcal{U}\left(\right[\mathrm{0},\mathrm{1}\left]\right)$ controls the shortscale 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 (Higham, 2001) with time discretization Δt≈15 (yr^{−1}) to simulate the resulting stochastic process on $[{T}_{min},\mathrm{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
The values ${W}_{i}\sim \mathcal{N}(\mathrm{0},\mathrm{\Delta}t)$ together with ν and τ and the initial value ${X}_{\mathrm{0}}\sim \mathcal{U}\left(\right[\mathrm{0},\mathrm{1}\left]\right)$ form the vector of auxiliary variables ϕ referred to above, which fully determines the process $({X}_{i}{)}_{i=\mathrm{0}}^{{T}_{\text{min}}/\mathrm{\Delta}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 $[{d}_{min},{d}_{\text{max}}{]}^{N}$ conditioned on the requirement that ${\sum}_{n=\mathrm{1}}^{N}{d}_{n}={H}_{\text{sc}}$, where H_{sc} (see Fig. 1) is the present height of the fault scarp, and d_{min} and d_{max} are faultdependent bounds on displacement sizes.
Finally, we assign prior distributions to the components of z: ${T}_{\text{init}}\sim \mathcal{U}\left(\right[\mathrm{12}\phantom{\rule{0.125em}{0ex}}\mathrm{000},\mathrm{20}\phantom{\rule{0.125em}{0ex}}\mathrm{000}\left]\right)$, ${\mathrm{\Psi}}_{\text{sp}}\sim \mathcal{N}(\mathrm{48.8},\mathrm{1.7})$, ${\mathrm{\Psi}}_{\mathit{\mu}}\sim \mathcal{N}(\mathrm{190},\mathrm{19})$, ${\mathrm{\Lambda}}_{\text{sp}}\sim \mathcal{U}\left(\right[\mathrm{180},\mathrm{220}\left]\right)$, and ${\mathrm{\Lambda}}_{\mathit{\mu}}\sim \mathcal{U}\left(\right[\mathrm{1300},\mathrm{1700}\left]\right)$, 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 ${\mathit{\rho}}_{\text{coll}}\sim \mathcal{U}\left(\right[\mathrm{1.2},\mathrm{1.8}\left]\right)$ for all case studies considered in this work.
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 X_{0}, 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 $[{T}_{\text{min}},\mathrm{0}]$ and the corresponding piecewise constant drift and volatility coefficients, we employ reversible jump MCMC (Green, 1995), which allows for the application of MCMC to variables whose state space contains subspaces of different dimensionality. To propose new values of the variables $({W}_{i}{)}_{i=\mathrm{1}}^{{T}_{\text{min}}/\mathrm{\Delta}t}$, which drive the process $({X}_{i}{)}_{i=\mathrm{1}}^{{T}_{\text{min}}/\mathrm{\Delta}t}$, we again use a global proposal, which redraws all values independently, as well as a Brownian bridgetype local proposal that redraws the values within a random subinterval of $[{T}_{\text{min}},\mathrm{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 resampled too, and if a proposal changes the number of switch points, the vectors of drift and volatility coefficients are resampled as well.
In this section we apply our Bayesian MCMC method to synthetic ^{36}Cl data. To generate these data, we drew d and s from the prior distributions described in Sect. 3.1 with ${d}_{min}=\mathrm{10}$ (cm) and d_{max}=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 $\mathbf{z}=({T}_{\text{init}}=\mathrm{19000}$, ${\mathrm{\Psi}}_{\text{sp}}=\mathrm{49.5},{\mathrm{\Psi}}_{\mathit{\mu}}=\mathrm{200}$, Λ_{sp}=195, ${\mathrm{\Lambda}}_{\mathit{\mu}}=\mathrm{1700},{\mathit{\rho}}_{\text{coll}}=\mathrm{1.6})$. The remaining parameters, which are not part of the inference problem, were chosen as H_{sc}=2705, H_{tr}=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 ^{36}Cl 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., ${T}_{\text{init}}\sim \mathcal{U}\left(\right[\mathrm{12}\phantom{\rule{0.125em}{0ex}}\mathrm{000},\mathrm{20}\phantom{\rule{0.125em}{0ex}}\mathrm{000}\left]\right)$, ${\mathrm{\Psi}}_{\text{sp}}\sim \mathcal{N}(\mathrm{48.8},\mathrm{1.7})$, ${\mathrm{\Psi}}_{\mathit{\mu}}\sim \mathcal{N}(\mathrm{190},\mathrm{19})$, ${\mathrm{\Lambda}}_{\text{sp}}\sim \mathcal{U}\left(\right[\mathrm{180},\mathrm{220}\left]\right)$, and ${\mathrm{\Lambda}}_{\mathit{\mu}}\sim \mathcal{U}\left(\right[\mathrm{1300},\mathrm{1700}\left]\right)$, ${\mathit{\rho}}_{\text{coll}}\sim \mathcal{U}\left(\right[\mathrm{1.2},\mathrm{1.8}\left]\right)$.
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 % burnin period. Among these samples were 843 735 distinct scenarios, whose repetitions correspond to rejected proposals and reflect their statistical weight. Finally, to accelerate postprocessing and to decrease memory consumption, we only saved each fifth scenario together with its number of repetitions (thinning).
Figure 3 shows that the medians of the ^{36}Cl 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 ^{36}Cl concentrations at the end of the LGM, T_{init}, as opposed to the alternative to start with zero concentrations or imposing ad hoc preexposure 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 ^{36}Cl samples involve timeconsuming 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 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, burnin, and thinning as before, which led to a maximal PSRF value of 1.02.
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 Michetti, 2004; 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 Celles, 1999; 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 longterm 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.
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 M_{w}∼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.
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 (Vittori, 2015). 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 Molin, 2014). 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 Molin, 2014). Minor damage to the Colosseum also occurred in an AD 1703 earthquake (Vittori, 2015).
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 ^{36}Cl analysis. ^{36}Cl 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 interevent times between sliprate change points and constant displacement sizes.
In our analysis, we use the same site parameter values and priors as in Sect. 4 except that we now use ${d}_{max}=\mathrm{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 burnin, 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 E52680 v2 at 2.80 GHz (20 cores) with MATLAB (2016a).
The sampled posterior ^{36}Cl 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 sliprate 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 highintensity 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
^{36}Cl 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 ^{36}Cl data were collected sparsely (15 samples), similar to the situation in our synthetic case with a reduced amount of data shown. The sitespecific parameters for the Frattura fault are α=25, β=53, and γ=28; fault scarp height H_{sc}=1570, trench depth H_{tr}=130; and the rock mean density ρ_{rock}=2.7. Here the maximum displacement size is again chosen conservatively as d_{max}=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 ${\mathit{\rho}}_{\text{coll}}\sim \mathcal{U}\left(\right[\mathrm{1.2},\mathrm{1.8}\left]\right)$.
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 M_{w}5.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, burnin 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 sliprate variability. The relatively shortlived 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 faultbased and timedependent seismic hazard models are based on the BPT distribution and require specification of the most recent earthquake time, the mean interevent time, and the coefficient of variation (CV) value (standard deviation of interevent times divided by the mean interevent time) to factor in variability in the interevent time (Pace et al., 2016). However, our results show that in addition to the variability in interevent times around a constant slip rate, faults show heightened activity and quiescence over time periods lasting a few millennia relative to the longerterm 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 sliprate 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 timevarying 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 sampledependent values of the parameters describing the behavior of the Brownian motion: the average time 1∕Λ between largescale changes in slip activity, the smallscale 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 stateoftheart methodology that accounts for largescale sliprate 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 d_{min}=10 and d_{min}=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 BEDBPT approach that uses the results of our MCMC algorithm but neglects different timescales of sliprate variability in future event time simulations. (Since typical historical earthquake records are too sparse to provide realistic estimates of interevent 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 interevent 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 d_{min}=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 sitespecific 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 longterm slip rates, and importantly the two timescales of sliprate variability (millennialscale heightened activity or quiescence, and interearthquake time variability).
However, a note of caution is that these results are probably only meaningful for nearfuture prediction. A physical basis that explains the cause of the sliprate 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 (Scholz, 2010; 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 ^{36}Cl results from many faults in the same region, is needed to elucidate such interaction.
This work provides a validation of the ^{36}Cl 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 ^{36}Cl modeling in Schlagenhauf et al. (2010), the Bayesian inference for ^{36}Cl 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.
The MATLAB code Bayesian Earthquake Dating (BED) v1.0.1 of our Bayesian MCMC method for ^{36}Cl earthquake dating and regional probabilistic seismic hazard assessment and the data required to reproduce our results are available as Supplement and at https://doi.org/10.5281/zenodo.1402093 (Beck et al., 2018). Future releases of BED can be found at https://github.com/beckjh/bed36Cl (last access: 22 August 2018; Beck et al., 2018).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1143832018supplement.
JB and SW developed and implemented the MCMC algorithm. All authors contributed to the underlying statistical modeling and the ^{36}Cl simulation itself. JB was the main author of the manuscript with SW authoring Sect. 3. GPR contributed with an indepth description of the physical problem with the corresponding figures, cowrote the introduction, collated the information on the historical accounts, and provided geological and seismological interpretation of the results.
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 36Cl 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/25840101. We thank the many
participants in these grants for discussions on earthquakes and
^{36}Cl, 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, https://doi.org/10.5281/zenodo.1402093, 2018.
Benedetti, L., Finkel, R., Papanastassiou, D., King, G., Armijo, R., Ryerson, F., Farber, D., and Flerit, F.: Postglacial slip history of the Sparta fault (Greece) determined by 36Cl cosmogenic dating: evidence for nonperiodic earthquakes, Geophys. Res. Lett., 29, 1–4, https://doi.org/10.1029/2001GL014510, 2002. a
Brynjarsdottir, J. and O'Hagan, A.: Learning about physical parameters: the importance of model discrepancy, Inverse Prob., 30, 114007, https://doi.org/10.1088/02665611/30/11/114007, 2014. a, b
Cavinato, G. and Celles, P. D.: Extensional basins in the tectonically bimodal central Apennines foldthrust 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, https://doi.org/10.1038/NGEO1991, 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.: Orogenscale 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 geologicoidraulico: 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 BC1997) and Mediterranean Area (760 BC1500), available at: http://storing.ingv.it/cfti4med/ (last access: 22 October 2018), 2007. a, b, c
Heisinger, B., Lal, D., Jull, A., Kubik, P., IvyOchs, 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 situproduced 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.: CRONUSEarth 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 nonplanar 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 seismichazard 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 moderatemagnitude earthquake significantly modify the regional seismic hazard? Hints from the L'Aquila 2009 normalfaulting 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 SeismicHazard 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 sliprate variations during crustalscale strain localisation, central Italy, Geophys. Res. Lett., 29, 1–4, https://doi.org/10.1029/2001GL013529, 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 Chlorine36 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 chlorine36 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 SimulationBased 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 situproduced 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 (precongress 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.: Nearfield fault slip of the 2016 Vettore M w 6.6 earthquake (Central Italy) measured using lowcost 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 chlorine36 in a bedrock fault scarp at Hebgen Lake, Science, 282, 1097–1099, 1998. a
 Abstract
 Introduction
 Simulation of ^{36}Cl concentrations along fault scarps
 Bayesian inference of displacement histories using MCMC
 Application to synthetic ^{36}Cl data
 Application to normal faults in the Italian Apennines
 Conclusions
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 References
 Supplement
 Abstract
 Introduction
 Simulation of ^{36}Cl concentrations along fault scarps
 Bayesian inference of displacement histories using MCMC
 Application to synthetic ^{36}Cl data
 Application to normal faults in the Italian Apennines
 Conclusions
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 References
 Supplement