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

Over the past 20 years, analyzing the abundance of the isotope chlorine-36 (³⁶Cl) has emerged as a popular tool for geologic dating. In particular, it has been observed that ³⁶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 ³⁶Cl accumulation for a set of fault-scarp rock samples, which are progressively exhumed during earthquakes, in order to infer displacement histories from ³⁶Cl measurements. While the physical models underlying such simulations have continuously been improved, the inverse problem of recovering displacement histories from ³⁶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 Markov-chain Monte Carlo algorithm for the inverse problem, and provides a validation of the ³⁶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.

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 (³⁶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 ³⁶Cl production in rock. The resulting characteristic ³⁶Cl concentration profiles along fault planes therefore provide information about the timing and intensity of past earthquakes.

A comprehensive mathematical model of ³⁶Cl production was provided in Gosse and Phillips (2001) and later formed the basis of a MATLAB code that computes ³⁶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 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 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 ³⁶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 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 ³⁶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 ³⁶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 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 ³⁶Cl 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 ³⁶Cl 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 Michetti, 2004; 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 (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 ³⁶Cl concentration in a fixed sample of limestone bedrock evolves according to

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}g\left(t\right)}{\mathrm{d}t}}=\mathit{\varphi }\left(t\right)-\mathit{\lambda }g\left(t\right),\end{array}$$

where g(t) is the ³⁶Cl concentration (atoms g⁻¹), ϕ(t) is the production rate (atoms g⁻¹ yr⁻¹), and λ is the decay rate of ³⁶Cl (yr⁻¹). The main ³⁶Cl production pathways in pure limestones are spallation of ⁴⁰Ca, muon capture by ⁴⁰Ca, and thermal neutron capture by ³⁵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 earthquake-induced changes to the surrounding environment. Roughly speaking, samples taken from the trench have experienced ³⁶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 ³⁶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 post-glacial alluvial fans and eroded gullies where exposure is also influenced by post-glacial erosion and sedimentation (see Fig. 1f)

Figure 1 Images 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 ³⁶Cl production (Lal, 1988; Phillips et al., 1996; Mitchell et al., 2001) have been assembled into a MATLAB code that computes ³⁶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 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 ³⁶Cl 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, T_init, before which erosion immediately eroded scarps formed by earthquakes, and after which erosion stopped completely.

The provided MATLAB code calculates ³⁶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⁻³);

  • colluvial wedge mean density, ρ_coll (g cm⁻³);

  • spallation production rate of ³⁶Cl by fast secondary neutrons at the surface from ⁴⁰Ca, Ψ_sp (atoms g⁻¹ yr⁻¹);

  • slow muon capture production rate of ³⁶Cl at the surface from ⁴⁰Ca, Ψ_μ (atoms g⁻¹ yr⁻¹);

  • neutron apparent attenuation length, Λ_sp, and muon apparent attenuation length, Λ_μ, for a horizontal unshielded surface (g cm⁻²);

  • 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⁻¹).

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 ³⁶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 ³⁶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

$$\begin{array}{}\text{(2)}& \mathbf{y}={\mathbf{y}}^{*}+\mathit{ϵ}=\mathbf{g}\left(\mathbf{x}\right)+\mathit{\delta }+\mathit{ϵ},\end{array}$$

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{ϵ}\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,

$$\begin{array}{}\text{(3)}& P_{\mathit{\theta }|\mathbf{y}}\propto P_{\mathbf{y}|\mathit{\theta }}{P}_{\mathit{\theta }},\end{array}$$

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

$$P_{\mathbf{y}|\mathit{\theta }}\sim \mathcal{N}\left(\mathbf{g}\right(\mathbf{x}),{\mathit{\sigma }}^{\mathrm{2}}+(\mathit{\rho }\mathbf{y}{)}^{\mathrm{2}}).$$

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,

