# Bayesian Parameter Estimation for Space and Time Interacting Earthquake Rupture Model Using Historical and Physics-Based Simulated Earthquake Catalogs

Luis Ceferino, Percy Galvez, Jean Paul Ampuero, Anne Kiremidjian, Gregory Deierlein, Juan Villegas-Lanza

<figure><img src="image_1.png" /></figure>

## To cite this version:

Luis Ceferino, Percy Galvez, Jean Paul Ampuero, Anne Kiremidjian, Gregory Deierlein, et al. Bayesian Parameter Estimation for Space and Time Interacting Earthquake Rupture Model Using Historical and Physics-Based Simulated Earthquake Catalogs. Bulletin of the Seismological Society of America, 2021, 111 (6), pp.3356-3373. 10.1785/0120210013. hal-03462495

HAL Id: hal-03462495

https://hal.science/hal-03462495v1

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Bayesian parameter estimation for space and time interacting
earthquake rupture model using historical and physics-based
simulated earthquake catalogs

Luis Ceferino¹, Percy Galvez²,³, Jean-Paul Ampuero⁴,⁵, Anne Kiremidjian¹, Gregory Deierlein¹, and Juan C. Villegas-Lanza⁶

¹Civil and Environmental Engineering Department, Stanford University

²King Abdullah University of Science and Technology

³ETH Zurich

⁴Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Goazur

⁵Seismological Laboratory, California Institute of Technology

⁶Instituto Geofísico del Perú

¹450 Serra Mall, Stanford, CA 94305

²Thuwal 23955, Saudi Arabia

³Rmistrasse 101, 8092 Zurich, Switzerland

⁴250 Rue Albert Einstein, 06560 Valbonne, France

⁵1200 E California Blvd, Pasadena, CA 91125, USA

⁶Calle Badajoz N 169 Urb. Mayorazgo IV Etapa, Lima, Perú

May 10, 2019

Corresponding author: Luis Ceferino. E-mail: lceferinor@gmail.com
Preprint submitted to the Bulletin of the Seismological Society of America
Abstract

This paper presents a robust parameter estimation technique for a probabilistic earthquake haz-
ard model that captures time and space interactions between earthquake mainshocks. The approach
addresses the existing limitations of parameter estimation techniques by developing a Bayesian for-
mulation and leveraging physics-based simulated synthetic catalogs to expand the limited datasets
of historical catalogs. The technique is based on a two-step Bayesian update that uses the syn-
thetic catalog to perform a first parameter estimation and then uses the historical catalog to further
calibrate the parameters. We applied this technique to analyze the occurrence of large-magnitude
interface earthquakes along 650 km of the central subduction zone in Peru, located offshore of Lima.
We built 2,000-years-long synthetic catalogs using quasi-dynamic earthquake cycle simulations based
on the rate-and-state friction law. The validity of the synthetic catalogs was verified by comparing
their annual magnitude exceedence rates to those of recorded seismicity and their predicted areas of
high interseismic coupling to those inferred from geodetic data. We show that when the Bayesian
update uses the combination of synthetic and historical data, instead of only the historical data, it
reduces the uncertainty of model parameter estimates by 45% on average. Further, our results show
that the time-dependent seismic hazard estimated with the both datasets is 40% smaller than the
one estimated with only the historical data.
# INTRODUCTION

Recent earthquake observations and improvements in the understanding of the mechanics of earthquake rupture processes have enabled the development of advanced probabilistic models of time-dependent earthquake hazards (Akinci et al., 2009; Field, 2015; Ceferino et al., 2018b). Even though these probabilistic models can capture complex spatio-temporal interactions between earthquakes, they still face two main challenges. First, there is a lack of robust statistical techniques to estimate the parameters of these complex models. Simplified parameter estimation techniques that are often used can lead to unreliable parameter estimates (Field, 2015; Ceferino et al., 2018b). Second, datasets, i.e., earthquake catalogs, are sparse and short, especially for earthquakes of large magnitude, which introduces large uncertainties in the parameter estimates. Thus, time-dependent seismic hazard analyses based on such complex earthquake models have low reliability and large uncertainties, particularly in fault systems with large earthquakes and long seismic gaps, i.e., long times without large ruptures (Hong and Goda, 2006). Here, we present a methodology that addresses these two limitations by utilizing robust statistical techniques and synthetic catalogs to expand sparse historical catalogs.

Robust statistical techniques for parameter estimation are mostly targeted to probabilistic models that only capture the time dependency of earthquake occurrence. For example, parameter estimation techniques based on maximum likelihood estimation (MLE) can be found for the Weibull, Lognormal, Gamma, and the Brownian passage-time (BPT) distributions (Hagiwara, 1974; Nishenko and Buland, 1987; Udias and Rice, 1975; Matthews et al., 2002). However, as the complexity of the probabilistic model increases, parameter estimation techniques become harder to formulate. For example, the Uniform California Earthquake Rupture Forecast version 3 (UCERF v3) time-dependent model utilizes a parameter estimation technique that is not consistent with the complexity of the model. UCERF v3 uses rupture data to find averages of mean rupture interevent times at multiple fault sections separately. However, UCERF v3 introduces spatial dependencies in earthquakes with overlapping rupture lengths, which requires a parameter estimation technique that evaluates interevent time data simultaneously instead of separately. As a result, UCERF v3 rupture simulations mismatch the assumed distributions of rupture interevent times on the fault (Field, 2015). Similarly, the model proposed by Ceferino et al. (2018b) presented a simple technique that estimates the means of rupture interevent times at different fault sections separately. Because the model incorporates spatial dependencies, such a technique can lead to unreliable parameter estimates (Ceferino et al., 2018a).

In addition to the lack of robust parameter estimation techniques, sparse earthquake data can significantly increase the uncertainties in earthquake occurrence modeling (Woessner and Wiemer, 2005). Most complete earthquake catalogs cover periods from a few decades to several hundreds of years depending on the region and earthquake magnitude of analysis. For example, catalogs including data with large earthquakes (i.e., Mw > 8.0) often cover only few hundred years. Because very large earthquakes have return periods larger than the catalogs span, these catalogs are often insufficient to characterize the recurrence time of large earthquakes. Data is inadequate even in regions with rich seismic catalogs. For example, the Tokai segment of the Nankai-Tokai subduction zone area in Japan, one of the Earth's best-identified set of fault segments, has a long catalog that covers 3,000 years because it was built using a combination of historical earthquakes, paleo-earthquakes and tsunami records. Yet the catalog contains only six large ruptures with magnitudes larger than 8.0, and such a small number of data points introduces large uncertainties in the interevent time distributions and in the earthquake recurrence modeling (Parsons et al., 2012).

To address these limitations, we propose a robust parameter estimation technique for the model proposed by Ceferino et al. (2018b) that leverages physics-based simulated synthetic catalogs to enlarge the historical earthquake occurrence dataset. The parameter estimation technique is based on a two-step Bayesian estimation that utilizes synthetic catalogs to perform an initial update of the parameters. It then uses these intermediary results as a prior distribution and performs a second update using the historical catalog to obtain the final distribution of the model parameters. The Bayesian formulation for each update is presented in this paper and is solved using a Markov Chain Monte Carlo (MCMC) algorithm. The physics-based synthetic catalog was developed using quasi-dynamic earthquake cycle modeling based on the rate-and-state friction law (Rubin and Ampuero, 2005; Ampuero and Rubin, 2008).

The remainder of the paper is structured as follows. First, a summary of the model proposed by Ceferino et al. (2018b) is provided for completeness. Then, the two-step Bayesian update is presented. Next, earthquake datasets for the subduction zone along the coast of Lima, Peru, are described, including the historical catalog and the physics-based simulated synthetic catalogs for large earthquakes. Finally,
the results of applying the parameter estimation technique are shown and discussed.

# PROBABILISTIC MODEL

Ceferino et al. (2018b, 2017) presented a probabilistic formulation for modeling time and space interactions between earthquake mainshocks. The model uses Brownian Passage Time (BPT) distributions to assess rupture interevent times in multiple fault segments to represent time interactions, and it uses spatial correlograms to represent spatial interactions of rupture occurrences. The model was proposed as an alternative to existing probabilistic models because (1) it captures time and space interactions of mainshocks, (2) it preserves the marginal distributions of interevent times after including the spatial interactions, and (3) it has an implicit physical interpretation compatible with earthquake rupture observations (Ceferino et al., 2018b). The following subsections briefly describe the earthquake rupture model for completeness.

## Fault and rupture representation

The model represents the contact surface between tectonic plates as a plane subdivided into smaller area sections. Figure 1 shows the fault geometry, the sections, and an earthquake rupture (in shaded areas) occurring at time t. N is the total number of sections in the fault, and t is a discrete time index with time steps of one year. A rupture at year t is represented by the rupture vector $X_t \in \{0, 1\}^N$. Each section has a corresponding element $X_t(j)$ of the vector $X_t$, representing its rupture state: $X_t(j) = 1$ if section j was involved in the rupture, or 0 otherwise. A vector $T_t \in \mathbb{N}^N$ contains the time since the last earthquake at year t for each section. Therefore, $T_{t+1}(j) = T_t(j) + 1$ if there is no rupture in the j-th section at time t (i.e., $X_t(j) = 0$), or $T_{t+1}(j) = 1$ otherwise. This is represented in the following equation:

