# Accounting for uncertain fault geometry in earthquake source inversions – I: theory and simplified application

Théa Ragon, Anthony Sladen, Mark Simons

► To cite this version:

Théa Ragon, Anthony Sladen, Mark Simons. Accounting for uncertain fault geometry in earthquake source inversions – I: theory and simplified application. Geophysical Journal International, 2018, 214 (2), pp.1174 - 1190. 10.1093/gji/ggy187. hal-01856119

HAL Id: hal-01856119

https://hal.science/hal-01856119

Submitted on 15 Oct 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Accounting for uncertain fault geometry in earthquake source inversions – I: theory and simplified application

Théa Ragon,<fnref n="1" /> Anthony Sladen<fnref n="1" /> and Mark Simons<fnref n="2" />

Accepted 2018 May 10. Received 2018 February 6; in original form 2017 December 6

## SUMMARY

The ill-posed nature of earthquake source estimation derives from several factors including the quality and quantity of available observations and the fidelity of our forward theory. Observational errors are usually accounted for in the inversion process. Epistemic errors, which stem from our simplified description of the forward problem, are rarely dealt with despite their potential to bias the estimate of a source model. In this study, we explore the impact of uncertainties related to the choice of a fault geometry in source inversion problems. The geometry of a fault structure is generally reduced to a set of parameters, such as position, strike and dip, for one or a few planar fault segments. While some of these parameters can be solved for, more often they are fixed to an uncertain value. We propose a practical framework to address this limitation by following a previously implemented method exploring the impact of uncertainties on the elastic properties of our models. We develop a sensitivity analysis to small perturbations of fault dip and position. The uncertainties of our fixed fault geometry are included in the inverse problem under the formulation of the misfit covariance matrix that combines both prediction and observation uncertainties. We validate this approach with the simplified case of a fault that extends infinitely along strike, using both Bayesian and optimization formulations of a static slip inversion. If epistemic errors are ignored, predictions are overconfident in the data and slip parameters are not reliably estimated. In contrast, inclusion of uncertainties in fault geometry allows us to infer a robust posterior slip model. Epistemic uncertainties can be many orders of magnitude larger than observational errors for great earthquakes ($M_w > 8$). Not accounting for uncertainties in fault geometry may partly explain observed shallow slip deficits for continental earthquakes. Similarly, ignoring the impact of epistemic errors can also bias estimates of near-surface slip and predictions of tsunamis induced by megathrust earthquakes.

**Key words:** Inverse theory; Probability distributions; Earthquake source observations.

## 1 INTRODUCTION

Imaging earthquake sources is a crucial step towards developing a better understanding of the physics of earthquake rupture. Such images allow us to improve our models of regional tectonics and potentially refine our assessment of seismic and tsunami hazards. Inferred source images, thus, need to be as accurate as possible to ensure a robust interpretation and should attempt to describe the limits of our understanding as well as possible.

Published models of the distribution of coseismic subsurface fault slip (a.k.a. slip models) for a given earthquake can be very discrepant, even models for well-recorded events such as the 1999 $M_w$ 7.6 Izmit earthquake (e.g. Duputel et al. 2014; Mai et al. 2016), the 2009 $M_w$ 6.3 I’Aquila earthquake (e.g. Yano et al. 2014; Volpe

et al. 2015) or the $M_w$ 9.0, 2011, Tohoku-Oki earthquake (we refer the reader to Razafindrakoto et al. 2015; Lay 2017, for an overview of selected published models). This variability relates to the inherently ill-posed nature of the inverse problem. The non-uniqueness of the inverted source model is affected by numerous factors, including data quality, quantity and processing (filtering, sampling, etc.) that controls the level of information available to reduce the uniqueness of the problem (e.g. Delouis et al. 2002; Lohman & Simons 2005; Sladen et al. 2010). Other factors affecting estimates of the source are related to our simplified description of the physics of the problem, which is inherently more difficult to quantify. Source parametrization and the inversion method will both affect our estimates of source characteristics. In particular, any simplification including the possible linearization of the inverse problem, the level

<footnote n="1">¹Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, IRD, Géoazur, 250 rue Albert Einstein, 06560 Valbonne, France. E-mail: ragon@geoazur.unice.fr</footnote>

<footnote n="2">²Seismological Laboratory, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA</footnote>
of complexity in the discretization of the source, smoothing and other types of regularization will impact the solutions of the inverse problem (e.g. Du et al. 1992; Beresnev 2003). Our imperfect forward problem can produce systematic errors and biases. Indeed, Earth structure (i.e. rheology or elastic properties) and fault geometry are always uncertain despite the fact that they are frequently fixed a priori in fault slip models. The impact of this prior choice is rarely investigated, yet, a small change can affect the Green's functions, which will, in turn, distort the link between the data space and the parameters and lead to potentially very different source models (e.g. Simons et al. 2002; Beresnev 2003; Hartzell et al. 2007; Razafindrakoto & Mai 2014; Gallovic et al. 2015; Diao et al. 2016).

While the quantity and quality of observations is ever-increasing, numerous procedures have been developed to enhance the reliability of inferred models or at least to investigate their variability. Validation and comparison procedures (e.g. Tichelaar & Ruff 1989; Gallovič & Ampuero 2015; Mai et al. 2016), Bayesian formulations (e.g. Minson et al. 2014) or posterior analysis (e.g. Piatanesi et al. 2007; Monelli & Mai 2008) allow for exploration of uncertainties. The inclusion of observation uncertainties with covariance matrices (e.g. Yagi & Fukahata 2008; Sudhaus & Jónsson 2009; Duputel et al. 2012) improves the estimates of model-parameter uncertainties. The most influential parameters on the forward problem are the assumed Earth structure and fault geometry (Das & Suhadolc 1996; Beresnev 2003). The impact of Earth structure uncertainty has previously been investigated by Duputel et al. (2014) and applied in Jolivet et al. (2014); Duputel et al. (2015); Jolivet et al. (2015); Gombert et al. (2017). The statistical impact of uncertainty of parameters describing fault geometry remains largely unexplored and is the central focus of this study.

Field studies as well as modelling analyses have documented faults complexity, with multiple segments that can only be well observed when they reach the surface (e.g. Segall & Pollard 1980; Okubo & Aki 1987; Peacock 1991; Walsh et al. 2003; Manighetti et al. 2015). The complexity of the surface trace may be similar at greater depths, as suggested by some observations (e.g. Aochi & Fukuyama 2002; Aochi & Madariaga 2003; Jolivet et al. 2014; Huang et al. 2017). Faults are generally non-planar, and the roughness of rupture surfaces can be observed at different scales (e.g. Power et al. 1987; Aochi & Madariaga 2003; Candela et al. 2012; Perrin et al. 2016; Zielke et al. 2017). Bends, jogs, relay zones, curvature and other complexities of fault structures clearly influence, or even drive, earthquake rupture processes (King & Náblek 1985; Zhang et al. 1991; Bouchon et al. 1998; Lee et al. 2006; Moreno et al. 2009; Wei et al. 2011; Moreno et al. 2012; Bletery et al. 2016b; Pizzi et al. 2017). Developing a realistic fault geometry is often a challenging component of finite source modelling (e.g. Aochi & Madariaga 2003).

Yet, even when incorporating increasingly complex structures in inversions, fault geometry will remain uncertain. Geometries determined from surface rupture, centroid moment tensor solutions, previous earthquakes, aftershocks distributions or tomography, are all inaccurate to varying extents. This uncertainty can result in a discrepancy between published models for a given earthquake. Even for particularly well-recorded and studied events, fault geometry parameters are poorly constrained. For instance, the published dip and strike for planar faults parameters for the 2009 l'Aquila event vary by 15° and 10°, respectively (Lavecchia et al. 2012). We frequently adopt fault models that are planar and thus smooth, sometimes multisegmented, with segments of fixed dip and fixed strike. Inaccurate fault dip and/or position distorts the Green's functions which will bias the inversion algorithm towards incorrect source models. In

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic view of the 2-D simplified toy model we use to explore the impact of including uncertainties in fault geometry on source estimation problems. The assumed geometry does not vary along strike. We investigate first the effect of accounting for inaccuracies on fault dip. Then, we analyze the impact of accounting for uncertainties in fault position.</figcaption></figure>

some cases, data allow identification of additional complexities in the fault geometry, such as curvature or segmentation.

In this study, we investigate the impact of uncertainty in fault geometry on inverted slip models. We propose a formalism based on small perturbation theory which limits its application to events for which the gross characteristics of the fault geometry are known. Indeed, we choose not to solve for the fault geometry parameters as it is computationally expensive and can be biased by uncertainties, as discussed further. Uncertainties on the *a priori* fixed fault parameters are specified using a covariance matrix of prediction errors following the approach developed by Duputel et al. (2014) for Earth structure. We present the corresponding theoretical approach and then describe implementations for both optimization and Bayesian formulations of the inverse problem. Through a simplified 2-D toy model, we explore the influence of accounting for uncertainties in both the fault dip and the surface trace of the fault (Fig. 1). We choose to investigate the simple case of a 2-D infinite fault to get a better understanding of the influence of basic fault geometry parameters, such as fault dip and position. A 2-D model avoids the 3-D complexity that will perturb the impact of a variation of fault dip or position on inferred models. The investigation of a realistic finite-fault inversion in 3-D is left for a subsequent study.

## 2 PREDICTION UNCERTAINTIES DUE TO INACCURACIES IN THE FORWARD MODEL

### 2.1 Uncertainties of the forward physical theory

We aim to infer a source model, **m**, based on observations, $d_{obs}$, and Green's functions, **G**, that relate the displacement response of the Earth's surface to a unit displacement on a portion of the assumed causative fault (i.e. a subfault or a fault patch). Estimated parameters can thus be impacted by both observational and forward modelling, or epistemic, errors.

Measurement errors can be instrumental or result from the measurement method and conditions, such as atmospheric or topographic bias on Global Positioning System (GPS) and InSAR observations. A measurement error will result in a discrepancy between
real surface displacement **d** and observations **d**<sub>obs</sub>. The observational uncertainty is generally accounted for in source estimation problems by assuming that the real surface displacement **d** follows a Gaussian distribution centred on the observations **d**<sub>obs</sub>. The discrepancy between observations and real surface displacement is characterized by a data covariance matrix **C**<sub>d</sub>. The diagonal components of **C**<sub>d</sub> will reflect the variance of independent data, whereas off-diagonal terms will add knowledge on the correlation of observational errors (e.g. Lohman & Simons 2005; Sudhaus & Jónsson 2009; Duputel et al. 2015; Bletery et al. 2016a). **C**<sub>d</sub> contains all the statistical information about the observations.

