# Bulletin of the Seismological Society of America

## Abstract

Inverting geophysical data has provided fundamental information about the behavior of earthquake rupture. However, inferring kinematic source model parameters for finite‐fault ruptures is an intrinsically underdetermined problem (the problem of nonuniqueness), because we are restricted to finite noisy observations. Although many studies use least‐squares techniques to make the finite‐fault problem tractable, these methods generally lack the ability to apply non‐Gaussian error analysis and the imposition of nonlinear constraints. However, the Bayesian approach can be employed to find a Gaussian or non‐Gaussian distribution of all probable model parameters, while utilizing nonlinear constraints. We present case studies to quantify the resolving power and associated uncertainties using only teleseismic body waves in a Bayesian framework to infer the slip history for a synthetic case and two earthquakes: the 2011 *M*_{w} 7.1 Van, east Turkey, earthquake and the 2010 *M*_{w} 7.2 El Mayor–Cucapah, Baja California, earthquake. In implementing the Bayesian method, we further present two distinct solutions to investigate the uncertainties by performing the inversion with and without velocity structure perturbations. We find that the posterior ensemble becomes broader when including velocity structure variability and introduces a spatial smearing of slip. Using the Bayesian framework solely on teleseismic body waves, we find rake is poorly constrained by the observations and rise time is poorly resolved when slip amplitude is low.

*Electronic Supplement:*Figures of histograms of slip and rise time, waveform comparisons between data and synthetics, and slip velocity along the fault plane for a synthetic case, as well as for the 2011 *M*_{w} 7.1 Van, east Turkey, earthquake, and the 2010 *M*_{w} 7.2 El Mayor–Cucapah, Baja California, earthquake.

## Introduction

Understanding earthquake rupture behavior is of fundamental importance for earthquake‐hazard analysis. Inferences on the evolution of slip on a fault are vital for scenario event modeling and the prediction of ground‐motion parameters through strong ground motion simulations. Current understanding of the behavior of earthquakes is primarily informed by inversions of geophysical data. These deductions are made possible through the utilization of observations of seismic ground motions, which have been commonly used in inverse problems for fault rupture history (Trifunac, 1974; Hartzell and Heaton, 1983; Olson and Anderson, 1988; Das and Kostrov, 1990, 1994; Hartzell *et al.*, 1991, 1996, 2007; Zeng and Anderson, 1996; Sekiguchi *et al.*, 2000; Graves and Wald, 2001; Olson and Apsel, 2001; Wald and Graves, 2001; Ji *et al.*, 2002; Liu and Archuleta, 2004; Custódio *et al.*, 2005; Liu, 2006; Piatanesi *et al.*, 2007; Liu *et al.*, 2008). However, underdetermined model parameters lead to a nonunique ill‐posed problem.

To make this problem tractable, regularization techniques have been employed to transform an ill‐conditioned inverse problem into a well‐conditioned one. These techniques include smoothing, moment minimization, and positivity constraints on the distribution of slip. Mendoza and Hartzell (2013) applied the L‐curve methodology from Tikhonov regularization to obtain an optimal smoothing weight for the finite‐fault teleseismic inversion problem. In this method, a trial‐and‐error search defines the trade‐off curve between fitting the data and satisfying the constraint and thereby defines the optimal constraint weight. Alternatively, Asano *et al.* (2005) employed Akaike’s Bayesian information criterion (Akaike, 1980) applied to the inversion of strong ground motion waveforms to select optimum constraint weights.

Although these types of regularization techniques seek a well‐posed inverse problem, application of the regularization constraints alters the solution space. How this solution space is altered can be fundamentally different given different regularization techniques, yielding contrasting results. However, Bayesian inference can be used to avoid these problems. Inversion schemes that employ regularization techniques most often conceive a single solution based on a single set of model parameters. Bayesian inference yields a set of potential models, using physical *a priori* constraints throughout the process, to form posterior probability density functions (PDFs) of the model parameters. The ability to infer a set of potential models through forward model evaluation eliminates the need for regularization techniques, while illuminating the uncertainties in the solution.

Though Bayesian theory has a long history (Bayes, 1763; Laplace, 1812; Jeffreys, 1931, 1939; Tarantola, 2005), its application to finite‐fault problems has historically been limited by computational power. However, increasing computational power has made tractable Bayesian analysis of larger scale geophysical problems. This has led to a wide range of Bayesian geophysical inverse problems, for example, from utilizing lunar seismic and gravity data (Khan *et al.*, 2007), to reservoir monitoring measurements (Xuan and Sava, 2010), to frequency‐domain electromagnetic data (Minsley, 2011). In earthquake seismology, specifically with reference to finite‐fault inversions, a multiplicity of different data sets has been incorporated into a Bayesian framework. With the inclusion of strong ground motion data (Monelli and Mai, 2008; Fan *et al.*, 2014); Global Positioning System (GPS) offsets, offshore seafloor geodesy, and tsunami records (Minson *et al.*, 2014); *W*‐phase waveforms (Dettmer *et al.*, 2014); and geodetic, tsunami, and strong ground motion records (Duputel *et al.*, 2015), Bayesian inference has been widely used to access the uncertainty question and yield solutions unaffected by regularization. To quantify such uncertainties in the Bayesian framework, Dettmer *et al.* (2014) presented 95% credibility intervals for slip magnitude and rupture velocity, as well as the marginal densities for the hypocenter location, rupture velocity, and moment magnitude. Duputel *et al.* (2015) quantified the uncertainties by presenting the posterior mean slip model with 95% confidence error ellipses, the PDFs for stress drop, rupture velocity, rise time, and seismic potency, and additionally highlighted the observational errors and prediction uncertainty using sensitivity plots for each of the data sets used. Despite the growing body of literature employing Bayesian methods for finite‐source inversions, we know of no work that has modeled only teleseismic body‐wave data within a Bayesian framework; thus, we seek to better understand their resolution and uncertainty.

In this study, we focus on uncertainty analysis from Bayesian inference using teleseismic body waves to infer the slip history for a synthetic case, as well as for the 2011 *M*_{w} 7.1 Van, east Turkey, earthquake and the 2010 *M*_{w} 7.2 El Mayor–Cucapah, Baja California, earthquake. Our objective in this study is to quantify the resolving power of teleseismic body waves in determining finite‐fault slip.

## Methodology

### Data