$$T_{t+1}(j) = (T_t(j))(1 - X_t(j)) + 1 \quad (1)$$

<figure><img src="image_3.png" /><figcaption>Figure 1: Fault subdivision into small sections. The shaded sections represent a rupture at year t. Extracted from Ceferino et al. (2018b)</figcaption></figure>

## Probabilistic model formulation

The rupture occurrence at year t conditioned on the time since the last earthquake, i.e., $X_t|T_t$, is modeled as a multivariate Bernoulli distribution:

$$X_t|T_t \sim \text{Multivariate Bernoulli}(p_t, \Sigma) \quad (2)$$

where $p_t$ is the vector containing the rupture occurrence probabilities of the N sections, and $\Sigma$ is the covariance matrix of size $N \times N$ containing rupture correlations between the sections. The vector's element $p_t(j)$ is a function of the time since the last rupture $T_t(j)$ at the j-th section. $p_t(j)$ can be estimated as

$$p_t(j) = P[X_t(j) = 1|T_t(j)] = P[T_t(j) + 1 \ge \tau_j | \tau_j > T_t(j)] \quad (3)$$

where $\tau_j$ is a random variable that represents the rupture interevent time of the j-th section. $\tau_j$ is modeled as a BPT distribution with parameters $\mu_j$ and $\alpha_j$, and therefore $p_t(j)$ can be estimated as
$$p_t(j) = \frac{(\Phi[u_1(T_t(j))] - \Phi[u_1(T_t(j) - 1)]) + e^{2/\alpha_j^2}(\Phi[-u_2(T_t(j))] - \Phi[-u_2(T_t(j) - 1)]) - (\Phi[u_1(T_t(j) - 1)] + e^{2/\alpha_j^2}\Phi[-u_2(T_t(j) - 1)]))}{1 - (\Phi[u_1(T_t(j) - 1)] + e^{2/\alpha_j^2}\Phi[-u_2(T_t(j) - 1)]))} \quad (4a)$$

$$u_1(t) = \alpha_j^{-1}[t^{1/2}\mu_j^{-1/2} - t^{-1/2}\mu_j^{1/2}] \quad (4b)$$

$$u_2(t) = \alpha_j^{-1}[t^{1/2}\mu_j^{-1/2} + t^{-1/2}\mu_j^{1/2}] \quad (4c)$$

The elements $\rho_{i,j}$ of the covariance matrix $\Sigma$ are defined by a spherical correlogram in order to capture earthquake nucleation and spatial rupture propagation effects. The correlogram is

$$\rho_{i,j} = \exp\left(-\left(\frac{\mathrm{dist}(i,j)}{\gamma}\right)^2\right) \quad (5)$$

The correlation decreases as function of the distance $\mathrm{dist}(i,j)$ between the sections $i$ and $j$, with a characteristic correlation length $\gamma$.

## Copulas method

Because correlated multivariate Bernoulli distributions cannot be written in closed-form, an approximation is used through the copula method (Jin et al., 2015). A vector $Z_t$ of normally distributed random variables is first defined, with a zero-valued mean vector and a covariance equal to the covariance matrix $\Sigma$ (i.e., from the correlogram in Equation 5). Then, $X_t(j)$ is obtained by evaluating whether $\Phi[Z_t(j)]$ is smaller than $p_t(j)$, as

$$X_t(j) = \mathbf{1}\{\Phi(Z_t(j)) < p_t(j)\} \quad (6)$$

where $\mathbf{1}\{\xi\}$ is an indicator function that equals 1 if the Boolean argument $\xi$ is true, or 0 otherwise.

## BAYESIAN PARAMETER ESTIMATION PROCEDURE

Ceferino et al. (2018b) described a simple parameter estimation method for the probabilistic model summarized thus far. The method uses Maximum Likelihood Estimation (MLE) at each individual fault section separately to estimate $\mu_j$ and $\alpha_j$. Then, the method calibrates $\gamma$ to match annual exceedance rates and average seismic moment releases in tectonic faults. This method was applied to the subduction zone along the coast of Lima, Peru, due to its implementation simplicity. However, Ceferino et al. (2018a) showed that this method can generate unreliable estimates of rupture occurrence, particularly, when earthquake rupture data is sparse.

In this paper, we propose a more robust parameter estimation technique. It uses a two-step Bayesian parameter update and has two advantages in comparison to the existing parameter estimation technique proposed by Ceferino et al. (2018b). First, in addition to the historical earthquake catalogs, the technique uses synthetic physics-based simulated earthquake catalogs to enlarge the earthquake dataset. Second, our technique uses a Bayesian approach that estimates all the parameters simultaneously rather than separately. As shown by Ceferino et al. (2018a), this feature gives more reliable parameter estimates, particularly, for the fault sections with few or no known earthquake ruptures.

The two-step parameter estimation technique first uses a physics-based simulated earthquake catalog to perform an initial Bayesian update of the parameters, and then uses the historical earthquake catalog to perform a second Bayesian update as shown in Figure 2. In the first Bayesian update, the prior distributions of the parameters are defined as independent, log-normally distributed random variables. The log-normal distribution parameters are chosen by a combination of expert opinion, existing common BPT parameter values for similar tectonic regions (Sykes and Menke, 2006; Field, 2015), or average historical earthquake interevent times in the region. Then, the technique updates the parameters of the model to a posterior distribution using the physics-based catalog. The posterior is estimated using a Monte Carlo Markov Chain (MCMC) technique as described later. In the second Bayesian update, the prior distributions of the parameters are also defined as independent, log-normally distributed random variables. The log-normal distribution parameters are estimated using MLE on the realizations sampled from the MCMC of the first Bayesian update. Finally, the final posterior of the parameters is estimated using the MCMC and the historical earthquake catalog. The following subsection describes the Bayesian update formulation and the MCMC technique.
<figure><img src="image_4.png" /><figcaption>Figure 2: Two-step Bayesian update using both a synthetic physic-based simulated earthquake catalog and a historical catalog.</figcaption></figure>

## Bayesian update

The posterior distribution of the parameter vector is estimated with the Bayesian update. There are 2 × N + 1 parameters in the model, two per each of the N sections of the fault (α_j and μ_j), and another defining the correlogram model (γ).

### Posterior

The Bayesian update for the posterior distribution of the parameters is formulated as

$$P(\alpha, \mu, \gamma | X) \propto P(X | \alpha, \mu, \gamma) P(\alpha, \mu, \gamma) \quad (7)$$

The posterior is proportional to the likelihood of observing the data, $P(X|\alpha, \mu, \gamma)$, times the prior distribution of the parameters, $P(\alpha, \mu, \gamma)$. X represents the collection of H years of rupture history in the earthquake catalog for all the fault sections: $X = \{X_1, X_2, \dots, X_H\}$. The vectors α and μ are the collections of α_j and μ_j, respectively, for all the fault sections.

### Prior

The prior distributions of α_j, μ_j and γ are assumed to be log-normal probability distributions. Because these parameters can only take positive values, a lognormal probability distribution is considered suitable for the analysis. Their marginal probability distributions are

$$\alpha_j \sim P[\alpha_j] : \text{Lognormal}(\mu_{\log \alpha_j}, \sigma_{\log \alpha_j}) \quad (8)$$

$$\mu_j \sim P[\mu_j] : \text{Lognormal}(\mu_{\log \mu_j}, \sigma_{\log \mu_j}) \quad (9)$$

$$\gamma \sim P[\gamma] : \text{Lognormal}(\mu_{\log \gamma}, \sigma_{\log \gamma}) \quad (10)$$

We assumed the parameters to be independent in the prior distribution, thus the joint probability distribution can be written as

$$P(\alpha, \mu, \gamma) = P[\gamma] \prod_{j=1}^{N} P[\mu_j] P[\alpha_j] \quad (11)$$

### Likelihood of observing the data

The likelihood of observing the rupture history, $P[X|\alpha, \mu, \gamma]$, can be computed using the following formulation. The rupture history X is initially conditioned by $T_1$, the number of years since the last earthquake at the starting year of the catalog. Given that $T_{t+1}$ is a deterministic function of $T_t$ and $X_t$ and there is a one-to-one correspondence between $T_{t+1}$ and $\{T_t, X_t\}$ (see Equation 1), then
$$P[X|\alpha, \mu, \gamma] = P_{\alpha, \mu, \gamma}[X_1, X_2, \dots, X_H|T_1] = P[X_1, T_2, X_2, \dots, T_H, X_H|T_1] \quad (12)$$

In addition, the set {$T_{t+1}, X_t$} is a Markov chain, because if {$T_{t+1}, X_t$} is conditioned on the last step {$T_t, X_{t-1}$}, it is independent of all the previous rupture history. This is because the set {$T_t, X_{t-1}$} contains all the required information to assess the rupture in the next year (i.e., whether a rupture occurred the previous year and the time since the last rupture in each fault section). Therefore

$$P[X|\alpha, \mu, \gamma] = P_{\alpha, \mu, \gamma}[T_2, X_1|T_1] \prod_{t=2}^{H} P_{\alpha, \mu, \gamma}[T_{t+1}, X_t|T_t, X_{t-1}] \quad (13)$$

