# Revisiting the 1992 Landers earthquake: a Bayesian exploration of co-seismic slip and off-fault damage

B. Gombert, Z. Duputel, R. Jolivet, C. Doubre, L. Rivera, M. Simons

► To cite this version:

B. Gombert, Z. Duputel, R. Jolivet, C. Doubre, L. Rivera, et al. Revisiting the 1992 Landers earthquake: a Bayesian exploration of co-seismic slip and off-fault damage. Geophysical Journal International, 2018, 212 (2), pp.839-852. 10.1093/gji/ggx455. hal-01639527

HAL Id: hal-01639527

https://hal.science/hal-01639527

Submitted on 5 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Revisiting the 1992 Landers earthquake: a Bayesian exploration of co-seismic slip and off-fault damage

B. Gombert,<fnref n="1" /> Z. Duputel,<fnref n="1" /> R. Jolivet,<fnref n="2" /> C. Doubre,<fnref n="1" /> L. Rivera<fnref n="1" /> and M. Simons<fnref n="3" />

Accepted 2017 October 17. Received 2017 October 16; in original form 2017 May 23

## SUMMARY

Existing models for the distribution of subsurface fault slip associated with the 1992 Landers, CA, earthquake ($M_w = 7.3$) show significant dissimilarities. In particular, they exhibit different amounts of slip at shallow depths (<5 km). These discrepancies can be primarily attributed to the ill-posed nature of the slip inversion problem and to the use of physically unjustifiable smoothing or regularization constraints. In this study, we propose a new coseismic model obtained from the joint inversion of multiple observations in a relatively unregularized and fully Bayesian framework. We use a comprehensive data set including GPS, terrestrial geodesy, multiple SAR interferograms and co-seismic offsets from correlation of aerial images. These observations provide dense coverage of both near- and far-field deformation. To limit the impact of modelling uncertainties, we develop a 3-D fault geometry designed from field observations, co-seismic offsets and the distribution of aftershocks. In addition, we account for uncertainty in the assumed elastic structure used to compute the Green's functions. Our solution includes the ensemble of all plausible models that are consistent with our prior information and fit the available observations within data and prediction uncertainties. Using near-fault high-resolution ground deformation measurements and the density of aftershocks, we investigate the properties of the damage zone and its impact on the inferred slip at depth. We attribute a part of the inferred slip deficit at shallow depth to our models not including the impact of a damage zone associated with a reduction of shear modulus in the vicinity of the fault.

**Key words:** Inverse theory; Probability distributions; Earthquake source observations; Fractures, faults, and high strain deformation zones.

## 1 INTRODUCTION

Following the 1979 Imperial Valley earthquake more than three decades ago (Olson & Apsel 1982; Hartzell & Heaton 1983), finite-fault source models have been routinely constructed after most significant earthquakes. Despite the increasing volume and quality of available geodetic and seismological data, we still observe a significant variability in inferred subsurface fault slip for a given event. Estimating the distribution of fault slip from surface deformation is fundamentally an ill-posed inverse problem with different models that can fit the data equally well. Therefore, different finite-fault models for the same earthquake often display significant dissimilarities. Over the past decade, there have been considerable efforts in the seismological community to study this problem and characterize the variability of the models (e.g. Mai et al. 2016). Furthermore, data and forward predictions are imperfect and the corresponding uncertainties are often difficult to account for. A standard approach to overcome the non-uniqueness of the solution relies on Tikhonov

regularization (e.g. Hansen 1998) involving minimization of first or second order spatial derivatives of the slip model to enforce smoothness of the slip distribution. However, various regularization strategies can affect the solution. The impact of different approaches to regularization, coupled with the lack of consideration of model uncertainties, can hamper our ability to draw clear conclusions about earthquake source processes.

Due to the availability of a comprehensive data set, many finite-fault models have been published for the 1992 $M_w = 7.3$ Landers earthquake (e.g. Murray et al. 1993; Cohee & Beroza 1994; Frey-mueller et al. 1994; Hudnut et al. 1994; Wald & Heaton 1994; Cotton & Campillo 1995; Fialko 2004b; Xu et al. 2016). Common patterns emerge in the inferred slip distributions including the fact that most of the slip occurred in the central section of the rupture (i.e. the Homestead Valley Fault). However, there are also clear inconsistencies. In particular, published studies have inferred shallow slip to vary between 30 per cent and 112 per cent of the slip inferred at 7 km depth. Since there is no indication of large inter- or

<footnote n="1">¹Institut de Physique du Globe de Strasbourg, UMR7516, Université de Strasbourg, EOST/CNRS, Strasbourg, France. E-mail: gombert@unistra.fr</footnote>

<footnote n="2">²Laboratoire de géologie, Département de Géosciences, École Normale Supérieure, PSL Research University, CNRS UMR Paris, France</footnote>

<footnote n="3">³Seismological Laboratory, Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USA</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** General overview of the area. (a) Tectonic context of Southern California. The dashed grey rectangle shows the extent of (b). The Landers earthquake surface rupture is plotted in red. The faults involved are part of the Eastern California Shear Zone (ECSZ). (b) Far-field observations used in this study. The thin black rectangles illustrate the InSAR track footprints. The ascending interferogram (Track 349) covers the time span between 26 May and 30 June 1992 and the descending interferogram (Track 399) between 24 April and 7 August. Topography is from the Space Radar Topographic Mission (SRTM) database.</figcaption></figure>

post-seismic slip at shallow depth (Shen et al. 1994; Savage & Svarc 1997; Fialko 2004a), the amount of the potential shallow co-seismic slip deficit has an impact on seismic risk assessment as this suggests that part of the accumulated strain is not released by the earthquake (Simons et al. 2002; Fialko et al. 2005). Simons et al. (2002) and Kaneko & Fialko (2011) suggested that such deficits might be an artefact due to inelastic response of the medium in the vicinity of the fault. Inelasticity would bias slip models where observations at short distances are modelled assuming elastic Green's functions. An apparent shallow slip deficit could also be caused by smoothing constraints and sparseness of near-fault data (Simons et al. 2002; Xu et al. 2016). Finally, unaccounted heterogeneities in the crust elastic properties can also result in a biased slip distribution at depth (Barbot et al. 2008). One way to evaluate these hypotheses is to derive all the models consistent with the available data without arbitrary regularization of the inverse problem and explore the potential mechanisms statistically.

We perform a Bayesian exploration of the 1992 Landers rupture to evaluate the population of plausible slip models given geodetic data and forward problem uncertainties. Our approach is exempt from any smoothing and allows us to assess the extent of any purported shallow slip deficit as constrained by available geodetic data. Using near-fault data, we also investigate the impact of lateral heterogeneities on the inferred slip distribution at depth.

## 2 DATA OVERVIEW

We use a large geodetic data set composed of GPS measurements at 82 sites, 23 trilateration measurements, 2 SAR interferograms and

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

14 optical correlation images. This combination of data provides good coverage in both the near- and far-fields.

### 2.1 GPS and trilateration data

We use 3-component observations from 82 GPS stations scattered across southern California (Hudnut et al. 1994) with a few stations in the vicinity of the fault (Figs 1 and 2). Observations of the vertical component of displacement are associated with significantly larger uncertainties than the horizontal components. In addition, a trilateration network covers the southern part of the rupture (Figs 1 and 2). We invert directly the horizontal relative line-length changes provided by Murray et al. (1993) instead of the pre-inverted displacement vectors of the trilateration stations. The GPS and trilateration data include up to a few months of inter-seismic and post-seismic deformation. However, the associated post-seismic displacements measured by GPS are expected to be less than ~10 cm, which is substantially smaller than the ~8 m of co-seismic displacement observed near the earthquake rupture. (Murray et al. 1993; Peltzer et al. 1998).