The teleseismic body‐wave data sets were obtained from the Incorporated Research Institutions for Seismology Data Management Center. The *P* waves used in the study are from the vertical component, whereas the *S* waves are from the transverse component. The records are first processed by removing the instrument response, over the frequency range of 0.0166 Hz (60 s) to 2 Hz (0.5 s). The instrument‐corrected data are then resampled to a time step of 0.1 s, and a band‐pass Butterworth filter from 0.0166 Hz (60 s) to 1 Hz (1 s) is applied. This filter is applied to aid in the evaluation of signal and noise present in the records and to help define an appropriate frequency band to use for the inversion. The data records are then transformed into the frequency domain, in which the Bayesian inversion is conducted. The frequency band used for the Bayesian inversion depends on the magnitude of the earthquake, with longer periods required for larger magnitude events.

The record length used in the inversion is chosen based on a consideration of the characteristics of the waveforms, the fault dimensions, and the assumed rupture velocity, to encompass the entire source process. Generally, the record length should not exceed the length of the synthetic seismogram for a subfault furthest from the hypocenter, accounting for the time it takes the wave front to propagate to that distance (Hartzell, 1989).

### Model Parameterization

Fault planes for our finite‐fault inversions are parameterized using fixed orientation (strike and dip), dimensions (along strike and down dip), and depth, embedded in a seismic velocity model. For a more complex rupture, like the El Mayor–Cucapah earthquake, we can consider multiple fault planes. The determination of the fault‐plane parameters is based on moment tensor solutions (e.g., U.S. Geological Survey, Global Centroid Moment Tensor [CMT], Southern California Seismic Network, and other published studies). The fault plane is discretized into a specified number of equal‐sized subfaults or patches. The size and total number of subfaults are chosen based on several considerations. With a linear formulation (Hartzell and Heaton, 1983), the fault plane is generally overparameterized using a relatively large number of subfaults and then applying smoothing constraints to adjacent slip. In this study, we wish to avoid smoothing constraints, and our choice of subfault size is driven by the time resolution between nearby subfaults and the frequency content of the data. There is also a practical limitation, based on computational effort, to the number of model parameters in a Bayesian style inversion. We will discuss this topic further in the Results and Discussions section.

Our finite‐fault spatiotemporal slip history is defined by the following set of parameters. There is a finite rupture duration to quantify the time it takes for a subfault to rupture, known as rise time. The rise time and its functional form define the source time function (STF) for each subfault. The direction of slip on a subfault defines the rake angle. The magnitude of slip for each subfault completes the parameterization of the rupture.

To find the spatiotemporal parameters, we first compute displacement waveform synthetics for two orthogonal unit steps in slip on a subfault known as the Green’s functions (Fig. 1). During the Bayesian process, two weights for the orthogonal Green’s functions and one for the STF duration are chosen to evaluate the goodness of fit between the observed and modeled data. In this study, we use a box‐car STF, in which the width of the box‐car is a parameter in the inversion. The two weights that scale the orthogonal Green’s functions define the slip amplitude and rake angle. In this study, we fix the rupture velocity to be a constant. This simplification is made because variations in rupture velocity can have a significant impact on the slip model (Lay *et al.*, 2010) and can lead to a model parameterization that greatly increases the computational effort. We choose to fix rupture velocity to highlight the uncertainties associated with slip and rise time.

The two orthogonal Green’s functions are computed such that they are ±45° from the inferred average rake. This rake estimate is generally from moment tensor analysis. The Green’s functions are calculated using the ak135 velocity structure (Kennett *et al.*, 1995).

### Bayesian Inference

#### Bayes Theorem

The Bayesian approach is used to estimate the posterior model distribution that defines the solution space. However, to reach the posterior model distribution we must define the prior PDF and data likelihood function, which are defined over both the model and data spaces. The prior and theoretical PDFs are then used to develop the posterior model distribution to infer the solution.

In a Bayesian inference framework, there are two main manifolds (or dimensions) that represent the model space and the data space . Within each respective space lie individual models **m**={*m*^{1},*m*^{2},…,*m*^{n}} and data **d**={*d*^{1},*d*^{2},…,*d*^{n}}. Given these manifolds, a new state of information about a model **m** can be achieved given data **d** using Bayes theorem: (1)(Tarantola, 2005), in which *P*(**d**|**m**) is the conditional PDF of the data **d** given a model **m**, is the model marginal PDF over the model space , is the data marginal over the space , and *P*(**m**|**d**) is the conditional PDF of the model **m** given the data **d**. Bayes theorem is the generic solution to probabilistic inverse problems and can be broken down into three fundamental parts: prior, theoretical, and posterior joint PDFs.

#### Prior

The prior joint PDF encapsulates both the *a priori* information on the models and the observed data, defined as (2)in which is the marginal prior probability of the data, is the marginal prior probability of the model, and *ρ*(**d**,**m**) is the resultant prior joint PDF defined over the manifold .

The *a priori* information on the models, or the marginal prior probability of the model , captures the prior knowledge of the model independent of all observed data. In this study, the *a priori* information on the models is represented by a truncated uniform distribution, meaning that all values are equally probable, which only reflects the range of possible values in each model **m**: (3)in which *m*_{maxi} and *m*_{mini} are the bounds chosen for the model parameter **m**, and *N* is the total number of parameters. In our parameterization, *N* equals three times the total number of subfaults in the finite‐fault problem, accounting for the two slip weights and rise‐time values for each subfault. The limits for both the slip and rise time are chosen based on previously published studies and the magnitude of the event. Table 1 highlights the bounds chosen for each of the cases presented in the Results and Discussions section.

With each observed value in the data there are associated uncertainties. These uncertainties in the data can be represented by a PDF, the prior marginal of the data space, . In this study, the state of information on the observed data is considered a Gaussian PDF given observed data, **d**_{obs}: (4)in which |·| represents the determinate, and *C*_{d} is the data covariance matrix.

#### Theoretical

The theoretical PDF captures the physical theory that relates the model with the data, defined as (5)in which *θ*(**d**|**m**) is the conditional PDF of the data given the model, is the marginal homogeneous PDF of the model, and Θ(**d**,**m**) is the resultant joint theoretical PDF.

In this study, the theoretical conditional PDF *θ*(**d**|**m**) links the slip and rise time (model parameters) to the teleseismic body waves (data). It is assumed to be a Gaussian PDF, (6)in which *G*(·) represents the forward operator, and *C*_{T} is the theory covariance matrix. The forward operator is applied to given models **m** to construct plausible data, or the forward problem.

#### The Posterior and Likelihood Function

The prior and theoretical states of information combine, through a conjunction, to form the *a posteriori* state of information, (7)in which *ρ*(**d**,**m**) is the prior joint PDF, Θ(**d**,**m**) is the theoretical PDF, *ν* is a normalizing constant (8)defined over the joint manifold , *μ*(**d**,**m**) is the homogeneous joint PDF (9)and *σ*(**d**,**m**) is the posterior joint PDF. To arrive at the *a posteriori* state of information about the model, or the solution, a posterior marginal PDF of the model can be constructed, (10)in which is the resultant *a posteriori* state of information on the model represented by a PDF, given the prior and theoretical states of information. This refined state of information on the model can be written as (11)in which *L*(**m**) is the likelihood function.