Finally, given the one-to-one correspondence between $T_{t+1}$ and {$T_t, X_t$} described previously, $T_{t+1}$ can be dropped from the conditional probabilities. $X_{t-1}$ can also be dropped from the conditional set {$T_t, X_{t-1}$}, because $T_t$, the time since the last rupture in each section, is the only information that is needed to evaluate the likelihood of $X_t$ (see Equation 2). Thus

$$P[X|\alpha, \mu, \gamma] = P_{\alpha, \mu, \gamma}[X_1|T_1] \prod_{t=2}^{H} P_{\alpha, \mu, \gamma}[X_t|T_t] \quad (14)$$

We evaluate $P_{\alpha,\mu,\gamma}[X_t|T_t]$ using the copula method described earlier. $P_{\alpha,\mu,\gamma}[X_t|T_t]$ is estimated as the probability that the vector $Z_t$ (from Equation 6) is within the region defined in Equation 15, where $Z_t(j)$ is the j-th element of the vector $Z_t$ that follows a normal distribution. This region is the intersection of the regions $A_j$, where $A_j$ extends over the set {$Z_t(j) \le \Phi^{-1}(p_t(j))$} if there is a rupture in the j-th section at year t, or extends over the set {$Z_t(j) > \Phi^{-1}(p_t(j))$} otherwise.

$$P[X_t|T_t] = P[\cap_{j=1}^N A_j], \text{ where } \begin{cases} A_j = \{Z_t(j) \le \Phi^{-1}(p_t(j))\} & \text{if } X_t(j) = 1, \\ A_j = \{Z_t(j) > \Phi^{-1}(p_t(j))\} & \text{otherwise} \end{cases} \quad (15)$$

## Posterior Estimation using Markov Chain Monte Carlo

As shown in Equation 7, the posterior distribution can only be partially evaluated by multiplying the prior distribution by the likelihood of observing the data. This product is proportional to the posterior distribution of the parameters. Thus, in order to evaluate the posterior with a closed-form equation, a normalizing factor must be estimated. However, because of the complexity and high dimensionality of the posterior distribution, it is not feasible to estimate such a normalizing factor by direct methods such as integration. Instead, in this paper we use the MCMC technique as it allows us to sample from and evaluate complex, high-dimensional distributions (Liu, 2004), avoiding the numerical challenges stemming from the highly dimensional integration.

We used the Metropolis-Hastings (MH) version of the MCMC algorithm because it is an approach designed to compute complex posterior distributions as in our case. Liu (2004) provides a detailed description of the MH algorithm. Briefly, the MH searches on the high-dimensional parameter vector space using a Markov Chain (MC) whose stationary distribution is the posterior distribution of the parameters. Algorithm 1 describes the procedure for generating samples from the posterior distribution. In order to generate a new parameter sample {$\alpha, \mu, \gamma$}_m, we estimate a sample candidate {$\alpha, \mu, \gamma$}^* using a random walk with an uncorrelated multivariate normal distribution

$$[\alpha, \mu, \gamma]^* \sim Q([\alpha, \mu, \gamma]|[\alpha, \mu, \gamma]_{m-1}) = \mathcal{N}([\alpha, \mu, \gamma]_{m-1}, \text{diag}([\sigma_{\alpha}^2, \sigma_{\mu}^2 \sigma_{\gamma}^2])) \quad (16)$$

In the random walk, the normal distribution has a mean that equals the values from parameter sample in the last step and has a covariance matrix that equals a diagonal matrix diag([σ_α^2, σ_μ^2 σ_γ^2]), i.e., the walk has fixed variances in all parameters at each step. Then, we estimate an acceptance rate A as

$$A = \frac{P([\alpha, \mu, \gamma]^*|X) \times Q([\alpha, \mu, \gamma]^*|[\alpha, \mu, \gamma]_{m-1})}{P([\alpha, \mu, \gamma]_{m-1}|X) \times Q([\alpha, \mu, \gamma]_{m-1}|[\alpha, \mu, \gamma]^*)} \quad (17)$$

Because the random walk has symmetrical probabilities, $Q([\alpha, \mu, \gamma]_{m-1}|[\alpha, \mu, \gamma]^*) = Q([\alpha, \mu, \gamma]^*|[\alpha, \mu, \gamma]_{m-1})$. Thus

$$A = \frac{P([\alpha, \mu, \gamma]^*|X)}{P([\alpha, \mu, \gamma]_{m-1}|X)} \quad (18)$$
Finally, the sample candidate is accepted with probability A. According to the MH properties, the samples will eventually converge to the posterior distribution of the parameters when the Markov chain reaches stationarity. Therefore, the first B samples belong to the burning period and are discarded. Only M-B samples are used to evaluate the posterior distribution.

<caption>**Algorithm 1** Metropolis-Hastings for Bayesian Parameter Estimation of Earthquake Rupture Model</caption>

1: **procedure** BAYESIAN PARAMETER UPDATE(X)
2:   **for** each $m$ in {1, 2, ..., $M$} **do**
3:     Sample $[\alpha, \mu, \gamma]^* \sim Q([\alpha, \mu, \gamma]|[\alpha, \mu, \gamma]_{m-1})$ as in Equation 16
4:     Calculate $A$ as in Equation 17
5:     Sample $U \sim \text{Uniform}(0,1)$
6:     **if** $U < A$ **then**
7:       $[\alpha, \mu, \gamma]_m = [\alpha, \mu, \gamma]^*$
8:       **else**
9:         $[\alpha, \mu, \gamma]_m = [\alpha, \mu, \gamma]_{m-1}$
10:      **end if**
11:    **end for**
12:  Use the last $M-B$ samples to estimate $P([\alpha, \mu, \gamma|X])$
13: **end procedure**

## MODEL SETUP AND EARTHQUAKE DATA

The parameter estimation technique was applied to the central portion of the subduction zone of Peru, where the oceanic Nazca plate is subducting beneath the continental South American plate at an average rate of 6 cm/yr (Kendrick et al., 2003). Figure 3 shows the seismotectonic setting of the region of analysis and the rupture areas of the four last large earthquakes that occurred in 1940, 1966, 1970, and 1974. According to historical reports and instrumental catalogs, the region has been very active and has generated more than ten earthquakes with magnitudes larger than 7.5 in the last 450 years (Dorbath et al., 1990; Villegas-Lanza et al., 2016; Ceferino et al., 2018a). This region has a size of 650 km along the strike direction, 200 km along the dip direction, and a dip angle of 15°, reaching depths of up to 50 km. The northern and southern boundaries of the selected seismic region correspond to the areas of significant fault creep inferred by Villegas-Lanza et al. (2016) using geodetic data.

The tectonic region was discretized into eight sections along the strike direction as shown by the dashed lines in Figure 3. Thus, in this application, the model has 17 parameters: one $\mu_j$ and $\alpha_j$ per each section, and $\gamma$ for the correlogram. The sections are labeled from 1 to 8 starting from south to north. Each section has approximately 81.3 km along the strike direction, whose rupture corresponds to an earthquake magnitude of approximately 7.5 according to the scaling law proposed by Strasser et al. (2010). Because the fault sections represent the minimum tectonic area unit, the model only captures earthquakes with magnitude equal to or larger than 7.5.

The following subsections describe the earthquake data used in the two-step parameter estimation in the study region.

### Historical Catalog

We used the historical catalog compiled by Ceferino et al. (2018b). It contains information on the earthquake rupture areas, locations, and magnitudes during the last 450 years. For earthquakes before 1940, the catalog uses the area, location, and magnitude estimates inferred by Dorbath et al. (1990) from isoseismal maps. For seismic events starting in 1940, the catalog uses earthquake estimates from aftershock information and other instrumental data (Kelleher, 1972; Kanamori, 1977; Langer and Spence, 1995; Chlieh et al., 2011). We included in the historical catalog two additional earthquakes that occurred in 1806 and 1828, respectively. Although several studies omit these two earthquakes, Silgado (1978) and Seiner (2011) showed that they were catastrophic events. Villegas-Lanza et al. (2016) estimated the magnitude, locations and sizes of these two events.

Next, we projected all rupture areas along the strike direction and discretized them to match the fault sections defined in Figure 3. Figure 4 shows these rupture projections in the Y axis and the earthquake occurrence year in the X axis.
<figure><img src="image_5.png" /></figure>

Figure 3: Central Peru subduction zone and its recent seismicity. The line parallel to the coastline shows the fault trench where the Nazca and South American plates converge. The black quadrilateral shows the region of study subdivided into eight sections, and the four enclosed areas show the rupture areas of the four most recent large earthquakes, with the year and magnitude indicated in each area. The arrows indicate the plate convergence rates in mm/year. Modified from Ceferino et al. (2018b)

## Synthetic Catalog Using Physics-based Simulation

