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

Revisiting the 1992 Landers earthquake: a Bayesian exploration of co-seismic slip and off-fault... Abstract Existing models for the distribution of subsurface fault slip associated with the 1992 Landers, CA, earthquake (Mw = 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. 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 Mw = 7.3 Landers earthquake (e.g. Murray et al. 1993; Cohee & Beroza 1994; Freymueller 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 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 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). Figure 1. View largeDownload slide 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. Figure 1. View largeDownload slide 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. Figure 2. View largeDownload slide 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 (Mw = 6.5) hypocentres from the Southern California Earthquake Center are indicated with a green and an orange star, respectively. Figure 2. View largeDownload slide 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 (Mw = 6.5) hypocentres from the Southern California Earthquake Center are indicated with a green and an orange star, respectively. 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 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 × 10 km2 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 Mw = 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). Figure 3. View largeDownload slide (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 conjugated 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’ Figure 3. View largeDownload slide (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 conjugated 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’ 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 ≥ 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 write the posterior PDF as:   \begin{eqnarray} p(\boldsymbol {\mathrm{m}}|\mathrm{{{\bf d}}_{obs}}) \!\propto p(\boldsymbol {\mathrm{m}})\exp \left[\!-\!\frac{1}{2}(\mathrm{{{\bf d}}_{obs}}\!-\!\boldsymbol {\mathrm{G}}\boldsymbol {\mathrm{m}})^\mathrm{T}{\bf C_{\chi }}^{\!-1}(\mathrm{{{\bf d}}_{obs}}\!-\!\boldsymbol {\mathrm{G}}\boldsymbol {\mathrm{m}})\right] \end{eqnarray} (1)where m is the model vector, p(m) is the prior distribution, dobs is the data vector, G is the Green functions matrix, and Cχ 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(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 Cp, which is added to the observation uncertainties matrix Cd to obtain the misfit covariance:   $$\mathbf {C_{\chi }}=\mathbf {C}_{\mathrm{d}}+\mathbf {C}_{\mathrm{p}}.$$ (2) We build Cp 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 marginal PDFs are nearly Gaussians (cf. Supporting Information text T2). Figure 4. View largeDownload slide 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). Figure 4. View largeDownload slide 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). 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 Cp 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 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.) Figure 5. View largeDownload slide 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. Figure 5. View largeDownload slide 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. Figure 6. View largeDownload slide 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. Figure 6. View largeDownload slide 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. Figure 7. View largeDownload slide 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). Figure 7. View largeDownload slide 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). 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 Mw = 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}$$ (3)where Δuik is the slip inferred in a patch of area Aik and width wk 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:   $${}\text{SSD} = 100\left(\frac{P_{k=3}-P_{k=1}}{P_{k=3}}\right).$$ (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. 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. Fault segment  Mean SSD  95 per cent conf. interval  Probability than SSD is greater or equal than...        0 per cent  25 per cent  50 per cent  Emerson and Camp Rock  2.6 per cent  −25.1 per cent to 33.5 per cent  62 per cent  3.1 per cent  0 per cent  Johnson Valley  25.4 per cent  −3.8 per cent to 57.6 per cent  94 per cent  58.2 per cent  <1 per cent  Homestead Valley  51.7 per cent  42.7 per cent to 61.9 per cent  100 per cent  97.0 per cent  67.5 per cent  All faults combined  40.9 per cent  35.2 per cent to 47.3 per cent  100 per cent  99.9 per cent  <1 per cent  All faults combined taking into account a compliant zone  29.6 per cent  14.32 per cent to 46.4 per cent  99.6 per cent  75.8 per cent  <1 per cent  Fault segment  Mean SSD  95 per cent conf. interval  Probability than SSD is greater or equal than...        0 per cent  25 per cent  50 per cent  Emerson and Camp Rock  2.6 per cent  −25.1 per cent to 33.5 per cent  62 per cent  3.1 per cent  0 per cent  Johnson Valley  25.4 per cent  −3.8 per cent to 57.6 per cent  94 per cent  58.2 per cent  <1 per cent  Homestead Valley  51.7 per cent  42.7 per cent to 61.9 per cent  100 per cent  97.0 per cent  67.5 per cent  All faults combined  40.9 per cent  35.2 per cent to 47.3 per cent  100 per cent  99.9 per cent  <1 per cent  All faults combined taking into account a compliant zone  29.6 per cent  14.32 per cent to 46.4 per cent  99.6 per cent  75.8 per cent  <1 per cent  View Large 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. 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 mest (Segall & Harris 1987; Ortega 2013):   $$\boldsymbol {\mathrm{m_{est}}}(\epsilon ) = \left(\boldsymbol {\mathrm{G}}^\mathrm{T}{\bf C_{\chi }}^{-1}\boldsymbol {\mathrm{G}}+(\epsilon \nabla ^2)^2\right)^{-1}\boldsymbol {\mathrm{G}}^\mathrm{T}{\bf C_{\chi }}\mathrm{{{\bf d}}_{obs}}$$ (5) where ∇2 is the Laplacian operator defined on fault slip surface coordinates, and ε is the damping parameter. As shown in Supporting Information Figs S15(c)–(h), the larger the damping ε, the smoother the solution. Supporting Information Fig. S15(a) shows that shallow slip deficit values vary widely as a function of ε, from 13 per cent to 57 per cent. Unsurprisingly, models with little regularization (e.g. ε ∼ 0.1) are quite consistent with our Bayesian solution, including in terms of shallow slip deficit. The choice of ε 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 Fig. S15b). Such a strong dependence on ε 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 8. View largeDownload slide 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. Figure 8. View largeDownload slide 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. The fault zone is often regarded as a highly deformed core surrounded 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 secondary faults reducing the elastic strain released on the main rupture 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 vertical 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 effective shear modulus contrast μ1/μ0 (where μ1 is the shear modulus of the fault zone while μ0 is the modulus of the surrounding crust). The compliant zone half-width and shear modulus ratio being typical Jeffreys parameters (Tarantola 2005), they are sampled in the 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. Figure 9. View largeDownload slide 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. Figure 9. View largeDownload slide 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. 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 ν(x) of seismicity where x is the fault normal distance, and assume a power law decay of the form   $$\nu (x) = \nu _0\left(1+\frac{x^2}{d^2}\right)^{-\gamma /2}$$ (6)where ν0 is the aftershock density at x = 0, d is the damage zone half-width and γ is the asymptotic roll-off of the seismicity away from the fault. Using a Metropolis inversion scheme, we then sample ν0, d, and γ given the seismicity density, ν(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 ∼ 800 m is larger than what is inferred from optical images (d ∼ 570 m), an inversion with a fixed d = 570 m also explains the data reasonably well (cf. Fig. 10b). Figure 10. View largeDownload slide 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. Figure 10. View largeDownload slide 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. 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 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. Figure 11. View largeDownload slide 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 11. View largeDownload slide 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 12. View largeDownload slide 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). Figure 12. View largeDownload slide 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). 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 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 inversion 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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Ben-Zion Y., Sammis C. G., 2003. Characterization of fault zones, in Seismic Motion, Lithospheric Structures, Earthquake and Volcanic Sources: The Keiiti 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. Google Scholar CrossRef Search ADS   Chester F.M., Evans J.P., Biegel R.L., 1993. Internal structure and weakening mechanisms of the San Andreas fault, J. geophys. Res. , 98( B1), 771– 786. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Fialko Y., 2004b. Probing the mechanical properties of seismically active crust with space geodesy: study of the coseismic deformation due to the 1992 Mw7.3 Landers (southern California) earthquake, J. geophys. Res. , 109( B3), B03307, doi:10.1029/2003JB002756. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  Lohman R.B., Simons M., 2005. Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling, Geochem. Geophys. Geosyst. , 6( 1), Q01007, doi:10.1029/2004GC000841. Google Scholar CrossRef Search ADS   Mai P.M. et al.  , 2016. The earthquake-source inversion validation (SIV) project, Seismol. Res. Lett. , 87( 3), 690– 708. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Murray M., Savage J., Lisowski M., Gross W., 1993. Coseismic displacements: 1992 Landers, California, earthquake, Geophys. Res. Lett. , 20( 7), 623– 626. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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, Geophys. J. Int. , 155( 3), 1021– 1041. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Rosen P.A., Hensley S., Peltzer G., Simons M., 2004. Updated repeat orbit interferometry package released, EOS, Trans. Am. geophys. Un. , 85( 5), 47. Google Scholar CrossRef Search ADS   Savage J.C., Svarc J.L., 1997. Postseismic deformation associated with the 1992 Mw=7.3 Landers earthquake, southern California, J. geophys. Res. , 102( B4), 7565– 7577. Google Scholar CrossRef Search ADS   Segall P., 2010. Earthquake and Volcano Deformation , Princeton University Press. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Segall P., Pollard D., 1980. Mechanics of discontinuous faults, J. geophys. Res. , 85( B8), 4337– 4350. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   Tarantola A., 2005. Inverse Problem Theory and Methods for Model Parameter Estimation , SIAM. Google Scholar CrossRef Search ADS   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 , 227, p. 255, eds Thomas M. Y., Mitchell T. M., Bhat H. S., John Wiley & Sons, Inc. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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 = CmGT(GCmGT + Cχ) − 1G · Cm 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 of which the two PDFs overlap is given on the top-left corner of each histograms. 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 φ and epicentral distance Δ. 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 φ and epicentral distance Δ. 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 ε (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-σ 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 aftershock density profile parameters ν0, d and γ. 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. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

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

, Volume 212 (2) – Feb 1, 2018
14 pages

/lp/ou_press/revisiting-the-1992-landers-earthquake-a-bayesian-exploration-of-co-9YigGvBBDu
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx455
Publisher site
See Article on Publisher Site

### Abstract

Abstract Existing models for the distribution of subsurface fault slip associated with the 1992 Landers, CA, earthquake (Mw = 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. 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 Mw = 7.3 Landers earthquake (e.g. Murray et al. 1993; Cohee & Beroza 1994; Freymueller 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 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 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). Figure 1. View largeDownload slide 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. Figure 1. View largeDownload slide 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. Figure 2. View largeDownload slide 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 (Mw = 6.5) hypocentres from the Southern California Earthquake Center are indicated with a green and an orange star, respectively. Figure 2. View largeDownload slide 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 (Mw = 6.5) hypocentres from the Southern California Earthquake Center are indicated with a green and an orange star, respectively. 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 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 × 10 km2 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 Mw = 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). Figure 3. View largeDownload slide (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 conjugated 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’ Figure 3. View largeDownload slide (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 conjugated 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’ 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 ≥ 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 write the posterior PDF as:   \begin{eqnarray} p(\boldsymbol {\mathrm{m}}|\mathrm{{{\bf d}}_{obs}}) \!\propto p(\boldsymbol {\mathrm{m}})\exp \left[\!-\!\frac{1}{2}(\mathrm{{{\bf d}}_{obs}}\!-\!\boldsymbol {\mathrm{G}}\boldsymbol {\mathrm{m}})^\mathrm{T}{\bf C_{\chi }}^{\!-1}(\mathrm{{{\bf d}}_{obs}}\!-\!\boldsymbol {\mathrm{G}}\boldsymbol {\mathrm{m}})\right] \end{eqnarray} (1)where m is the model vector, p(m) is the prior distribution, dobs is the data vector, G is the Green functions matrix, and Cχ 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(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 Cp, which is added to the observation uncertainties matrix Cd to obtain the misfit covariance:   $$\mathbf {C_{\chi }}=\mathbf {C}_{\mathrm{d}}+\mathbf {C}_{\mathrm{p}}.$$ (2) We build Cp 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 marginal PDFs are nearly Gaussians (cf. Supporting Information text T2). Figure 4. View largeDownload slide 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). Figure 4. View largeDownload slide 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). 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 Cp 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 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.) Figure 5. View largeDownload slide 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. Figure 5. View largeDownload slide 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. Figure 6. View largeDownload slide 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. Figure 6. View largeDownload slide 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. Figure 7. View largeDownload slide 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). Figure 7. View largeDownload slide 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). 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 Mw = 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}$$ (3)where Δuik is the slip inferred in a patch of area Aik and width wk 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:   $${}\text{SSD} = 100\left(\frac{P_{k=3}-P_{k=1}}{P_{k=3}}\right).$$ (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. 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. Fault segment  Mean SSD  95 per cent conf. interval  Probability than SSD is greater or equal than...        0 per cent  25 per cent  50 per cent  Emerson and Camp Rock  2.6 per cent  −25.1 per cent to 33.5 per cent  62 per cent  3.1 per cent  0 per cent  Johnson Valley  25.4 per cent  −3.8 per cent to 57.6 per cent  94 per cent  58.2 per cent  <1 per cent  Homestead Valley  51.7 per cent  42.7 per cent to 61.9 per cent  100 per cent  97.0 per cent  67.5 per cent  All faults combined  40.9 per cent  35.2 per cent to 47.3 per cent  100 per cent  99.9 per cent  <1 per cent  All faults combined taking into account a compliant zone  29.6 per cent  14.32 per cent to 46.4 per cent  99.6 per cent  75.8 per cent  <1 per cent  Fault segment  Mean SSD  95 per cent conf. interval  Probability than SSD is greater or equal than...        0 per cent  25 per cent  50 per cent  Emerson and Camp Rock  2.6 per cent  −25.1 per cent to 33.5 per cent  62 per cent  3.1 per cent  0 per cent  Johnson Valley  25.4 per cent  −3.8 per cent to 57.6 per cent  94 per cent  58.2 per cent  <1 per cent  Homestead Valley  51.7 per cent  42.7 per cent to 61.9 per cent  100 per cent  97.0 per cent  67.5 per cent  All faults combined  40.9 per cent  35.2 per cent to 47.3 per cent  100 per cent  99.9 per cent  <1 per cent  All faults combined taking into account a compliant zone  29.6 per cent  14.32 per cent to 46.4 per cent  99.6 per cent  75.8 per cent  <1 per cent  View Large 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. 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 mest (Segall & Harris 1987; Ortega 2013):   $$\boldsymbol {\mathrm{m_{est}}}(\epsilon ) = \left(\boldsymbol {\mathrm{G}}^\mathrm{T}{\bf C_{\chi }}^{-1}\boldsymbol {\mathrm{G}}+(\epsilon \nabla ^2)^2\right)^{-1}\boldsymbol {\mathrm{G}}^\mathrm{T}{\bf C_{\chi }}\mathrm{{{\bf d}}_{obs}}$$ (5) where ∇2 is the Laplacian operator defined on fault slip surface coordinates, and ε is the damping parameter. As shown in Supporting Information Figs S15(c)–(h), the larger the damping ε, the smoother the solution. Supporting Information Fig. S15(a) shows that shallow slip deficit values vary widely as a function of ε, from 13 per cent to 57 per cent. Unsurprisingly, models with little regularization (e.g. ε ∼ 0.1) are quite consistent with our Bayesian solution, including in terms of shallow slip deficit. The choice of ε 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 Fig. S15b). Such a strong dependence on ε 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 8. View largeDownload slide 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. Figure 8. View largeDownload slide 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. The fault zone is often regarded as a highly deformed core surrounded 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 secondary faults reducing the elastic strain released on the main rupture 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 vertical 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 effective shear modulus contrast μ1/μ0 (where μ1 is the shear modulus of the fault zone while μ0 is the modulus of the surrounding crust). The compliant zone half-width and shear modulus ratio being typical Jeffreys parameters (Tarantola 2005), they are sampled in the 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. Figure 9. View largeDownload slide 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. Figure 9. View largeDownload slide 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. 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 ν(x) of seismicity where x is the fault normal distance, and assume a power law decay of the form   $$\nu (x) = \nu _0\left(1+\frac{x^2}{d^2}\right)^{-\gamma /2}$$ (6)where ν0 is the aftershock density at x = 0, d is the damage zone half-width and γ is the asymptotic roll-off of the seismicity away from the fault. Using a Metropolis inversion scheme, we then sample ν0, d, and γ given the seismicity density, ν(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 ∼ 800 m is larger than what is inferred from optical images (d ∼ 570 m), an inversion with a fixed d = 570 m also explains the data reasonably well (cf. Fig. 10b). Figure 10. View largeDownload slide 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. Figure 10. View largeDownload slide 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. 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 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. Figure 11. View largeDownload slide 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 11. View largeDownload slide 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 12. View largeDownload slide 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). Figure 12. View largeDownload slide 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). 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 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 inversion 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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Ben-Zion Y., Sammis C. G., 2003. Characterization of fault zones, in Seismic Motion, Lithospheric Structures, Earthquake and Volcanic Sources: The Keiiti 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. Google Scholar CrossRef Search ADS   Chester F.M., Evans J.P., Biegel R.L., 1993. Internal structure and weakening mechanisms of the San Andreas fault, J. geophys. Res. , 98( B1), 771– 786. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Fialko Y., 2004b. Probing the mechanical properties of seismically active crust with space geodesy: study of the coseismic deformation due to the 1992 Mw7.3 Landers (southern California) earthquake, J. geophys. Res. , 109( B3), B03307, doi:10.1029/2003JB002756. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  Lohman R.B., Simons M., 2005. Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling, Geochem. Geophys. Geosyst. , 6( 1), Q01007, doi:10.1029/2004GC000841. Google Scholar CrossRef Search ADS   Mai P.M. et al.  , 2016. The earthquake-source inversion validation (SIV) project, Seismol. Res. Lett. , 87( 3), 690– 708. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Murray M., Savage J., Lisowski M., Gross W., 1993. Coseismic displacements: 1992 Landers, California, earthquake, Geophys. Res. Lett. , 20( 7), 623– 626. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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, Geophys. J. Int. , 155( 3), 1021– 1041. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Rosen P.A., Hensley S., Peltzer G., Simons M., 2004. Updated repeat orbit interferometry package released, EOS, Trans. Am. geophys. Un. , 85( 5), 47. Google Scholar CrossRef Search ADS   Savage J.C., Svarc J.L., 1997. Postseismic deformation associated with the 1992 Mw=7.3 Landers earthquake, southern California, J. geophys. Res. , 102( B4), 7565– 7577. Google Scholar CrossRef Search ADS   Segall P., 2010. Earthquake and Volcano Deformation , Princeton University Press. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Segall P., Pollard D., 1980. Mechanics of discontinuous faults, J. geophys. Res. , 85( B8), 4337– 4350. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   Tarantola A., 2005. Inverse Problem Theory and Methods for Model Parameter Estimation , SIAM. Google Scholar CrossRef Search ADS   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 , 227, p. 255, eds Thomas M. Y., Mitchell T. M., Bhat H. S., John Wiley & Sons, Inc. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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 = CmGT(GCmGT + Cχ) − 1G · Cm 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 of which the two PDFs overlap is given on the top-left corner of each histograms. 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 φ and epicentral distance Δ. 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 φ and epicentral distance Δ. 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 ε (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-σ 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 aftershock density profile parameters ν0, d and γ. 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. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Feb 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations

Abstract access only

18 million full-text articles

Print

20 pages / month

PDF Discount

20% off