The likelihood function describes how well a model **m** explains data **d**. Using the prior marginal of the data (equation 4) and the theoretical (equation 6), the likelihood function used is constructed as (12)Equation (12) highlights the convolution of two Gaussian probability densities and can be simplified (e.g., Tarantola, 2005) to (13)in which *C*_{D} is the covariance that encapsulates the uncertainties in the data *C*_{d} and the uncertainties in the theory *C*_{T}, given as *C*_{D}=*C*_{T}+*C*_{d}.

#### Covariance Matrix

The covariance of a multidimensional Gaussian, synonymous to variance in 1D, describes both the uncertainties (diagonal elements) and correlations (off‐diagonal elements), (14)in which *f*(**x**) is the PDF on which to define the covariance, is the mean of that PDF, and is the resulting covariance at . The covariance needed for Bayesian inference, shown in the likelihood function (equation 13), includes the covariance matrix *C*_{D}, which combines the associated uncertainties and correlations in the theory and data. In this study, we adopt a diagonal form for the covariance matrix. This simplification invokes the uncertainties, saving the computational expense of updating a covariance matrix.

The data covariance *C*_{d} describes the uncertainties associated with the observations or the teleseismic body waves. To capture the uncertainties in the seismograms, the pre‐event noise was used to determine the fractional difference between the peak amplitude of the pre‐event noise and the peak amplitude of the *P* waves, (15)in which represents the uncertainty of the data for a specific station, *M*_{E} represents the maximum amplitude of the *P* wave present in the waveform, and *M*_{P} represents the maximum amplitude of the pre‐event noise.