### 2.2 InSAR data

We use two SAR interferograms computed from pre- and post-earthquake acquisitions on both ascending and descending tracks of the ERS satellite (Fig. 1b). Interferograms are computed using the ROI_PAC software (Rosen et al. 2004). We downsample the unwrapped interferograms using a recursive quad-tree algorithm (Simons et al. 2002; Lohman & Simons 2005) to reduce the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Near-field observations. Lines are coloured according to length changes in the trilateration network. The optical correlation mosaic is plotted around the fault trace from Sieh *et al.* (1993). Main shock and Big Bear aftershock ($M_w = 6.5$) hypocentres from the Southern California Earthquake Center are indicated with a green and an orange star, respectively.</figcaption></figure>

number of observation points. The final downsampled ascending and descending interferograms contain 730 and 663 pixels, respectively. Downsampled observations, predictions, and residuals are shown in Supporting Information Fig. S1. Using the procedure described by Jolivet *et al.* (2014) for each InSAR scene, we estimate an empirical data covariance function, which statistically represents atmospheric noise. We find standard deviations of 3.5 cm and 0.9 cm for the descending and ascending tracks, respectively. The correlation length is 11 km for both images. Covariance functions are shown in Supporting Information Fig. S2. While the second image of the interferogram on the ascending track was acquired only 2 d after the main shock, the interferogram on the descending track includes more than one month of post-seismic deformation.

## 2.3 Optical correlation images

We use optical correlation images of the ground displacement from Ayoub *et al.* (2009). Maps of ground displacement are made using 14 pairs of aerial photographs acquired before and after the earthquake. Cross-correlation is performed to derive horizontal co-seismic displacements in the vicinity of the fault. Pre-earthquake photographs were acquired during the summer 1989 while post-earthquake were acquired during the autumn 1995. The footprint of each pair is slightly less than $10 \times 10 \text{ km}^2$ and the data set covers almost the entire surface rupture of the fault (Figs 2 and 7a). Because of their near-field coverage, optical data can finely constrain shallow slip in our models. However, as pointed out by Kaneko & Fialko (2011), near-fault observations may include inelastic effects that can bias slip estimates assuming linear elasticity. To avoid such artefact, we remove any near-fault pixels within 300 m of the fault. This cut-off length is in agreement with measurements by Milliner *et al.* (2015) showing that off-fault deformation is generally limited to a narrow

zone around the fault (with an average half-width smaller than 80 m). Removing data in the vicinity of the fault also reduces the impact of modelling errors due to fault parametrization. Indeed, the assumption of constant slip in fault patches and the discretization of the fault trace (every ~1.5 km) induce artefacts in the predicted deformation field very close to the fault (see Supporting Information text T1 and Fig. S3). In addition, using the same technique as for InSAR data in Section 2.2, each image is downsampled and data covariance is estimated using empirical covariograms. The resulting standard deviation is typically around 30 cm and the correlation length ranges from 300 m to 1 km. Most of the post-seismic deformation is included in the timespan separating the two acquisitions (Fialko 2004a). However, as mentioned by Milliner *et al.* (2015), the detection threshold of optical image correlation is about 10 cm, suggesting that ~15 cm of near-field post-seismic deformation lie in the uncertainties of the measurement.

# 3 PROBABILISTIC SLIP INVERSION

## 3.1 Model parametrization

While most previous studies used relatively simplified linear geometries, our fault parametrization shown in Fig. 3 consists of nine segments following the surface rupture trace. The three main segments are the Johnson Valley, Homestead Valley, and Emerson and Camp Rock faults (Sieh *et al.* 1993). Those three segments are linked by two small junctions and completed by the small Galway Lake Fault in the northern part of the rupture. In addition, we parametrize two antithetic faults on the eastern side of the Emerson segment. These two faults were not directly mapped by Sieh *et al.* (1993) but have been previously incorporated as linear segments by Fialko (2004b) from the distribution of aftershocks. In the present study, the northern antithetic segment is refined as a curved fault from the detailed analysis of InSAR ground deformation profiles along with the Hauksson *et al.* (2012) relocated earthquake catalogue (see Fig. 3). Finally, we use an additional fault corresponding to the $M_w = 6.5$ Big Bear aftershock, which orientation is derived from the Hauksson *et al.* (2012) catalogue. Consistent with Fialko (2004b), faults segments are assumed to be vertical and to extend down to 15 km. Although this depth is roughly in agreement with the maximum depth of aftershock, we cannot exclude a more complex geometry at depth as often reported when multiple fault segments interact (Segall & Pollard 1980). To evaluate the effect of such complexities, we propose an alternative geometry in which shallow parallel branches merge on a single deeper segment. This geometry is similar to a flower structure that can be observed in some strike-slip faults (e.g. Zigone *et al.* 2015).

For both assumed fault geometries, each segment is discretized in four rows of subfaults extending down to 1.5 km, 4.5 km, 9.0 km, and 15.0 km depth. The size of each subfault is designed to have an acceptable resolution at depth (resolution $R \ge 0.8$ as defined in the Supporting Information for the strike-slip component, see Supporting Information Fig. S4). This strategy ensures small posterior model uncertainty but more importantly, it enables good convergence of the Bayesian sampling algorithm used for the inversion.

## 3.2 Bayesian sampling

We use a Bayesian approach to obtain the full posterior probability density function (PDF) of the slip distribution given the observations and uncertainties. According to the Bayes–Laplace theorem, we
<figure><img src="image_5.png" /><figcaption>Figure 3. (a) Surface trace of the parametrized fault segments. Each segment is plotted as a thick black line. 1. Emerson and Camp Rock Faults, 2. Homestead Valley Fault, 3. Johnson Valley Fault, 4. Northern conjugate Fault, 5. Galway Lake Fault, 6. Southern conjugate Fault, 7. Emerson-Homestead Valley junction, 8. Kickapoo Fault, 9. Big Bear Fault. Blue dots represent aftershock locations from Hauksson et al. (2012). Dashed white rectangle shows the extent of (b). (b) Surface trace of the northern conjugate segment (dashed line). Rectangles show the position of the profiles shown in (c) and (d). Background colour represents the InSAR ascending track LOS displacement pattern. (c,d) InSAR data profiles A-A' and B-B'</figcaption></figure>

write the posterior PDF as:

$$p(\mathbf{m}|\mathbf{d}_{\text{obs}}) \propto p(\mathbf{m}) \exp \left[ -\frac{1}{2} (\mathbf{d}_{\text{obs}} - \mathbf{G}\mathbf{m})^{\mathrm{T}} \mathbf{C}_{\chi}^{-1} (\mathbf{d}_{\text{obs}} - \mathbf{G}\mathbf{m}) \right] \quad (1)$$

where $\mathbf{m}$ is the model vector, $p(\mathbf{m})$ is the prior distribution, $\mathbf{d}_{\text{obs}}$ is the data vector, $\mathbf{G}$ is the Green functions matrix, and $\mathbf{C}_{\chi}$ is the misfit covariance describing both data and forward prediction uncertainties. We compute the Green's functions for a semi-infinite stratified elastic medium using the EDKS software (Zhu & Rivera 2002)