$$P_{\mathbf{y}|\mathit{\theta }}^{\text{(l)}}\sim \mathcal{N}\left(\mathbf{g}\right(\mathit{\theta }),{\mathit{\kappa }}_{\text{l}}({\mathit{\sigma }}^{\mathrm{2}}+\left(\mathit{\rho }\mathbf{y}{)}^{\mathrm{2}}\right)),\mathrm{1}={\mathit{\kappa }}_{\mathrm{1}}<\mathrm{\dots }<{\mathit{\kappa }}_{\text{l}}<\mathrm{\dots }<{\mathit{\kappa }}_{\text{L}},$$

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, 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 2 Sample path of stochastic process used for earthquake event time generation. Sample path generated with X₀=0.16, $T_{\text{min}}=-\mathrm{30}$ kyr, $t_{\mathrm{1}}=-\mathrm{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 t₁ is represented by a red vertical line.

Figure 3 Measured and sampled ³⁶Cl concentrations for the synthetic case with 146 ³⁶Cl 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, $\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 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 (Λ 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³,10⁴]. 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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left(a\right(t),b(t{)}^{\mathrm{2}}):=(\mathrm{1}/{\mathit{\nu }}_j,{\mathit{\tau }}^{\mathrm{2}}/{\mathit{\nu }}_j)\text{for}{t}_j\le t\le t_{j+\mathrm{1}}\\ \text{(4)}& {\displaystyle }& {\displaystyle }(\text{with}{t}_{\mathrm{0}}:=-T_{min}\text{and}{t}_{J+\mathrm{1}}:=\mathrm{0}),\end{array}$$

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 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 (Higham, 2001) with time discretization Δt≈15 (yr⁻¹) 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

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\widehat{X}}_{i+\mathrm{1}}& :=X_i+a\left(i\mathrm{\Delta }t\right)\mathrm{\Delta }t+b\left(j\mathrm{\Delta }t\right)W_{i+\mathrm{1}},\\ X_{i+\mathrm{1}}& :=\left\{\begin{array}{l}{\widehat{X}}_{i+\mathrm{1}}\text{if }{\widehat{X}}_{i+\mathrm{1}}<\mathrm{1},\\ \mathrm{0}\text{otherwise.}\end{array}\right.\end{array}\end{array}$$

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 fault-dependent bounds on displacement sizes.

Figure 4 The synthetic case with 146 ³⁶Cl 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 5 Gelman–Rubin diagnostic for the synthetic case with 146 ³⁶Cl measurements

Finally, we assign prior distributions to the components of z: $T_{\text{init}}\sim \mathcal{U}\left(\right[\mathrm{12}\mathrm{000},\mathrm{20}\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.

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

Figure 7 The synthetic case with 16 ³⁶Cl 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 X₀, 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 bridge-type 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 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 8 Location 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 ³⁶Cl data

In this section we apply our Bayesian MCMC method to synthetic ³⁶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 ³⁶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}\mathrm{000},\mathrm{20}\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 % 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 9 Damage to the Colosseum in Rome.

Figure 3 shows that the medians of the ³⁶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 ³⁶Cl concentrations at the end of the LGM, T_init, 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 ³⁶Cl 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 10 Fiamignano fault: measured and sampled ³⁶Cl concentrations.

Figure 11 Fiamignano 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 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 long-term rates of vertical motion across the faults are in the order of 0.1 cm yr⁻¹.

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 12 Fiamignano 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 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.

Figure 13 Frattura fault: measured and sampled ³⁶Cl concentrations.

Figure 14 Frattura 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 15 Frattura 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 (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 ³⁶Cl analysis. ³⁶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 inter-event times between slip-rate change points and constant displacement sizes.

Figure 16 Posteriors of next earthquake time for Fiamignano (a) and Frattura (b) according to BED with d_min=10 and d_min=50, and according to a hybrid BED-BPT approach (with d_min=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 $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 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 ³⁶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 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⁻¹, we find rapid slip of 1–1.4 cm yr⁻¹ 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

³⁶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 ³⁶Cl 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 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_w5.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⁻¹) and quiescence (<0.1 cm yr⁻¹) 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 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 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 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 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 (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 ³⁶Cl results from many faults in the same region, is needed to elucidate such interaction.

6 Conclusions

This work provides a validation of the ³⁶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 ³⁶Cl modeling in Schlagenhauf et al. (2010), the Bayesian inference for ³⁶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.