Uncertainties due to imperfect forward modelling will also influence inferred source model (e.g. Duputel et al. 2014). But, contrary to observational uncertainties, such prediction or epistemic errors are often not accounted for. When studying large events, epistemic uncertainties may become larger (scaling approximately with slip amplitude) while observation uncertainties will remain constant (independent of slip amplitude). In contrast, the effect of uncertainties in the forward theory will be negligible for small events with poor data: observational errors will dominate the error budget in such cases.

The physics of the forward model is subject to numerous simplifications such as material properties (e.g. rheology), fault geometry or regional characteristics (e.g. topography). Inaccuracy in one of these parameters or an oversimplification in the forward model will corrupt the link between forward predictions **d**<sub>pred</sub> and real surface displacement **d**<sub>obs</sub>. Suppose we assume a set of parameters of the forward model Ψ<sub>prior</sub> that are the most realistic a priori. The prediction error thus manifests our imperfect knowledge of **d**<sub>pred</sub> = **G**(Ψ<sub>prior</sub>, **m**) for a particular resulting source model **m** and the assumed set of non-inverted model parameters Ψ<sub>prior</sub>. We assume that the real surface displacement **d** follows a Gaussian distribution centred on the predictions **d**<sub>pred</sub> with a prediction covariance matrix **C**<sub>p</sub>(**m**) that depends on the resulting source model.

The discrepancies between observations **d**<sub>obs</sub> and forward predictions **d**<sub>pred</sub> = **G**(Ψ<sub>prior</sub>, **m**) can thus be determined by a misfit function of the form

$$ \chi(\mathbf{m}) = \frac{1}{2} [\mathbf{d}_{\text{obs}} - \mathbf{G}(\Psi_{\text{prior}}, \mathbf{m})]^T \cdot \mathbf{C}_{\chi}^{-1} \cdot [\mathbf{d}_{\text{obs}} - \mathbf{G}(\Psi_{\text{prior}}, \mathbf{m})], \quad (1) $$

where **C**<sub>χ</sub>(**m**) is the misfit covariance matrix defined as (Tarantola 2005; Minson et al. 2013, 2014; Duputel et al. 2014)

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

The misfit function χ(**m**) is defined as a criterion to select the preferred source model and is independent of the chosen inversion method. The misfit function is minimized in an optimization approach or is part of the likelihood function of Bayesian analysis. While the formalism of prediction uncertainty modelling has already been explored for shear modulus by Duputel et al. (2014), we expand it here to include epistemic uncertainties due to imperfect knowledge of fault geometry.

## 2.2 Theoretical development to account for inaccuracies in fault geometry

Fault models are often described by a fixed dip, a fixed strike and a fixed surface position, although occasionally some of these parameters have been solved for in simultaneous inversions of fault parameters and slip distribution (e.g. Fukuda et al. 2010; Sun et al. 2011Tinti et al. 2016; Liu et al. 2017) or via sequential estimation (e.g. Huang et al. 2017). The inversion of the fault geometry

independent of slip should be interpreted carefully because of the trade-off between assumed fault morphology and slip behaviour.

Mirroring the approach described in Duputel et al. (2014), we can compute a prediction covariance matrix **C**<sub>p</sub> accounting for the uncertainty of fault geometry parameters. Suppose we choose an uncertain and presumably inaccurate set of parameters Ψ<sub>prior</sub> for given generic fault properties Ψ. We can explore the prediction uncertainty by assuming that the predictions **d**<sub>pred</sub> = **G**(Ψ, **m**) for the generic fault properties Ψ can be approximated by linearized perturbations of the predictions **G**(Ψ<sub>prior</sub>, **m**) for the fixed set of fault geometry parameters Ψ<sub>prior</sub> and the unknown source parameters **m**:

$$ \mathbf{G}(\Psi, \mathbf{m}) \approx \mathbf{G}(\Psi_{\text{prior}}, \mathbf{m}) + \mathbf{K}_{\Psi}(\Psi_{\text{prior}}, \mathbf{m}) \cdot (\Psi - \Psi_{\text{prior}}), \quad (3) $$

where the matrix **K**<sub>Ψ</sub>(Ψ<sub>prior</sub>, **m**) is the sensitivity kernel of the predictions with respect to fault geometry parameters. For most earthquakes, fault geometry is approximately known and the chosen parameters Ψ<sub>prior</sub> are reasonable approximations of the true properties Ψ. Even if the predictions **G**(Ψ, **m**) are nonlinear with respect to fault geometry parameters, the limited variability between the true Ψ and the chosen Ψ<sub>prior</sub> makes the perturbation Ψ − Ψ<sub>prior</sub> small enough that the problem may be linearized.

**C**<sub>p</sub> is calculated a priori, and as the predictions **d**<sub>pred</sub> = **G**(Ψ<sub>prior</sub>, **m**<sub>prior</sub>) depend both on the assumed fault parameters Ψ<sub>prior</sub> and on the assumed source model **m**<sub>prior</sub>, any variation in the source model **m**<sub>prior</sub> will affect **C**<sub>p</sub>(Ψ<sub>prior</sub>, **m**<sub>prior</sub>) for a given fault geometry Ψ<sub>prior</sub>. Once the source model **m**<sub>prior</sub> and fault parameters Ψ<sub>prior</sub> are selected, the source inversion problem being linear, we have **d**<sub>pred</sub> = **G**(Ψ<sub>prior</sub>) · **m**<sub>prior</sub>. Following the development of Duputel et al. (2014), the prediction matrix describing uncertainty in fault geometry parameters can be written as

$$ \mathbf{C}_p = \mathbf{K}_{\Psi} \cdot \mathbf{C}_{\Psi} \cdot \mathbf{K}_{\Psi}^T, \quad (4) $$

where

$$ \mathbf{K}_{\Psi} = \mathbf{K}_{\Psi}^{G} \cdot \mathbf{m}_{\text{prior}}, \quad (5) $$

$$ (\mathbf{K}_{\Psi}^{G})_{ijk}(\Psi_{\text{prior}}) = \frac{\partial G_{ik}}{\partial \Psi_{j}}(\Psi_{\text{prior}}), \quad (6) $$

and **C**<sub>Ψ</sub> is the standard deviation of the a priori distribution of parameters Ψ. We can then pre-compute the sensitivity kernels **K**<sub>Ψ</sub><sup>G</sup> without any dependence on the assumed model **m**<sub>prior</sub>. Note that we apply this approach to a static slip inversion. The slip parameters we solve for are thus included in the vector **m** while fixed fault geometry parameters correspond to Ψ<sub>prior</sub> in the above-mentioned formulas.

## 2.3 Implementation of C<sub>p</sub> in a source inversion problem

The prediction covariance matrix **C**<sub>p</sub> can be included in the inversion following two different approaches. In the first, **C**<sub>p</sub> is calculated a priori and included in the inversion process. The alternative approach that we use for our Bayesian formulation of the problem is to update **C**<sub>p</sub> with interim models at each step of the tempered inversion (see later discussion in Section 2.3.2).

### 2.3.1 Calculation of C<sub>p</sub> a priori

In this case, the data covariance matrix **C**<sub>d</sub> is replaced by the misfit covariance matrix **C**<sub>χ</sub> = **C**<sub>d</sub> + **C**<sub>p</sub>. Once the sensitivity kernels **K**<sub>Ψ</sub><sup>G</sup> are pre-computed, the **C**<sub>p</sub> value relies upon the choice of a priori uncertainty on chosen parameters **C**<sub>Ψ</sub> and of the assumed model
m prior . Practically, m prior can be selected among a range of a priori source models, such as the centroid moment tensor solution or a resulting model solved without including C p , that is, with C χ = C d .

2.3.2 Implementation of C$_{p}$ in our Bayesian approach

In this study, most of the presented results are solved through a Bayesian sampling approach that relies on the AlTar package, a rewrite of the code CATMIP (Minson et al. 2013). AlTar combines the Metropolis algorithm with a tempering process to realize an iterative sampling of the solution space of the source models. A large number of samples are tested in parallel at each transitional step. Additionally, a resampling is performed at the end of each tempering step to replace less probable models. The probability of each sample to be selected depends on its likelihood. The ability of each model parameter to solve the source problem is evaluated through

$$
f(\mathbf{m}, \beta_i) \propto p(\mathbf{m}) \cdot \exp[-\beta_i \cdot \chi(\mathbf{m})], \quad (7)
$$

with **m** being the current sample and *p*(**m**) the prior information on this sample, **χ**(**m**) being the misfit function, *i* corresponding to each iteration and the tempering parameter **β** evolving dynamically from 0 to 1 to improve the efficiency of the parameter space exploration (Minson *et al.* 2013). Practically, this implementation of AlTar allows us to recompute **C**<sub>p</sub> at each tempering step by choosing the mean of tested samples as the assumed model **m**<sub>prior</sub>. The prediction covariance is thus updated following the evolution of the target source model. We choose a uniform model of consistent magnitude as the initial model. In the case where our target model is uniform, we choose a uniform zero slip model as the initial model.

3 A SIMPLIFIED 2-D MODEL:
ACCOUNTING FOR INACCURACIES IN
THE FAULT DIP

To highlight the effect of uncertainty in fault dip, we consider a synthetic 2-D case (Fig. 1). We assume a fault with a 20 km downdip extent that extends infinitely along strike. The fault is discretized along dip into subfaults of 1 km. The true dip of the fault is 55° for our dip-slip case, and of 80° for the strike-slip case. For both cases, we compute the corresponding synthetic observations using the expressions of surface displacement in a homogeneous half-space, modified after Segall (2010) to be dip-dependent (see equations in the Appendix). These synthetic observations are computed for 100 data points at the surface, spaced every kilometre if within a distance of 30 km from the fault, and every 2 km if at greater distances (i.e. from 30 to 70 km away from the fault trace). This uneven distribution of data points is driven by the rapid decay of the surface deformation with distance and the better resolution of fault slip by observations close to the fault (a few fault widths in this case, e.g. Lohman & Simons 2005). An uncorrelated Gaussian noise of 7 mm is added to the observations to simulate measurement errors. The choice of a 7 mm noise amplitude derives from values of uncertainties associated with coseismic GPS offsets of continental earthquakes, that are in average of 2–5 per cent the amount of displacement (for instance, for the 2009 Mw 6.3 l’Aquila, 2015 Mw 7.8 Pedernales or 2015 Mw 7.8 Gorkha earthquakes). The effect of spatial correlation of the noise is explored in Section 2.3.