To sample the model space we use AlTar, a parallel Markov Chain Monte Carlo (MCMC) algorithm based on the CATMIP formalism (Minson et al. 2013). Using multiple MCMC chains in parallel, AlTar initially samples the prior PDF, $p(\mathbf{m})$, and then slowly increases the information brought by the data until it samples the posterior PDF. The implementation benefits from the use of high efficiency Graphic Processing Units (GPUs), allowing us to run more than 500 000 chains in parallel. Our final solution consists of an ensemble of models that are statistically distributed according to the posterior PDF. No spatial smoothing constraint is used in this procedure. We adopt different priors for the two different slip directions. The strike-slip component prior is a uniform PDF between -1 m and 30 m, hence promoting right-lateral faulting. The dip-slip prior is a Gaussian PDF centred on 0 m with a standard deviation of 5 m.

## 3.3 Model prediction uncertainties

Accounting for uncertainties in our forward predictions uncertainties is crucial since they correspond to one of the largest sources of variability between published slip models. Moreover, these uncertainties are important in our Bayesian framework as we do not use smoothing regularization. The model prediction uncertainties are described by the matrix $\mathbf{C}_p$, which is added to the observation uncertainties matrix $\mathbf{C}_d$ to obtain the misfit covariance:

$$\mathbf{C}_{\chi} = \mathbf{C}_{d} + \mathbf{C}_{p}. \quad (2)$$

We build $\mathbf{C}_p$ using the approach of Duputel et al. (2014) to account for uncertainties in the elastic model used to compute the Green's functions. The layered elastic model used in this study is derived from the Southern California Earthquake Center 3-D velocity model (Kohler et al. 2003). Uncertainties on the elastic parameters are inferred by comparing different models in the source region along with the distribution of 3-D velocity models from Kohler et al. (2003), as shown in Supporting Information Fig. S5.

## 3.4 Probabilistic slip model

Using our Bayesian framework, we generate 500 000 models representing our posterior information on slip distribution given available geodetic data. To interpret this ensemble, we need to extract a representative model and the corresponding uncertainties. In Fig. 4, we show the posterior mean model (i.e. the average of all sampled models) along with 95 per cent confidence ellipses. A more detailed view is available in Supporting Information Fig. S6. The posterior mean model is a common choice as the Bayesian approach encourages one to think in terms of an ensemble solution instead of one single model. However, as shown in Supporting Information Fig. S7, other models can also be depicted such as the maximum a posteriori model (i.e. the mode of the posterior distribution) or the best fitting model (i.e. the sample in our population having the maximum posterior value). In our case, the maximum a posteriori model is insignificantly different to the posterior mean model since most
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Posterior mean co-seismic slip model. The colour of each subfault patch indicates the slip amplitude. Arrows and their associated 95 per cent confidence ellipse indicate the slip direction and uncertainty. The bottom left inset shows the potency normalized by patch row width as a function of depth. PDFs of shallow slip deficit (SSD) are presented for the entire fault system and for individual fault segments. Vertical lines on the same plots indicate the SSD of two published models (Cotton & Campillo 1995; Fialko 2004b).</figcaption></figure>

marginal PDFs are nearly Gaussians (cf. Supporting Information text T2).

The results in Fig. 4 are based on vertical fault segments. They can be compared with the solution in Supporting Information Fig. S8 obtained assuming a more complicated flower parametrization introduced in Section 3.1. Despite different fault dips, the inferred slip distributions are fairly similar in both geometries, showing the lack of sensitivity to the parametrization at depth. Although posterior PDFs of both geometries generally overlap in fault patches with large slip, we still observe significant differences as shown in Supporting Information Fig. S8. This suggests that modelling uncertainties included in $C_p$ are still underestimated as we only incorporate

Earth model uncertainties and neglect errors in the fault parametrization. In the following, we focus on the results obtained using vertical fault segments.

As expected, we observe predominately strike-slip motion along the entire fault system. Most of the slip concentrates along the central and northern parts of the rupture, with a peak amplitude of ~11 m. These features are to first order comparable to previous results, although published models have lower peak slip amplitudes (e.g. Cohee & Beroza 1994; Fialko 2004b; Xu et al. 2016). This difference is probably due to smoothing imposed in previous studies that decreases the maximum slip amplitude. The two small junctions (shown in Supporting Information Fig. S6) show relatively
<figure><img src="image_7.png" /><figcaption>Figure 5. Model performance for GPS and trilateration data. (a) GPS observations (blue) and predictions (red) with their 1σ error ellipses. (b) Length changes residuals for the posterior mean model.</figcaption></figure>

large slip at depth, although they are associated with significant posterior uncertainties. In addition, these estimates are associated with significant along-dip correlation of slip amplitudes (cf. Supporting Information Fig. S9).

The model predictions reproduce the observations reasonably well. The performance of the models for GPS and trilateration data is presented in Fig. 5 with associated posterior uncertainties. Posterior mean InSAR predictions and residuals are shown in Fig. 6 in high-resolution, and decimated in Supporting Information Fig. S1. In high-resolution, we observe some moderate residuals in the vicinity of the fault, mainly due to the finiteness of the fault patches. Some larger wavelength residuals are visible on the southern part of the descending track. We suspect that this signal originates from post-seismic deformation (Fialko 2004a) as the second pass of this track is 5 weeks after the main shock. Finally, our model explains reasonably well the optical correlation images despite large uncertainties associated with this data set (Fig. 7). We also computed an equivalent moment tensor and centroid location and tested it against long period seismological observations (details are provided in Supporting Information text T3 and Figs S10–S12.)

## 3.5 Shallow slip deficit

A shallow slip deficit is commonly observed for large strike-slip earthquakes (Simons et al. 2002; Fialko et al. 2005). Although, in a simple linear elastic model, a uniform slip distribution at depth is expected when averaged over many seismic cycles (Tse & Rice 1986), this deficit does not seem to be recovered by either inter-seismic creep or post-seismic deformation (Fialko 2004a). Some exceptions with no detectable shallow slip deficit have nonetheless been documented such as the 2013 M<sub>w</sub> = 7.7 Balochistan earthquakes (Jolivet et al. 2014; Vallage et al. 2015).

Although a shallow slip deficit is observed in most published models of the Landers earthquake, there is a large variability in the actual amount of shallow slip deficit between different inversion

results. To investigate this, we compute the normalized potency as a function of depth:

$$P_k = \frac{\sum_i \Delta u_{ik} \times A_{ik}}{w_k} \quad (3)$$

where Δu<sub>ik</sub> is the slip inferred in a patch of area A<sub>ik</sub> and width w<sub>k</sub> located in the kth row and at an along-strike position i. This formulation allows us to avoid any bias due to the increase of patch size with depth. As shown in Fig. 4, we find a maximum potency on the 3rd row of patches (i.e. between 4.5 and 9 km depth, consistent with Simons et al. (2002)) that is nearly 1.7 times larger than surface estimates (i.e. at depth between 0 and 1.5 km). To highlight this for individual fault segments, we define the percentage of shallow slip deficit (SSD) as:

$$SSD = 100 \left( \frac{P_{k=3} - P_{k=1}}{P_{k=3}} \right). \quad (4)$$