We generated the synthetic catalog of seismic events using a physics-based model for earthquake cycles. In this model, the tectonic load resulting from the convergence between plates and fault slip is balanced by a fault friction resistance governed by the rate-and-state friction law (Dieterich, 1979; Tullis, 1988; Marone and Kilgoret, 1993). The black solid line in Figure 3 represents the fault interface as a 2D frictional plane embedded in a 3D elastic medium dipping 15°. As a result of the tectonic load, some regions on the fault become unstable, nucleate and trigger earthquakes of different magnitudes (Rubin and Ampuero, 2005). We used the implementation in the software QDYN to carry out the physics-based modeling (Luo et al., 2017). To limit computational cost, this physics-based model necessarily contains simplifying assumptions and approximations, like the use of rate-and-state friction even during fast sliding stages. However, it maintains the advantage of tracking the stress history on the fault, accounting for the quasi-static elastic interactions between asperities mediated by co-seismic stress transfer and intervening transient creep, and utilizing it to compute the evolution of the future earthquakes.

### Fault asperities

We represent the main fault asperities that produce earthquakes by assigning heterogeneous fault friction properties defining four velocity-weakening (VW) regions on the fault as shown in Figure 5. These VW regions are locked in the interseismic period between earthquakes. Their shapes were idealized as ellipses with length $L_{asp}$ along the strike direction and width $W_{asp}$ along the dip direction. They were located according to the rupture areas of the last four large earthquakes shown in Figure 3. We assume that these previous large rupture areas revealed the main asperities.

### Input data for the physics-based model

We set the reference steady-state slip velocity equal to 6 cm/year, which is the convergence rate of the tectonic plates. Other input data were defined according to commonly used values in previous studies in subduction zone regions (Liu and Rice (2007)). The values are shown in Table 1. We assume that these values are constant through the entire fault, except for the state-effect frictional coefficient $b$. The
<figure><img src="image_6.png" /><figcaption>Figure 4: Along-strike distribution of historical earthquake ruptures through time. The rupture lengths were discretized to match the fault sections.</figcaption></figure>

sign of the difference between *b* and the direct-effect frictional coefficient *a* determines whether a region has VW behavior (*b* − *a* > 0) or velocity-strengthening (VS) behavior (*b* − *a* < 0). The VS regions tend to creep and do not nucleate earthquakes, but the areas of partial interseismic locking generally extend out of the VW areas into the VS areas. In order to have smoothness, we applied an exponential spatial transition of *b* values between the VW and VS regions.

<table><caption>Table 1: Input data summary for the physics-based model.</caption><thead><tr><th>Input Data Description</th><th>Value</th></tr></thead><tbody><tr><td>Shear modulus (G)</td><td>32.5 GPa</td></tr><tr><td>Elastic modulus (E)</td><td>32.5 GPa</td></tr><tr><td>Shear wave speed (V<sub>s</sub>)</td><td>3000 m/s</td></tr><tr><td>Initial effective normal stress (σ)</td><td>100 MPa</td></tr><tr><td>Reference frictional coefficient (f)</td><td>0.6</td></tr><tr><td>Direct-effect frictional coefficient (a)</td><td>see Table 2</td></tr><tr><td>State-effect frictional coefficient in VS (b)</td><td>see Table 2</td></tr><tr><td>State-effect frictional coefficient in VW (b)</td><td>see Table 2</td></tr><tr><td>Characteristic slip distance (D<sub>c</sub>)</td><td>see Table 2</td></tr><tr><td>Reference steady-state slip velocity (ν<sub>ss</sub>)</td><td>6 cm/year</td></tr><tr><td>Fault length along the strike direction (L)</td><td>650 km</td></tr><tr><td>Fault width along the dip direction (W)</td><td>200 km</td></tr><tr><td>Dip angle (φ)</td><td>15°</td></tr><tr><td>Grid size along the strike direction (dl)</td><td>634 m</td></tr><tr><td>Grid size along the dip direction (dw)</td><td>809 m</td></tr></tbody></table>

We tested multiple values of *a*, *b*, and characteristic slip distance $D_c$, under the following constraints. The parameter $b-a$ in the VW regions controls the stress drop during earthquakes. It directly influences the earthquake magnitude and frequency distribution in the simulation. The ratio $D_c/b$ controls the required mesh size in the simulation, as described in the next section. These parameters were varied in order to ensure that the results capture the main features of the previous historical seismicity and satisfy computational constraints. Table 2 shows the three cases that best represented the seismicity in the region while capturing different possible earthquake behaviors.

<table><caption>Table 2: Values of *a*, *b* and $D_c$ for the three selected models.</caption><thead><tr><th>Input Data</th><th>Case A</th><th>Case B</th><th>Case C</th></tr></thead><tbody><tr><td>a</td><td>0.0014</td><td>0.0028</td><td>0.0042</td></tr><tr><td>b (VS)</td><td>0.0008</td><td>0.0017</td><td>0.0025</td></tr><tr><td>b (VW)</td><td>0.0028</td><td>0.0056</td><td>0.0084</td></tr><tr><td>D<sub>c</sub></td><td>2.3 cm</td><td>4.6 cm</td><td>6.8 cm</td></tr></tbody></table>
<figure><img src="image_7.png" /></figure>

Figure 5: Distribution of asperity regions in the physics-based model. The shaded, rectangular area represents the extent of the 2D frictional area. Its size is 650 km along the strike direction and 200 km along the dip direction. The four VW asperity regions with high interseismic coupling are represented by ellipses with length $L_{asp}$ along the strike direction and $W_{asp}$ along the dip direction and are shown in dark shades. The VS areas are shown in a white shade. The grey scale indicates the b/a ratio of the rate-and-state friction law.

## Computational constraints

The physics-based model is computationally demanding both in terms of memory usage and execution time. A large stiffness matrix $K$ stores the stresses due to unitary slips at every fault node. The size of $K$ is $n \times n$, where $n$ is the number of fault nodes. The stress at each time step is updated by multiplying $K$ with the vector $V$, which contains the slip velocities at all fault nodes. The software QDYN uses the Message Passing Interface (MPI) and Fast Fourier Transform (FFT) methods to compute the matrix multiplication efficiently. We utilized 256 processors of the Stanford Sherlock computer cluster. The computational time increases significantly with reduced grid sizes. The model running time was approximately seven days for the grid size shown in Table 1 for 2,000 years of simulation. In order to have a sound physics-based model, the model parameters in this application were set up according to the two following constraints.