Using these 100 synthetic data points, we then estimate the depth
distribution of slip with AlTar, still assuming a homogeneous elas-
tic half-space, but with an incorrect fault geometry: the dip is 5°

different from the true value, that is 50° instead of 55° for the dip-
slip scenario and 75° instead of 80° for the strike-slip scenario. We
account for uncertainty of the fault dip following the formulation
of C<sub>p</sub> defined in Section 3.2.

3.1 Calculation of C<sub>p</sub>

Here, we want to account for the uncertainty in just the fault dip. To calculate C<sub>p</sub> following eq. (4), we need to pre-compute the sensitivity kernel K<sub>dip</sub> of Green's functions with respect to the dip parameter, and to choose a covariance C<sub>dip</sub> that will describe the uncertainty of the fault dip.

Suppose the assumed, but incorrect dip parameter is d prior . For each of the ten dip integer values d i on the range [ d prior − 5 ◦ : d prior + 5 ◦ ], we compute the matrix G d i of Green’s functions. Regarding the small perturbations of the Green’s functions on this restricted dip interval, we linearize the evolution of each element of matrix G d i . Given

$$
\mathbf{K}_{\text{dip}} = \frac{\partial \mathbf{G}_{\text{d}}}{\partial \mathbf{d}}, \qquad (8)
$$

each element of matrix **K**<sub>dip</sub> is thus equal to the slope of the linear regression of **G**<sub>d</sub>. The evolution and linear regression of **G**<sub>d</sub> is shown in Fig. 2 for six elements of the matrix. For most of **G**<sub>d</sub> components, the linearization well approximates the evolution of the Green’s functions (Figs 2b–d,f,g). When there is a critical point on this dip interval (Fig. 2e), the linearization is not a good approximation of **G**<sub>d</sub> but rather an average of the variation of the Green’s functions.

We choose a covariance C dip that corresponds to a standard de-
viation of 5° from the assumed value (50° or 75° in our case). C p is
recalculated at each step of our inversion. The posterior C p (i.e. the
C p calculated at the last tempering iteration) is represented in Fig. 3.
Epistemic uncertainties are generally higher near the fault (Fig. 3),
and roughly proportional to the amplitude of observations as pro-
posed by Minson et al. (2013). C p values also follow the variation
with distance to the fault of the difference between dip-dependent
surface displacements (Fig. 4). For instance, C p for the vertical dis-
placement reproduces this variation: left of the fault (upper left-hand
part of Fig. 3b), the distance between vertical displacement from a
50° and a 55° dipping fault increases towards the fault (Fig. 4b).
Right of the fault (lower right-hand part of Fig. 3b), the distance
between the 50° and 55° curves is low near the fault, tends to zero
around 10 km from the fault, increases around 20 km and decreases
to zero in far field (Fig. 4b). C p values will vary with the assumed
dip, and this variation will be conditioned by the discrepancy be-
tween assumed and true dips, but also by their own value (80° and
85° or 30° and 35°). C p values are also 10 times larger for predic-
tions in the strike-perpendicular direction than in other directions
(Fig. 3c), as surface displacement is dominant in this direction for a
55° dipping fault (Fig. 4). The difference between surface displace-
ments due to faults with different dips (assumed or true one) thus
explains the scaling and asymmetries of C p (between left and right
side of the fault).

3.2 Influence of prediction uncertainty on the inferred slip distribution

We now evaluate the effect of accounting for prediction uncertainty
in the inversion by comparing slip models inferred with or without
C<sub>p</sub>. The synthetic data are computed for a uniform slip model of
1 m from 20 km downdip to the surface (grey vertical bars in Fig. 5).
Vertical surface displacement for 1 m of dip slip, fault dip of 50°

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

Evolution of the Green's functions with the fault dip, between 45° and 55°

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Variation of six components of the matrix of the Green's functions for an infinitely long fault along strike with a dip varying between 45° and 55°. The Green's function components are vertical surface displacement (a) components for a distance of -3, 5 and 48 km to the fault [respectively (b) and (e), (c) and (f) and (d) and (g)]. Components are shown for the third downdip parameter [shallow part of the fault, plots (b), (c) and (d)] or the 13th parameter [deeper part of the fault, plots (e), (f) and (g)]. Green's functions variation is shown in blue, whereas linearized Green's functions are shown in orange. The slope of linearized Green's functions is taken as component of the sensitivity kernel $K_{\text{dip}}$.</figcaption></figure>

We use a uniform prior distribution $p(\mathbf{m}) = U(-0.5 \text{ m}, 5 \text{ m})$ (uniform implies that all values are considered equally likely with no *a priori* knowledge) and include 7 mm of observational uncertainty in $\mathbf{C}_d$.

The posterior mean model as well as the marginal probability density functions (PDFs) for each fault patch are presented in Fig. 5. The inversion results have been solved neglecting $\mathbf{C}_p$ (Figs 5a and c) or including $\mathbf{C}_p$ (Figs 5b and d), respectively, for the dip-slip and the strike-slip scenarios. PDFs reproduce the number of samples in the last iteration that account for a particular value of a model parameter. A narrower posterior PDF means that the parameter is more accurately estimated as a large set of samples agree on a very limited range of values.

The introduction of $\mathbf{C}_p$ significantly improves the posterior model estimations. Both in the dip-slip and strike-slip cases, the inferred model accounting for $\mathbf{C}_p$ is a far better approximation of the target model than the one without $\mathbf{C}_p$, with this improvement found at all depths. For the shallow slip distribution, posterior PDFs are very narrow but yet their mean is offset from the target model if we ignore epistemic uncertainties (Figs 5a and c). Accounting for $\mathbf{C}_p$, posterior PDFs for shallow patches are broader but their mean

matches the target slip model. $\mathbf{C}_p$ has an even stronger influence on the deeper parts of the fault. When $\mathbf{C}_p$ is not accounted for, the error of the mean model distribution is about 1 m for dip-slip and up to 70 cm for strike-slip. For the dip slip case, the deepest subfault totally excludes the target value from the interval of uncertainty (Fig. 5a). Again, including $\mathbf{C}_p$ gives posterior distributions that are broader (see patch 14–15 for dip-slip in Fig. 5b) but, more interestingly, centred around a mean model agreeing almost exactly with the target model. Accounting for $\mathbf{C}_p$, offsets between target and inferred models are less than 10 cm, and generally range less than 5 cm.

In Fig. 6, we compare the synthetic observations, the surface displacements induced by the target model and the predictions of the mean model with or without including $\mathbf{C}_p$ (the last three with an incorrect dip). As expected, the predictions are slightly different from the observations because they were not produced assuming the same fault geometry (dip = 50° and 75° used for the inversion, and dip = 55° and 80° used to create the data). The posterior uncertainty on predictions is large as following the uncertainty on the parameters that generally reaches 70 cm. Interestingly, we note that the RMS is smaller for the worst slip
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Matrix $C_p$ quantifying the influence of fault dip uncertainty on predictions in the case of a strike-slip fault (a) or dip-slip fault: (b) and (c) for vertical and strike perpendicular displacements, respectively. Near centre $C_p$ values correspond to near-fault data points, whereas near edge values (upper left and bottom right) correspond to far-field predictions (respectively in the hanging wall and foot wall). $C_p$ is a covariance matrix of the predictions, and thus the variance of predictions for a variation of dip corresponds to diagonal values. The assumed fault is 20 km wide, with a 55° dip and is infinite along strike. These figures are computed using a uniform $m_{prior}$ slip value of 1 m.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Influence of dip variation on surface displacement in function of the distance to the fault, for an infinite fault along strike with a 1 m uniform slip from 20 km downdip to the surface.</figcaption></figure>

model: the one with no $C_p$ has a RMS of 5.5 mm for horizontal displacement and 3.5 mm for vertical, while the almost perfectly solved slip distribution with $C_p$ has a RMS of 6.4 mm for horizontal and 4.4 mm for vertical displacement. A similar behaviour is described in Duputel *et al.* (2014) and is consistent with allowing for correlated misfits. The misfit function can be rewritten (from eq. 1) $\chi(\mathbf{m}) = \frac{1}{2}(\mathbf{R}\mathbf{d}_{\text{obs}} - \mathbf{R}\mathbf{d}_{\text{pred}})^T \cdot (\mathbf{R}\mathbf{d}_{\text{obs}} - \mathbf{R}\mathbf{d}_{\text{pred}})$, where $\mathbf{R}$ is the Cholesky decomposition of $C_\chi^{-1} = \mathbf{R}^T \cdot \mathbf{R}$. The standardized data and predictions, respectively $\mathbf{R}\mathbf{d}_{\text{obs}}$ and $\mathbf{R}\mathbf{d}_{\text{pred}}$, thus differ in

function of which $C_\chi$ is used (including $C_p$ or not). When including $C_p$, the offset between standardized data and predictions from the target model is two times larger than when no $C_p$ is used (see Fig. S1 of the Supporting Information for a comparison of standardized data and predictions with or without $C_p$). In other words, the inclusion of $C_p$ increases the impact of the fault dip on surface displacements, and makes the fit of observations harder to reach.

Thus, the model with no $C_p$ fits the observations better but is not in agreement with the target model. In contrast, the model accounting for $C_p$ does not overfit components of the misfit associated with potential model deficiencies. Similarly, when we assume a strictly positive prior ($p(\mathbf{m}) = U(0 \text{ m}, 5 \text{ m})$ instead of $p(\mathbf{m}) = U(-0.5 \text{ m}, 5 \text{ m})$), the model with $C_p$ better approaches the target model (Fig. S2, Supporting Information). These conclusions also hold when the fault dip is off by 10° from the correct fault dip (Fig. S3, Supporting Information). When assuming a fault dipping 45°, the target model is less well approached than with a 50° dipping fault but the inferred model is still a good approximation. We also explore the impact of including in an inverse problem for a non-uniform target model. As for the uniform target model case, the differences between posterior PDFs and mean model distribution with or without including are compelling (Figs S4 and S5, Supporting Information). We also investigate the impact of the assumed standard deviation. If $C_\psi$ is undervalued, inferred model is farther from the target model than with a right standard deviation (Fig. S6, Supporting Information). Additionally, if the observational errors are underestimated, assuming undervalued standard deviation can lead to an inferred model at odd with the target model (Fig. S7, Supporting Information). And even if the standard deviation $C_\psi$ is not undervalued, only the assumption of correct observational errors allows to infer a perfect estimation of the target model (compare Figs S6 and S7, Supporting Information). Similarly, the effect of spatially correlated noise in the observations is negligible if the noise amplitude is not underestimated (Fig. S8, Supporting Information).
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of inversion results with and without neglecting $C_p$. The uniform target slip model is shown as a grey vertical bar. The posterior PDFs are shown in (a) and (c) when the prediction uncertainty is neglected and in (b) and (d) if $C_p$ is included in the inversion, respectively, for dip-slip and strike-slip behaviour. The posterior mean model is indicated with a vertical grey dotted line. The offset between target model and posterior mean is displayed with a colour scale saturated at 50 cm.</figcaption></figure>