According to this definition, SSD > 0 indicates some amount of shallow slip deficit while SSD ≤ 0 means that potency is equal or larger at the surface than at depth (i.e. no shallow slip deficit). The posterior distribution of SSD is shown in Fig. 4 for the three main fault segments and the overall rupture. Results and probability estimates are also summarized in Table 1.

Although the overall rupture depicts a shallow slip deficit of about 41 per cent, we find different behaviours for different fault segments. We observe the smallest deficit along the Emerson and Camp Rock segment where the probability of shallow slip deficit is only 62%. The Johnson Valley fault is more likely to present a shallow deficit, but the SSD is relatively moderate (SSD~25 per cent). The largest deficit is measured for the Homestead Valley fault where the mean SSD is 52 per cent with a probability close to 1 that the deficit is larger than 25 per cent. The remaining fault segments are either too small, with too large uncertainties or did not slip enough to contribute significantly to the overall rupture estimate.
<figure><img src="image_8.png" /><figcaption>Figure 6. Model performance for InSAR. (a, d) InSAR observations. (b, e) Predictions for the posterior mean model. (c, f) InSAR residuals of the descending (top) and ascending (bottom) tracks.</figcaption></figure>

## 4 DISCUSSION

As pointed out in Section 1, previously published models differ, in particular regarding the amount of shallow slip deficit. A detailed comparison between our solution and previous models is provided in Fig. 4 and Supporting Information Fig. S13. The SSD values for previously published models extend from 70 per cent (i.e. a large shallow slip deficit) to −12 per cent (shallow slip exceeds slip at 7 km depth). Our slip deficit is thus smaller than some models (e.g. Zeng & Anderson 2000) but larger than others (e.g. Cohee & Beroza 1994; Wald & Heaton 1994; Cotton & Campillo 1995; Hernandez et al. 1999). Overall, there is a fairly good agreement with the model of Fialko (2004b) which closely matches our estimate of slip deficit. Unlike most of these previous models, our inversion includes near-field optical images which give a solid constraint on slip along the shallow part of the fault, hence improving our estimates of SSD. This is presented in Supporting Information Fig. S14 showing slip posterior uncertainties obtained with and without incorporating optical images, illustrating their significance in our inversion.

To assess the impact of smoothing constraints on the shallow slip deficit, we also performed damped least squares inversions incorporating a second-order Tikhonov regularization minimizing the roughness of the slip model $m_{est}$ (Segall & Harris 1987; Ortega 2013):

$$ m_{est}(\epsilon) = (G^T C_{\chi}^{-1} G + (\epsilon \nabla^2)^2)^{-1} G^T C_{\chi} d_{obs} \quad (5) $$

where $\nabla^2$ is the Laplacian operator defined on fault slip surface coordinates, and $\epsilon$ is the damping parameter. As shown in Supporting Information Figs S15(c)–(h), the larger the damping $\epsilon$, the smoother the solution. Supporting Information Fig. S15(a) shows that shallow slip deficit values vary widely as a function of $\epsilon$, from 13 per cent to 57 per cent. Unsurprisingly, models with little regularization (e.g. $\epsilon \sim 0.1$) are quite consistent with our Bayesian solution, including in terms of shallow slip deficit. The choice of $\epsilon$ is to a large extent arbitrary. However, we still note large variations of the SSD by selecting a few models localized around the corner of the L-curve (Supporting Information
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Model performance for optical image correlation data. (a) Observations. (b) Predictions for the posterior mean model. (c) Residuals. Positive displacements are towards the northwest (see arrow in the legend).</figcaption></figure>

<table><caption>**Table 1.** Shallow slip deficit estimated for different fault segments and for the whole rupture. A zero or negative SSD means that there is no deficit. An SSD value of 50 per cent means that there is twice more slip at depth than at the surface.</caption><thead><tr><th rowspan="2">Fault segment</th><th rowspan="2">Mean SSD</th><th rowspan="2">95 per cent conf. interval</th><th colspan="3">Probability than SSD is greater or equal than...</th></tr><tr><th>0 per cent</th><th>25 per cent</th><th>50 per cent</th></tr></thead><tbody><tr><td>Emerson and Camp Rock</td><td>2.6 per cent</td><td>-25.1 per cent to 33.5 per cent</td><td>62 per cent</td><td>3.1 per cent</td><td>0 per cent</td></tr><tr><td>Johnson Valley</td><td>25.4 per cent</td><td>-3.8 per cent to 57.6 per cent</td><td>94 per cent</td><td>58.2 per cent</td><td>&lt;1 per cent</td></tr><tr><td>Homestead Valley</td><td>51.7 per cent</td><td>42.7 per cent to 61.9 per cent</td><td>100 per cent</td><td>97.0 per cent</td><td>67.5 per cent</td></tr><tr><td>All faults combined</td><td>40.9 per cent</td><td>35.2 per cent to 47.3 per cent</td><td>100 per cent</td><td>99.9 per cent</td><td>&lt;1 per cent</td></tr><tr><td>All faults combined taking into account a compliant zone</td><td>29.6 per cent</td><td>14.32 per cent to 46.4 per cent</td><td>99.6 per cent</td><td>75.8 per cent</td><td>&lt;1 per cent</td></tr></tbody></table>

Fig. S15b). Such a strong dependence on $\epsilon$ complicates any interpretation of the results of smoothed models in terms of shallow slip deficit. Of course, other factors can possibly impact the inferred slip distribution such as the choice of fault geometry or the data sets included in the inversion. As shown in Supporting Information Table S2, we do not see any clear direct relationship between used data sets and the inferred SSD. For example, both Fialko (2004b) and Xu et al. (2016) used observations similar to ours but with different estimates of the SSD. Such variability does not seem to be explained by the assumed fault parametrization since both studies used a complex geometry similar to the one we use (cf. Supporting Information Table S2). Another example is Cohee & Beroza (1994) and Zeng & Anderson (2000) that are based on similar fault planes and data sets but with different SSD estimates. Inversion results can be affected by other parameters such as fault discretization, data weighting, and elastic structure (whose uncertainty is accounted for in the present study). A better understanding of the variability of previous models would require extensive tests using different geometries, data sets, and weighting schemes, which is beyond the scope of this study.

Different artefacts affecting co-seismic slip models are often proposed to explain the shallow slip deficit inferred for large strike-slip earthquakes. One of them is the inelastic strain in the vicinity of the fault that is usually unaccounted in finite-fault inversions (e.g. Simons et al. 2002; Fialko et al. 2005). Such inelastic response can indeed bias slip inversions that are based on elastic Green's functions and artificially decrease the amount of slip at shallow depth (Kaneko & Fialko 2011). However, as reported by Milliner et al. (2015), inelastic strain for the 1992 Landers earthquake is limited to a relatively narrow region around the fault (e.g. within ~65 m of the fault trace in Fig. 8c). To avoid any strong bias due to our elastic assumption and reduce modelling errors due to fault discretization at shallow depth, we have removed displacement data within a minimum distance of 300 m from the fault trace (see Section 2.3). This procedure is roughly equivalent to localizing the inelastic contribution of the strain field onto an idealized fault plane (Dahlen & Tromp 1998). Although removing near-fault pixels should reduce artefacts due to inelastic effects, unaccounted lateral heterogeneities due to accumulated damage around the fault can also have a significant impact on surface deformation patterns and by extension on the inverted slip distribution (Barbot et al. 2008).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Modelling of near-field deformation data. (a) Overall view of optical correlation data. The profile shown in (c) is localized with a black line. (b) Close up view of near-field data. Grey rectangle indicate the location of the profile shown in (c). (c) Comparison between observed displacement (in red) and the stochastic predictions (in grey). Black arrows labelled F1 and F2 in (b) and (c) highlight two small secondary ruptures visible in the data. These small ruptures are incorporated in our modelling approach assuming two vertical dislocations. Data inside the black brackets are not used in the inversion of the full 3-D slip distribution presented in Fig. 4 to reduce the impact of inelastic effects in the vicinity of the main rupture.</figcaption></figure>