First, the grid sizes $dl$ and $dw$ along the strike and dip directions, respectively, were set such that they were at least three times smaller than the cohesive zone length (Day et al., 2005; Lapusta and Liu, 2009). The cohesive zone length is comparable to $L_b = \frac{G' D_c}{\sigma b}$, where $G'$ is the shear modulus for anti-plane rupture, and equal to the shear modulus $G$ divided by one minus the Poisson's ratio for in-plane rupture. The values of $G$, $D_c$, $\sigma$, and $b$ are shown in Table 1. Because $b$ is larger in the VW regions, these regions impose heavier constraints on the grid size than the VS regions. We chose the grid sizes such that $dl/L_b \approx 0.3$ and $dw/L_b \approx 0.25$. Either decreasing the grid sizes or increasing $L_b$ would have resulted in smaller ratios. However, computational constraints did not allow to further decrease the grid sizes, and the second model constraint described in the next paragraph did not allow to increase $L_b$. Implementations that use the hierarchical matrix method could enable simulations with smaller grid size (Bradley, 2014).

Second, the dimensions of the VW regions, $L_{asp}$ and $W_{asp}$, were set such that they were significantly larger than the nucleation length $L_\infty = (b/b-a)^2(L_b/\pi)$ derived by Ampuero and Rubin (2008). Here $L_{asp}/L_\infty$ and $W_{asp}/L_\infty$ range from 35 to 80 for the four VW regions. Smaller ratios resulted in simulations with very large earthquake ruptures breaking all asperities in each event. Because this behavior is unsuitable for generating a catalog with multiple earthquake magnitudes, the ratios were not further
reduced, and thus, $L_b$ could not be increased.

## Sensitivity to the asperity size

The sensitivity to different asperity sizes was analyzed by varying the sizes of the four VW regions. As described earlier, the four VW regions were idealized as ellipses with size and location corresponding to the last four large earthquakes that occurred in the region. Initially, the ellipses' centers were located at the rupture areas' centroids. Then, the ellipses' locations and sizes were varied until the simulated interseismic coupling (ISC) areas matched the results inferred from recent GPS measurements of crustal motions (Villegas-Lanza et al., 2016) as will be described later.

In this paper, three representative cases with different VW geometries are described in Table 3. Because each VW region was assigned according to the location of a previous rupture, they were labeled according to their respective earthquake occurrence year. In the three cases, the VW regions have the same centers. The ellipses' diameters of cases 2 and 3 are 90% and 80% of the diameters of case 1, respectively, to show the effect of varying asperity sizes. Figure 5 shows the extent of the VW regions for case 2.

<table><caption>Table 3: Three representative asperity sizes tested in the model calibration. Lengths are provided in km.</caption><thead><tr><th rowspan="2">Asperity</th><th colspan="2">Center</th><th colspan="2">Case 1</th><th colspan="2">Case 2</th><th colspan="2">Case 3</th></tr><tr><th>lat</th><th>lon</th><th>L<sub>asp</sub></th><th>W<sub>asp</sub></th><th>L<sub>asp</sub></th><th>W<sub>asp</sub></th><th>L<sub>asp</sub></th><th>W<sub>asp</sub></th></tr></thead><tbody><tr><td>1940</td><td>-77.22°</td><td>-11.42°</td><td>255</td><td>115</td><td>229</td><td>104</td><td>204</td><td>92</td></tr><tr><td>1966</td><td>-79.16°</td><td>-10.57°</td><td>150</td><td>135</td><td>135</td><td>122</td><td>120</td><td>108</td></tr><tr><td>1974</td><td>-78.20°</td><td>-12.51°</td><td>280</td><td>118</td><td>252</td><td>107</td><td>224</td><td>107</td></tr><tr><td>2007</td><td>-76.79°</td><td>-13.91°</td><td>133</td><td>128</td><td>119</td><td>116</td><td>106</td><td>102</td></tr></tbody></table>

## Discussion on the resulting physics-based synthetic catalog

In this subsection, we present the results of the physics-based simulation for the cases A-1, A-2, A-3, B-1, B-2, B-3, C-1, C-2 and C-3 obtained from varying the values of *a*, *b*, and $D_c$ and the asperity sizes according to Tables 2 and 3. We selected these cases to represent a broad set of feasible earthquake behaviors in the physics-based simulation. The physics-based simulation outputs instantaneous slip rate throughout the fault over the 2,000 years of simulation. During nucleation the slip rate increases significantly and reaches values above 0.1 m/s during the coseismic phase. We define the beginning of an earthquake when the slip rate at any point exceeds 1.6 cm/s. This threshold value is 10% of the slip rate at which inertial effects (radiation damping) start to dominate (Rubin and Ampuero, 2005). While high slip rates propagate through the fault, the total earthquake slip accumulates in the rupture area. The earthquake is considered to be finished once all points have slip rates below 1.6 cm/s. The final earthquake area is estimated as the area where the final earthquake slip exceeded 10% of the maximum earthquake slip. For example, Figure 6 shows an earthquake rupture area from the simulation in case C-2. In this case the earthquake rupture area is 1,533 km², and the magnitude is 8.4. This earthquake occurs in the VW area 2007, where the 2007 Pisco earthquake occurred. This earthquake occurs at year 656 in the simulation, and its average stress drop is 2.8 MPa.

In general, the earthquakes in our simulation have longer durations than historical earthquakes. For example, the magnitude 8.4 earthquake shown in Figure 6 has a duration of about 4 minutes, and other earthquakes of similar magnitude reach durations of 7 minutes. Such long durations stem from the quasi-dynamic approximation in the physics-based model, which combines quasi-static elasticity with radiation damping. Despite this limitation of the model, we obtained reasonable rupture areas and slips, which are the key parameters to compute the earthquake magnitudes in the synthetic catalog.

We also observed that some sequences of simulated earthquakes occurred within very short timespans. For example, in the simulation B-2, a magnitude 7.7 earthquake is followed by a magnitude 8.5 earthquake within 6 minutes. Considering the quasi-dynamic approximation of the model, it is not possible to conclude whether these events are two different earthquakes or just a single larger and longer earthquake.

Because the focus of this paper is to leverage the synthetic simulations to reduce parameter uncertainty in our time-dependent hazard model, rather than answering fundamental questions arising from physics-based simulations, we considered two extreme cases to represent the spectrum of possible earthquake durations. In the first approach, we assumed that all seismic events in the simulation are separate earthquakes. In the second approach, we grouped together the seismic events that have intersecting
<figure><img src="image_8.png" /><figcaption>Figure 6: One example of simulated earthquake slip from case C-2 (29.7). The event occurs at year 656 in the simulation and has a magnitude of 8.4.</figcaption></figure>

rupture areas and that start within two weeks from each other. Figure 7 shows the synthetic catalogs for the 2000 years of simulation for the second approach. In these models, seismic events can break multiple asperities at once, causing earthquakes with a wide range of magnitudes. We observed the following two main features from the earthquake catalogs in the sensitivity analysis.

First, earthquakes in cases A, B, and C have different shear stress drops. From the simulation outputs, we defined stress drops as the difference between the shear stress at the end and start of each earthquake, spatially-averaged within its rupture area. The stress drops were controlled by the values of b, a, and σ, roughly consistent with the expectation in rate-and-state friction models that stress drops are approximately proportional to σ(b − a) (Ampuero and Rubin (2008)). According to the values in Table 2, the stress drops in case A should be half those in case B and a quarter of those in case C. In cases A-1, A-2 and A-3, the median drops were 0.58, 0.65, and 0.64 MPa, respectively, with some earthquake drops as low as 0.38 MPa and as high as 1.12 MPa. In cases B-1, B-2 and B-3, the median drops were 1.2, 1.4, and 1.5 MPa, respectively, with some earthquake drops as low as 0.8 MPa and as high as 2.1 MPa. In cases C-1, C-2 and C-3, the median drops were all about 2 MPa, with some earthquake drops as low as 1.1 MPa and as high as 3.3 MPa. After analyzing 860 historical earthquakes worldwide, Allmann and Shearer (2009) determined that in subduction zones the median earthquake stress drop is about 3 MPa and the standard deviation is 0.2 MPa. Though the median stress drops in the simulations fall in the lower tail of the stress drop distribution of real earthquakes, especially in cases B and C, we still considered that these simulations are plausible representations of the seismicity in the region. The quasi-dynamic approximation in the physics-based model can induce a bias towards lower stress drops, because it does not include the extra stress drop due to the dynamic overshooting effect of seismic waves (Madariaga, 1979).

Second, earthquakes in cases 1, 2 and 3 have different moment magnitude distributions. By construc-
tion, the asperities in case 1 are larger and closer to each other than those in case 2, which in turn are
larger and closer than in case 3. As a result, earthquakes in case 1 involve multiple asperities more often
than in cases 2 and 3, generating earthquakes with larger magnitudes. For large magnitudes, the annual
magnitude exceedance rates are the highest in case 1 and the lowest in case 3 (Figure 8).

Figure 8 also shows the comparison of the exceedance rates from the synthetic and historical catalogs. The agreement is better for magnitudes lower than 8.0 than for larger magnitudes. Such non-uniform matching is a desirable feature in synthetic catalogs to better represent the uncertainty in earthquake occurrences because the rates derived from the historical catalog also exhibit increasing uncertainty at larger magnitudes (Figure 8). A historical catalog spanning thousands of years would be required to
<figure><img src="image_9.png" /><figcaption>Figure 7: Physics-based synthetic earthquake catalogs for nine cases. Earthquake rupture lengths in the along-strike direction are shown as gray vertical lines at their time of occurrence. For each event, a symbol indicates which VW regions broke (see legend). The red vertical lines in Cases B-1, B-2 and B-3 indicate snapshot times selected for the comparison of ISC ratios in Figure 9.</figcaption></figure>

obtain rate estimates with high statistical confidence. However, in our study region the historical catalog spans only 450 years. It contains only two earthquakes that can be used to estimate the exceedance rates for magnitudes higher than 8.2. Additionally, there is large uncertainty on the magnitude estimation. Earthquakes with magnitudes larger than 8.0 were observed before instrumentation was available, thus their magnitudes were inferred through a combination of subjective estimates of isoseismals and highly uncertain empirical formulations (Dorbath et al., 1990). Thus, the differences between the rates for large magnitude earthquakes in the historical and synthetic catalogs are a strength of the physics-based simulation rather than a drawback, as they can provide new information on the earthquake rupture occurrences that are not observed in the reduced timespans of historical catalogs.

Finally, we verified that the snapshots of ISC in the synthetic catalog match the high-coupling regions inferred from GPS data. The ISC ratio can be estimated at any point of the fault as $(\nu - \nu_{ss})/\nu_{ss}$, where $\nu$ is the instantaneous slip rate and $\nu_{ss}$ is the steady-state (long-term) slip rate imposed by the plate convergence (from Table 1). When the ISC ratio equals 0 there is total creep, and when the ISC ratio equals 1 there is full locking. Villegas-Lanza et al. (2016) utilized multiple GPS velocity measurements collected between 2008 and 2013 along the South American plate to infer spatially varying ISC through an inversion procedure. Figure 9a shows one of the resulting ISC maps that best fitted the GPS data in Villegas-Lanza et al. (2016). ISC ratios are high in a stripe roughly parallel to the fault strike direction along most of the fault. The stripe intersects the VW regions 1966, 1974 and 2007 of our model; however, it does not completely cover all of them. Moreover, the GPS-based ISC shows high ratios in the central
<figure><img src="image_10.png" /></figure>

Figure 8: Annual magnitude exceedance rates from the synthetic, physics-based catalog and from the historical catalog. Eight of the nine synthetic catalogs are shown in light lines, and case B-2 is highlighted in a different color. Case B-2 will be used in the next section to perform the two-step Bayesian parameter estimation. The bell curves show that there are higher uncertainties for larger earthquakes in both the magnitude axis and the rate axis.

east part of the fault where the VW region of 1940 is located. Because we selected VW regions based on previous earthquakes, these intersections show that where previous large earthquakes occurred, currently there is high ISC.

The ISC varies over time as the fault undergoes multiple complex processes of interaction between seismic and aseismic slip, and it heavily depends on the preceding earthquake cycles. Consequently, we seek earthquake cycles in the synthetic catalog that are similar to ones preceding the GPS data collection (before 2008). For example, before the data collection, the 1966 earthquake ruptured the northern sections of the study area, the 1940 and 1974 earthquakes the central sections, and the recent 2007 earthquake ruptured the southern sections (Figure 3). We selected snapshots in the physics-based simulations within the first five years after a large rupture in the southern sections that were preceded by ruptures in the central and northern sections in previous decades. We observed that the high-coupling regions have strong correspondence with the predefined VW regions, but extend beyond (Figure 9). In the simulations that have different shear stress drops but identical VW regions, e.g., cases A-2, B-2 and C-2, the high ISC regions are similar. In contrast, in the cases with similar stress drops but different asperity sizes the high-coupling regions were different. We found that larger VW regions have larger areas with high coupling. For example, Figures 9b, 9c and 9d show the ISC ratio snapshots at years 1214.2, 1501.8 and 1710.4 for cases B-1, B-2 and B-3, respectively.

The comparison between the physics-based simulations illustrates that in case B-1 the regions with high ISC ratios are more extensive than in case B-2, and similarly the regions in case B-2 are more extensive than in B-3. Additionally, the simulated ISC ratios have high ratios along the west part of the study area due to a boundary effect in the physics-based simulation. Models B-1 and B-2 cover most of the regions with high GPS-based ISC. However, case B-1 has regions with high ISC ratios that extend over larger areas than the GPS-based ISC. Case B-3 has four pockets with high ISC ratios that fall mostly within the GPS-based high ISC regions, but that cover significantly less area than the GPS-based ISC high ratios. Because cases B-1 and B-3 have significantly larger and smaller areas than the GPS-based ISC map, we selected case B-2 as the best representation of the ISC coupling. Case B-2 does not fully show a stripe with high ISC as found in the GPS-based ISC, whose stripe has ratios larger 0.8 (Figure 9a). Yet, case B-2 is close to representing this stripe because the northern VW region of 1966 is connected to the central VW region of 1974 with a high ISC coupling bridge with ratios larger than 0.6, and the distance between the high ISC areas of this central VW region and the southern VW region of 2007 is relatively small.

## PARAMETER ESTIMATION APPLICATION

We utilize the synthetic catalog in case B-2 together with the historical catalog to perform the two-step Bayesian update shown in Figure 2. A similar approach using multiple-step Bayesian updates can incorporate the historical catalog and all the synthetic catalogs shown in Figure 7 into the parameter estimation; however, because the goal of this paper is to show how to combine historical and synthetic catalogs, we use only the catalog in case B-2 to provide the proof of concept for the parameter estimation.

As described earlier, we use independent log-normal distributions to model the prior distributions of
<figure><img src="image_11.png" /></figure>

Figure 9: Comparison of observed and modeled ISC ratios. (a) ISC ratio inferred from GPS data collected from 2008 to 2013 (Villegas-Lanza et al. (2016)'s model 4.5). (b–d) ISC ratios from three physics-based simulations, estimated after a sequence of large earthquakes similar to the real one before 2008. The dashed ellipses show the boundaries of the VW asperity, defined by $b-a = 0$. High ISC areas can extend beyond the asperity boundaries.
<table><caption>Table 4: Double update of earthquake rupture parameters</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Prior</th><th colspan="2">1st Update</th><th colspan="2">2nd Update</th></tr><tr><th>Median</th><th>Std. Dev.</th><th>Median</th><th>Std. Dev.</th><th>Median</th><th>Std. Dev.</th></tr></thead><tbody><tr><td>μ<sub>1</sub></td><td>175</td><td>56</td><td>159</td><td>16</td><td>169</td><td>17</td></tr><tr><td>μ<sub>2</sub></td><td>175</td><td>56</td><td>159</td><td>16</td><td>153</td><td>14</td></tr><tr><td>μ<sub>3</sub></td><td>175</td><td>56</td><td>189</td><td>22</td><td>184</td><td>18</td></tr><tr><td>μ<sub>4</sub></td><td>175</td><td>56</td><td>169</td><td>21</td><td>149</td><td>17</td></tr><tr><td>μ<sub>5</sub></td><td>175</td><td>56</td><td>169</td><td>19</td><td>151</td><td>15</td></tr><tr><td>μ<sub>6</sub></td><td>175</td><td>56</td><td>277</td><td>64</td><td>228</td><td>41</td></tr><tr><td>μ<sub>7</sub></td><td>175</td><td>56</td><td>177</td><td>23</td><td>190</td><td>21</td></tr><tr><td>μ<sub>8</sub></td><td>175</td><td>56</td><td>184</td><td>27</td><td>168</td><td>24</td></tr><tr><td>α<sub>1</sub></td><td>0.70</td><td>0.22</td><td>0.43</td><td>0.09</td><td>0.54</td><td>0.10</td></tr><tr><td>α<sub>2</sub></td><td>0.70</td><td>0.22</td><td>0.46</td><td>0.09</td><td>0.63</td><td>0.09</td></tr><tr><td>α<sub>3</sub></td><td>0.70</td><td>0.22</td><td>0.53</td><td>0.10</td><td>0.63</td><td>0.09</td></tr><tr><td>α<sub>4</sub></td><td>0.70</td><td>0.22</td><td>0.54</td><td>0.11</td><td>0.66</td><td>0.09</td></tr><tr><td>α<sub>5</sub></td><td>0.70</td><td>0.22</td><td>0.47</td><td>0.09</td><td>0.56</td><td>0.07</td></tr><tr><td>α<sub>6</sub></td><td>0.70</td><td>0.22</td><td>0.72</td><td>0.20</td><td>0.95</td><td>0.17</td></tr><tr><td>α<sub>7</sub></td><td>0.70</td><td>0.22</td><td>0.50</td><td>0.14</td><td>0.56</td><td>0.11</td></tr><tr><td>α<sub>8</sub></td><td>0.70</td><td>0.22</td><td>0.50</td><td>0.12</td><td>0.80</td><td>0.12</td></tr><tr><td>γ</td><td>375</td><td>120</td><td>456</td><td>53</td><td>407</td><td>36</td></tr></tbody></table>

the parameters. For interevent time mean $\mu_j$, the prior median was set to 175 years, which represents the average rupture time in the historical catalog. The prior logarithmic standard deviation was set to 0.3 in all the fault sections. For the interevent time coefficient of variation $\alpha_j$, the prior median was set to 0.7, a value that has been used extensively in BPT distributions and earthquake hazard models (Field et al., 2015). The prior logarithmic standard deviation was set to 0.3 in all the fault sections. For the correlation length $\gamma$, the prior median was set to 375 km, which equals the calibrated values found in previous analyses of the seismic region (Ceferino et al., 2018a). The prior logarithmic standard deviation was also set to 0.3. Table 4 lists all the parameters of the prior distributions, and Figure 10a shows the probability density function (pdf) of the joint prior distribution for the mean interevent times at sections 4 and 5, $\mu_4$ and $\mu_5$. The standard deviations in Table 4 were estimated using both the median and the logarithmic standard deviations of the lognormal prior distributions.

In the first update, we use the synthetic catalog B-2 that was computed based on the two-week duration window. Before the update, we preprocessed the catalog by removing two seismic events that we considered to be aftershocks because they occurred in the same regions and shortly after larger earthquakes. Additionally, because the probabilistic model has one-year time steps, we adjusted the occurrence times of earthquakes in the synthetic catalog to be at least one year apart. For example, if two different earthquakes occurred during the same year, we assigned them to two consecutive year indexes instead of to the same index. After prepocessing, we ran the first Bayesian update using Algorithm 1. We obtained 10,000 samples of the parameters and calibrated the random walk to achieve a 25% acceptance rate. This rate falls within the desirable range for good sample mixing, i.e., effective exploration of the high-probability regions in the high-dimensional parameter space (Chib and Greenberg (1995); Robert (2014)). We deemed that the first 900 samples were part of the burning period because they still belonged to low-probability regions of the parameter space. Table 4 shows the list of means and standard deviations estimated with the samples of all parameters, and Figure 10b shows the 2D joint distributions of the parameters $\mu_4$ and $\mu_5$ after the update, which was estimated using a Gaussian-kernel pdf estimator on the samples. The comparison between the prior and the first update in both the table and the plot shows that the first update is effective at reducing the uncertainty in the parameters. For example, the standard deviation of $\mu_4$, $\alpha_4$, and $\gamma$ decreased by 62, 53, and 56%, respectively. Additionally, the first update is also effective at correcting the initial biases in the prior distribution as the means of the parameters change by incorporating the information in the synthetic catalog data.

In the second update, we use the historical catalog. To follow Algorithm 1 in the second update, we fitted independent lognormal distributions to the samples of the first update in order to use them as prior distributions. Similarly to the first update, we obtained 10,000 samples and calibrated the random walk to achieve 25% acceptance rate. Table 4 shows the list of medians and standard deviations estimated with the samples from the second update for all the parameters. Figure 10c shows the 2D joint distribution of
<figure><img src="image_12.png" /><figcaption>Figure 10: Two-step parameter update for interevent means in sections 4 and 5, $\mu_4$ and $\mu_5$. We show the prior distribution based on expert opinion and previous studies in (a), the distribution update after incorporating the physics-based catalog in (b), and the distribution update after including the historical catalog in (c). The shades show the value of a Gaussian-kernel pdf estimator.</figcaption></figure>

the parameters $\mu_4$ and $\mu_5$ calculated with the Gaussian-kernel pdf estimator. Both Table 4 and Figure 10 indicate that the second update further reduces the uncertainty in the parameter estimates but to a lesser degree than the first update. For example, whereas in the first update the standard deviation of $\mu_4$, $\alpha_4$, and $\gamma$ decreased by 62, 53, and 56%, respectively, in the second update the standard deviation only decreased by 19, 10, and 32%, respectively. We observed that the second update is more effective at further calibrating the parameter estimates as the means of the parameters still improve their values by incorporating the historical catalog data.

To show how effective the two-step Bayesian update is at calibrating the probabilistic model, we compare the two-step update using both the historical and the synthetic catalogs with a one-step update using only the historical catalog. Through Algorithm 1, we obtained 10,000 samples of the parameter posterior using the prior distribution described in Table 4 and the historical catalog data using a one-step update. Figures 11a, 11b and 11c show the comparison between the prior, two-step and one-step updates for $\mu_4$, $\alpha_4$ and $\gamma$.

The one-step update shifted the median of $\mu_4$, $\alpha_4$ and $\gamma$ from 175 years, 0.7 and 375 km in the prior to 121 years, 0.71 and 390 km, respectively, whereas the two-step update, which included the additional synthetic data, shifted these medians to 149 years, 0.66 and 407 km, respectively. Similarly, the two-step update resulted in larger medians of $\mu_j$ and smaller medians of $\alpha_j$ in most sections. Although these relative differences in the medians were slight on average, some sections had significant differences, particularly the central sections. For example, whereas the medians of $\mu_j$ with the two-step update are only 5% higher than with the one-step update on average, in sections 4 and 5 the medians are up to 23% (Figure 11a) and 15% higher. The medians of $\mu_j$ in the two-step update, which incorporates both historical and B-2 data, tend to be higher than the ones in the one-step update, which incorporates historical data and not B-2 data, because sections break less recurrently in the B-2 data, especially in the central sections (Figure 7e). Looking to the entire fault, the occurrence rate is 0.027/yr in the historical data (i.e., 12 earthquakes in 450 years) and 0.020/yr. in the B-2 data (i.e., 39 earthquakes in 2000 years of simulation) for magnitudes larger than 7.5 in the fault. However, turning to section 4, the occurrence rate in the historical data (i.e., 6 earthquakes in 450 years) is twice the one in data B-2 (i.e., 13 earthquakes in 2000 years), thus B-2 data generate larger interevent times $\mu_j$ in the two-step update.

Both the one-step and the two-step updates reduce the standard deviation of the parameters. However, the two-step update is more effective at the reduction because it incorporates the additional information from the synthetic catalog. For example, the one-step update reduced the standard deviations of $\mu_4$, $\alpha_4$ and $\gamma$ from 56 years, 0.22 and 120 km in the prior to 32 years, 0.22 and 55 km, respectively, whereas the two-step reduced these standard deviations to 17 years, 0.09, and 36 km, respectively. On average, the two-step update provided an additional 45% in the uncertainty reduction of the parameters to the reduction provided by the one-step update, with some cases with additional reductions as high as 75% for $\mu_2$.

Finally, we evaluate how the discrepancies between the two-step and one-step update propagate to the annual magnitude exceedance rates and the time-dependent seismic hazard in the region. To compute the magnitude exceedance rates, we used a similar procedure as Ceferino et al. (2018b) with model parameters corresponding to the medians of the estimates, which were obtained previously with the two-step update and one-step update. Figure 12 shows the annual magnitude exceedance rates from
<figure><img src="image_13.png" /><figcaption>Figure 11: Comparison from synthetic catalog information in parameter estimation for $\mu_4$, $\alpha_4$ and $\gamma$. H: historical. S: synthetic.</figcaption></figure>

the synthetic and historical catalogs and from the models with the one-step and two-step updates. The resulting exceedance rates are a piecewise constant function because of the discrete nature of the model. The exceedance rates varies only at magnitude values that correspond to an integer number of fault sections rupturing. The two-step update resulted in slightly lower exceedance rates than the one-step update. For example, the exceedance annual rate of an earthquake with magnitude larger than 8.0 is 0.012 (i.e., return period of 82 years) and 0.013 (i.e., return period of 77 years) with the two- and one-step updates, respectively. The slightly lower rates are due to the slightly larger values of $\mu_j$ from the two-step update. Though the comparison of $\mu_j$ medians across sections showed that there are some localized high variations in the central sections between the one- and two-step update, the magnitude exceedance rates are not able to capture them because they measure the aggregated magnitude occurrence in all the fault rather than at individual sections.

<figure><img src="image_14.png" /><figcaption>Figure 12: Annual magnitude exceedance rates from the historical and synthetic (case 2-B) catalogs and from the model calibrated with a one-step (only H.: historical) and with a two-step (both S. & H.: synthetic and historical) updates.</figcaption></figure>

Next, we estimated the time-dependent seismic hazard as the likelihood of exceeding peak ground accelerations (PGA) of 0.4 g in the next 30 years due to earthquakes with magnitudes larger than 7.5. We also utilized a similar procedure as Ceferino et al. (2018b) and estimates corresponding to the medians of the parameters with the two-step and one-step updates. Figures 13a and 13b show the seismic hazard for one-step update and two-step update, respectively, and Figure 13c shows the ratio between both estimations. These plots indicate that the seismic hazard with the two-step update is 40% less than the hazard with one-step update in most of the region of analysis. These larger ratios are prevalent in the central region close to section 5, where the hazard has the highest values. The one-step update yields larger time-dependent hazard estimates because, unlike the two-step update, it has larger $\mu_j$ as a consequence of only incorporating the historical catalog and not the synthetic catalog. The high ratios are less prevalent in the southern region, where the recent 2007 earthquake occurred. The large relative differences in the south of the study area are not significant in absolute terms because the hazard is small there, as shown in Figures 13a and 13b.

# CONCLUSIONS

This paper presents a robust parameter estimation technique for a time- and space-interaction model of main-shock earthquake occurrence. Existing complex models that capture these interactions often use
<figure><img src="image_15.png" /><figcaption>Figure 13: Estimations of time-dependent seismic hazard as the probability of PGA exceeding 0.4 g in the next 30 years using the rupture model fitted with (a) only historical data and (b) a combination of synthetic and historical data. (c) Ratio of these two hazard estimations.</figcaption></figure>

simplified techniques that lead to unreliable parameter estimates. Additionally, these complex models utilize limited earthquake datasets that introduce high uncertainties in the parameter estimates, especially for large earthquake magnitudes. This paper addresses both issues by proposing a robust Bayesian formulation for parameter estimation and leveraging physics-based earthquake simulation to enlarge the existing datasets, which are currently limited.

Two steps compose the proposed parameter estimation: a first Bayesian parameter update utilizing a physics-based synthetic catalog followed by a second Bayesian update utilizing the historical catalog. We developed and presented the algorithm for each Bayesian update by (1) constructing formulations to evaluate both the prior distribution and the likelihood of observing earthquake data for the time- and space-interaction model of main-shock earthquake occurrence developed by Ceferino et al. (2018b), and (2) utilizing MH MCMC to find the posterior distribution of the parameters.

We applied the proposed parameter estimation technique to study the occurrence of earthquakes with magnitude larger than 7.5 in the subduction zone off-shore of Lima, Peru. We utilized a historical catalog for large-magnitude earthquakes containing 450 years of data. Additionally, we built multiple synthetic catalogs containing 2,000 years of data for the region utilizing physics-based simulation and high-performance computing. The simulations were based on a quasi-dynamic approximation and on the rate-and-state friction law. We performed a sensitivity analysis on the physics-based model by varying key input data and the properties of asperity regions based on past earthquakes. We also validated the synthetic catalogs by verifying that they are able to reproduce feasible earthquake cycles and annual magnitude exceedance rates compared to existing historical catalogs. We verified that the synthetic catalogs have high ISC in the same regions as the ones inferred from recent GPS measurements.

After building the synthetic catalogs, we utilized one of them in combination with the historical catalog to perform the two-step Bayesian update. A comparison with a one-step Bayesian update using only the historical catalog showed that the two-step update is a more effective technique for parameter estimation because it leverages the synthetic catalog. Our results show that, unlike the one-step update, the two-step update strategically uses the 2000-year synthetic catalog to perform a first approximation to the parameter medians and to reduce the initial high uncertainty in the parameters. After this update, it further calibrates the estimates by utilizing the 450-year historical catalog. As a result, the two-step update gives estimates with better calibration and lower uncertainties than only using the historical data. A comparison with a one-step update using only historical data shows that the incorporation of the synthetic data reduced the average uncertainty in the parameters by an additional 45%.

Finally, we propagated the effectiveness of the two-step update to the annual magnitude exceedance rates and the time-dependent seismic hazard. We showed that the annual exceedance rates only vary slightly when the one-step update is compared to the two-step update. The return period for earthquakes with magnitudes larger than 8.0 in the fault is close to 80 years for both approaches. In terms of the time-dependent hazard, we calculated the probability of exceeding a PGA of 0.4 g during the next 30
years using the medians of the parameters from both the two-step update and the one step-update with
only historical data. The incorporation of synthetic data in the parameter estimation reduced the hazard
estimates by more than 40% in most of the region of analysis; however, the variation can be different in
other tectonic plates according to the specific seismicity properties of the synthetic catalog.
DATA AND RESOURCES

The QDYN software is open-source and is available at https://github.com/ydluo/qdyn.
# 1 Acknowledgments

We thank Stanford University and the Stanford Research Computing Center for providing computational resources and financial support. We acknowledge Pablo Heresi from Stanford University for insightful discussions on the uncertainty of historical catalogs. We acknowledge the support by the Shah Family Fellowship and the Jhon A. Blume Fellowship from the Civil Engineering Department at Stanford University. We acknowledge the support by the French government through the UCAJEDI Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01.
## 2 References

### References

Akinci, A., Galadini, F., Pantosti, D., Petersen, M., Malagnini, L., and Perkins, D. (2009). Effect of time dependence on probabilistic seismic-hazard maps and deaggregation for the central Apennines, Italy. *Bulletin of the Seismological Society of America*, 99(2 A):585-610.

Allmann, B. P. and Shearer, P. M. (2009). Global variations of stress drop for moderate to large earthquakes. *Journal of Geophysical Research: Solid Earth*, 114(1):1-22.

Ampuero, J.-P. and Rubin, A. M. (2008). Earthquake nucleation on rate and state faults. *Aging and slip laws. Journal of Geophysical Research*, 113(B01302):1-21.

Bradley, A. M. (2014). Software for Efficient Static Dislocation-Traction Calculations in Fault Simulators. *Seismological Research Letters*, 85(6):1358-1365.

Ceferino, L., Kiremidjian, A., and Deierlein, G. (2017). Space and time interaction modeling of earthquake rupture occurrence. In *12th International Conference on Structural Safety and Reliability*, pages 694-703, Vienna, Austria.

Ceferino, L., Kiremidjian, A., and Deierlein, G. (2018a). Parameter Estimation Methods for Modeling of Time and Space Interactions of Earthquake Rupture. In *16th European Conference on Earthquake Engineering*, Thessaloniki.

Ceferino, L., Kiremidjian, A., and Deierlein, G. (2018b). Probabilistic space- and time-interaction modeling of mainshock earthquake rupture occurrence. *Bulletin of the Seismological Society of America*, In review.

Chib, S. and Greenberg, E. (1995). Understanding the Metropolis-Hastings Algorithm. *The American Statistician*, 49(4):327-335.

Chlieh, M., Perfettini, H., Tavera, H., Avouac, J. P., Remy, D., Nocquet, J. M., Rolandone, F., Bondoux, F., Gabalda, G., and Bonvalot, S. (2011). Interseismic coupling and seismic potential along the Central Andes subduction zone. *Journal of Geophysical Research: Solid Earth*, 116(12):1-21.

Day, S., Dalguer, L., Lapusta, N., and Liu, Y. (2005). Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture. *Journal of Geophysical Research*, 110(B12307).

Dieterich, J. H. (1979). Modeling of rock friction: 1. Experimental results and constitutive equations. *Journal of geophysical research*, 84(9):2161-2168.

Dorbath, L., Cisternas, A., and Dorbath, C. (1990). Assessment of the size of large and great historical earthquakes in Peru. *Bulletin of the Seismological Society of America*, 80(3):551-576.

Field, E. H. (2015). Computing Elastic-Rebound-Motivated Earthquake Probabilities in Unsegmented Fault Models: A New Methodology Supported by Physics-Based Simulators. *Bulletin of the Seismological Society of America*, 105(2A):544-559.

Field, E. H., Biasi, G. P., Bird, P., Dawson, T. E., Felzer, K. R., Jackson, D. D., Johnson, K. M., Jordan, T. H., Madden, C., Michael, A. J., Milner, K. R., Page, M. T., Parsons, T., Powers, P. M., Shaw, B. E., Thatcher, W. R., Weldon, R. J., and Zeng, Y. (2015). Long-Term Time-Dependent Probabilities for the Third Uniform California Earthquake Rupture Forecast (UCERF3). *Bulletin of the Seismological Society of America*, 105(2A):511-543.

Hagiwara, Y. (1974). Probability of earthquake occurrence as obtained from a Weibull distribution analysis of crustal strain. *Tectonophysics*, 23(3):313-318.

Hong, H. P. and Goda, K. (2006). A comparison of seismic-hazard and risk deaggregation. *Bulletin of the Seismological Society of America*, 96(6):2021-2039.

Jin, R., Wang, S., Yan, F., and Zhu, J. (2015). Generating Spatial Correlated Binary Data Through a Copulas Method. 3(4):206-212.
Kanamori, H. (1977). The energy release in great earthquakes. *Journal of Geophysical Research*, 82(20):2981.

Kelleher, J. a. (1972). Rupture zones of large South American earthquakes and some predictions. *Journal of Geophysical Research*, 77(11):2087.

Kendrick, E., Bevis, M., Smalley, R., Brooks, B., Vargas, R. B., Lauría, E., and Fortes, L. P. S. (2003). The Nazca-South America Euler vector and its rate of change. *Journal of South American Earth Sciences*, 16(2):125-131.

Langer, C. and Spence, W. (1995). The 1974 Peru Earthquake Series. *Bulletin of the Seismological Society of America*, 85(3):665-687.

Lapusta, N. and Liu, Y. (2009). Three-dimensional boundary integral modeling of spontaneous earthquake sequences and aseismic slip. *Journal of Geophysical Research: Solid Earth*, 114(9):1-25.

Liu, J. S. (2004). *Monte Carlo Strategies in Scientific Computing*. Springer Series in Statistics. Springer New York, New York, NY.

Liu, Y. and Rice, J. R. (2007). Spontaneous and triggered aseismic deformation transients in a subduction fault model. *Journal of Geophysical Research: Solid Earth*, 112(9):1-23.

Luo, Y., Ampuero, J. P., Galvez, P., van den Ende, M., and Idini, B. (2017). QDYN: a Quasi-DYNamic earthquake simulator (v1.1).

Madariaga, R. (1979). On the Relation Between Seismic Moment and Stress Drop in the Presence of Stress and Strenght Heterogeneity. *Journal of Geophysical Research*, 84(B5).

Marone, C. and Kilgoret, B. (1993). (1993) Scaling of the critical slip distance for seismic faulting with shear strain in fault zones [Marone & Kilgore 1993 Nature].pdf. 362(April):618-621.

Matthews, M. V., Ellsworth, W. L., and Reasenberg, P. a. (2002). A Brownian model for recurrent earthquakes. *Bulletin of the Seismological Society of America*, 92(6):2233-2250.

Nishenko, S. and Buland, R. (1987). A generic recurrence interval distribution for earthquake forecasting. *Bulletin of the Seismological Society of America*, 77(4):1382-1399.

Parsons, T., Console, R., Falcone, G., Murru, M., and Yamashina, K. (2012). Comparison of characteristic and Gutenberg-Richter models for time-dependent M 7.9 earthquake probability in the Nankai-Tokai subduction zone, Japan. *Geophysical Journal International*, 190(3):1673-1688.

Robert, C. P. (2014). The Metropolis-Hastings algorithm. *Wiley StatsRef: Statistics Reference Online*, pages 1-15.

Rubin, A. M. and Ampuero, J. P. (2005). Earthquake nucleation on (aging) rate and state faults. *Journal of Geophysical Research: Solid Earth*, 110(11):1-24.

Seiner, L. (2011). *Historia de los sismos en el Perú. Catálogo: Siglos XVIII-XIX*. Universidad de Lima. Fondo Editorial, Lima, Peru.

Silgado, E. (1978). Historia de los Sismos más Notables Ocurridos en el Perú (1513-1974). In Instituto de Geología y Minería, editor, *Serie C. Geodinámica e Ingeniería Geológica. Boletín No. 3*, Lima, Peru.

Strasser, F. O., Arango, M., and Bommer, J. J. (2010). Scaling of the Source Dimensions of Interface and Intraslab Subduction-zone Earthquakes with Moment Magnitude. *Seismological Research Letters*, 81(6):951-954.

Sykes, L. R. and Menke, W. (2006). Repeat times of large earthquakes: Implications for earthquake mechanics and long-term prediction. *Bulletin of the Seismological Society of America*, 96(5):1569-1596.

Tullis, T. E. (1988). *Rock friction constitutive behavior from laboratory experiments and its implications for an earthquake prediction field monitoring program*, volume 126.

Udias, A. and Rice, J. (1975). Statistical analysis of microearthquake activity near San Andreas geophysical observatory, Hollister, California. *Bulletin of the Seismological Society of America*, 65(4):809-827.
Villegas-Lanza, J. C., Chlieh, M., Cavalié, O., Tavera, H., Baby, P., Chire-Chira, J., and Nocquet, J.-M. (2016). Active tectonics of Peru: Heterogeneous interseismic coupling along the Nazca megathrust, rigid motion of the Peruvian Sliver, and Subandean shortening accomodation. *Journal of Geophysical Research : Solid Earth*, pages 1–24.

Woessner, J. and Wiemer, S. (2005). Assessing the quality of earthquake catalogues: Estimating the magnitude of completeness and its uncertainty. *Bulletin of the Seismological Society of America*, 95(2):684–698.