## 4 A SIMPLIFIED 2-D MODEL: ACCOUNTING FOR INACCURACIES ON THE FAULT SURFACE TRACE LOCATION

Given the importance of accounting for dip uncertainty, we examine the influence of uncertainties on the fault location within the same simplified approach (Fig. 1). The three main parameters allowing to describe a planar fault geometry are the fault dip, already investigated, and the fault strike and position. The variation of fault strike is inherently a 3-D problem. The impact of a priori inaccuracies on strike parameter will thus be dependent on (1) the considered earthquake and (2) the data distribution. In contrast, the variation of fault tip location can be investigated in a 2-D synthetic case.

We still consider the 2-D synthetic case of a fault with a 20 km long downdip extent that prolongs infinitely along strike, discretized along dip into subfaults of 1 km and with a true dip of 55° or 80°, respectively, for dip-slip or strike-slip cases. This time, we invert for a fault location shifted 2 km from the original one, with dip parameters set to their true value. To avoid the unrealistic case of having ground displacement measurements incompatible with the fault trace, we do not take into account the nearest faults observations (i.e. data points adjacent to the true fault, adjacent to the assumed shifted fault and in between observations). Five observation points are thus removed, between -1 and 3 km offset from the true fault position. Given the zero offset synthetic observations, we account for uncertainty of the fault location according to the formulation of $C_p$ defined in Section 2.3.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparisons between synthetic data and predictions with and without accounting for $C_p$. We assume the slip on the fault to be uniform and equal to 1 m. The location of the fault surface rupture is shown with a grey vertical line. Dip-slip behaviour (assumed fault dip of 50°) generates strike-perpendicular surface displacement (a) and vertical surface displacement (b), whereas strike-slip (assumed fault dip of 75°) behaviour produces only strike-parallel displacement (c). The data points (i.e. generated by a slip on the true fault) are shown in blue. The displacements in khaki green are produced by the 1 m uniform target model assuming an incorrect fault dip. The predictions shown in orange and green are calculated from the posterior mean model, respectively, without and with $C_p$ included. The orange and green vertical bars correspond to the $1\sigma$ posterior uncertainty on the predictions, respectively, without and with $C_p$ included.</figcaption></figure>

## 4.1 Calculation of $C_p$

We pre-compute the sensitivity kernels $K_{\text{shift}}$ following the same framework as for the dip parameter (Section 3.1). Assuming an incorrect shift of the fault $s_{\text{prior}}$ and a range of uncertainty $e_s$ of 3 km, we compute the matrices $G_{s_i}$ Green's functions for each fault offset in the range $[s_{\text{prior}} - e_s, s_{\text{prior}} + e_s]$. Each element of $K_{\text{shift}}$ is estimated by the slope of the linear regression of $G_s$. We choose an incorrect shift of the fault equal to 2 km and a covariance $C_{\text{shift}}$ corresponding to a standard deviation of 2 km.

## 4.2 Influence of prediction uncertainty on inferred slip distributions

Synthetic observations are calculated from a non-shifted fault and for a uniform slip model of 1 m as shown in Fig. 7. We use a uniform prior distribution $p(\mathbf{m}) = U(-0.5 \text{ m}, 5 \text{ m})$ and we include 7 mm of measurement errors in $C_d$. The posterior mean model as well as the PDFs for each fault patch are presented in Fig. 7. We run the inversion with and without $C_p$ (Figs 7b and d or Figs 7a and c, respectively) respectively for dip-slip and strike-slip cases.

The inclusion of uncertainty on fault position improves the estimation of the slip parameters for both the dip-slip and strike-slip scenarios. However, the improvements brought by the inclusion of $C_p$ are not the same as when the dip is uncertain. When the fault location is shifted, the shallowest subfaults (first 7 km) are the less well resolved and, if we ignore $C_p$, the mean model can be off from the target model (1 m homogenous slip) by more than 3 m (Figs 7a and c). For some shallow subfaults, the target value is totally excluded from the interval of probable parameters (Figs 7a and c). Even when introducing $C_p$, the recovery of the target model is not perfect for the shallowest subfaults for the dip-slip case, with up to 20 cm offset between target and mean models (Fig. 7b). Estimates of slip on deeper subfaults is improved when accounting for $C_p$.

We explain the poor resolution at shallow depth by the compound effect of (1) the lack of near-fault data, and (2) the greater influence of an incorrect fault position on the near-field data. These data are mainly sensitive to shallow slip. Indeed, the largest difference between data and model predictions is near the surface trace of the fault (Fig. 8). In Fig. 8, we compare the synthetic observations (calculated with a zero offset fault) and the surface displacements induced by the target model and the predictions with or without including $C_p$ (the last three with a 2 km shifted fault). As expected, the target model predictions and the observations are different as produced by faults of various position. As in the case of an uncertain dip angle, the RMS of residuals between observations and predictions are higher when $C_p$ is included (RMS of 6.1, 1.3 and 3.3 mm, respectively, for strike-perpendicular, strike-parallel and up displacement) than when $C_p$ is neglected (RMS of 5.9, 1.0 and 3.4 mm, respectively, for strike-perpendicular, strike-parallel and up displacement). Thus, the model inferred without $C_p$ better fits the observations while being unable to reach the target slip distribution. In contrast, the model including $C_p$ solves the problem while being less confident on the accuracy of the observations. The introduction of fault position uncertainty in the inverse problem allows to increase the resolution of shallow slip distribution. These conclusions also hold for a non-uniform slip model (Figs S9 and S10, Supporting Information).

The simplified toy models presented here, thus, illustrate the relevance of accounting for fault position inaccuracies in the inverse problem. The formulation of the prediction covariance matrix $C_p$ allows us to reliably estimate the posterior distribution of source model parameters inverted with an incorrect fault trace location. In particular, the inclusion of $C_p$ allows more realistic estimates of shallow slip when there is uncertainty about the exact surface extent of a fault rupture.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison of inversion results with and without incorporating $C_p$. The uniform target slip model is shown as a grey vertical bar. The posterior PDFs are shown in (a) and (c) when the prediction uncertainty is neglected and in (b) and (d) if $C_p$ is included in the inversion, respectively, for dip-slip and strike-slip behaviour. The posterior mean model is indicated with a vertical grey dotted line. The offset between target model and posterior mean is displayed with a colour scale saturated at 50 cm.</figcaption></figure>

# 5 DISCUSSION

## 5.1 Influence of uncertainties in fault geometry

When there is ambiguity in the fault geometry, which will always be the case to some extent, our synthetic tests reveal that accounting for this uncertainty allows one to match the target slip distribution with minimal errors. Not accounting for these uncertainties can induce artefacts in the inferred slip distributions. The inclusion of $C_p$ improves the robustness of the solution.

$C_p$ is not limited to probabilistic algorithms and can also be accounted for when treating the inverse problem with an optimization approach. To ensure that the proposed $C_p$ formalism works for both Bayesian and optimization approaches, we perform similar tests as

in the Bayesian framework (Section 3.1). The tests and results are detailed in Supporting Information. We compute $C_p$ *a priori* and use the result of the optimization neglecting $C_p$ as initial model (see Section 2.3.1). We perform a non-negative least squares inversion and include spatial smoothing following the approach of Tarantola (2005). The choice of regularization parameters is explained in Fig. S11 (Supporting Information). The inferred parameters and their posterior uncertainty are presented in Fig. S12 (Supporting Information). Similar to the Bayesian tests, the ability of the optimization algorithm to approach the target slip model depends on the inclusion of $C_p$. If the assumed fault geometry is incorrect and related uncertainties not accounted for (Fig. S12a, Supporting Information), the inferred model is at odds with the target model. Conversely, the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparisons between synthetic data and predictions with and without accounting for $C_p$. We assume the slip on the fault to be uniform and equal to 1 m. The correct location of the fault surface rupture is shown with a grey vertical line, whereas the assumed fault position is shown by a dotted grey vertical line. Dip-slip behaviour generates a strike-perpendicular surface displacement (a) and vertical surface displacement (b), whereas strike-slip behaviour produces only strike-parallel displacement (c). The data points (i.e. generated by a slip on the true fault) are shown in blue. The displacements in khaki green are produced by the 1 m uniform target model assuming an incorrect fault dip. The predictions shown in orange and green are calculated from the posterior mean model, respectively, without and with $C_p$ included. The orange and green vertical bars correspond to the $1\sigma$ posterior uncertainty on the predictions, respectively, without and with $C_p$ included. The five missing data points around the fault's position have been removed to simulate a real case.</figcaption></figure>

inclusion of $C_p$ in the optimization allows the quasi-perfect approximation of the target model (Fig. S12b, Supporting Information).

When neglecting $C_p$, whether it be with an optimization or a Bayesian approach, predictions of inferred source models approach the observations as much as possible ignoring the limitations of the forward model. In this case, we place too much confidence or weight in our residuals given our limited ability to predict the observations. Instead, with the introduction of $C_p$, we allow for a larger misfit between observations and predictions. The result of the inversion might seem counter-intuitive because the inferred model approaches the target model but with an increase in misfit. This apparent paradox illustrates the aim of $C_p$, which adjusts the confidence in each observation such that the inversion algorithm is guided by information that is less dependent on inaccurate forward theory. The end result of including $C_p$ is a higher standard deviation of the posterior slip distribution but a more reliable estimate of its maximum likelihood values. If the assumed geometry is incorrect by only a few degrees or a couple of kilometres, as in our simple 2-D tests, observations can be well fit by a slip model at odds with the true model.

The issue of overconfidence brings into question the rigour of solving for some fault geometry parameters before or during the inversion of the slip distribution if $C_p$ is not included. If only one or two fault geometry parameters (e.g. strike and dip) are solved, while another (fault position) is assumed correct when it is not, the resulting source model and fault parameters are possibly all incorrect. Hence, to infer a robust solution, the uncertainty of all parameters that are not part of the inversion should be quantified.

Thus, if the studied fault is of particularly complex morphology, one would want to include the uncertainty related to this complexity in the inversion. Indeed, faults are frequently modelled as planar but