The fault zone is often regarded as a highly deformed core sur-
rounded by a more or less broad damage zone of reduced stiffness
(e.g. Chester *et al.* 1993; Ben-Zion & Sammis 2003; Dor *et al.* 2006;
Mitchell & Faulkner 2009). The damage zone consists of cracks
and microfractures in the host rock and can be associated with sec-
ondary faults reducing the elastic strain released on the main rup-
ture interface (Chester & Chester 1998; Dieterich & Smith 2009).
Such secondary cracks have been reported around the Landers fault
system (McGill & Rubin 1999). An example is given in Fig. 8,
showing two secondary ruptures (labelled F1 and F2) visible in
optical correlation images near the Emerson Valley fault. Such off-

fault ruptures are not accounted for in our slip model presented in
Fig. 4.

To investigate the properties of the damage zone and secondary
ruptures, we analyse a profile across the fault using simple verti-
cal elastic screw dislocations embedded in a compliant fault zone
(Segall 2010). Using a Metropolis algorithm, we invert for the slip
distribution on each fault, a compliant zone half-width and an effec-
tive shear modulus contrast μ₁/μ₀ (where μ₁ is the shear modulus
of the fault zone while μ₀ is the modulus of the surrounding crust).
The compliant zone half-width and shear modulus ratio being typ-
ical Jeffreys parameters (Tarantola 2005), they are sampled in the
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Posterior joint probability distribution of the compliant zone half-width and shear modulus ratio. Dots are model samples that are coloured according to the PDF value. Blue histograms are marginal PDFs for both parameters.</figcaption></figure>

logarithmic domain. To avoid any effect of off-fault inelasticity, we remove the data within 65 m of the fault, which is consistent with fault-width measurements by Milliner et al. (2015) at this location. The results presented in Fig. 8(c) indicate very shallow secondary ruptures with 32 ± 8 cm and 36 ± 5 cm of slip down to 84 ± 30 m and 180 ± 40 m respectively for faults F1 and F2. Although such slip amplitudes are not negligible, these off-fault dislocations are relatively shallow and thus represent only 3.3 per cent of the total seismic slip inferred from the surface down to 0.5 km. Of course, these measurements are only valid locally since the properties of secondary faults might vary significantly along the main rupture (Lewis & Ben-Zion 2010; Milliner et al. 2015; Thomas et al. 2017).

The results shown in Fig. 9 highlight the existence of a ~1.1 km wide compliant zone around this part of the fault. Although there is some correlation between the compliant zone width and rigidity, our solution indicates that shear modulus can be reduced by as much as a factor ~5 within the damage zone (i.e. a shear modulus ratio of ~0.2). This estimate is consistent with measurements from guided seismic waves (Li et al. 1994, 2007; Peng et al. 2003) that indicate shear modulus ratios between 0.1 and 0.4, corresponding to 80 per cent of our models. On the other hand, these studies suggest relatively small damage zone widths of a few hundred meters, which is narrower than our estimates.

Using the aftershock catalogue of Hauksson et al. (2012), we compare our estimates with the distribution of seismicity around the main fault, which is another indicator of distributed damage in the host rock (Amitrano 2006; Powers & Jordan 2010). As shown in Fig. 10(a), we select two profiles across the main rupture surrounding the southern antithetic fault to avoid any bias due to events located on this segment. Following Powers & Jordan (2010), we

compute the horizontal density $v(x)$ of seismicity where $x$ is the fault normal distance, and assume a power law decay of the form

$$v(x) = v_0 \left( 1 + \frac{x^2}{d^2} \right)^{-\gamma/2} \qquad (6)$$

where $v_0$ is the aftershock density at $x = 0$, $d$ is the damage zone half-width and $\gamma$ is the asymptotic roll-off of the seismicity away from the fault. Using a Metropolis inversion scheme, we then sample $v_0$, $d$, and $\gamma$ given the seismicity density, $v(x)$. Comparison between observations and stochastic predictions are shown in Fig. 10(b) and the full posterior PDFs for the 3 parameters are shown in Supporting Information Fig. S16. Although the posterior mean damage-zone half-width $d \sim 800$ m is larger than what is inferred from optical images ($d \sim 570$ m), an inversion with a fixed $d = 570$ m also explains the data reasonably well (cf. Fig. 10b).

To estimate the impact of the damage zone on the inverted slip distribution, we also invert the fault-parallel displacement profile of Fig. 8(c) without a compliant zone and after removing the data within 300 m of the fault (i.e. the same way it is done in our main slip inversion). The posterior PDFs of shallow slip and stochastic predictions with and without accounting for the damage zone are shown in Fig. 11. Although far field deformation is well-predicted in both inversions, predictions neglecting a compliant zone fail to reproduce near-fault observations and underestimate slip at shallow depth. On average, accounting for the compliant zone increases shallow slip by a factor of 1.2. On the other hand, neglecting lateral shear modulus heterogeneities will systematically lead to smaller slip (with a probability of 98 per cent). To roughly estimate the effect of the damage zone, we can empirically correct the surface mean slip of the Landers rupture by factors drawn from posterior PDFs with and without accounting for the compliant zone. Results
<figure><img src="image_12.png" /></figure>

**Figure 10.** Distribution of seismicity across the fault. (a) Our parametrized fault trace is indicated with thick black lines. Blue dots are aftershock epicentres from Hauksson *et al.* (2012). Grey rectangles illustrate the location of profiles used for the seismicity density analysis. (b) Seismicity density as a function of fault normal distance. Densities are computed over the two stacked profiles using 100 m wide distance bins. Black circles are resulting event density measurements used in the power-law inversion. Red circles are observations not included in the inversion since they correspond to events located at distance larger than ~2 km that may be partly linked to the southern antithetic fault segment. The 1σ error bars were obtained by computing the standard deviation of density in each bin from 1000 random catalogues generated according to event location uncertainties.

presented in Fig. 12 and Table 1, indicate that this significantly reduces the overall shallow slip deficit from 41 per cent to 27 per cent. These results should, however, be considered with caution, as the damage behaviour can vary significantly along the fault (Lewis & Ben-Zion 2010). We tried to conduct similar experiments in other locations on the fault but did not obtain reliable constraints on the compliant zone parameters (see e.g. Supporting Information Figs S17 and S18). Even if damage properties can widely vary along the fault, such structures will necessarily impact slip estimated at shallow depth, thereby reducing the inferred shallow deficit.

## 5 CONCLUSION

We used an extensive geodetic data set, careful uncertainty estimates and a realistic fault geometry to produce a stochastic finite-fault model of the Landers earthquake. Our Bayesian approach to the inversion has two main advantages: (1) the solution is not biased by any kind of smoothing and (2) posterior parameter uncertainties are available and provide valuable information on the validity of the

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