The covariance *C*_{T} describes the uncertainty in the theory, which can include incorrect source geometry or hypocenter location and poorly chosen Earth structure. In this study, we chose to highlight the effects of the Earth structure, specifically the velocity model, on the *a posteriori* state of information on the model. To assess this uncertainty, 10 random crustal model perturbations were drawn from a Gaussian distribution with a 5% standard deviation in wave velocity and layer thickness from the reference *P*‐ and *S*‐wave velocity model (Fig. 2; Razafindrakoto and Mai, 2014). For each perturbation we maintain the ratio 1.60≤*V*_{P}/*V*_{S}≤1.90. Table 2 gives the reference velocity model values.

### CATMIP: Algorithm for Bayesian Inference

Our Bayesian inference employs an implementation of a parallel Markov chain Monte Carlo (MCMC) algorithm: Cascading Adaptive Transitional Metropolis in Parallel (CATMIP; Minson *et al.*, 2013). The Metropolis algorithm (Metropolis *et al.*, 1953; Hastings, 1970) uses a Markov chain (random walk) to draw samples of a target PDF. A Markov chain is a stochastic sequence that transitions between different states, in which the probability of each state is dependent on the previous state. The utilization of a single Markov chain consequently is a nonparallelized algorithm. However, CATMIP can employ thousands of Markov chains simultaneously, allowing it to be developed in a parallel framework and permitting solutions of higher dimensional modeling problems that once were computationally intractable. To ensure higher efficiency for higher dimensional spaces, CATMIP also uses transitioning as well as resampling (Ching and Chen, 2007).

CATMIP is a generic Bayesian MCMC sampler that employees transitioning, which shares several traits with simulated annealing optimization. Algorithms that employ annealing start off from an initial hot state and cool to a final cold state. In this particular implementation, the initial state is the prior PDF that cools to the final posterior PDF. To ensure that a proper cold state has been reached, CATMIP samples from intermediate PDFs during cooling steps, or transitional stages, which are controlled by an annealing parameter. During transitioning, resampling is conducted in which less probable models from the previous transitioning step are replaced with more probable models. This allows CATMIP to relocate models that are trapped in areas of low probability to become seeds for new Markov chains in regions of higher probability. For a complete discussion of CATMIP, see Minson *et al.* (2013)

To ensure that the posterior PDF is stable, the total number of samples is regulated by the total number of model parameters, in which an increase in dimensionality requires an increase in samples. In this study, the values chosen for the number of Markov chains and the length of each chain are chosen purely on the basis of repetition of synthetic cases. The values are chosen such that the solution stabilized and did not change significantly with further sampling of the model space.

### Uncertainty Analysis

To investigate the effect of variability in seismic velocity models, the Bayesian inversion is run 11 times, for a reference model and 10 velocity model perturbations. For each of the velocity models, the inversion is run using 1000 Markov chains with 200,000 steps per Markov chain, in which 1000 Markov chains with 200,000 steps was found to stabilize the solution after repetition of different combinations of Markov chains and steps. Once all 11 velocity models have been run, the resultant posterior PDFs are combined, giving a total of 11,000 estimates for each of the model parameters: weights for the two orthogonal slip components and rise time for each subfault. Then, the most probable solution is given by the mean of this stacked posterior PDF.

To highlight the single velocity model, or the reference model, the inversion is also run using 11,000 Markov chains with 200,000 steps per Markov chain, equating to an equivalent number of Markov chains as for the multiple velocity models.

After performing multiple inversions for all three fault configurations (synthetic case, 2011 Van earthquake, and the 2010 El Mayor–Cucapah earthquake), it was found that the rake for each subfault was an unstable parameter, varying unrealistically between adjacent subfaults. To overcome this difficulty, we constructed the rake such that each subfault could vary ±10° from the average inferred rake (Table 3). The variance of ±10° allows for a significant variation in rake that is typically not exceeded in moment tensor estimates. As an alternative, the inversion could be conducted using a single Green’s function; however, there is a clear bias in this parameterization, because the inferred rake is assumed to be correct with no uncertainty. There is no constraint placed on the moment or any other constraint on the two slip components or the rise time.

## Results and Discussions

Here, we present the Bayesian inversion results for three different fault configurations, including a synthetic test, the 2011 *M*_{w} 7.1 Van, east Turkey, earthquake, and the 2010 *M*_{w} 7.2 El Mayor–Cucapah, Baja California, earthquake. For each of these cases, we present two distinct solutions for comparison, one utilizing a single given velocity model and another to show the effect of uncertainty in velocity model. The solutions presented here are obtained by calculating the mean value from the distribution of slip or rise time for each subfault. Rise‐time values in the plots corresponding to values of slip 20 cm and less are shown in gray (Figs. 4, 8, and 12).

### Synthetic Tests

We first consider a synthetic test case, using teleseismic *P* waves, for a thrust fault with 36 idealized stations, all at a distance of 60°, uniformly spaced in azimuth every 10°. The synthetic model consists of a single 120 km (along strike) × 16 km (along dip) fault plane with 4 km×4 km subfaults (120 total subfaults). We construct the synthetic model to have four major asperities, containing a total of 16 subfaults with slip of 500 cm per subfault (Fig. 3a), a rake of 60°, a rise time of 3 s (Fig. 4a), and a constant rupture velocity of 2.5 km/s. The slip history for each subfault is constructed from two orthogonal components of slip, at 105° and 15° (±45° from the true rake). The inversion is conducted using a frequency range from 0.0167 Hz (60 s) to 1 Hz (1 s).

Although our synthetic data are noise free, we construct a data covariance *C*_{d} (Minson *et al.*, 2013) using equation (15) that is representative of the quietest stations for both the 2011 Van, east Turkey, and 2010 El Mayor–Cucapah, Baja California, data sets. The data covariance consists of only diagonal elements with each element consisting of an uncertainty of 0.02 for all observations.

#### Synthetic: Single Velocity Model with No Angle Constraint

Figures 3b and 4b present the posterior mean solution given a single velocity model for slip and rise time, respectively. The solution is obtained using no angle constraint, meaning that the weights found in the inversion for the two orthogonal Green’s function have open range, and that the resulting rake can vary from 15° to 105°. In this particular scenario, the change in rake between adjacent subfaults is reasonable, with no sudden jumps in rake. However, given real data, the behavior of the rake is unreasonable, as we show below.

#### Synthetic: Single Velocity Model with No Angle Constraint and Addition of SH Waves

The solution obtained in Figures 3b and 4b only considers *P* waves. Figures 3c and 4c present the posterior mean solution when an additional 36 idealized teleseismic *SH* waves are added. The additional stations reside at the same locations as the 36 stations for *P* waves. The solution is obtained using a single velocity model with no angle constraint. In this particular implementation, the addition of the *SH* waves does not help resolve the rake; because the behavior of the rake is more erratic, evidenced by the approximately 90° jump in rake near the hypocenter. To be consistent throughout the remaining synthetic scenarios, we only implement teleseismic *P* waves with a ±10° variance allowed from the 60° rake.

#### Synthetic: Single Velocity Model

With the implementation of a single velocity model with a ±10° angle constraint on the inferred rake, the slip distribution (Fig. 3d) and rise‐time distribution (Fig. 4d) highlight characteristics similar to the true model (Fig. 3a). The distribution of slip is well resolved except near the hypocenter, where there is a smoothing effect. The poor resolution near the hypocenter is due to the small arrival‐time difference at teleseismic distances between subfaults at similar but opposite directions from the hypocenter. However, as the rupture propagates further, the difference in arrival times increases, allowing the slips to be differentiated from one another. The rise times, as shown in Figure 4d, are also not as well resolved near the hypocenter and tend to be more consistent with the true model further from the hypocenter. In addition, subfaults with longer rise times and smaller slips make negligible contribution to the total synthetic. The overall rake obtained is consistent with the rake of the synthetic model, given the angle constraint. The mean of the posterior model PDF also does a good job recovering the moment of 3.98×10^{19} N·m (*M*_{w} 7.0) compared with the correct moment of 3.75×10^{19} N·m (*M*_{w} 7.0).

#### Synthetic: Multiple Velocity Models

The implementation of multiple velocity models is presented in Figures 3e and 4e for slip and rise time, respectively. The results highlight similar characteristics with (Figs. 3d and 4d) and without (Figs. 3e and 4e) velocity model variability. Here, we see a similar lack of resolution near the hypocenter as a consequence of the small arrival‐time difference between certain subfaults. However, in this scenario with velocity model variability, there is somewhat greater smoothing of slip on different regions of the fault. With velocity model variations, changes in timing of phases cause subfaults to acquire slip that are not near the active subfaults of the true model (Fig. 3a). With velocity model variability, the moment increases slightly to 4.11×10^{19} N·m (*M*_{w} 7.0).

#### Synthetic: Uncertainties

To investigate the uncertainties, Figure 5a presents scatter plots of rise time versus slip for the four numbered subfaults (subfaults 6, 51, 66, and 91), in Figure 3d and 3e, for a single velocity model (shown in blue) and multiple velocity models (shown in red), respectively. Each scatter plot represented in Figure 5a shows the posterior ensemble of 11,000 possible solutions for each subfault. The posterior ensemble is no longer confined to a single region but contains multiple areas. These represent the different velocity models used in the inversion, because each velocity model creates its own distribution of solutions. In general, the scatter plots for the single and multiple velocity models show an independence between slip and rise time; however, there is a slight linear trend on some subfaults, such as in subfault number 66, which suggest that slip and rise time are codependent. The rise times are generally better constrained than slip, except for low‐slip areas of the fault plane. This conclusion is based on the composite of the multiple velocity model solutions, in which changes in the velocity model have a much larger effect on slip than on rise time.

To highlight the uncertainties between the single and multiple velocity models, a comparison of the model slip histograms for the four numbered subfaults in Figure 3d and 3e is conducted (see Ⓔ the electronic supplement to this article). The distribution of slip for the single velocity model shows highly peaked and well‐defined distributions relative to the PDFs with velocity model variability. With multiple velocity models, the distribution of slip for each subfault is greater with accompanying larger standard deviations, on the order of 3.3 times greater on average.

Table 4 highlights the values for the mean slip and standard deviations for the subfaults compared. Taking subfault number 6 as an example, the mean slip value corresponding to a single velocity model is 468 cm, with an accompanying standard deviation of 49 cm, compared with the multiple velocity model solution of 403 cm, with an accompanying standard deviation of 250 cm. In general, the uncertainties increase with multiple velocity models, altering the solution space by having a wider range of slip values in the posterior (Razafindrakoto and Mai, 2014).

Given velocity model variability, certain distributions are best represented by a multimodal distribution, not a unimodal one, as highlighted by the scatter plots. This phenomenon can be explained by two effects: the timing difference caused by using different velocities and the use of a finite number of velocity models. Each of the 11 velocity models contains different *P*‐ and *S*‐wave velocities, which causes timing differences between solutions. This timing difference causes slip to reside on different subfaults.

The phenomenon of the multimodal distribution is further explained by considering an additional 10 velocity models (see Ⓔ the electronic supplement to this article). When more velocity models are used, the histograms start to fill in. As an example, consider the histograms for subfault number 6 for 11 and 21 velocity models. There is also a resulting small change in the mean and standard deviation. From this point forward, we opt to consider velocity model variability using a total of 11 models, due to computational expenses.

A comparison of the model rise‐time histograms for the solution containing a single velocity model (Fig. 4d) and multiple velocity models (Fig. 4e) can be seen in Ⓔ the electronic supplement to this article. The standard deviations corresponding to the rise times show similar characteristics to the slip standard deviations, with an increase in standard deviation when considering multiple velocity models, increasing on average by a factor of 1.3. Table 4 shows the rise‐time mean values and associated standard deviations for the four subfaults. With the inclusion of multiple velocity models, the rise‐time distributions are also best represented by a multimodal distribution.

#### Synthetic: Waveform Fit

The calculated waveforms from the true model (Figs. 3a and 4a), the inversion results for the mean posterior model from the reference velocity model inversion (Figs. 3d and 4d), and the inversion including velocity model variability (Figs. 3e and 4e) are shown in Ⓔ the electronic supplement to this article. The waveforms with and without velocity model variability show no visible deviations from the true data (see Table 5 for L2 norms). This result points out the limitations in resolution of teleseismic body waves given the differences in these fault models, as well as illustrating how the inverse problem is underdetermined.

### 2011 Van Earthquake

Now, we turn to the inversion of real earthquake data and the 2011 *M*_{w} 7.1 Van, east Turkey, earthquake with 55 teleseismic *P*‐wave records azimuthally dispersed with distances from 31° to 94°, as shown in Figure 6a. The rupture is modeled on a single plane that is 100 km (along strike) × 40 km (along dip), having 5 km×5 km subfaults (160 total subfaults), with a strike of 246°, a dip of 38°, and a rupture velocity of 2.5 km/s (Mendoza and Hartzell, 2013). The strike and dip values used in the inversion refer to the Global CMT solution (see Data and Resources). Konca (2015) used teleseismic *P* and *SH* waves, GPS displacements, strong‐motion data, and high‐rate GPS station waveforms for the 2011 *M*_{w} 7.1 Van earthquake to resolve the uncertainty in fault geometry by performing a grid search using different combinations of strike and dip. Their analysis showed that the teleseismic and GPS data were unable to definitively constrain the fault geometry, demonstrating that different combinations of strike and dip had little impact on the solutions. We opt to use the strike and dip given by the Global CMT. Given the magnitude of the 2011 Van earthquake, a range of frequencies from 0.0167 Hz (60 s) to 1 Hz (1 s) is used in the inversion. The *a priori* values (Table 1) chosen for the 2011 Van earthquake come from the magnitude of the event as well as published studies (Hayes, 2011; Fielding *et al.*, 2013; Mendoza and Hartzell, 2013; Utkucu, 2013; Konca, 2015).

#### 2011 Van Earthquake: Single Velocity Model with No Angle Constraint

The posterior mean solution, for both the slip distribution (Fig. 7a) and rise‐time distribution (Fig. 8a), is obtained using a single velocity model with no angle constraint. With no constraint, the rake is allowed to vary from 15° to 105°. In this scenario, as seen in Figure 7a, the rake has a sudden change of 90° from adjacent subfaults. This behavior is unrealistic and highlights the lack of rake resolution and the requirement for an angle constraint.

#### 2011 Van Earthquake: Single Velocity Model with No Angle Constraint and Addition of *SH* Waves

In an attempt to help constrain the rake, 18 teleseismic *SH* waveforms (Fig. 6a) are added together with the *P* waves. The posterior solution for the slip distribution (Fig. 7b) and rise time (Fig. 8b) is obtained with the use of a single velocity model and no angle constraint. The rake shows similar results to the solution with only *P* waves. Given that the rake still has poor resolution, we conclude that the addition of *SH* waves does not help constrain the finite‐fault problem in this application of the Bayesian inference. From this point forward, we only consider scenarios using teleseismic *P* waves with an angle constraint applied, allowing the rake to vary ±10° from the inferred rake of the earthquake.

#### 2011 Van Earthquake: Single Velocity Model

Utilizing a single velocity model and the angle constraint of ±10°, Figures 7c and 8c show the slip and rise‐time distribution, respectively. Using a single crustal model, the peak slip achieved is 411 cm near the hypocenter and a maximum slip of 320 cm at the surface. Given the dip of the fault of 38° and rake of 60°, the surface slip is equivalent to a vertical displacement of 160 cm and a horizontal displacement of 171 cm. These values are greater than those from field reports (Doğan and Karakaş, 2013), in which the maximum vertical and horizontal offsets were found to be 150 and 90 cm, respectively. However, the surface measurements may not be an accurate representation of slip over the depth range of the subfault, and our estimate of slip is dependent on the velocity model. One or two subfaults at the boundaries of the model show larger slip. These subfaults only contribute to the late portion of the teleseismic records, in which they are tapered for the inversion. Arrivals around this time are probably not source related. The posterior mean solution gives a moment of 5.24×10^{19} N·m (*M*_{w} 7.1) for a shear modulus of 29.3 GPa, based on the velocity model. The average slip velocities, found by dividing the slip by the rise time for each subfault (Fig. 9a), show a general trend of higher slip velocities near the hypocenter that decrease radially away from the hypocenter. This suggests an energetic initiation followed by decreasing radiation of high frequencies as the rupture propagates away from the hypocenter.

#### 2011 Van Earthquake: Multiple Velocity Models

The slip distribution (Fig. 7d) and rise‐time distribution (Fig. 8d) obtained using multiple velocity models show similar characteristics to that of the single velocity model (Figs. 7c and 8c). With uncertainty in velocity model, the peak slip decreases from 411 to 361 cm, and the peak surface slip has decreased from 320 to 180 cm, equivalent to 90 and 96 cm of vertical and horizontal offsets, respectively. The patch of higher slip amplitude near the hypocenter obtained in the single velocity model solution (Fig. 7c) decreases in amplitude. The most prominent consequence of including velocity model uncertainty is the activation of slip on more subfaults, as we saw in the synthetic case. As observed in Figures 7d and 8d, compared with the results of a single velocity model (Figs. 7c and 8c), a smoothing effect causes the peak slip values to lower by distributing the slip among more subfaults. The moment rises slightly from 5.24×10^{19} N·m (*M*_{w} 7.1) to 5.70×10^{19} N·m (*M*_{w} 7.1), using the same shear modulus of 29.3 GPa. The rise times, as shown in Figure 8d, show similar patterns to that of using a single velocity model (Fig. 8c), starting at lower values near the hypocenter and increasing away from the hypocenter. The slip velocities (Fig. 9b) show a similar distribution to a single velocity model (Fig. 9a).

#### 2011 Van Earthquake: Solution Comparison

Figure 10 presents a comparison of four solutions for slip distributions: Figure 10a shows the single velocity model solution, Figure 10b shows the multiple velocity model solution, Figure 10c presents the solution obtained by Mendoza and Hartzell (2013), and Figure 10d presents the solution obtained by Hayes (2011). Mendoza and Hartzell (2013) used teleseismic *P* waves to obtain a slip distribution using the linear least‐squares inversion scheme of Hartzell and Heaton (1983). Hayes (2011) incorporated teleseismic *P* waves, *S* waves, and long‐period surface waves, using the finite‐fault algorithm of Ji *et al.* (2002), which is a simulated annealing regularized inversion in the wavelet domain, to show a main area of slip around the hypocenter with a peak slip of 400 cm. All four solutions show similar characteristics, highlighting a major patch of slip near the hypocenter as well as surface slip. Mendoza and Hartzell (2013) obtained a moment of 5.3×10^{19} N·m, similar to our estimate. The Global CMT estimate of 6.3×10^{19} N·m is somewhat larger, probably due to the incorporation of longer period data.

#### 2011 Van Earthquake: Uncertainties

Figure 5b presents the posterior ensemble of rise time versus slip as scatter plots for the four numbered subfaults (subfaults 35, 57, 86, and 131), in Figure 7c,d, for a single velocity model and multiple velocity models. Figure 5b highlights the same attributes as shown in the synthetic case (Fig. 5a), because each velocity model tends to create its own solution space. The linear relationship between rise time and slip is also depicted in certain subfaults. Subfault number 57, when using multiple velocity models, shows that rise times corresponding to lower slip values are poorly resolved.

The uncertainties associated with the four numbered subfaults for a single velocity model and multiple velocity models are shown as histograms in Ⓔ the electronic supplement to this article. The distribution of slip for both a single velocity model and multiple velocity models shares the same characteristics as discussed in the synthetic case. With the implementation of velocity model uncertainty, the slip distribution becomes broader, with an increase in standard deviation by a factor of 4.7. Each subfault compared has different mean values and accompanying standard deviations, as shown in Table 4.

The associated rise‐time uncertainties of the four numbered subfaults are additionally shown in Ⓔ the electronic supplement to this article for a single velocity model and multiple velocity models. Rise‐time standard deviations increase by an average factor of 2.8 when considering uncertainty in the velocity model. The rise‐time mean values and associated standard deviations for the subfaults compared are given in Table 4.

#### 2011 Van Earthquake: Waveform Fit

Comparison of the 55 teleseismic *P*‐wave records with the waveforms calculated from the slip distribution of the single velocity model (Figs. 7c and 8c) and the waveforms calculated from the slip distribution using velocity model variability (Figs. 7d and 8d) can be seen in Ⓔ the electronic supplement to this article. Both solutions fit the data well, showing a small misfit between the data and synthetics for the two solutions (see Table 5 for L2 norms). In general, the synthetics tend to match each other for the first 10 s and then slightly deviate from one another. The ability of both solutions to represent the data well highlights the limitation in resolution of teleseismic body waves. Although the solution is altered given uncertainty in velocity model, the waveform fits for both solutions are comparable.

### 2010 El Mayor–Cucapah Earthquake

We now turn to the last, and more complicated, scenario of the 2010 *M*_{w} 7.2 El Mayor–Cucapah, Baja California, earthquake. Wei *et al.* (2011) showed that the earthquake initially ruptured as a normal fault lasting ∼15 s with an *M*_{w} 6.3, then evolving into an *M*_{w} 7.2 earthquake on the main strike‐slip fault. Wei *et al.* (2011) used GPS, Interferometric Synthetic Aperture Radar, subpixel correlation, synthetic aperture radar, and teleseismic body waves in the simulated annealing regularized inversion of Ji *et al.* (2002) for both the sub‐ and main event of the 2010 El Mayor–Cucapah earthquake. Initial modeling was unable to achieve good waveform fits when using a single source.

#### 2010 El Mayor–Cucapah Earthquake: Subevent

Because the earthquake rupture initiated as a normal fault of *M*_{w} 6.3, our modeling approach is to initially invert only the first 15 s of 47 teleseismic *P*‐wave records azimuthally dispersed in distance from 31° to 93°, as shown in Figure 6b. Because of the relatively small magnitude of the initial rupture, spectral values of 0.067 Hz (15 s) to 1 Hz (1 s) are chosen for the inversion. The initial rupture is considered to occur on a single plane of 32 km (along strike) × 16 km (along dip), having 4 km×4 km subfaults (32 total subfaults), with a strike of 335°, a dip of 45°, and rupture velocity of 2.0 km/s (Wei *et al.*, 2011). The *a priori* values (Table 1) are chosen based on the magnitude of the event and from Wei *et al.* (2011).

#### Subevent: Single Velocity Model

The posterior mean solution for slip and rise‐time distribution are shown in Figures 11a and 12a, respectively, for the subevent of the 2010 El Mayor–Cucapah, Baja California, earthquake using a ±10° constraint on the inferred rake. Using a single velocity model, the peak slip near the hypocenter is 135 cm, which is close to the peak value found by Wei *et al.* (2011) of 160 cm, also near the hypocenter. Wei *et al.* (2011) found higher amplitude slip near the hypocenter and on the southeast side of the fault. We also find high slip near the hypocenter with a more bilateral distribution. The posterior mean solution obtained a seismic moment of 3.86×10^{18} N·m (*M*_{w} 6.3), consistent with that of the 6.3 moment magnitude found by Wei *et al.* (2011). The rise times (Fig. 12a) show a similar pattern to the Van earthquake results with lower rise times near the hypocenter and generally increasing rise times radially away from the hypocenter. The slip velocity distribution (see Ⓔ the electronic supplement to this article) shows similar characteristics to that of the 2011 Van earthquake, suggesting an energetic initiation.

#### Subevent: Multiple Velocity Models

With the utilization of multiple velocity models, Figures 11b and 12b show the posterior mean solution for slip and rise time, respectively. The results are similar to that of a single velocity model (Figs. 11a and 12a). The seismic moment slightly rises from 3.68×10^{18} N·m (*M*_{w} 6.3) to 4.16×10^{18} N·m (*M*_{w} 6.4). The rise times obtained (Fig. 12b) using velocity variability show a very similar pattern as with a single earth structure. The slip velocity distribution (see Ⓔ the electronic supplement to this article) obtained from multiple velocity models highlights similar patterns to that of a single velocity model.

#### Subevent: Uncertainties

The posterior ensemble of the solution for rise time and slip are plotted against one another (Fig. 5c) for the single velocity model solution (shown in blue) and multiple velocity model solution (shown in red), for the four numbered subfaults (5, 14, 19, and 26) in Figure 11a and 11b. Figure 5c shows similar attributes to that of the synthetic case (Fig. 5a) and the 2011 Van earthquake (Fig. 5b), highlighting multiple solution spaces when using multiple velocity models. Figure 5c additionally shows that there are broader distributions when using multiple velocity models. The subevent of the 2010 El Mayor–Cucapah earthquake further shows the lack of resolution of rise times for small slip, as seen in subfaults 14 and 26.

The model slip histograms (see Ⓔ the electronic supplement to this article) compare the four numbered subfaults for the solution utilizing a single velocity model (Fig. 11a) and variability in velocity model (Fig. 11b). With the implementation of velocity model uncertainty, the slip distributions become broader, with an increase in standard deviations by a factor of 2.2. Table 4 shows the mean values and standard deviations for the subfaults compared.

The model rise‐time histograms (see Ⓔ the electronic supplement to this article) compare the four numbered subfaults for a single velocity model (Fig. 12a) and multiple velocity models (Fig. 12b). Table 4 shows the mean values and standard deviations of rise time associated with the subfaults. Standard deviations increase by a factor of 1.3 when considering velocity model variability.

#### Subevent: Waveform Fit

The comparison of the 47 teleseismic *P*‐wave records with the forward‐modeled waveforms calculated from the slip distribution of the single velocity model (Figs. 11a and 12a) and the forward‐modeled slip distribution using velocity model variability (Figs. 11b and 12b) are shown in Ⓔ the electronic supplement to this article. In general, waveforms for both solutions are very similar (see Table 5 for L2 norms). Both slip distributions yield comparable waveform fits to that of Wei *et al.* (2011).

#### 2010 El Mayor–Cucapah Earthquake: Main Event

The initial rupture of *M*_{w} 6.3 was followed by the main strike‐slip *M*_{w} 7.2 event. The inversion of the main event of the El Mayor–Cucapah earthquake uses 40 teleseismic *P*‐wave records azimuthally dispersed at distances from 31° to 94°, as shown in Figure 6b, with spectral values chosen for the inversion from 0.0167 Hz (60 s) to 1 Hz (1 s). *P*‐wave modeling of the main event waveforms starts where the initial subevent waveforms end. Following the work of Wei *et al.* (2011), the main event hypocenter is the same as the subevent hypocenter. The geometry of the main event’s fault plane is simplified to be a single plane of 120 km (along strike) × 16 km (along dip), having 4 km×4 km subfaults (120 total subfaults), with a strike of 313°, a dip of 88°, and a rupture velocity of 2.0 km/s. The Green’s functions for each subfault are ±45° from the inferred rake of −180°. The *a priori* values (Table 1) are chosen based on the magnitude of the event as well as on previous published studies (Wei *et al.*, 2011; Mendoza and Hartzell, 2013).

#### Main Event: Single Velocity Model

The posterior solution obtained by implementing a single velocity model is shown in Figures 11c and 12c for slip and rise‐time distribution, respectively, for the main event of the 2010 El Mayor–Cucapah earthquake. Using a single velocity model, the peak slip inferred is 983 cm near the hypocenter, which is higher than the peak value found by Wei *et al.* (2011) of 700 cm but consistent with the value of Mendoza and Hartzell (2013). Wei *et al.* (2011) shows higher amplitude slip near the hypocenter with another strong source 30 km to the northwest and a third source 50 km to the southeast. The Bayesian solution shows a major source at the hypocenter with another source 35 km to the southeast. There tend to be more nonzero slip subfaults in the solution of Wei *et al.* (2011), which could be a consequence of their Laplacian smoothing. The differences with Wei *et al.* (2011) could also be due to their more complex fault parameterization. They subdivided the main event of the El Mayor–Cucapah earthquake into three distinct fault planes with different strikes and dips. Mendoza and Hartzell (2013) obtained a slip distribution using an even simpler fault geometry, consisting of a single fault plane for both the initial subevent and the main event. The *P* waves for the El Mayor–Cucapah earthquake are long and complicated, requiring that they be tapered at a point at which direct source arrivals end. Our solution shows a few subfaults at the far end of the fault with larger slip. These subfaults only contribute to the extreme ends of the *P* waveforms, near where they are tapered for the inversion, and may not be source related. The posterior mean solution (Fig. 11c) obtained using a single velocity model gives a seismic moment of 8.79×10^{19} N·m (*M*_{w} 7.3), which is slightly lower than that of the 9.90×10^{19} N·m (*M*_{w} 7.3) found by Wei *et al.* (2011). The rise‐time distribution (Fig. 12c) shows a similar pattern to the subevent. The slip velocity distribution (see Ⓔ the electronic supplement to this article) is characterized by the most energetic source at the hypocenter, with laterally decreasing slip velocities away from the hypocenter.

#### Main Event: Multiple Velocity Model

Figures 11d and 12d present the posterior mean solution, given velocity model variability in velocity and layer thickness for slip and rise time, respectively. The slip distribution shares similar characteristics with a single velocity model (Fig. 11c). The peak slip of 984 cm is nearly identical. The most prominent consequence of including velocity‐model uncertainty is the smoothing effect, causing slip to be distributed among a greater number of subfaults. The seismic moment increases from 8.79×10^{19} N·m (*M*_{w} 7.3) to 9.35×10^{19} N·m (*M*_{w} 7.3), which is more consistent with the finding of Wei *et al.* (2011). The rise‐time (Fig. 12d) and slip velocity distributions (see Ⓔ the electronic supplement to this article), given velocity model uncertainty, are similar to those of the single velocity solution.

#### Main Event: Uncertainties

Figure 5d presents scatter plots of rise times versus slip for the four numbered subfaults (27, 64, 67, and 73), in Figure 11c and 11d, for a single velocity model and multiple velocity models, respectively. Figure 5d shows further similarities to that of the previous cases, highlighting the slight linear trend between slip and rise time with multiple regions representing the multiple velocity models.

The associated uncertainties for the subfaults for the single velocity solution (Fig. 11c) and multiple velocity solution (Fig. 11d) are also plotted as histograms (see Ⓔ the electronic supplement to this article). The histograms show similar characteristics to the synthetic case, the 2011 Van earthquake, and the subevent of the 2010 El Mayor–Cucapah earthquake, showing broader distributions when considering multiple velocity models. With the implementation of velocity model variability, the uncertainty in slip increases by a factor of 3.9. Table 4 shows the mean values and standard deviations for the subfaults compared.

The model rise‐time histograms for a single velocity model and multiple velocity models are shown in Ⓔ the electronic supplement to this article, highlighting the same attributes presented in the slip distributions. Given multiple velocity models, the distribution of rise time broadens. Table 4 shows the mean values and corresponding standard deviations for the rise time of the subfaults compared. Introducing velocity model variability causes uncertainty to increase by a factor of 1.3.

#### Main Event: Waveform Fit

The comparison of the 40 teleseismic *P*‐wave records with the forward‐modeled waveforms calculated from the slip distribution of the single velocity model (Figs. 11c and 12c) and the forward‐modeled slip distribution using velocity model variability (Figs. 11d and 12d) are shown in Ⓔ the electronic supplement to this article. Given the complexity of the event, waveforms for both solutions do not fit the data as well as Wei *et al.* (2011; see Table 5 for L2 norms). To improve the model misfit, a more complex fault parameterization may be required.

## Conclusions

We conducted a finite‐fault analysis of the resolving powers and associated uncertainties of teleseismic body waves in a Bayesian framework using a synthetic case, the 2011 Van, east Turkey, earthquake, and the 2010 El Mayor–Cucapah, Baja California, earthquake. The Bayesian framework (CATMIP) allows us to refrain from using regularization techniques, dismissing smoothing, moment minimization, and positivity constraints on the solution, all of which are fairly standard in other common source inversion techniques. For each scenario, we present two distinct solutions to explore some of the epistemic uncertainty. The first solution utilizes a single velocity model, whereas the second solution utilizes 11 velocity models, varying the *P*‐ and *S*‐wave velocity and layer thickness.

The synthetic case shows degraded resolution near the hypocenter, with and without velocity model variability, as a consequence of small arrival‐time differences at teleseismic stations between subfaults, resulting in a spatial smearing. The rake is reasonably well resolved when no angle constraint is applied for the synthetic case, as is the moment. The addition of *SH* waves did not help constrain the rake. With the inclusion of multiple velocity models, the standard deviations for slip values increase by a factor of 3.3, increasing the uncertainty of the solution while sampling a wider range of values for each subfault. For all solutions (single velocity model with and without angle constraint and multiple velocity models), the waveform fits to the data are nearly identical.

The 2011 Van earthquake shows a spatial smoothing in slip when including velocity model variability, which increases the total number of nonzero slip subfaults. The rake is unresolvable when no angle constraint is applied when using *P* or *P* and *SH* waves. Standard deviations for the slip increase by a factor of 4.7 when considering velocity model variability compared with using a single velocity model. In general, the waveforms representing both solutions fit the data well, highlighting qualitatively small misfits between the data and solution waveforms. Both solutions show similar slip distributions to those presented by Hayes (2011) and Mendoza and Hartzell (2013), each using a different inversion scheme.

The 2010 El Mayor–Cucapah earthquake is the most complicated scenario, because the event is divided into an initial subevent and a main event. The subevent accounts for the first 15 s of record and contains relatively simplistic waveforms. Solutions with and without velocity model uncertainty show similar attributes to those of Wei *et al.* (2011). With the implementation of multiple velocity models, the standard deviations increase on average by a factor of 2.2, when compared with using a single velocity model. Given the simplistic nature of the subevent, both solutions with and without velocity model variability explain the data equivalently. The overall waveform fits are similar to that found by Wei *et al.* (2011).

The main source of the El Mayor–Cucapah earthquake highlighted similar characteristics to the previous scenarios, obtaining a slip distribution containing a spatial smearing effect when considering multiple velocity models. The slip distribution, with and without velocity model variability and an angle constraint, shows significant differences from the slip distributions found by Wei *et al.* (2011) and Mendoza and Hartzell (2013). However, Wei *et al.* (2011) utilized a more complicated fault geometry, employing three fault planes with different strikes and dips for the main event. Mendoza and Hartzell (2013) used a simpler fault model with a single fault plane for the entire event. When implementing multiple velocity models, the standard deviation increases by a factor of 3.9, in comparison with using a single velocity model. The synthetic waveforms, with and without multiple velocity models, do not fit the data as well as Wei *et al.* (2011), who used multiple data sets to help constrain the complexity of the earthquake.

Although each scenario shows unique attributes of the resolution and uncertainties, there are similar characteristic behaviors throughout all scenarios. Given a fixed rupture velocity with a specific fault geometry and parameterization, the small timing difference between subfaults around the hypocenter can induce a spatial smearing and degraded resolution. This is an issue with teleseismic data that can be removed if near‐field records are available. Velocity model variability introduces a greater number of nonzero subfaults, also creating a spatial smoothing effect. A multimodal distribution is best representative of the PDFs of slip and rise times, for certain subfaults, when considering velocity model variability. However, this is an artifact of the timing differences caused by the velocity models and the model parameterization, as well as by only considering a few velocity models. When a change in velocity structure is introduced, the timing difference can shift in which slip resides to adjacent subfaults. Rake is generally not well resolved without an angle constraint, except given a simple source; as such, we usually allow the rake to diverge ±10° from the inferred rake of the earthquake. The total seismic moment for each event is well resolved. The rise times corresponding to lower slip values are poorly resolved when considering either single or multiple velocity models. For the two earthquakes considered, there is a tendency for shorter rise times and large slip velocities near the hypocenter. The slip distributions with and without velocity model variability contain no dramatic differences, resulting, generally, in equivalent waveforms. The standard deviations of slip and rise time increase when considering uncertainty in the velocity model. The slip distributions obtained from teleseismic body waves in the Bayesian inference, excluding the main event of the 2010 El Mayor–Cucapah earthquake, were comparable to other studies using different and multiple data sets with different inversion schemes. The utilization of both teleseismic *P* and *SH* waves provided no additional constraints to help resolve the rake, possibly reflecting the difficulty of obtaining clean *SH* waveforms.

To foster the knowledge of the resolving power of teleseismic body waves in a Bayesian framework with associated uncertainties, further exploration should be conducted to see the effect on the posterior, including the effect of using different functional forms for STFs, choosing different model parameterizations, and utilizing more complicated fault geometries. With a Bayesian framework, the number of model parameters can be restricted compared with more conventional approaches, due to the computational effort involved. Even with a few hundred parameters, a supercomputer is desirable. However, significant insight is gained from an examination of the model parameter posterior PDFs that is not obtained from other methods.

## Data and Resources

All data used in this study were obtained through Incorporated Research Institutions for Seismology (IRIS) and can be downloaded from the IRIS Data Management Center at www.iris.edu (last accessed July 2016). Source mechanisms are from the Global Centroid Moment Tensor (CMT) Project database (www.globalcmt.org, last accessed July 2016). The algorithm for the Bayesian inversion is published in *Geophysical Journal International* (Minson *et al.*, 2013). The teleseismic body‐wave synthetics code is available at https://github.com/usgs/finite-fault (last accessed July 2016). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant Number ACI‐1053575. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing high performance computing resources (www.tacc.utexas.edu, last accessed July 2016).

## Acknowledgments

We would like to acknowledge Carlos Mendoza for helpful discussions; and thoughtful reviews from Martin Mai and an anonymous reviewer.

- Manuscript received 23 August 2016.