are a lot more complex in reality: curved along dip, bent or segmented along strike. The $C_p$ matrix is designed to include all kinds of geometrical complexities, the only limitation being our ability to compute the Green's Functions accordingly. For a fault segmented along strike, the covariance between azimuths of fault segments can be included within the non-diagonal terms of the covariance matrix $C_\psi$. For instance, for a two-segments planar fault, $C_\psi$ will be of size (2,2). Then, the variation of the Green's functions with the strike of each segment will be introduced with two sensitivity kernels $K_\psi$, one for each segment. For a curved fault, the relation between geometrical parameters (dip of subfaults or strike of subfaults) can also be introduced within $C_\psi$. Using $C_p$, random geometrical complexities (bends, segments) could be included in the inversion process to simulate the impact of an unknown fault roughness. Yet, the uncertainties in fault dip, strike or position, which are the first-order parameters controlling the fault geometry, are assumed to have the largest contribution on uncertainties of the fault geometry. Thus, including uncertainties of the main fault geometry parameters may be sufficient to acknowledge (or overwhelm) most of uncertainties due to geometrical complexities.

To calculate the sensitivity kernels, we have linearized the evolution of Green's functions with respect to fault geometry, the perturbation of predictions being small enough on a restricted fault geometry parameter interval. We show that this approximation is validated in the case of synthetics tests. For highly uncertain fault geometries, this linearization may not be justified.

Note that, in the simplified case of 1 m dip slip, epistemic uncertainties related to fault dip (reaching maximal values of ~$7.10^{-2}$ m for strike-perpendicular displacement, Fig. S13c of the Supporting Information) can be of more than $10^3$ times higher than observational uncertainties ($\sigma^2$, ~$10^{-5}$ m), when we assume a 5° standard
deviation on our preferred fault dip. If we consider the case of great earthquakes with an average of 10 m of dip slip, such as the $M_w$ 9.0, 2011, Tohoku-Oki earthquake, observational errors will remain on the order of $10^{-6}-10^{-3}$ m (around the mm standard deviation for GPS, cm for InSAR). In contrast, epistemic uncertainties could reach 1 m, as they scale with the square of the slip. Epistemic uncertainties could thus be three to six orders of magnitude larger than observational errors for great earthquakes, assuming that only one fault geometry parameter is inaccurate. This ratio is up to $10^3$ times greater than for our simplified 2-D case, showing that uncertainties in fault geometry would be of critical importance when studying great ($M_w$ > 8) and well observed earthquakes.

Looking more closely at the details of our 2-D synthetic tests, we observe that the influence of uncertainty on the fault geometry is different depending on which parameter is analyzed. The uncertainty in fault dip has a stronger influence on the deepest parts of the fault, which can be explained by the increase with depth of the offset between true and assumed fault structures (the location of the fault surface rupture remains the same). Conversely, the uncertainty on fault position mainly impacts the fit of near-fault observations and the estimation of shallow slip. Hence, the impact of uncertainty in the fault geometry will depend on the unknown slip distribution.

The result of the inverse problem will also depend on the distribution of surface displacement data. As illustrated in Fig. 4, the influence of dip variation on surface displacement varies with the distance to the fault. Therefore, if observations are concentrated in areas where the surface displacement has a strong gradient (i.e. where observations allow to better resolve slip on the fault), the impact of dip uncertainty on source estimation will increase. Near-fault observations will be the most affected by an incorrect fault position. If near-fault data are unavailable, the weight of uncertainties will thus be reduced. Thus, the impact of uncertainty in fault geometry will also depend on the data distribution. The link between epistemic uncertainties and surface displacement amplitude is particularly striking when looking at $C_p$ values (Fig. 3). In the case of a subduction event with observations very far to the fault, the impact of uncertainties in the fault geometry thus could be lesser than for a well-instrumented continental earthquake.

The impact of an inaccurate fault position might have contributed to the debate on a possible shallow slip deficit for continental earthquakes, that is, a systematic discrepancy between the amount of slip in the near surface and on the rest of the fault (Simons et al. 2002; Kaneko & Fialko 2011). The shallow slip deficit could be explained by combined effects: a lack of near-fault observations (Xu et al. 2016), inelastic deformation (Simons et al. 2002; Kaneko & Fialko 2011; Gombert et al. 2017) and/or incorrect forward modelling of Earth elastic properties (Simons et al. 2002; Xu et al. 2016). But the shallow slip deficit studies are, by happenstance, usually based on strike-slip earthquakes with long and complex surface fault trace (Xu et al. 2016). Thus, even the construction of an elaborated fault morphology remains a simplification of the reality. For the $M_w$ 7.2 2010 El Mayor Cucapah earthquake, Mexico, the assumed fault geometry can be offset by more than 2 km in the case of a two-segment fault model (Wei et al. 2011) or by more than 1 km for a six-segment fault model (Xu et al. 2016). Yet, in our toy models, a shift of the fault of 2 km induces a significant bias on the estimate of the shallow slip distribution. Part of the observed shallow slip deficit could, thus, result from inaccurate models of the fault geometry. The introduction of $C_p$ in such studies could provide more insight on the importance of the shallow slip deficit.

Earthquake-induced tsunamis are primarily controlled by the amount of seafloor deformation. Therefore, tsunami amplitude is

highly sensitive to the distribution and amplitude of shallow slip (e.g. Geist & Dmowska 1999; Geist et al. 2006) and variations on fault dip (Tanioka & Satake 1996; Goda et al. 2014; Bleryet et al. 2015). Slip and geometry of the ruptured fault can thus be considered as major sources of uncertainty for predicting tsunami hazards and observations (Goda et al. 2014). Unfortunately, up dip ruptured fault geometry is poorly constrained for megathrust earthquakes (Lay et al. 2005; Moore et al. 2007). We showed that inferred slip can be biased if the assumed fault structure is inaccurate. Uncertainty in fault geometry and its impact on shallow slip could thus also influence near-surface slip estimates of megathrust events and thus have potential implications for tsunami predictions.

## 5.2 Accounting for uncertainties in both fault geometry and Earth structure

Fault geometry parameters are not the only parameters that are held fixed in source estimation problems. Earth structure is often chosen a priori despite being an important source of epistemic uncertainty (e.g. Beresnev 2003; Duputel et al. 2014). And just as with fault geometry, Earth structure is almost never known perfectly. With the same toy model, we investigate the relative influence of uncertainty on both fault geometry and Earth structure. To generate our synthetic data, we assume a pure vertical and strike-slip fault with a 20 km long down-dip extent, that prolongs infinitely along strike, and which is discretized every kilometre along dip. This fault is embedded in a 30 km thick low compliance layer of shear modulus $\mu_1$, overlying an elastic half-space of shear modulus $\mu_2$. We define the shear modulus of both layers so that the ratio of the moduli ($\mu_2/\mu_1$) is equal to 1.3. Synthetic observations are calculated for 100 data points irregularly distributed at the surface, using analytical expressions of surface displacements for 1 m of uniform slip on a strike-slip fault within a layer over a half-space (Segall 2010).

For the inversion, we assume a fault with a dip of 85° and calculate our Green's functions for an homogeneous elastic half-space of shear modulus $\mu_2$. To account for uncertainties on both fault geometry and Earth model, $C_p^{\text{fault}}$ and $C_p^{\text{earth}}$ are calculated a priori, following the approach defined in Section 2.3 and in Duputel et al. (2014). The initial model is chosen as the result of an inversion not accounting for any prediction uncertainty. $C_p^{\text{fault}}$ or $C_p^{\text{earth}}$ or both are then added to a 7 mm data covariance matrix $C_d$, and the misfit covariance matrix $C_\chi = C_d+(C_p^{\text{fault}})+(C_p^{\text{earth}})$ replaces $C_d$ in the inversion process. We then estimate the slip distribution with our Bayesian approach (approximately 60 000 samples are explored at each tempering step) with a uniform prior distribution $p(\mathbf{m}) = U(-0.5 \text{ m}, 5 \text{ m})$.

The posterior mean model as well as the PDFs for each fault patch are presented in Fig. 9. The inversion results have been solved either neglecting $C_p^{\text{fault}}$ and $C_p^{\text{earth}}$ (Fig. 9a), including $C_p^{\text{fault}}$ only (Fig. 9b), $C_p^{\text{earth}}$ only (Fig. 9c) or both (Fig. 9d). As expected, with no $C_p$ included, the inversion recovers the uniform 1 m target model poorly (Fig. 9a). The inclusion of only one type of uncertainty does not improve the fit to the target model either (Figs 9b and c): in both cases, the model misfit does not exceed 1 m but is also rarely below 25 cm. Only the introduction of both types of uncertainty allows us to recover the target model, with a model misfit below 10 cm for most of subfaults [difference reach 20 cm (20 per cent) on the two deepest fault segments]. The predictions are very close to the observations whether $C_p$ is included or not. The RMS on the surface displacement data is lower when one or no $C_p$ term is
STRIKE SLIP, 1D marginal PDFs

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Comparison of inversion results with and without incorporating the prediction uncertainty. The uniform target slip model is shown as a grey vertical bar. The true fault is vertical and embedded in a 30 km thick shallow low compliance layer overlying an elastic half-space. The assumed fault is of 85° dip embedded in an elastic half-space. The posterior PDFs are shown in (a) when the prediction uncertainty is neglected and in (b) and (c) if only $C_p^{\text{fault}}$ or $C_p^{\text{earth}}$ is included in the inversion. In (d), both uncertainties on fault geometry and Earth structure are included in the source problem. The posterior mean model is indicated with a vertical grey dotted line. The offset between target model and posterior mean is displayed with a colour scale saturated at 50 cm.</figcaption></figure>

included (2.0 mm). When both $C_p^{\text{fault}}$ and $C_p^{\text{earth}}$ are accounted for, data residuals are higher (RMS of 2.7 mm).

The relative importance of inaccuracies on fault geometry and Earth structure will be variable between an earthquake and another. Yet, in our toy model, realistic errors on fault geometry and rheology have an equivalent influence on inferred models. This simple approach thus highlights the need to account for the uncertainties associated with all parameters of the forward model that are not inverted for.

## 6 CONCLUSION

Inversion method, parametrization of the problem and simplifications of the forward physics are many factors responsible for the ill-posed nature of the inverse problem. Among these factors, the assumed Earth structure and fault geometry are simplifications of the reality and have a major impact on the source analysis. Indeed, fixed values of fault geometry parameters or Earth elastic properties are inherently unrealistic (as not reflecting the roughness of the fault and of the Earth), but more importantly are uncertain. The uncertainties of a source inversion problem are defined in the misfit covariance matrix $C_X$. This covariance matrix is involved in the
inversion process as part of the misfit function, which is used as a criterion to select the preferred source model in any inversion method. Although observational errors $C_d$ are generally included in $C_\chi$, epistemic uncertainties that are related to the forward problem are often ignored. However, we show that uncertainties in fault geometry could be $10^3-10^6$ times higher than observational errors ($\sigma^2$) for well-observed great earthquakes ($M_w > 8$).