**Figure 11.** Comparison between shallow slip posterior PDFs assuming an homogeneous half-space (in blue) and accounting for a damage zone of reduced stiffness (in green). The inset shows stochastic predictions for both inversions. Observations are plotted as a thin black line. Blue results are inferred without the data inside the brackets in Fig. 8(c) and green results without the data inside the red brackets at ±65 m

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

**Figure 12.** Overall shallow slip deficit (SSD). The black PDF indicate the SSD for the overall rupture presented in Fig. 4. The purple PDF is the SSD corrected from the effect of the damage zone with reduced stiffness. Blue and green vertical lines are the SSD for two published models (Cotton & Campillo 1995; Fialko 2004b).

model. The predictions from our solution agree well with various observations.

Consistent with previous studies, our solution indicates a substantial shallow slip deficit that is particularly pronounced for the Homestead Valley Fault. We argue that part of this deficit results from unmodelled lateral heterogeneities in shear modulus, corresponding to a damage zone surrounding the fault. Using high resolution optical correlation images, we highlight a ~1 km wide damage zone on the Emerson Valley Fault responsible for an apparent reduction in shallow slip by a factor ~1.2. Our results also show the presence of secondary ruptures with significant slip amplitudes at shallow depth. By reducing the elastic strain on the main fault, these features also contribute to the apparent slip deficit budget.

Although we do not include data in the immediate vicinity of the fault where inelastic behaviour is commonly observed, we cannot
rule out that some wide plastic deformation is included in our in-
version and participates in the observed deficit. Following the same
procedure, other near-field displacement data of large strike-slip
earthquakes could provide new insights on fault zone properties
and their link to co-seismic slip distribution.

ACKNOWLEDGEMENTS

We are grateful to J.P. Avouac and F. Ayoub for providing the optical
images correlation observations. We also thank Y. Fialko for sending
us his co-seismic slip model. This study contributed from fruitful
discussions with Sarah Minson, Lijun Zhu, Michael Aivazis and
Gilles Peltzer. Some GPUs used for this research were donated by
the NVIDIA Corporation. This work was supported by the Initiative
d'Excellence (IDEX) funding framework (Université de Strasbourg)
and the CNRS PICS program (Zacharie Duputel). This work was
also funded by NSF grant 1447107 awarded to Mark Simons. We
thank the Editor, Martin Mai, and two anonymous reviewers for their
constructive comments, which helped improve this manuscript.

REFERENCES

Amitrano, D., 2006. Rupture by damage accumulation in rocks, *Int. J. Fract.*, **139**(3), 369–381.

Ayoub, F., Leprince, S. & Avouac, J.-P., 2009. Co-registration and correlation of aerial photographs for ground deformation measurements, *ISPRS J. Photogramm. Remote Sens.*, **64**(6), 551–560.

Barbot, S., Fialko, Y. & Sandwell, D., 2008. Effect of a compliant fault zone on the inferred earthquake slip distribution, *J. geophys. Res.*, **113**(B6), B06404, doi:10.1029/2007JB005256.

Ben-Zion, Y. & Sammis, C. G., 2003. Characterization of fault zones, in *Seis-
mic Motion, Lithospheric Structures, Earthquake and Volcanic Sources:
The Keiti'i Aki Volume*, pp. 677–715, Springer.

Chester, F.M. & Chester, J.S., 1998. Ultracataclasite structure and friction processes of the Punchbowl fault, San Andreas system, California, *Tectonophysics*, **295**(1), 199–221.

Chester, F.M., Evans, J.P. & Biegel, R.L., 1993. Internal structure and weak-
ening mechanisms of the San Andreas fault, *J. geophys. Res.*, **98**(B1),
771–786.

Cohee, B.P. & Beroza, G.C., 1994. Slip distribution of the 1992 Landers earthquake and its implications for earthquake source mechanics, *Bull. seism. Soc. Am.*, **84**(3), 692–712.

Cotton, F. & Campillo, M., 1995. Frequency domain inversion of strong motions: application to the 1992 Landers earthquake, *J. geophys. Res.*, **100**(B3), 3961–3975.

Dahlen, F. & Tromp, J., 1998. *Theoretical Global Seismology*. Princeton University Press.

Dieterich, J.H. & Smith, D.E., 2009. Nonplanar faults: Mechanics of slip and off-fault damage, *Pure appl. Geophys.*, **166**(10-11), 1799-1815.

Dor, O., Ben-Zion, Y., Rockwell, T.K. & Brune, J., 2006. Pulverized rocks in the Mojave section of the San Andreas Fault Zone, *Earth planet. Sci. Lett.*, **245**(3), 642–654.

Duputel, Z., Agram, P.S., Simons, M., Minson, S.E. & Beck, J.L., 2014. Accounting for prediction uncertainty when inferring subsurface fault slip, *Geophys. J. Int.*, **197**(1), 464–482.

Fialko, Y., 2004a. Evidence of fluid-filled upper crust from observations of postseismic deformation due to the 1992 Mw7.3 Landers earthquake, *J. geophys. Res.*, **109**(B8), B08401, doi:10.1029/2004JB002985.

Fialko, Y., 2004b. Probing the mechanical properties of seismically active crust with space geodesy: study of the coseismic deformation due to the 1992 Mw 7.3 Landers (southern California) earthquake, *J. geophys. Res.*, **109**(B3), B03307, doi:10.1029/2003JB002756.

Fialko, Y., Sandwell, D., Simons, M. & Rosen, P., 2005. Three-dimensional deformation caused by the Bam, Iran, earthquake and the origin of shallow slip deficit, *Nature*, **435**(7040), 295–299.

Freymueller, J., King, N. & Segall, P., 1994. The co-seismic slip distribution of the Landers earthquake, *Bull. seism. Soc. Am.*, **84**(3), 646–659.

Hansen, P.C., 1998. *Rank-deficient and Discrete Ill-posed Problems: Numerical Aspects of Linear Inversion*, SIAM.

Hartzell, S.H. & Heaton, T.H., 1983. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake, *Bull. seism. Soc. Am.*, **73**(6A), 1553–1583.

Hauksson, E., Yang, W. & Shearer, P.M., 2012. Waveform relocated earthquake catalog for southern California (1981 to June 2011), *Bull. seism. Soc. Am.*, **102**(5), 2239–2244.

Hernandez, B., Cotton, F. & Campillo, M., 1999. Contribution of radar interferometry to a two-step inversion of the kinematic process of the 1992 Landers earthquake, *J. geophys. Res.*, **104**(B6), 13 083–13 099.

Hudnut, K.W. et al., 1994. Co-seismic displacements of the 1992 Landers earthquake sequence, *Bull. seism. Soc. Am.*, **84**(3), 625–645.

Jolivet, R. et al., 2014. The 2013 Mw 7.7 Balochistan earthquake: seismic potential of an accretionary wedge, *Bull. seism. Soc. Am.*, **104**(2), 1020–1030.

Kaneko, Y. & Fialko, Y., 2011. Shallow slip deficit due to large strike-slip earthquakes in dynamic rupture simulations with elasto-plastic off-fault response, *Geophys. J. Int.*, **186**(3), 1389–1403.

Kohler, M., Magistrale, H. & Clayton, R., 2003. Mantle heterogeneities and the SCEC reference three-dimensional seismic velocity model version 3, *Bull. seism. Soc. Am.*, **93**(2), 757–774.

Lewis, M.A. & Ben-Zion, Y., 2010. Diversity of fault zone damage and trapping structures in the Parkfield section of the San Andreas Fault from comprehensive analysis of near fault seismograms, *Geophys. J. Int.*, **183**(3), 1579–1595.

Li, H., Zhu, L. & Yang, H., 2007. High-resolution structures of the Landers fault zone inferred from aftershock waveform data, *Geophys. J. Int.*, **171**(3), 1295–1307.

Li, Y.-G., Vidale, J.E., Aki, K., Marone, C.J. & Lee, W.H., 1994. Fine structure of the Landers fault zone: segmentation and the rupture process, *Science*, **265**, 367–367.

Lohman, R.B. & Simons, M., 2005. Some thoughts on the use of In-
SAR data to constrain models of surface deformation: Noise structure
and data downsampling, *Geochem. Geophys. Geosyst.*, **6**(1), Q01007,
doi:10.1029/2004GC000841.

Mai, P.M. et al., 2016. The earthquake-source inversion validation (SIV) project, *Seismol. Res. Lett.*, **87**(3), 690–708.

McGill, S.F. & Rubin, C.M., 1999. Surficial slip distribution on the central Emerson fault during the June 28, 1992, Landers earthquake, California, *J. geophys. Res.*, **104**(B3), 4811–4833.

Milliner, C.W., Dolan, J.F., Hollingsworth, J., Leprince, S., Ayoub, F. & Sammis, C.G., 2015. Quantifying near-field and off-fault deformation patterns of the 1992 Mw 7.3 Landers earthquake, *Geochem. Geophys. Geosyst.*, **16**(5), 1577–1598.

Minson, S., Simons, M. & Beck, J., 2013. Bayesian inversion for finite fault earthquake source models I—Theory and algorithm, *Geophys. J. Int.*, **194**(3), 1701–1726.

Mitchell, T. & Faulkner, D., 2009. The nature and origin of off-fault damage surrounding strike-slip fault zones with a wide range of displacements: a field study from the Atacama fault system, northern Chile, *J. Struct. Geol.*, **31**(8), 802–816.

Murray, M., Savage, J., Lisowski, M. & Gross, W., 1993. Coseismic displacements: 1992 Landers, California, earthquake, *Geophys. Res. Lett.*, **20**(7), 623–626.

Olson, A.H. & Apsel, R.J., 1982. Finite faults and inverse theory with applications to the 1979 Imperial Valley earthquake, *Bull. seism. Soc. Am.*, **72**(6A), 1969–2001.

Ortega, F., 2013. Aseismic Deformation in Subduction Megathrusts: Central Andes and North-East Japan, *PhD thesis*, California Institute of Technology.

Peltzer, G., Rosen, P., Rogez, F. & Hudnut, K., 1998. Poroelastic rebound along the Landers 1992 earthquake surface rupture, *J. geophys. Res.*, **103**(B12), 30 131–30 145.

Peng, Z., Ben-Zion, Y., Michael, A.J. & Zhu, L., 2003. Quantitative analysis of seismic fault zone waves in the rupture zone of the 1992 Landers,
California, earthquake: evidence for a shallow trapping structure, *Geo-phys. J. Int.*, **155**(3), 1021–1041.

Powers, P.M. & Jordan, T.H., 2010. Distribution of seismicity across strike-slip faults in California, *J. geophys. Res.*, **115**(B5), B05305, doi:10.1029/2008JB006234.

Rosen, P.A., Hensley, S., Peltzer, G. & Simons, M., 2004. Updated repeat orbit interferometry package released, *EOS, Trans. Am. geophys. Un.*, **85**(5), 47.

Savage, J.C. & Svarc, J.L., 1997. Postseismic deformation associated with the 1992 M<sub>w</sub>=7.3 Landers earthquake, southern California, *J. geophys. Res.*, **102**(B4), 7565–7577.

Segall, P., 2010. *Earthquake and Volcano Deformation*, Princeton University Press.

Segall, P. & Harris, R., 1987. Earthquake deformation cycle on the San Andreas Fault near Parkfield, California, *J. geophys. Res.*, **92**(B10), 10 511–10 525.

Segall, P. & Pollard, D., 1980. Mechanics of discontinuous faults, *J. geophys. Res.*, **85**(B8), 4337–4350.

Shen, Z.-K., Jackson, D.D., Feng, Y., Cline, M., Kim, M., Fang, P. & Bock, Y., 1994. Postseismic deformation following the Landers earthquake, California, 28 June 1992, *Bull. seism. Soc. Am.*, **84**(3), 780–791.

Sieh, K. et al., 1993. Near-field investigations of the Landers earthquake sequence, April to July 1992, *Science*, **260**, 171–176.

Simons, M., Fialko, Y. & Rivera, L., 2002. Coseismic deformation from the 1999 Mw 7.1 Hector Mine, California, earthquake as inferred from InSAR and GPS observations, *Bull. seism. Soc. Am.*, **92**(4), 1390–1402.

Tarantola, A., 2005. *Inverse Problem Theory and Methods for Model Parameter Estimation*, SIAM.

Thomas, M.Y., Bhat, H.S. & Klinger, Y., 2017. Effect of Brittle Off-Fault Damage on Earthquake Rupture Dynamics, in *Fault Zone Dynamic Processes: Evolution of Fault Properties During Seismic Rupture*, vol. 227, p. 255, eds Thomas, M. Y., Mitchell, T. M. & Bhat, H. S., John Wiley & Sons, Inc.

Tse, S.T. & Rice, J.R., 1986. Crustal earthquake instability in relation to the depth variation of frictional slip properties, *J. geophys. Res.*, **91**(B9), 9452–9472.

Vallage, A., Klinger, Y., Grandin, R., Bhat, H. & Pierrot-Deseilligny, M., 2015. Inelastic surface deformation during the 2013 Mw 7.7 Balochistan, Pakistan, earthquake, *Geology*, **43**(12), 1079–1082.

Wald, D.J. & Heaton, T.H., 1994. Spatial and temporal distribution of slip for the 1992 Landers, California, earthquake, *Bull. seism. Soc. Am.*, **84**(3), 668–691.

Xu, X., Tong, X., Sandwell, D.T., Milliner, C.W., Dolan, J.F., Hollingsworth, J., Leprince, S. & Ayoub, F., 2016. Refining the shallow slip deficit, *Geophys. J. Int.*, **204**(3), 1867–1886.

Zeng, Y. & Anderson, J., 2000. *Evaluation of numerical procedures for simulating near-fault long-period ground motions using Zeng method, Report 2000/01 to the PEER Utilities Program*, available at http://peer.berkeley.edu.

Zhu, L. & Rivera, L.A., 2002. A note on the dynamic and static displacements from a point source in multilayered media, *Geophys. J. Int.*, **148**(3), 619–627.

Zigone, D., Ben-Zion, Y., Campillo, M. & Roux, P., 2015. Seismic tomography of the Southern California plate boundary region from noise-based Rayleigh and Love waves, *Pure appl. Geophys.*, **172**(5), 1007–1032.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Supporting text T1:** Modelling errors due to fault parameterization

**Supporting Text T2:** Maximum a posteriori model

**Supporting text T3:** Model prediction for long period seismological data

**Figure S1.** Decimated InSAR images. (a, d) Decimated InSAR observations used in the inversion. (b, e) Predictions for the posterior mean model. (c, f) InSAR residuals of the descending (top) and ascending (bottom) tracks.