To account for uncertainties in the fault geometry, we propose a formalism based on small perturbation theory that limits its application to events for which the gross characteristics of the fault geometry are known. This formalism is designed for any inversion method, and is tested here for optimization and Bayesian approaches. A prediction covariance matrix $C_p$ is quantified from the sensitivity of forward predictions with respect to the variation of a particular fault geometry parameter. Then, we assume an evolved misfit covariance matrix $C_\chi = C_d + C_p$ in the inversion process, regardless of which sampling method is chosen. In practice, we try our formalism for uncertainties on the fault dip and fault position, a parametrization motivated by simplicity and the large range of possible fault configurations which can then be addressed. The mere inclusion of uncertainty in the fault dip and position reproduces the first-order inaccuracies when defining a particular fault structure. Accounting for related uncertainties in the inversion process, whatever the inversion method, allows to reliably estimate the posterior distribution of source model parameters. In particular, we find that uncertainty on the fault dip and position have complementary influence, and allow to improve the reliability of slip estimates both on the shallow and deep parts of the fault. When introducing $C_p$ in the inversion process, we prevent the predictions to be overconfident in the data. This modification of the inverse problem leads to a robust posterior source model. Finally, we investigate the joint impact of uncertainties in the Earth structure and fault geometry. Our results emphasize the need to account for both sources of uncertainty to obtain a robust source model.

Uncertainty in fault geometry could thus have biased shallow slip estimates for numerous earthquakes. Inaccuracy of assumed up dip fault structures could be another contributor to the observed shallow slip deficit for continental earthquakes. The poorly constrained morphology of ruptured faults for megathrust events could also have a high impact on near-surface slip estimates and potential implications for the tsunami hazard assessment.

The results of this study are accurate for our 2-D simplified application of an infinitely long fault. Although the discrepancy between assumed and true fault geometries is similar to what can be observed for real events, the striking impact of uncertainties in fault geometry may be enhanced by our simplified approach. Nevertheless, accounting for uncertainty of the forward problem allows us to infer source models with a level of precision that is not reachable otherwise.

## ACKNOWLEDGEMENTS