**Figure S2.** Empirical covariance functions for the InSAR observations: 1D empirical covariance functions and the associated best-fit exponential function for the descending (left) and ascending (right) tracks. For each image, we compute the empirical covariance as a function of the distance between pixels and then fit an exponential function to these covariances (Jolivet et al., 2012). This exponential function is then used to build the data covariance matrix used in the inversion.

**Figure S3.** Effect of geometry on forward modelling. (a) Forward model predictions for the one of the optical images mosaic imposing 4.1m of slip on a shallow fault with 300m long patches. (b) same as (a) but with a broader geometry (1.5km-long patches). (c) and (d) Difference between (a) and (b).

**Figure S4.** Problem resolution. For each slip component, we compute the Resolution matrix as $R = C_m G^T (G C_m G^T + C_\chi)^{-1} G \cdot C_m$ is a diagonal matrix constructed from our model a priori distribution standard deviation. The diagonal values are plotted on the fault. The closer to 1, the better is the resolution of the parameter.

**Figure S5.** Different models variability of the P-wave, S-wave, and density as a function of depth in the Landers area. Grey lines are model values of the 3D Community Velocity Model (CVM, Kohler et al. 2003) available at http://scedc.caltech.edu/research-tools/3d-velocity.html (last accessed January 2016). The dashed black line represents the averaged CVM value for this area. A layered version used in this study for Green's function [GF] calculations is plotted as a solid black line. Models from Cotton & Campillo (1993), Wald & Heaton (1994), Hauksson (1993), and Jones and Helmberger (1998) are plotted as solid green, dashed green, red, and blue lines, respectively. Grey histograms are the probability density function representing our confidence level on the elastic properties, as used to build the model prediction error. Histograms are derived from the averaged CVM assuming a Gaussian distribution.

**Figure S6.** Posterior mean co-seismic slip model. The color of each subfault patch indicates the slip amplitude. Arrows and their associated 95% confidence ellipse indicate the slip direction and uncertainty.

**Figure S7.** Comparison between posterior mean, maximum a posteriori and best fitting models. (a) Maximum a posteriori coseismic slip model. It is built by considering the maximum of each marginal PDF (cf., supporting text T2). The 10 patches where the slip is the most important are labelled in purple. (b) Best-fitting model sample. This model represents the sample in our population having the maximum posterior value (cf., supporting text T2). The colour of each subfault patch indicates the slip amplitude. Arrows and their associated 95% confidence ellipse indicate the slip direction and uncertainty. (c) Boxplot of the strike-slip within the 10 patches labelled in (a). Horizontal red lines show posterior mean values (Figure 5 in main text). Horizontal blue lines show the maximum a posteriori model (a), and horizontal green lines show the best fitting sample (b). Notice that the best fitting sample is a poor estimate of the MAP.

**Figure S8.** The subfigure at the centre is the posterior mean coseismic slip models for an alternative “flower” geometry. The colour of each subfault patch indicates the slip amplitude. Arrows and their associated 95% confidence ellipse indicate the slip direction and uncertainty. The patches that slip the most in the vertical geometry are numbered from 1 to 10. We show the PDF of SSD as a black line on the bottom-right insert. The magenta line illustrates the SSD value when corrected from a compliant zone. On the same plots are
represented the SSD for two published models, Cotton & Campillo (1995) and Fialko (2004b). The histograms on the sides show the strikeslip PDF of the 10 patches that are labelled on the finite-fault model for both the vertical and flower geometries. The percentage of which the two PDFs overlap is given on the top-left corner of each histogram.

**Figure S9.** Posterior covariance of two along-dip patches of the Homestead Valley segment and the Homestead Valley Camp Rock segments junction. (a) Homestead Valley and (b) Homestead Valley Camp Rock junction posterior mean coseismic slip. In each one of the two segments, the across-patch correlation is computed for the two coloured patches. (c) Joint posterior PDF of the strike-slip component of the two coloured patches in (a), also labelled 1 and 2. (d) Same as (c), but for the two coloured patches in (b) (labelled 3 and 4). For both (c) and (d), dots are model samples that are coloured according to the PDF value. Blue histograms are marginal PDFs for both parameters.

**Figure S10.** Red dots indicate the posterior ensemble of centroid locations derived from our solution. The red focal mechanism is the moment tensor computed from our posterior mean model. The blue and yellow focal mechanisms come from the Global CMT (Ekström et al., 2012) and W-Phase (Duputel et al. 2012) catalogues, respectively.

**Figure S11.** Broadband seismograms (black line) and synthetics computed from the posterior mean model moment tensor (red line) are plotted for 5 stations along with their locations. On each map, the blue star and the red dot indicate the hypocenter and station locations, respectively. For each trace is indicated the station azimuth $\varphi$ and epicentral distance $\Delta$.

**Figure S12.** Figure S9 continued. Broadband seismograms (black line) and synthetics computed from the posterior mean model moment tensor (red line) are plotted for 5 stations along with their locations. On each map, the blue star and the red dot indicate the hypocenter and station locations, respectively. For each trace is indicated the station azimuth $\varphi$ and epicentral distance $\Delta$.

**Figure S13.** Comparison of SSD values. The thick black line is the probability density of SSD values for this study. Vertical coloured lines represent the SSD values of 6 published models.

**Figure S14.** Posterior slip uncertainties for (a) the solution obtained by inverting all datasets and (b) the solution obtained by inverting all the datasets minus the optical correlation mosaic.

**Figure S15.** Impact of smoothing constraints on the shallow slip deficit (SSD). (a) SSD value of models obtained by a least-square inversion as a function of the damping parameter $\epsilon$ (see equation 5 in the main text). Red dots indicate the models shown in (c) to (g). The horizontal dashed line marks the mean SSD value of our stochastic solution, and the grey shaded area represent the $1-\sigma$ deviation. (b) L-curve of the regularized models. Dots colour indicates the damping value. The red rectangle shows the extent of the top-right inset. The position of the models (c-h) is indicated with their damping value. (c)-(h) Least-square models for six different damping values.

**Figure S16.** Results of the Metropolis sampling of the after-shock density profile parameters $v_0$, $d$ and $\gamma$. 1D plots are posterior marginal PDFs and 2D plots are posterior joint PDFs. On the 2D histograms dots are model samples that are coloured according to the PDF value. Hot colours indicate region of high-probability.

**Figure S17.** Modeling of Near-field deformation data. (a) Localization of the profiles in the optical correlation observations. (b) Close up view of near-field data. Grey rectangle indicate the location of the inverted profile. (c) Comparison between observed displacement (in red) and the stochastic predictions (in grey). Data inside the black brackets are not used in the inversion of the full 3D slip distribution presented in Figure 6 of the main manuscript.

**Figure S18.** Posterior joint probability distribution of the compliant zone half-width and shear modulus ratio for the profile presented in Figure S12. Dots are model samples that are coloured according to the PDF value. Blue histograms are marginal PDFs for both parameters. Red lines are the prior information used in the sampling.

**Table S1.** Source parameters for 3 solutions. Wphase (Duputel et al., 2012), Global CMT (Ekström, 2012), and this study.

**Table S2.** Summary of fault geometries and datasets used in this study and in previously published models, and associated SSD values. In the InSAR column, the number in parentheses is the number of interferograms used.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.