We thank Martin Mai, František Gallovič and Kaj Johnson for their thorough and constructive reviews. We thank Quentin Bletery, Baptiste Gombert, Zacharie Duputel and Romain Jolivet for helpful discussions. The Bayesian simulations were performed on the HPC-Regional Center ROMEO (https://romeo.univ-reims.fr) of the University of Reims Champagne-Ardenne (France). We thank Michael Aivazis for his meticulous help for the installation. The Classic Slip Inversion (CSI) Python library created by Romain Jolivet was used

to build inputs for the Bayesian algorithm and to perform the optimization of the inverse problem. Figures were generated with the Matplotlib Python library. This study was partly supported by the French National Research Agency (ANR) EPOST project ANR-14-CE03-0002. TR is supported by a fellowship from the French Ministry of Research and Higher Education. MS was supported by U.S. National Science Foundation grant 1447107.

## REFERENCES

* Aochi, H. & Fukuyama, E., 2002. Three-dimensional nonplanar simulation of the 1992 Landers earthquake, *J. geophys. Res.*, **107**(B2), ESE 4–1.

* Aochi, H. & Madariaga, R., 2003. The 1999 Izmit, Turkey, earthquake: nonplanar fault structure, dynamic rupture process, and strong ground motion, *Bull. seism. Soc. Am.*, **93**(3), 1249–1266.

* Beresnev, I.A., 2003. Uncertainties in finite-fault slip inversions: to what extent to believe? (a critical review), *Bull. seism. Soc. Am.*, **93**(6), 2445–2458.

* Bletery, Q., Sladen, A., Delouis, B. & Mattéo, L., 2015. Quantification of Tsunami Bathymetry effect on finite fault slip inversion, *Pure appl. Geophys.*, **172**(12), 3655–3670.

* Bletery, Q., Sladen, A., Jiang, J. & Simons, M., 2016a. A Bayesian source model for the 2004 great Sumatra–Andaman earthquake, *J. geophys. Res.*, **121**(7), 5116–5135.

* Bletery, Q., Thomas, A.M., Rempel, A.W., Karlstrom, L., Sladen, A. & Barros, L.D., 2016b. Mega-earthquakes rupture flat megathrusts, *Science*, **354**(6315), 1027–1031.

* Bouchon, M., Campillo, M. & Cotton, F., 1998. Stress field associated with the rupture of the 1992 Landers, California, earthquake and its implications concerning the fault strength at the onset of the earthquake, *J. geophys. Res.*, **103**(B9), 21 091–21 097.

* Candela, T., Renard, F., Klinger, Y., Mair, K., Schmittbuhl, J. & Brodsky, E.E., 2012. Roughness of fault surfaces over nine decades of length scales, *J. geophys. Res.*, **117**(B8), B08409, doi:10.1029/2011JB009041.

* Das, S. & Suhadolc, P., 1996. On the inverse problem for earthquake rupture: the Haskell-type source model, *J. geophys. Res.*, **101**(B3), 5725–5738.

* Delouis, B., Giardini, D., Lundgren, P. & Salichon, J., 2002. Joint inversion of INSAR, GPS, teleseismic, and strong-motion data for the spatial and temporal distribution of earthquake slip: application to the 1999 Izmit Mainshock, *Bull. seism. Soc. Am.*, **92**(1), 278–299.

* Diao, F., Wang, R., Aochi, H., Walter, T.R., Zhang, Y., Zheng, Y. & Xiong, X., 2016. Rapid kinematic finite-fault inversion for an $M_w 7+$ scenario earthquake in the Marmara Sea: an uncertainty study, *Geophys. J. Int.*, **204**(2), 813–824.

* 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.

* Duputel, Z., Rivera, L., Fukahata, Y. & Kanamori, H., 2012. Uncertainty estimations for seismic source inversions, *Geophys. J. Int.*, **190**(2), 1243–1256.

* Duputel, Z. et al., 2015. The Iquique earthquake sequence of April 2014: Bayesian modeling accounting for prediction uncertainty, *Geophys. Res. Lett.*, **42**(19), 7949–7957.

* Du, Y., Aydin, A. & Segall, P., 1992. Comparison of various inversion techniques as applied to the determination of a geophysical deformation model for the 1983 Borah Peak earthquake, *Bull. seism. Soc. Am.*, **82**(4), 1840–1866.

* Fukuda, Jun'ichi & Johnson, Kaj M., 2010. Mixed linear-non-linear inversion of crustal deformation data: Bayesian inference of model, weighting and regularization parameters, *Geophys. J. Int.*, **181**(3), 1441–1458.

* Gallovič, F., Imperatori, W. & Mai, P.M., 2015. Effects of three-dimensional crustal structure and smoothing constraint on earthquake slip inversions: case study of the Mw6.3 2009 L'Aquila earthquake, *J. geophys. Res.*, **120**(1), 428–449.

* Gallovič, F. & Ampuero, J.-P., 2015. A new strategy to compare inverted rupture models exploiting the eigenstructure of the inverse problem, *Seismol. Res. Lett.*, **86**(6), 1679–1689.
Geist, E.L., Bilek, S.L., Arcas, D. & Titov, V.V., 2006. Differences in tsunami generation between the December 26, 2004 and March 28, 2005 Sumatra earthquakes, *Earth Planets Space*, **58**(2), 185–193.

Geist, E.L. & Dmowska, R., 1999. Local tsunamis and distributed slip at the source, *Pure appl. Geophys.*, **154**, 485–512.

Goda, K., Mai, P.M., Yasuda, T. & Mori, N., 2014. Sensitivity of tsunami wave profiles and inundation simulations to earthquake slip and fault geometry for the 2011 Tohoku earthquake, *Earth Planets Space*, **66**, 105, doi:10.1186/1880-5981-66-105.

Gombert, B., Duputel, Z., Jolivet, R., Doubre, C., Rivera, L. & Simons, M., 2018. Revisiting the 1992 Landers earthquake: a Bayesian exploration of co-seismic slip and off-fault damage, *Geophys. J. Int.*, **212**, 2, 839–852

Hartzell, S., Liu, P., Mendoza, C., Ji, C. & Larson, K.M., 2007. Stability and uncertainty of finite-fault slip inversions: application to the 2004 Parkfield, California, Earthquake, *Bull. seism. Soc. Am.*, **97**(6), 1911–1934.

Huang, M.-H., Fielding, E.J., Dickinson, H., Sun, J., Gonzalez-Ortega, J.A., Freed, A.M. & Bürgmann, R., 2017. Fault geometry inversion and slip distribution of the 2010 M<sub>w</sub> 7.2 El Mayor-Cucapah earthquake from geodetic data, *J. geophys. Res.*, **122**(1), 607–621.

Jolivet, R., Simons, M., Agram, P.S., Duputel, Z. & Shen, Z.-K., 2015. Aseismic slip and seismogenic coupling along the central San Andreas Fault, *Geophys. Res. Lett.*, **42**(2), 297–306.

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

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

King, G. & Náblek, J., 1985. Role of fault bends in the initiation and termination of earthquake rupture, *Science*, **228**(4702), 984–987.

Lavecchia, G., Ferrarini, F., Brozzetti, F., Nardis, R.D., Boncio, P. & Chialaluce, L., 2012. From surface geology to aftershock analysis: constraints on the geometry of the L'Aquila 2009 seismogenic fault system, *Ital. J. Geosci.*, **131**(3), 330–347.

Lay, T., 2018. A review of the rupture characteristics of the 2011 Tohoku-Oki Mw 9.1 earthquake, *Tectonophysics*, **733**, 4–36.

Lay, T. et al., 2005. The great Sumatra–Andaman earthquake of 26 December 2004, *Science*, **308**(5725), 1127–1133.

Lee, S.-J., Ma, K.-F. & Chen, H.-W., 2006. Effects of fault geometry and slip style on near-fault static displacements caused by the 1999 Chi-Chi, Taiwan earthquake, *Earth planet. Sci. Lett.*, **241**(1-2), 336–350.

Liu, C., Zheng, Y., Xie, Z. & Xiong, X., 2017. Rupture features of the 2016 M<sub>w</sub> 6.2 Norcia earthquake and its possible relationship with strong seismic hazards, *Geophys. Res. Lett.*, **44**, 1320–1328.

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.

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

Manighetti, I., Caulet, C., De Barros, L., Perrin, C., Cappa, F. & Gaudemer, Y., 2015. Generic along-strike segmentation of afar normal faults, East Africa: implications on fault growth and stress heterogeneity on seismogenic fault planes, *Geochem. Geophys. Geosyst.*, **16**(2), 443–467.

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

Minson, S.E. et al., 2014. Bayesian inversion for finite fault earthquake source models - II: the 2011 great Tohoku-Oki, Japan earthquake, *Geophys. J. Int.*, **198**(2), 922–940.

Monelli, D. & Mai, P.M., 2008. Bayesian inference of kinematic earthquake rupture parameters through fitting of strong motion data, *Geophys. J. Int.*, **173**(1), 220–232.

Moore, G.F., Bangs, N.L., Taira, A., Kuramoto, S., Pangborn, E. & Tobin, H.J., 2007. Three-dimensional splay fault geometry and implications for tsunami generation, *Science*, **318**(5853), 1128–1131.

Moreno, M. et al., 2012. Toward understanding tectonic control on the M<sub>w</sub> 8.8 2010 Maule Chile earthquake, *Earth planet. Sci. Lett.*, **321-322**, 152–165.

Moreno, M.S., Bolte, J., Klotz, J. & Melnick, D., 2009. Impact of megathrust geometry on inversion of coseismic slip from geodetic data: application to the 1960 Chile earthquake, *Geophys. Res. Lett.*, **36**(16), L16310, doi:10.1029/2009GL039276.

Okubo, P.G. & Aki, K., 1987. Fractal geometry in the San Andreas fault system, *J. geophys. Res.*, **92**(B1), 345–355.

Peacock, D.C.P., 1991. Displacements and segment linkage in strike-slip fault zones, *J. Struct. Geol.*, **13**(9), 1025–1035.

Perrin, C., Manighetti, I., Ampuero, J.-P., Cappa, F. & Gaudemer, Y., 2016. Location of largest earthquake slip and fast rupture controlled by along-strike change in fault structural maturity due to fault growth, *J. geophys. Res.*, **121**(5), 3666–3685.

Piatanesi, A., Cirella, A., Spudich, P. & Cocco, M., 2007. A global search inversion for earthquake kinematic rupture history: application to the 2000 western Tottori, Japan earthquake, *J. geophys. Res.*, **112**(B7), B07314, doi:10.1029/2006JB004821.

Pizzi, A., Di Domenica, A., Gallovic, F., Luzi, L. & Puglia, R., 2017. Fault segmentation as constraint to the occurrence of the main shocks of the 2016 central Italy seismic sequence, *Tectonics*, **36**(11), 2370–2387.

Power, W.L., Tullis, T.E., Brown, S.R., Boitnott, G.N. & Scholz, C.H., 1987. Roughness of natural fault surfaces, *Geophys. Res. Lett.*, **14**(1), 29–32.

Razafindrakoto, H.N.T. & Mai, P.M., 2014. Uncertainty in earthquake source imaging due to variations in source time function and earth structure, *Bull. seism. Soc. Am.*, **104**(2), 855–874.

Razafindrakoto, H.N.T., Mai, P.M., Genton, M.G., Zhang, L. & Thingbaijam, K.K.S., 2015. Quantifying variability in earthquake rupture models using multidimensional scaling: application to the 2011 Tohoku earthquake, *Geophys. J. Int.*, **202**(1), 17–40.

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

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

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.

Sladen, A. et al., 2010. Source model of the 2007 Mw 8.0 Pisco, Peru earthquake: implications for seismogenic behavior of subduction megathrusts, *J. geophys. Res.*, **115**(B2), B02405, doi:10.1029/2009JB006429.

Sudhaus, H. & Jónsson, S., 2009. Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: application to the June 2000 Kleifarvatn earthquake, Iceland, *Geophys. J. Int.*, **176**(2), 389–404.

Sun, J., Johnson, K.M., Cao, Z., Shen, Z., Bürgmann, R. & Xu, X., 2011. Mechanical constraints on inversion of coseismic geodetic data for fault slip and geometry: Example from InSAR observation of the 6 October 2008 Mw 6.3 Dangxiong-Yangyi (Tibet) earthquake, *J. geophys. Res.*, **116**, B01406.

Tanioka, Y. & Satake, K., 1996. Tsunami generation by horizontal displacement of ocean bottom, *Geophys. Res. Lett.*, **23**(8), 861–864.

Tarantola, A., 2005. *Inverse Problem Theory and Methods for Model Parameter Estimation*, Philadelphia: Society for Industrial and Applied Mathematics.

Tichelaar, B.W. & Ruff, L.J., 1989. How good are our best models? Jackknifing, bootstrapping, and earthquake depth, *EOS, Trans. Am. Geophys. Un.*, **70**(20), 593–606.

Tinti, E., Scognamiglio, L., Michelini, A. & Cocco, M., 2016. Slip heterogeneity and directivity of the M<sub>L</sub> 6.0, 2016, Amatrice earthquake estimated with rapid finite-fault inversion, *Geophys. Res. Lett.*, **43**, 10 745–10 752.

Volpe, M., Atzori, S., Piersanti, A. & Melini, D., 2015. The 2009 L'Aquila earthquake coseismic rupture: open issues and new insights from 3D finite element inversion of GPS, InSAR and strong motion data, *Ann. Geophys.*, **58**(2), S0221, doi:10.4401/ag.6711.
Walsh, J.J., Bailey, W.R., Childs, C., Nicol, A. & Bonson, C.G., 2003. Formation of segmented normal faults: a 3-D perspective, *J. Struct. Geol.*, **25**(8), 1251–1262.

Wei, S. et al., 2011. Superficial simplicity of the 2010 El Mayor-Cucapah earthquake of Baja California in Mexico, *Nat. Geosci.*, **4**(9), 615–618.

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

Yagi, Y. & Fukahata, Y., 2008. Importance of covariance components in inversion analyses of densely sampled observed data: an application to waveform data inversion for seismic source processes, *Geophys. J. Int.*, **175**(1), 215–221.

Yano, T.E., Shao, G., Liu, Q., Ji, C. & Archuleta, R.J., 2014. Coseismic and potential early afterslip distribution of the 2009 Mw 6.3 L'Aquila, Italy earthquake, *Geophys. J. Int.*, **199**(1), 23–40.

Zhang, P., Slemmons, D.B. & Mao, F., 1991. Geometric pattern, rupture termination and fault segmentation of the Dixie Valley-Pleasant Valley active normal fault system, Nevada, U.S.A., *J. Struct. Geol.*, **13**(2), 165–176.

Zielke, O., Galis, M. & Mai, P.M., 2017. Fault roughness and strength heterogeneity control earthquake size and stress drop, *Geophys. Res. Lett.*, **44**, 777–783.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** Comparisons between standardized data and predictions with and without accounting for $C_p$, for a fault with dip-slip motion and assuming an incorrect dip angle. The location of the fault surface rupture is shown with a grey vertical line. Standardized data and predictions without accounting for $C_p$ (i.e. $C_χ = C_d$) are shown in (a) and (b), respectively, for strike perpendicular and vertical surface displacements. Standardized data and predictions accounting for $C_p$ (i.e. $C_χ = C_p + C_d$) are shown in (c) and (d). The standardized observations $Rd_{obs}$ (with $C_χ^{-1} = R^T R$) are shown in blue. The standardized displacements in khaki green are produced by the 1 m uniform target model assuming an incorrect fault dip. The standardized predictions $Rd_{pred}$ shown in orange and green are calculated from the posterior mean model, respectively, without and with $C_p$ included (and thus with a different R matrix).

**Figure S2.** Comparison of inversion results with and without neglecting the prediction uncertainty, for a dip-slip scenario with varying assumed prior. The uniform target slip model is shown as a grey vertical bar. In (a) and (b), the assumed prior allows for 50 cm of up dip slip, whereas in (c) and (d) the prior is strictly positive. The posterior PDFs are shown in (a) and (c) when the prediction uncertainty is neglected and in (b) and (d) if $C_p$ is included in the inversion. The posterior mean model is indicated with a vertical grey dotted line. The offset between target model and posterior mean is displayed with a colour scale saturated at 50 cm.

**Figure S3.** Comparison of inversion results with and without neglecting the prediction uncertainty, for a dip-slip scenario. The uniform target slip model is shown as a grey vertical bar. The correct fault is of 55° dip, whereas the assumed fault for inversion is of 50° dip in (a) and (b) and of 45° dip in (c) and (d). The posterior PDFs are shown in (a) and (c) when the prediction uncertainty is neglected and in (b) and (d) if $C_p$ is included in the inversion. The posterior mean model is indicated with a vertical grey dotted line. The offset between target model and posterior mean is displayed with a colour scale saturated at 50 cm.

**Figure S4.** Comparison of inversion results with and without neglecting the prediction uncertainty. The non-uniform target slip model is shown as a grey vertical bar. The posterior PDFs are

shown in (a) and (c) when the prediction uncertainty is neglected and in (b) and (d) if $C_p$ is included in the inversion, respectively, for dip-slip and strike-slip behaviour. The posterior mean model is indicated with a vertical grey dotted line. The offset between target model and posterior mean is displayed with a colour scale saturated at 50 cm.

**Figure S5.** Comparisons between synthetic data and predictions with and without accounting for $C_p$. We assume the slip on the fault to be non-uniform. The location of the fault surface rupture is shown with a grey vertical line. Dip-slip behaviour (assumed fault dip of 50°) generates cross strike surface displacement (a) and vertical surface displacement (b), whereas strikeslip (assumed fault dip of 75°) behaviour produces only strike parallel displacement (c). The data points (i.e. generated by a slip on the true fault) are shown in blue. The displacements in khaki green are produced by the 1 m uniform target model assuming an incorrect fault dip. The predictions shown in orange and green are calculated from the posterior mean model, respectively, without and with $C_p$ included.

**Figure S6.** Comparison of inversion results with or without neglecting for the prediction uncertainty for a dip-slip scenario. The uniform target slip model is shown as a grey vertical bar. The posterior PDFs are shown in (a) when the prediction uncertainty is neglected. When $C_p$ is included in the inversion, the assumed standard deviation varies from 1°, 5° to 10° respectively, in (b), (c) and (d). The posterior mean model is indicated with a vertical grey dotted line. The offset between target model and posterior mean is displayed with a colour scale saturated at 50 cm.

**Figure S7.** Comparison of inversion results with or without neglecting for the prediction uncertainty for a dip-slip scenario. The uniform target slip model is shown as a grey vertical bar. The posterior PDFs are shown in (a) when the prediction uncertainty is neglected. When $C_p$ is included in the inversion, the assumed standard deviation varies from 1°, 5° to 10°, respectively, in (b), (c) and (d). In this case, the observational errors are assumed of 4 mm, and are thus underestimated as the added noise is of 7 mm. The posterior mean model is indicated with a vertical grey dotted line. The offset between target model and posterior mean is displayed with a colour scale saturated at 50 cm.

**Figure S8.** Comparison of inversion results with or without neglecting for the prediction uncertainty for a dip-slip scenario. The uniform target slip model is shown as a grey vertical bar. The posterior PDFs are shown in (a) and (c) when the prediction uncertainty is neglected and in (b) and (d) if $C_p$ is included. In (a) and (b), the noise added to the data has an amplitude of 7 mm. In (c) and (d), the noise has an amplitude of 7 mm and is correlated spatially. The posterior mean model is indicated with a vertical grey dotted line. The offset between target model and posterior mean is displayed with a colour scale saturated at 50 cm.

**Figure S9.** Comparison of inversion results with and without neglecting the prediction uncertainty. The non-uniform target slip model is shown as a grey vertical bar. The offset between correct and assumed faults is of 2 km. The posterior PDFs are shown in (a) and (c) when the prediction uncertainty is neglected and in (b) and (d) if $C_p$ is included in the inversion, respectively, for dip-slip and strike-slip behaviour. The posterior mean model is indicated with a vertical grey dotted line. The offset between target model and posterior mean is displayed with a colour scale saturated at 50 cm.

**Figure S10.** Comparisons between synthetic data and predictions with and without accounting for $C_p$. We assume the slip on the fault to be non-uniform. The correct location of the fault surface rupture is shown with a grey vertical line, whereas the assumed fault position is shown by a dotted grey vertical line. Dip-slip behaviour
generates cross strike surface displacement (a) and vertical surface displacement (b), whereas strike-slip behaviour produces only strike parallel displacement (c). The data points (i.e. generated by a slip on the true fault) are shown in blue. The displacements in khaki green are produced by the 1 m uniform target model assuming an incorrect fault dip. The predictions shown in orange and green are calculated from the posterior mean model, respectively, without and with C<sub>p</sub> included. The five missing data points around the faults position have been removed to simulate a real case.

**Figure S11.** Misfit (see eq. 1, with C<sub>χ</sub>(m) = C<sub>d</sub>) as a function of maximum slip amplitude (indication of model roughness) for different optimizations made assuming correlation lengths varying from 2 to 26 km. We finally assume for our optimizations the value of 14 km. The a priori standard deviation of model parameters σ has been chosen from the posterior standard deviation of model parameters of a Bayesian optimization (see Fig. 5). Our assumed σ is of 50 cm. The scaling factor λ<sub>0</sub> is usually chosen as average distance between subfaults, and is thus assumed here of 1 km.

**Figure S12.** Comparison of results for a positive least squares optimization, with and without incorporating C<sub>p</sub>. The uniform target slip model is shown as a grey vertical bar. The inversion results are shown in (a) and (c) when the prediction uncertainty is neglected and in (b) and (d) if C<sub>p</sub> is included in the inversion, respectively for a 20 subfaults or a 8 subfaults fault parametrization. The inferred model corresponds to the mean of the posterior gaussian distribution, the standard deviation being the uncertainty on inferred parameters C<sub>m</sub> = (G<sup>T</sup> · C<sub>χ</sub><sup>-1</sup> · G)<sup>-1</sup>. The posterior mean model is indicated with a grey vertical dotted line, or with a colored dotted line if the posterior distribution is almost uniform (i.e the posterior uncertainty is high). The offset between target and resulting model is displayed with a colorscale saturated at 50 cm.

**Figure S13:** C<sub>p</sub> due to imprecise fault dip in a simplified 2D application. The assumed fault does not vary along strike, is 20 km large and is dipping of 35° while the correct fault is dipping 30°. The assumed model m<sub>prior</sub> is uniform at 1 m. The prediction covariance matrix C<sub>p</sub> is shown for strike parallel, vertical and cross strike surface displacement respectively in (a), (b) and (c).

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.

# APPENDIX: DIP-DEPENDENT EXPRESSIONS OF SURFACE DISPLACEMENT FOR AN INFINITE FAULT

We consider a simple 2-D model of a fault that extends infinitely along strike. For a strike-slip case, the displacement on the fault is restricted to the x<sub>3</sub> direction (i.e. along strike direction). For a dip-slip case, the displacements are restricted to the x<sub>1</sub>, x<sub>2</sub> plane (respectively perpendicular to strike and vertical directions). Each surface displacement is calculated for a particular data point located at x<sub>1</sub> along the strike-perpendicular direction x<sub>1</sub>. These expressions of surface displacement in an homogeneous half-space are modified after the fixed-dip equations of Segall (2010) to become dip-dependent.

## A1 For a uniform slip distribution

The fault has a width w, a dip δ and the slip is constant over the fault with a value s.

### A1.1 Strike-slip dislocation

The free surface displacement (i.e. at x<sub>2</sub> = 0) of a vertical fault (δ = 90°) with a strike-slip surface breaking dislocation is (Segall 2010)

$$u_3(x_2 = 0) = -\frac{s}{\pi} \left[ \tan^{-1} \frac{x_1}{d_1} - \frac{\pi}{2} \operatorname{sgn}(x_1) \right], \quad (A1)$$

with x<sub>1</sub> being the distance from the dislocation along the x<sub>1</sub> axis and d<sub>1</sub> being the depth end of dislocation.

For a non-vertical fault, we develop from eq. (A1). For a fault dipping of δ, the distance from the dislocation thus becomes x<sub>1</sub> − w cos δ and the depth end of dislocation becomes . The surface displacement is thus

$$u_3(x_2 = 0) = -\frac{s}{\pi} \left[ \tan^{-1} \xi - \frac{\pi}{2} \operatorname{sgn}(x_1) \right], \quad (A2)$$

with

$$\xi = \frac{x_1 - w \cos \delta}{w \sin \delta}. \quad (A3)$$

### A1.2 Dip-slip dislocation

From Segall (2010), the surface displacement due to a surface breaking dislocation of a fault dipping of δ are (u<sub>2</sub> being the vertical displacement and u<sub>1</sub> the displacement perpendicular to the strike direction)

$$u_2(x_2 = 0) = \frac{s}{\pi} \left[ \sin \delta \left( \tan^{-1} \xi - \frac{\pi}{2} \operatorname{sgn}(x_1) \right) + \frac{\cos \delta + \xi \sin \delta}{1 + \xi^2} \right], \quad (A4)$$

$$u_1(x_2 = 0) = -\frac{s}{\pi} \left[ \cos \delta \left( \tan^{-1} \xi - \frac{\pi}{2} \operatorname{sgn}(x_1) \right) + \frac{\sin \delta - \xi \cos \delta}{1 + \xi^2} \right], \quad (A5)$$

with

$$\xi = \frac{x_1 - w \cos \delta}{w \sin \delta}. \quad (A6)$$

## A2 For a non-uniform slip distribution

The fault has a width w, a dip δ and is divided into N subfaults. On each subfault, the slip has a different value{s<sub>n</sub>}<sub>n=0:N</sub>.

The surface displacement due to a non-uniform slip distribution can be calculated as the sum of the surface displacements produced by the uniform slip on each of the N subfaults. The surface displacement produced by an umpteenth subfault n can be calculated using uniform slip and surface breaking dislocation equations. Indeed, the surface displacement due to subfault n corresponds to the surface displacement due to slip from surface to subfault n minus the same slip from surface to subfault (n − 1).
## A2.1 Strike-slip dislocation

The free surface displacement (i.e. at $x_2 = 0$) for a strike-slip case is thus (from eq. A2)

$$u_3(x_2 = 0) = \sum_{n=0}^{N-1} -\frac{s_n}{\pi} [\tan^{-1} \xi_{n+1} - \tan^{-1} \xi_n] \quad (A7)$$

as

$$\lim_{n \to 0} \tan^{-1} \frac{Nx_1 - nw \cos \delta}{nw \sin \delta} = \frac{\pi}{2} \operatorname{sgn}(x_1), \quad (A8)$$

with

$$\xi_n = \frac{Nx_1 - nw \cos \delta}{nw \sin \delta}. \quad (A9)$$

## A2.2 Dip-slip dislocation

The surface displacement in the vertical direction for a dip-slip case is (from eq. A4)

$$\begin{aligned} u_2(x_2 = 0) = & \frac{s_0}{\pi} \left[ \sin \delta \left( \tan^{-1} \xi_1 - \frac{\pi}{2} \operatorname{sgn}(x_1) \right) \right. \\ & \left. + \frac{\cos \delta + \xi_1 \sin \delta}{1 + \xi_1^2} \right] \\ & + \sum_{n=1}^{N-1} -\frac{s_n}{\pi} \left[ \sin \delta \tan^{-1} \xi_n + \frac{\cos \delta + \xi_n \sin \delta}{1 + \xi_n^2} \right] \\ & + \frac{s_n}{\pi} \left[ \sin \delta \tan^{-1} \xi_{n+1} + \frac{\cos \delta + \xi_{n+1} \sin \delta}{1 + \xi_{n+1}^2} \right]. \end{aligned} \quad (A10)$$

The surface displacement in the strike-perpendicular direction for a dip-slip case is (from eq. A5)

$$\begin{aligned} u_1(x_2 = 0) = & \frac{s_0}{\pi} \left[ \cos \delta \left( \tan^{-1} \xi_1 - \frac{\pi}{2} \operatorname{sgn}(x_1) \right) \right. \\ & \left. + \frac{\sin \delta - \xi_1 \cos \delta}{1 + \xi_1^2} \right] \\ & + \sum_{n=1}^{N-1} \frac{s_n}{\pi} \left[ \cos \delta \tan^{-1} \xi_n + \frac{\sin \delta - \xi_n \cos \delta}{1 + \xi_n^2} \right] \\ & - \frac{s_n}{\pi} \left[ \cos \delta \tan^{-1} \xi_{n+1} + \frac{\sin \delta - \xi_{n+1} \cos \delta}{1 + \xi_{n+1}^2} \right], \end{aligned} \quad (A11)$$

with

$$\xi_n = \frac{Nx_1 - nw \cos \delta}{nw \sin \delta}. \quad (A12)$$