Accounting For Prediction Uncertainty When Inferring Subsurface Fault Slip

Geophysical Journal International
Geophys. J. Int. (2014) doi: 10.1093/gji/ggt517

Geophysical Journal International Advance Access published January 24, 2014
Accounting for prediction uncertainty when inferring subsurface

fault slip
Zacharie Duputel,1,∗ Piyush S. Agram,1,† Mark Simons,1 Sarah E. Minson1

and James L. Beck2
1 Seismological Laboratory, California Institute of Technology, 1200 E California Blvd., Pasadena, CA 91125-2100, USA. E-mail: zacharie.duputel@unistra.fr
2 Division of Engineering and Applied Science, California Institute of Technology, 1200 E California Blvd., Pasadena, CA 91125-2100, USA
Accepted 2013 December 20. Received 2013 December 20; in original form 2013 August 1
SUMMARY
This study lays the groundwork for a new generation of earthquake source models based on
Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

a general formalism that rigorously quantifies and incorporates the impact of uncertainties in
fault slip inverse problems. We distinguish two sources of uncertainty when considering the
discrepancy between data and forward model predictions. The first class of error is induced
by imperfect measurements and is often referred to as observational error. The second source
of uncertainty is generally neglected and corresponds to the prediction error, that is the
uncertainty due to imperfect forward modelling. Yet the prediction error can be shown to
scale approximately with the size of earthquakes and thus can dwarf the observational error,
particularly for large events. Both sources of uncertainty can be formulated using the misfit
covariance matrix, Cχ , which combines a covariance matrix for observation errors, Cd and a
covariance matrix for prediction errors, Cp , associated with inaccurate model predictions. We
develop a physically based stochastic forward model to treat the model prediction uncertainty
and show how Cp can be constructed to explicitly account for some of the inaccuracies
in the earth model. Based on a first-order perturbation approach, our formalism relates Cp
to uncertainties on the elastic parameters of different regions (e.g. crust, mantle, etc.). We
demonstrate the importance of including Cp using a simple example of an infinite strike-slip
GJI Seismology
fault in the quasi-static approximation. In this toy model, we treat only uncertainties in the
1-D depth distribution of the shear modulus. We discuss how this can be extended to general
3-D cases and applied to other parameters (e.g. fault geometry) using our formalism for Cp .
The improved modelling of Cp is expected to lead to more reliable images of the earthquake
rupture, that are more resistant to overfitting of data and include more realistic estimates of
uncertainty on inferred model parameters.
Key words: Inverse theory; Probability distributions; Earthquake source observations.
2007), earthquake energy budget analysis (e.g. Kanamori & Rivera

1 I N T RO D U C T I O N
2006) and Coulomb stress transfer calculations (e.g. King 2007).
Inferring earthquake source models is an essential ingredient in ef- Source estimates are also useful for rapid assessment and response
forts to understand the physics of seismic rupture phenomena and to seismic disasters when they occur (e.g. Duputel et al. 2012b).
the relationship of an earthquake with its tectonic and geodynamic However, source inversion algorithms usually do not include real-
environment. As such, the earthquake source model is not only an istic error analyses and their results are generally not accompanied
end into itself but serves as input into a variety of other related ap- by reliable estimates of uncertainty. These limitations reduce the
plications such as studies of fault zone rheology (e.g. Rice & Cocco utility of inferred rupture models and associated by-products. Fur-
thermore, uncertainty in both data and model predictions can cause
current source models to be significantly biased due to overfitting
∗ Now at: Institut de Physique du Globe de Strasbourg, UdS and EOST/CNRS of seismic and geodetic observations.
UMR 7516, France.
Descriptions of earthquake sources come in various flavours
† Now at: Jet Propulsion Laboratory, California Institute of Technology, depending on the nature of the data (e.g. seismological, geode-
Pasadena, California 91109, USA. tic, geological), the observation scale (i.e. regional or global), the

C The Authors 2014. Published by Oxford University Press on behalf of The Royal Astronomical Society. 1
2 Z. Duputel et al.

Figure 1. An illustration of variability of kinematic earthquake source models. Results of finite-source rupture modelling obtained by different research
groups are presented for the 1999 Izmit earthquake (Yagi & Kikuchi 2000; Bouchon et al. 2002; Delouis et al. 2002; Sekiguchi & Iwata 2002). The origin of
XY-coordinates is set at the epicentre location. These source models are available through the SRCMOD source inversion database (Mai 2012).
parametrization of the source (e.g. linear versus non-linear) and the One approach to estimate the uncertainty in source parameters
approach used to infer relevant parameters (e.g. Yabuki & Matsu’ura for a given earthquake is to compare fault slip models obtained
1992; Wald & Heaton 1994; Ji et al. 2002; Minson et al. 2013). The by various research groups using different inversion approaches
reliability of any source inversion depends on many factors includ- (Mai 2012). Fig. 1 shows selected kinematic rupture models for the
ing the size and complexity of the event, the amount and quality 1999 Izmit earthquake. Although these models are generally de-
of data, the way in which data sample the source region and, while rived from similar data sets, there is a large variability in inversion
usually disregarded, uncertainties in our forward models (i.e. our results. The 1999 Izmit earthquake is not an isolated case. For many
model predictions). events, such as the 1992 Landers or 2001 Arequipa earthquakes,
The last decade has seen considerable improvements in the small differences in modelling techniques and data lead to striking
fidelity of forward modelling capability (e.g. Komatitsch & differences in inferred slip models (Wald & Heaton 1994; Hernan-
Vilotte 1998; Williams et al. 2005) and a substantial expan- dez et al. 1999; Pritchard et al. 2007). When different methodologies
sion of geophysical observations including broad-band data from yield different results for the same event, it is not obvious how any
dense seismic networks (e.g. USArray, http://www.usarray.org; conclusion about the rupture process can be drawn.
Geonet, http://geonet.org.nz; CENC, http://www.csndmc.ac.cn; F- This study focus on theoretical and algorithmic developments
net, Okada et al. 2004), continuously recording geodetic positioning needed for the next generation of finite-fault earthquake source
data from permanent GPS installations (e.g. the Plate Boundary models by providing a general formalism to explicitly quantify the
Observatory, http://pbo.unavco.org; Geonet, http://geonet.org.nz, impact of uncertainties in our forward models and to rigorously
Taiwan GPS Network, Yu et al. 1997) and spatially synoptic geode- incorporate such uncertainties in large ill-posed source inversion
tic imaging data from orbiting radar and optical satellites (e.g. Si- problems. We stress the importance of using a stochastic forward
mons & Rosen 2007). Despite this progress in forward modelling modelling approach in this process. It allows us to describe a proba-
and data acquisition, one of the biggest obstacles to significant bility distribution of predictions for a given source model, contrary
progress in earthquake source modelling arises from imperfect pre- to a deterministic approach that provides a single set of (poten-
dictions of geodetic and seismic data due to uncertainties in (or tially inaccurate) predictions. This idea of incorporating stochastic
imperfect knowledge of) the Earth structure—whose impact is gen- (probabilistic) models in the inverse problem is not new and was in-
erally ignored. Indeed, for large earthquakes and even aseismic troduced in geophysics around 1980, notably by Tarantola & Valette
processes, our ability to measure ground motions frequently far ex- (1982). More recently, Yagi & Fukahata (2011) used such a formal-
ceeds our ability to model them. As discussed latter in Section 2 for ism and proposed a stochastic forward model based on adding Gaus-
linear elastic deformation, the prediction errors due to earth model sian noise to the unattenuated 1-D teleseismic Green’s functions.
inaccuracies scale with the fault slip. This aspect is particularly This Gaussian noise is characterized by a covariance matrix that is
important since large events with large amounts of slip will mag- partially specified a priori. Minson et al. (2013) also presented a
nify earth model errors in contrast to small earthquakes for which Gaussian model for the uncertain prediction error in the forward
measurement errors are dominant. Besides the necessity to continue modelling, taking a diagonal covariance matrix with variances that
improving the accuracy and efficiency of forward calculations, one scale with the square of observed amplitudes. In the two approaches,
of the main challenges today is thus to develop an accurate stochas- the scale factor that controls the prediction-error variances is incor-
tic model that better describes modelling uncertainty in predicting porated in the model parameters to be inverted. Based on these early
geodetic and seismic data. studies, we develop here a new formulation exploiting more of the
Accounting for model prediction uncertainty 3
physics of the forward problem to improve the modelling of the for the measurement process is thus given by
prediction covariance matrix. This general formalism can be used
p(d∗ |d) = N (d∗ |d, Cd )
for various problems (e.g. earthquake or volcanic source inversions

based on seismic or geodetic data) and relates input uncertainties in 1 1
the earth model or source geometry to the corresponding distribu- = exp − (d∗ − d)T C−1
d (d ∗
− d) , (2)
(2π ) N |Cd | 2
tion of predictions. As recognized by Yagi & Fukahata (2011) and
Minson et al. (2013), a physically based stochastic model must also where p(d∗ |d) is the probability (density) for getting the measured
account for the dependence of the prediction uncertainty upon the value d∗ when the uncertain physical quantity being measured has
slip model. the value d. The PDF p(d∗ |d) and the associated measurement co-
We begin by developing the concept of the misfit covariance ma- variance matrix Cd depend of course on the nature of measurement
trix as used in inversions of slip on subsurface faults. This matrix is and on the type of instrument used. A common model is to take
the sum of a covariance matrix for observations (often assumed in- independent observational errors (i.e. diagonal Cd ). However, for
dependent) and a covariance matrix for prediction errors associated observations like InSAR or seismic data, off-diagonal components
with inaccurate model predictions (often entirely ignored). We then should be included in Cd to allow correlation of measurement errors
describe how a physically informed prediction covariance matrix between neighbouring data samples (e.g. Lohman & Simons 2005;
can be obtained. In particular, we consider the effect of uncertain- Fukahata & Wright 2008; Duputel et al. 2012a).
ties in the earth model. Our description is based on a Bayesian The second source of uncertainty corresponds to our imperfect
formulation of the inverse problem but our formalism can also be knowledge of d for a given source model m, which comes from the
used in optimization methods. Although our approach is general prediction error due to imperfect forward modelling, also referred
and can be used for various seismic and geodetic data sets, we ex- to as epistemic error. For earthquake source modelling problems,
plore here the advantage of including more structured reasonable this component includes but is not limited to, lack of fidelity in the

prediction covariance matrices by using simple quasi-static models. fault geometry, oversimplifications of the mechanical earth model
We also discuss how to account for other sources of prediction un- and approximations made when calculating the Earth’s response to
certainty such inaccurate as fault geometries. Given the increased an applied force.
computational complexity, we leave a similar development of the Let g(, m) be a deterministic model for the forward predictions
prediction covariance for kinematic modelling to a future study. for a source model m. In addition to the source model, the forward
model depends on a set of uncertain properties parametrized by a
vector that is not solved for (e.g. earth model elastic properties,
fault geometry). We are uncertain about what value to take for the
2 O N T H E I M P O RTA N C E O F parameters . Suppose denotes the most plausible value a priori.
P R E D I C T I O N U N C E RTA I N T Y I N Then, for a given source model m, we obtain the corresponding
S O U RC E I N V E R S I O N P RO B L E M S m). We define the uncertain prediction error
prediction dpred = g(,
by
Let dobs = [d1obs , d2obs , . . . , d Nobs ]T be the set of N field observations
used in the source inversion process. Let also define m as the set m),
= d − dpred = d − g(, (3)
of source model parameters in the M-dimensional model space M.
which is taken as a stochastic variable associated with a probability
In the source inversion process, we want to use dobs to learn about
density p(|m), that describes the uncertainty in the actual physical
m. As in many inverse problems, error models play a central role
quantity d (i.e. displacement), given dpred . Using a similar argument
in the formulation of uncertainties and our knowledge of m will
for as for e above, we choose a maximum entropy distribution
strongly depend on the information they provide. In our problem,
subject to a zero mean and a covariance matrix Cp (m) to get the
we can separate two sources of error that we can account for using
Gaussian distribution p(|m) = N [|0, Cp (m)]. The correspond-
appropriate stochastic models. These two fundamental probability
ing stochastic forward model for the predictions is then given by
models are discussed in the following.
the conditional Gaussian PDF
The first source of uncertainty comes from the observational
1
error that is induced by the measurement process. When measuring p(d|m) = N d|g(, m), Cp (m) =

an N-dimensional quantity d, measurements are associated with an (2π ) Cp (m)
N
uncertain error e given by

∗
1
T −1

e = d − d. (1) × exp − d − g(, m) Cp (m) d − g(, m) .
2
In this equation, d∗ is a stochastic variable representing uncertain (4)
measurement and d is the actual displacement value. The set of field
observations dobs = [d1obs , d2obs , . . . , d Nobs ]T can be seen as a single Even though the prediction error is generally neglected (i.e. Cp = 0),
realization of d∗ . In other words, dobs is a fixed vector corresponding we know that its contribution can be comparable or even larger than
to actual measured values and d∗ is a stochastic vector representing measurement errors, in particular for large earthquakes. This can be
uncertainty on those field observations. We will assume that positive easily understood if we consider, for example, a linear formulation
and negative errors of equal magnitude are equally plausible and of the forward problem g(, m) = G() · m. In this case, we can
so take e to have zero mean. We also choose a covariance matrix write
Cd for e. Following the Principle of Maximum Entropy (Jaynes
= G() − G() · m = (, ) · m, (5)
1983, 2003), the probability density function (PDF) that assumes
the least additional information about e under these conditions is a where the matrix describes the uncertainty in G due to the un-
Gaussian probability density p(e|d) = N (e|0, Cd ) with covariance certainty in . In this equation, the level of the prediction error
matrix Cd and zero mean (cf. Bishop 2006). Our stochastic model modelled by Cp scales with the magnitude of the source model
4 Z. Duputel et al.
m (this dependence is also discussed in Yagi & Fukahata 2011). matrix structure based on an explicit treatment of uncertainties in
The observational uncertainties, on the other hand, are independent the predictions.
of the model parameters and are essentially controlled by the na-
ture and quality of the measurements. For large earthquakes, the
contribution of Cd is thus frequently negligible compared to Cp . 3 A STOCHASTIC MODEL FOR THE
At this point, we therefore have (1) a stochastic model p(d∗ |d) P R E D I C T I O N U N C E RTA I N T Y
associated with Cd in eq. (2) describing the measurement uncer-
tainty and (2) a stochastic forward model p(d|m) associated with The development of a covariance matrix for the predictions (Cp ) is
Cp in eq. (4) describing the prediction uncertainty. Using Bayesian important regardless of the particular approach one uses to invert for
source inversion, our goal here is to combine the available infor- source model parameters. Indeed, in most source inversion problems
mation from observations and prior information about the model (e.g. Hartzell & Heaton 1983; Delouis et al. 2000; Ji et al. 2002;
parameters and forward modelling to construct a posterior distri- Simons et al. 2011; Minson et al. 2013, 2014), the discrepancies
between data dobs and forward predictions g(, m) for a source
bution for the source model parameters. To do so, we use Bayes’
theorem to get the posterior PDF p(m|dobs ) over the model space model m are quantified by defining a least-squares misfit function
(Bayes 1763): of the form
1
m) T · C−1

p(m|dobs ) = κ p(dobs |m) p(m), (6) χ (m) = dobs − g(, χ · dobs − g(, m) . (11)
2
where κ is a normalization constant and p(dobs |m) is a likelihood In a parameter optimization process, χ (m) is minimized while
function: in a Bayesian formulation the likelihood function is given by

p(dobs |m) = η(m)exp ( − χ (m)). Therefore, whatever the source

p(dobs |m) = p(dobs |d) p(d|m) dd.

(7) estimation method, a central role is played by the misfit covariance
D pred
matrix, Cχ = Cd + Cp . As we discussed previously, Cp is the domi-
Strictly speaking, p(dobs |m) is a marginal or integrated likelihood in nant term for large earthquakes and the construction of an improved
which d has been marginalized out over the prediction space Dpred , prediction covariance matrix can lead to improved source imaging
which is the space of all conceivable forward model predictions. techniques. If one adopts a regularized least-squares approach for
Eq. (6) states that the posterior probability density of m given the finite-fault models using a Tikhonov regularization (e.g. minimum
measurements dobs is proportional to the product of a PDF p(m) moment, minimum roughness), then the particular choice of the
describing the prior information on m and the likelihood function penalty parameter depends on the effective information content
p(dobs |m), which gives a measure of how well the model m explains of the data. This dependence on Cχ is true regardless of whether
the data dobs . A similar result can be obtained using the concept one chooses the amount of damping by some L-curve approach,
of conjunction of states of information as introduced by Tarantola cross-validation or Bayesian model class selection. Similarly, if one
& Valette (1982). In eq. (7), p(dobs |d) comes from substituting chooses to adopt a fully Bayesian approach, then having the appro-
d∗ = dobs in the probability model p(d∗ |d). It describes the likelihood priate covariance matrix can control the extent to which one does
of having observed dobs if the actual displacement was d. Using the or does not overfit the observations.
stochastic models developed in eqs (2) and (4) in eqs (6) and (7),
we can then show (cf. Appendix A)

1 m) T
3.1 Prediction uncertainty due to inaccuracies in the
p(m|dobs ) = η(m) p(m) exp − dobs − g(, forward model
2
The derivation proposed here is developing ideas proposed by

× Cχ (m)−1 dobs − g(, m) , (8) Duputel et al. (2012a), who demonstrated that centroid location
uncertainty in point-source moment tensor inversions can be ac-
where Cχ (m) is the misfit covariance matrix defined as counted for by designing a prediction covariance matrix Cp based
on the relevant physics of the forward problem. In this study, we
Cχ (m) = Cd + Cp (m) (9) consider the effects of properties of the forward model which are
and η(m) is a normalization factor: not solved for but can significantly affect the predictions (e.g. Earth
−1/2 −1 elastic properties, fault geometry).
η(m) = κ (2π )−N /2 |Cd |−1/2 Cp (m) C + Cp (m)−1 −1/2 . Hereafter, we derive the statistics of the prediction uncertainty
d
(10) by assuming that g(, m) for given generic properties of the
forward model is well approximated by linearized perturbations of
It is important to note that the covariance Cχ plays a central role in m) for the a priori set of parameters :
our predictions g(,
the inversion process. First of all, Cχ controls the shape of the poste-
g(, m) ≈ g(, m) · ( − ),
m) + K (, (12)
rior probability density p(m|dobs ) and therefore affects the solution
of the inverse problem. Secondly, Cχ is the only term in eq. (8) that m) is the so-called sensitivity kernel of the
where the matrix K (,
describes the statistics of measurement (Cd ) and model prediction predictions with respect to :
errors (Cp ). Among these two contributors, Cd can easily be more
readily taken into account because the statistics of observational un- m) = ∂gi
(K )i j (, (, m). (13)
certainties are generally well characterized. Cp , on the other hand, ∂ j
is usually neglected even if its contribution can be larger than Cd Under the first-order approximation of eq. (12), we can write the
as discussed above. In this paper, we advocate improved modelling prediction error as
of the prediction-error uncertainty in any source inversion prob-
lem by developing a stochastic forward model using a covariance m) ≈ K (,
= g(, m) − g(, m) · δ. (14)
∂d
Using eq. (12) we can write Jacobian rule p() = p(d|m)| ∂ |, model perturbations (e.g. δln = δ / , Tromp et al. 2005). The
∂d
where | ∂ | is the Jacobian determinant of the transformation prediction covariance matrix associated with such uncertainty in the
d = g() and we write the prediction covariance as earth model is then given by

Cp = K · C · K T , (20)
Cp (m) = m)
g(, m) − g(,
where the matrix K is the sensitivity kernel of the predictions with

m) T p() d,
× g(, m) − g(, (15) respect to the earth model parameters:
where p() is the prior probability density describing the uncer- m) = ∂gi (,
(K )i j (, m). (21)
tainty in the generic properties . We assume here that p() is a ∂ ln j
Gaussian distribution: In Section 4, we consider a simple 1-D case for which j represents
C ),
p() = N (|, (16) the shear modulus μj in the jth layer of the tabular elastic model
=
μ used to compute the predictions gi (μ, m). In the 3-D case,
which corresponds to the least informative PDF that is adequate for j = μj can represent the shear modulus in the jth region of the
given a priori parameters and a covariance matrix C defined earth model (cf. Appendix B).
as

C = ( − )( T p() d.

− ) (17) 3.3 Practical implementation in a Bayesian framework
If more information is available about , one can of course choose In this study, we use a Bayesian sampling algorithm called Cascad-
another more informative form of p(). By plugging eqs (14) and ing Adaptive Transitional Metropolis In Parallel (CATMIP), which

(17) into eq. (15), it follows that allows sampling in very high dimensional problems in a parallel
computing framework (Minson et al. 2013). CATMIP combines the
Cp = K · C · KT , (18) Metropolis algorithm with elements of simulated annealing and ge-
where we drop the variables and m for clarity. Using the for- netic algorithms to dynamically optimize the algorithms efficiency
mulation developed above, we can do much better than predict- as it runs. As proposed initially by Beck & Au (2002) and Ching
ing a set of observations for a given source model. We can now & Chen (2007), the CATMIP algorithm samples from a series of
construct a probability distribution for the predictions p(d|m) = intermediate PDFs:
m), Cp ], that is, a stochastic forward model based on
N [d|g(,
f (m, βi ) ∝ p(m) p(dobs |m)β ,
more of the physics of the problem and reflecting uncertainties in
properties which are not inverted but can affect the predictions. i = 1, . . . , B
0 = β0 < β1 < β2 < . . . < β B = 1. (22)
3.2 Accounting for inaccuracies in the earth model
In this procedure, we start at β 0 = 0 to sample the prior and slowly
Among the different sources of uncertainty, the approximate na- increase β i through several transitional steps. When β i = β B = 1,
ture of the chosen earth model is of great interest as it is one of we obtain an accurate sampling of the posterior probability den-
the largest contributors to forward modelling errors (Savage 1987; sity p(m|dobs ). As proposed in Ching & Chen (2007) and Beck &
Hjörleifsdóttir & Ekström 2010; Yagi & Fukahata 2011). We explic- Zuev (2013), we optimize the transitional process by choosing each
itly treat this aspect in constructing a stochastic forward model for d. δβ i = β i + 1 − β i adaptively rather than using a pre-set schedule
In some extreme cases, an earth model can be very different from for the sequence of β i . More specifically, each δβ i is chosen such
the actual elastic structure and there will be no simple relationship that the effective sample size between f(m, β i ) and f(m, β i + 1 ) is
between the corresponding predictions and the actual displacement about 50 per cent (Beck & Zuev 2013). Between each transitional
values. In such situations, results of our inversion would be unreli- step, a large number of Markov chains run in parallel. Each chain
able and we should focus our efforts on obtaining a more accurate is governed by the Metropolis algorithm in which the probability
earth model. However, in most applications, elastic models are rea- of acceptance of one candidate sample is determined by comparing
sonable approximations of the true underlying elastic structure of the value of the intermediate PDF with that of the current sample
the Earth and we may still be able to solve the source inversion (Minson et al. 2013).
problem by designing an appropriate misfit covariance model, Cχ . Our derivation is consistent with ideas proposed initially by
As in previous section, we can derive the statistics of the pre- Tarantola & Valette (1982), who argued that modelling uncertain-
diction uncertainty by assuming that g(, m) for given elastic pa- ties can be incorporated into the inversion framework by adding a
rameters is well approximated by linearized perturbations of our term to the data covariance matrix that represents the uncertainties
predictions g(, m) for the a priori earth model . The elastic in the physical theory (i.e. Cp ). However, we differ from Tarantola
parameters in being strictly positive, we use the formulation de- & Valette (1982) by recognizing that our prediction covariance ma-
scribed in Section 3.1 with = ln . We therefore assume here that trix Cp should also depend on the earthquake source model and
p() is a log-normal distribution which corresponds to the least- is not just a constant matrix, that is, changing the magnitude and
informative PDF (i.e. the maximum entropy PDF) that is adequate distribution of fault slip will change Cp for a given elastic model.
for a Jeffrey’s parameter (Tarantola 2005): For practical implementation, there are different ways of dealing
C ), with the dependence of Cχ = Cd + Cp upon the source model m.
p(ln ) = N (ln | ln , (19)
For example, one can calculate Cp (mprior ) using an a priori source
and C , respectively, the mean and covariance of ln .
with ln model, mprior , such as a centroid-moment-tensor solution or a pre-
This choice of a log-normal distribution is also justified by the liminary finite-fault model and assume that Cp is constant. In con-
fact that modern tomography techniques are often based on relative trast, we propose here to update the prediction covariance during the
6 Z. Duputel et al.
source inversion process. More precisely, in our implementation of 4.1 Calculation of Cp in two dimensions
CATMIP, the sample mean m at each transitional step is used as
In this application example, we want to take into account the pre-
a new model to re-compute Cp . Therefore, we assume that Cp does
diction error due to uncertainty in the 1-D shear modulus structure
not vary significantly in the neighbourhood of m, which ensure
( = μ). In this case, the prediction covariance Cp defined in eq.
that the likelihood term p(dobs |m) in eq. (22) is Gaussian for a given
(20) can be rewritten as
value of β. For computational efficiency, the slip inversion being a
linear problem d = Gm, we can pre-calculate the sensitivity kernels
for each Green’s functions in the matrix G Cp = Kμ · Cμ · KμT . (25)
∂ G ik
(K G )i jk = , (23)
∂ log j
From this equation, we know that in order to obtain Cp , we need
such that to estimate the shear modulus sensitivity kernel Kμ and to choose
an appropriate covariance Cμ describing the uncertainty on μ. For
K = K G · m. (24) practical implementation, we discretized the earth model into 50
small layers from depth of 0 to 5H, where H is the thickness of the
In practice, this approach can be used for static and seismic data. shallow layer. We then compute the sensitivity of the predictions
Although K G can be pre-computed, its calculation remains a with respect to the shear modulus in each layer using the first-order
challenging problem. We propose here to use the perturbation the- perturbations introduced initially by Du et al. (1994) for infinite
ory which has been extensively employed in seismic tomography strike-slip faults. The calculation of Kμ in the quasi-static case is
through the Born approximation (e.g. Marquering et al. 1998; detailed in Appendix B.
Tromp et al. 2005; Virieux & Operto 2009) and has been intro- In this simple implementation of Cp , the covariance matrix Cμ

duced for quasi-static problems by Du et al. (1994) and Cervelli is used to describe the uncertainty and correlations of the shear
et al. (2002). In the motivational example presented in Section 4, modulus for each of the 50 layers used to discretize the earth model.
we use analytical solutions that are available in the quasi-static case Various covariance structures in Cμ can designed depending on
for infinite strike-slip faults embedded in a layered earth model (Du the amount of information available about the earth model. Two
et al. 1994). For non-infinite faults, one can increase the efficiency different forms of Cμ are presented in Figs 2(a) and (c). To illustrate
of the sensitivity kernel calculation by adopting an adjoint formu- the corresponding uncertainty on μ, Figs 2(b) and (d) show 1000
lation. This approach allows a reduction in the computational cost stochastic earth model realizations that are drawn for each form of
when the number of stations is small compared to the number of Cμ (cf. eq. 16). It should be noted that static predictions are only
elements and is popular in various fields such as 3-D seismic tomog- affected by gradients in elastic parameters and are not sensitive to
raphy (Tarantola 1984, 1988; Tromp et al. 2005), inverse problems absolute values in μ. Therefore, in this implementation, we take
in elasticity (Bonnet & Constantinescu 2005) and meteorological the half-space as the reference and allow the shear modulus to vary
studies (Talagrand & Courtier 2007). To further increase the com- in the shallow layers. In Figs 2(a)–(b), we assume that the shallow
putational tractability of the sensitivity kernel calculations, we can layer thickness is known but that the shear modulus contrast with the
also consider a limited number of tectonically parametrized regions homogeneous half-space is uncertain. In this case, we have a sharp
where the uncertainties on elastic parameters are prescribed. Indeed, transition in Cμ between the shallow layer and the homogeneous
we are not interested in assigning variable uncertainties on a fine half-space. In Figs 2(c)–(d), we consider the layer thickness as
mesh but rather describing the earth model uncertainty in a limited uncertain and we assume a smooth transition using a decaying
number of tectonic regions (e.g. crust, mantle, etc.). exponential correlation function for depth larger than 0.8H.
Figs 3(a) and (c) show the corresponding covariance matrices Cp
obtained for a simple uniform unit slip distribution between 0 to
0.9H where H is the thickness of the shallow layer. Figs 3(b) and
4 A P P L I C AT I O N T O G E O D E T I C D ATA :
(d) present 1000 stochastic prediction realizations for both versions
A M O T I VAT I O N A L E X A M P L E
of Cp . Note that Cp is not very different for the two forms of Cμ ,
In this synthetic application, we study the static surface displace- even though the prediction uncertainty for distant stations is slightly
ments due to a vertical strike-slip fault and solve for the depth larger if we include uncertainty in the layer thickness. Several inter-
distribution of slip given uncertainties in the depth distribution of esting features can be identified in the prediction covariance matrix
shear modulus (μ) in the elastic medium. With this purpose, we structure. From the diagonal elements of Cp , we note first an in-
consider a simple model of an infinite 2-D strike-slip fault embed- crease of the prediction uncertainty as one approaches the fault,
ded in a shallow low compliance layer of thickness H overlying an in agreement with the idea proposed by Minson et al. (2013) that
elastic half-space (cf. red line on Figs 2b and d). We define μ2 , mismodelling error should be roughly proportional to the data am-
the shear modulus in the half-space and μ1 , the shear modulus in plitude. However, we also note a drop of prediction uncertainty for
the shallow layer so that μ2 /μ1 = 1.4. Given synthetic data at 100 observation points very close to the fault, since they provide direct
observation points regularly sampled from the surface deformation slip measurements that have little sensitivity to the shear modulus.
of this 2-D model, we infer the depth distribution of slip assuming a The off-diagonal covariance components show a clear anticorrela-
homogeneous elastic half-space and accounting for the uncertainty tion with values on the opposite side of the fault for the form of
in μ using the formulation of Cp described in Section 3.2. For sim- Cμ considered here. If the predicted data happen to be anomalously
plicity, we assume no uncertainty on the Poisson’s ratio although large on one side, the prediction must also be anomalously large but
this can easily be implemented using our formalism of Cp . The 2-D of opposite sign on the other side of the fault. The prediction covari-
calculation of Cp , the fault parametrization as well as practical im- ance Cp thus provides essential information about data weighting
plications of using a non-constant prediction covariance structure and data correlation that should be taken into account in the inverse
are addressed in the following sections. problem.

Figure 2. Uncertainty in the 1-D shear modulus structure μ. The shear modulus uncertainty is presented in (a)–(b) assuming a fixed thickness H for the shallow
layer and in (c)–(d) considering the layer thickness as uncertain. For both cases, we present the shear modulus covariance matrix Cμ in (a) and (c). Using eq.
(16) to sample for each form of Cμ , we show 1000 stochastic earth model realizations in (b) and (d), where individual stochastic realizations are plotted in
black and the actual earth model is presented in red. We define μ1 , the shear modulus in the shallow layers and μ2 , the shear modulus in the half-space.
4.2 A side comment on fault parametrization Np -dimensional Dirichlet distribution produces sets of Np stochastic
positive numbers that sum to a given constrained value (cf. Min-
In this section, we investigate the model resolution in order to design son et al. 2013). The marginal of a Dirichlet distribution is a beta
a proper parametrization of the fault. We aim to derive a fault distribution of the form
discretization that allows us to accurately reflect the slip distribution
and account for the resolving power of available measurements.
β(m i |1, N p − 1) ∝ (1 − m i ) N p −2 , (26)
Furthermore, to allow a natural comparison between the inferred
model and the true slip distribution, we need to understand the
fundamental spatial resolution of the estimated slip model. where we assumed unit concentration parameters (Bishop 2006).
Fig. 4 shows the inversion results obtained for a target slip dis- The best-fitting Beta distributions shown in black in Fig. 4(a) in-
tribution with 1-m uniform slip (i.e. m = 1 m) at depths between 0 dicate that optimum values for Np range between 2 and 4. This
and 0.9H. We use p(m) = U(−0.3, 20) M as the prior information suggests that the sum of slip on Np ∼ 3 neighbouring patches can
on the slip parameters m (i.e. a uniform probability distribution be resolved while it is poorly constrained on individual subfaults.
from −0.3 to 20 m in M-dimensions). As described above, the data Fig. 5 shows some model samples chosen randomly near the mean of
computed for a layered half-space with μ2 /μ1 = 1.4 (cf. red line the posterior distribution. For these models, the distribution of slip
in Fig. 2 d) is inverted using Green’s functions for a homogeneous shows strong oscillations over adjacent patches while the average
half-space including our formulation of Cp (using the shear modulus slip over group of subfaults are consistent with the target model. The
structure Cμ shown in Fig. 2c). The results obtained if we neglect best-fitting Dirichlet distributions along with this checker-boarding
Cp are presented in Fig. S1 of the Supporting Information and the of neighbouring subfaults indicate that the slip can be well resolved
advantages of including Cp in the inversion are discussed latter (see over Np ∼ 3 neighbouring patches, that is, at a scale of about 0.2H.
Section 4.3). Histograms in Fig. 4(b) present marginal PDFs similar to Fig. 4(a)
Histograms in Fig. 4(a) present the marginal PDFs obtained if the with the fault discretized into 32 fault patches. In this case, the best-
fault is discretized into 16 fault patches. Apart for the shallowest fitting Dirichlet distribution indicates checkerboarding over Np ∼
subfault, the shape of these marginals clearly suggests a multivari- 6 patches which corresponds to a resolution scale of about 0.2H,
ate Dirichlet distribution of slip over group of adjacent patches. An consistent with results obtained using the 16 patches discretization.
8 Z. Duputel et al.

Figure 3. Prediction uncertainty in a 2-D quasi-static case. We assume a simple infinite strike-slip fault with unit slip from the free surface to 0.9H, where H
is the thickness of the shallow layer (cf. red line in Fig. 2). Panels (a) and (b) illustrate the prediction uncertainty statistics if we assume a fixed thickness for
the shallow layer (cf. Figs 2a–b). Panels (c) and (d) present the predictive stochastic model obtained if we consider some uncertainty on the layer thickness (cf.
Figs 2c–d). For both cases, we present the corresponding prediction covariance matrices in (a) and (c). For each stochastic model, we show 1000 stochastic
predictions in (b) and (d) in black assuming a homogeneous half-space. The predicted data for the actual layered earth model is shown in red and the predictions
for a homogeneous half-space are shown as yellow dashed lines.
There are several different possibilities to take into account this layered half-space presented in red in Fig 2(d) and add 5 mm of un-
limited resolution. A first possibility is to increase the size of correlated Gaussian observational noise. The resulting data vector
patches, therefore using a coarse fault discretization. This possi- is presented in red in Figs 6(d)–(f). Once again, the source inversion
bility has already been explored by Pritchard et al. (2002) and is performed assuming a homogeneous half-space. We therefore in-
Barnhart & Lohman (2010) who proposed to use variable patch clude the two classes of errors discussed in Section 2—errors on
sizes depending on the resolution scale on the fault. On the other the measurements and in the predictions. As described before, we
hand, one should assign small patch sizes to enhance the accuracy of use a simple uniform prior p(m) = U(−0.3, 20) M on the slip distri-
the forward modelling in order to minimize parametrization errors bution, a fault discretized into 16 patches, and a smoothing window
due to the assumption of constant slip on elements that produces over three patches. The measurement covariance matrix Cd is diag-
sharp discontinuities. A better practice is thus to use a discretization onal with standard deviation of 5 mm. For the calculation of Cp , we
smaller than the actual resolution scale and eventually to filter the consider the first layer thickness as uncertain (i.e. Cμ in Fig. 2c) and
slip distribution a posteriori using a smoothing or averaging length use a prediction covariance structure similar to the one presented in
comparable to the resolution scale. Therefore, we prefer here to use Fig. 2.
16 fault patches and to account for the model resolution scale us- The comparison of the posterior model distribution mean with
ing local averaging rather than using coarse discretization. Fig. 4(c) and without neglecting the prediction uncertainty are shown, re-
shows marginal distributions of filtered model samples using an spectively, on Figs 6(b) and (c). The corresponding 1-D and 2-D
arithmetic mean over a sliding window of Np = 3 patches. These re- marginal posterior PDFs for each fault patch are also presented
sults show a Gaussian-like distributions that are well centred around in Fig. 7 [each p(mi |dobs )] and Fig. 8 [each p(mi , mj |dobs )] based
the target slip value (i.e. m = 1 m), which confirms the possibility on nearly 700 000 samples of the slip vector m. We note signif-
to resolve slip over three to four patches. icant discrepancies in the inversion results depending on whether
the prediction covariance matrix Cp is included or ignored. If Cp is
neglected, the mean of the distribution shown in Fig. 6(b) is very
different from the target model in Fig. 6(a). Estimated slip values
4.3 Comparison of inversion results with and without are larger than the target model for depths between 0.2 and 0.4H
neglecting the prediction uncertainty and significantly lower than the target slip value at larger depth. We
To assess the impact of including Cp in the estimation process, we note also that the uncertainty on the slip distribution is clearly un-
compare inversion results with and without neglecting the predic- derestimated: the marginal PDFs in Figs 7(a) and 8(lower left) show
tion uncertainty. We assume a non-uniform target slip distribution very narrow peaks at large depth that are clearly shifted with respect
from 0 to 0.9H presented in Fig. 6(a). We compute the data for the to the target slip values. On the other hand, if Cp is included, we

Figure 4. Source estimation for a constant unit slip distribution. Raw marginal PDFs for the slip in each fault patch are shown in (a) if the fault is discretized
into 16 patches and in (b) if the fault is discretized into 32 patches. The target slip model is indicated as dashed black lines and the mean of the distribution is
shown in blue. The number of parameters Np of the best-fitting marginal Dirichlet distribution (black curves) is indicated for each patch. The marginal PDFs
for 16 patches after using a moving average over three neighbouring patches is shown in (c).
This correlation among adjacent patches will not affect inversion

results as long as the slip distribution is relatively smooth as in
Fig. 6(a). It can however be problematic if the target model includes
short-scale variations of slip. To explore this possibility, Fig. 9 shows
inversion results for a target model including an abrupt drop of 1 m
in the slip distribution at depth 0.56H (cf. dashed black lines). As
in the case of the smoother slip model in Fig. 7 , the introduction of
Cp significantly improves the posterior model estimates. However,
in Fig. 9(b), we note that the mean model (blue lines) is relatively
smooth and does not reproduce the slip step at a depth of 0.56H.
The step in slip cannot be recovered because in this example we can
only resolve features over a scale corresponding to three neighbour-
ing patches (cf. Section 4.2). Consequently, the mean model shows
better agreement with the filtered target model shown in white in
Fig. 9(b), which corresponds to a three-neighbouring-patches mov-
ing average of the initial target model. In any case, note that both
raw and filtered target models are within the error bounds defined
Figure 5. Examples of model samples chosen near the mean of the posterior by the posterior marginal distribution histograms. The two simple
distribution. The fault from 0 to 0.9H is discretized into 16 fault patches. models presented here clearly illustrate the importance of using an
The depicted models are the four samples closest to the mean of the posterior accurate stochastic model for the prediction uncertainty in order to
distribution in a CATMIP simulation of 8192 samples. reliably infer the posterior distribution of source model parameters.
Using our formulation for Cp , we can obtain more realistic estimates
of the posterior uncertainty but also improve the slip distribution
obtain much broader posterior distributions centred around a mean model.
model that is in agreement with the target slip model. This increase
of estimated posterior uncertainty is particularly visible on the 2-D
4.4 Dependence of prediction uncertainty upon the source
marginal posterior PDFs in Fig. 8. In this figure, we also note a
model
slight correlation between neighbouring patches at large depth (i.e.
for patches P 5–P 16), which is certainly related to the averaging of An important finding of this work is the dependence of the prediction
slip on neighbouring patches as described in the previous section. error upon the source model parameters m. Two critical questions
10 Z. Duputel et al.

Figure 6. Comparison of inversion results with and without neglecting the prediction uncertainty. Panel (a) is the target slip model. The mean of the posterior
distributions are shown in black (b) when the prediction uncertainty is neglected and (c) when Cp is included in the estimation process. The continuous red
line corresponds to the target slip model. Comparisons between data and predictions for each model in a homogeneous half-space are shown in (d)–(f). The
data are presented in red. The predictions shown in yellow are calculated for each model presented in black in (a)–(c). Black lines correspond to a set of 1000
stochastic realizations drawn from the posterior predictive PDF defined in eq. (27). The residuals between data and predictions are shown in (g)–(i). The yellow
line is the difference between the data and the predictions presented, respectively, in red and yellow in (d)–(f). Black lines are the residuals between data and
the 1000 stochastic predictions shown in (d)–(f).
need to be addressed. First, how is it possible to account for the model distribution. The prediction uncertainty covariance matrix
coupling between Cp and the source model? Secondly, what is the Cp is therefore considered here as a by-product of the inversion.
variation of Cp as a function of m? We explored different possibilities to design Cp at β 0 = 0: (1)
As discussed in Section 3.3, we propose here to account for the calculating Cp from the mean of the distribution m at β 0 = 0,
coupling between Cp and m by updating the covariance at each (2) computing Cp for a uniform unitary slip distribution as a func-
transitional step (i.e. increase of β in eq. 22) using the mean of the tion of depth, (3) calculating Cp from the solution of a prior source

Figure 7. 1-D marginal posterior PDFs for each patch as a function of depth. The marginal probability density histograms are shown in green (a) when the
prediction uncertainty is neglected and (b) when the prediction uncertainty is taken into account by including Cp in the inversion problem. The target slip
model is indicated as dashed black lines and the mean of the distribution is shown in blue.
inversion in which Cp is neglected and (4) using a preliminary form

5 DISCUSSION
of Cp whose diagonal elements are proportional to the observations
as suggested by Minson et al. (2014). Fig. 10 presents the evolu- The previous synthetic application clearly illustrates the importance
tion of the prediction covariance in the sampling process using the of incorporating inaccuracies of the earth model in the prediction
mean of the distribution m at β 0 = 0 (cf. Fig. 11). This evolution uncertainty in source estimation problems. When neglecting Cp ,
of Cp corresponds to the estimation process for the target model the posterior mean model shows strong oscillations with a slip
presented in Figs 6–8. Similar figures using approaches (2), (3) and distribution that is overestimated between 0.2 and 0.4H and under-
(4) are shown in Figs S2–S10. For all four cases, we note that Cp estimated at larger depth. We observe similar features for a constant
converges properly to the covariance structure estimated from the unit slip distribution in Fig. S1(b). Note that the negative slip at
target model and that the posterior distributions are almost iden- larger depth is due to the lower bound of the prior distribution
tical. Less transitional steps are however necessary when using a p(m) = U(−0.3, 20) M used to obtain these results. Using strictly
diagonal Cp proportional to the observations at β 0 = 0 (cf. Figs S9 positive constraints leads to similar slip distributions as shown in
and S10). Fig. S11. These oscillatory artefacts when Cp is neglected are clearly
The mean models used to compute Cp are presented in Fig. 11. related to errors in the Green’s functions that are calculated for
The transitional process is clearly illustrated by the evolution of the an homogeneous half-space instead of a layered medium. To ex-
model sample mean shown in this figure. The mean models at the plain these artefacts, Fig. 12 shows the result of a simple linear
first few transitional steps are largely controlled by the prior PDF least-square inversion without positivity constraints (i.e. a purely
[i.e. uniform probability distribution p(m) = U(−0.3, 20) M ] and Gaussian case without prior information). In this example, we con-
then slowly converge to the target model as more weight is given to sider noise-free data but still use incorrect Green’s functions (i.e.
the likelihood (i.e. as β increases). The corresponding covariance assuming an homogeneous half-space instead of the actual layered
matrix structures do not change significantly at each update although half-space with μ2 /μ1 = 1.4). As shown in Fig. 12(a), if we assume
the overall magnitude of Cp is quite different for small and large shallow slip on a fault with large a depth extent, the inferred least-
values of β. Note that an additional sampling step is performed squares solution has strong artefacts at depth due to inaccuracies
after reaching β = 1 in order to resample the posterior PDF with a in the predictions based on a homogeneous elastic half-space. The
prediction covariance estimated from the sample mean at β B = 1. model including Cp in Fig. 12(b) shows similar artefacts but has
However, only minor differences are visible after the final update error bars indicating that deep slip cannot be resolved given the
of Cp . uncertainty in μ. The inferred slip model is actually identical to

Figure 8. 2-D marginal posterior probability densities for all possible pairs of fault patches. Fault patches are numbered increasingly as a function of depth
(e.g. P 16 is the deepest fault patch). The 2-D histograms are shown (bottom left) when Cp is neglected and (top right) when the prediction uncertainty is
accounted for. Dashed black lines indicate the target slip value for each fault patch.
the series of dislocations derived by Savage (1987) and presented prediction uncertainty for the simple model presented in Fig. 6. To
in Fig. 12(c). This slip distribution produces surface displacements do so, we calculate the posterior predictive distribution p(d|dobs ) us-
that are exactly identical to the input target slip model in a layered ing the stochastic forward model p(d|m) in eq. (4) and the posterior
half-space. Of course, in common source inversion practices, faults model distribution p(m|dobs ) in eq. (8):
are parametrized with a limited depth extent and the prediction er-
ror cannot be perfectly mapped in the distribution of slip. Fig. 13 p(d|dobs ) = p(d|m) p(m|dobs ) dm. (27)
illustrates such case with a constant target slip model discretized M
in four slices down to 0.9H. We note strong oscillations when Cp This equation can be obtained directly from the total probability
is neglected. Although there is no obvious analytical solution such theorem (e.g. Bishop 2006) and describes the posterior variability
as Savage (1987) in this case, these oscillations are clearly due to on the predictions d due to modelling error [i.e. p(d|m)] and pos-
inaccuracies of half-space predictions since this is the only source terior uncertainty on the slip distribution [i.e. p(m|dobs )]. The PDF
of uncertainty considered here. On the other hand, the model in- p(d|dobs ) can also be seen as the posterior information on the dis-
cluding Cp is able to properly resolve the input target model since placement field. This posterior predictive PDF is illustrated in black
error in the Green’s function cannot be perfectly reproduced by slip in Figs 6(e)–(f) by showing 1000 stochastic realizations drawn from
oscillations. p(d|dobs ). In Figs 6(h)–(i), we show the corresponding residuals af-
To get more insight on the improvement of inversion results using ter subtracting the data vector from each predictive realization. For
our formulation of Cp , we propose here to estimate the posterior comparison, Figs 6(d) and (g) also show the posterior predictive

Figure 9. 1-D marginal posterior PDFs for a target model including a slip step of 1 m at 0.56H. The marginal probability density histograms are shown in
green (a) when the prediction uncertainty is neglected and (b) when the prediction uncertainty is taken into account by including Cp in the inversion problem.
The mean of the posterior model distribution is shown in blue. The raw target slip model is indicated as dashed black lines and the locally averaged target slip
model using a moving window over three adjacent patches is shown in white.
Figure 10. Evolution of the prediction covariance Cp at each transitional step. The values of β in eq. (22) are specified on top of each subfigure.

Figure 11. Slip models used for the calculation of the prediction covariance matrix Cp . These slip models correspond to the mean of the model sample
distributions at each transitional step. The values of β in eq. (22) are specified on top of each subfigures.
Figure 12. An illustration of the effects of prediction uncertainty on inver-

sion results for an extend fault at depth. The target slip model, indicated as
dashed black lines, is used to generate noise-free data in a layered half-space. Figure 13. An illustration of the importance of the prediction covariance
Simple least-square inversion results for an extended fault in a homogeneous Cp . The target slip model, indicated as dashed black lines, is used to generate
half-space are presented in (a) when Cp is neglected and (b) when Cp is in- noise-free data in a layered half-space. Simple least-square inversion results
cluded in the inversion. The calculation of Cp used in (b) is based on the for an extended fault in a homogeneous half-space are presented in (a) when
maximum a posteriori (MAP) slip model shown in (a). These two MAP Cp is neglected and (b) when Cp is included in the inversion. The calculation
models are very similar to the theoretically derived half-space equivalent of Cp used in (b) is based on the maximum a posteriori slip model shown
slip distribution by Savage (1987) presented in (c). in (a).
PDF for the target slip model, assuming no uncertainty on the slip use the stochastic forward model based on Cp for the same g(m) (i.e.
distribution (i.e. p(m|dobs ) = δ(m − mtrue ) where mtrue is the target p(d|m) = N [d|g(m), Cp ]). The posterior predictive uncertainty is
model). Figs 6(e) and (h) show the data fit if the prediction uncer- significantly larger if Cp is taken into account. Even if the posterior
tainty is neglected [i.e. p(d|m) = δ(d − g(m)), where g(m) are the variability of the source model is incorporated when Cp is neglected,
predictions for a homogeneous half-space]. In Figs 6(f) and (i), we we note that the posterior predictive uncertainty is negligible
compared to the variability when Cp is included. Moreover, when should be spent to improve forward modelling capabilities. The
comparing data with predictions for a homogeneous half-space, the ultimate goal would be to allow for updating the uncertainty in all
data misfit is clearly smaller if the prediction uncertainty is ne- model parameters (i.e. m, , φ) by sampling from their posterior
glected. This small data residual is due to data overfitting resulting PDF but this requires greatly increasing the computational speed for
from the use of a deterministic forward model. By neglecting Cp , the forward modelling, either by computer hardware and/or faster
too much information is conferred to imperfect forward predictions, algorithms (e.g. building surrogate (meta-) models of the forward
leading to spurious posterior distributions which favour source mod- model using machine learning methods; Bishop 2006).
els that closely explain observations. The differences between in-
version results in Fig. 6(b) and the target model in Fig. 6(a) are thus
certainly due to overfitting of observations using predictions for an
6 C O N C LU S I O N
inaccurate earth model (i.e. a homogeneous half-space instead of
a layered medium). On the other hand, using a stochastic forward This study improves the modelling of the misfit covariance ma-
problem, we gain an extra flexibility by allowing some variability trix as used in inversions for the distribution of slip on subsurface
in the predictions, which allows us to select models having larger faults. The misfit covariance, Cχ , is a combination of the observa-
and correlated data residuals. Thus, the impact of including the pre- tional covariance matrix, Cd , and the modelling covariance matrix,
diction covariance Cp is not just a better analysis of the posterior Cp . The latter class of uncertainty is often entirely ignored even
uncertainty. In fact, the use of Cp also improves source parame- though prediction errors scale with the size of earthquakes and
ter estimates by using physically based relative weights between are thus generally larger than the observation uncertainty for large
measurements and by preventing overfitting of observations. events. Furthermore, prediction errors can induce important cor-
In this study, we have mainly focused on the impact of inaccu- relation between the observation points that should be taken into
racies in the earth model, which significantly contribute to forward account in source inversion problems. This work provides a gen-

modelling errors. We recognize, however, that other contributors eral formalism to explicitly quantify the impact of uncertainties
can have profound impact on the predictions. In particular, we know in our forward models and to account for such prediction error in
that approximation of the true source geometry by a fault surface source estimation problems. In particular, we describe a physically

φ based on prior information can be a major source of uncertainty based stochastic model for the prediction uncertainty that allows
in source estimation problems. A judicious approach is to incorpo- for inaccuracies in the earth model by using sensitivity of the data
rate φ into the source model parameters m for which the data are predictions to first-order perturbations of elastic properties. More-
inverted. However, this approach may be difficult to implement in over, to account for the dependence of the prediction error upon the
practice since we would have to recalculate the Green’s functions source model parameters m, the covariance matrix Cp is updated as
for each explored model, a computationally challenging prospect. the model evolves in the source inversion process. We also discuss
Alternatively, if the imperfect fault surface is not grossly wrong, the possibility to account for other sources of uncertainty such as
we can use our formulation of Cp in order to account for the un- imperfections in the fault geometry.
certainty in the fault geometry. One possibility is to design a C For a simple 2-D synthetic data inversion, we demonstrated the
that incorporates a high variability of elastic properties in a narrow importance and feasibility of implementing Cp in source inversion
zone around the fault. Another approach is to take into account a problems. The results indicate two main advantages of using an
variability of φ in the inversion problem. The latter possibility is improved covariance matrix structure. First, it improves the poste-
formulated in the following. rior error description for the source model parameters. Secondly, it
Let us define a parameter = [ln T , φ T ]T incorporating all improves the source model characterization itself notably because
information about the earth model and fault surface φ (which using a more structured covariance matrix prevents overfitting of the
are assumed independent). As in Section 3.1, we assume that pre- observations. Although our implementation is based on a Bayesian
dictions g(, m) are modelled as first-order perturbations of our formulation of the inversion problem, it should be noted that incor-
m). In this case, the prediction covariance can be
predictions g(, porating Cp can be beneficial regardless of the particular inversion
written as approach used. Indeed, the prediction covariance matrix offers a
natural way to specify the relative information content found in dif-
Cp = K · C · K T + Kφ · Cφ · KφT . (28)
ferent observations. Furthermore, if different data types are used,
The prediction uncertainty due to the fault geometry can therefore Cp provides a physical basis for the relative weighting between
be included by augmenting the prediction covariance matrix with a disparate data sets.
new term Kφ · Cφ · KφT , where Cφ is the covariance matrix on fault The advances proposed here can thus enable production of the
m) is the corresponding sensitivity
geometry parameters and Kφ (, next generation of source models that are more resistant to over-
kernel fitting of data, provide a physical basis for the relative weighting
between disparate data sets, and include more realistic description
m) = ∂gi (,
(Kφ )i j (, m). (29) of uncertainty in the inferred source model parameters. These im-
∂φ j
provements in the model prediction uncertainty require developing
In this equation, the vector φ can, for example, include averages of computationally tractable approaches to estimate the sensitivity of
dip and strike values over the entire (or portions of) the fault. This the predicted geodetic and seismic observations of large earthquakes
approach can be extended to other sources of modelling uncertainty to perturbations in material properties assumed in any given earth
and we should be able to define a gallery of prediction covariance model. Much of the underlying theory and many of the tools al-
matrices to be included in the source inversion problem. However, ready exist in the domain of seismic tomography and could thus
this formulation is not appropriate if prediction uncertainties are be exploited to obtain more reliable images of the earthquake rup-
too large since there will be no simple relationship between the ture phenomena. While this work is motivated by specific goals
corresponding predictions g(, m) and g(, m). In such situations, related to the study of large earthquakes, these techniques can also
the results of our inversion would be unreliable and more efforts be applied to a broad range of inverse problems in geophysics and
earthquake engineering such as volcano monitoring and earthquake Duputel, Z., Rivera, L., Fukahata, Y. & Kanamori, H., 2012a. Uncertainty
early warning. estimations for seismic source inversions, Geophys. J. Int., 190(2), 1243–
1256.
Duputel, Z., Rivera, L., Kanamori, H. & Hayes, G., 2012b. W phase source
AC K N OW L E D G E M E N T S inversion for moderate to large earthquakes (1990–2010), Geophys. J.
Int., 189(2), 1125–1147.
We thank Bernard Valette and Yukitoshi Fukahata for their helpful
Fukahata, Y. & Wright, T.J., 2008. A non-linear geodetic data inversion
reviews. We have benefited from discussions with Jean-Paul Am- using ABIC for slip distribution on a fault with an unknown dip angle,
puero, Jean Virieux, Romain Brossier, Victor Tsai, Romain Jolivet Geophys. J. Int., 173(2), 353–364.
and Luis Rivera. MS acknowledges sabbatical support from the Hartzell, S.H. & Heaton, T.H., 1983. Inversion of strong ground motion
Univ. Joseph Fourier (Grenoble, France) during which this project and teleseismic waveform data for the fault rupture history of the 1979
started. PSA was supported by the Keck Institute of Space Studies Imperial Valley, California, earthquake, Bull. seism. Soc. Am., 73(6),
Postdoctoral Fellowship. Our Bayesian sampling algorithm is based 1553–1583.
on the ALTAR implementation of CATMIP developed by Michael Hernandez, B., Cotton, F. & Campillo, M., 1999. Contribution of radar
Aivasis and Hailiang Zhang. This work made use of the Matplotlib interferometry to a two-step inversion of the kinematic process of the
python library created by John D. Hunter. P. Agram was supported 1992 Landers earthquake, J. geophys. Res., 104(B6), 13 083–13 099.
Hjörleifsdóttir, V. & Ekström, G., 2010. Effects of three-dimensional Earth
by the Keck Institute of Space Studies Postdoctoral Fellowship. This
structure on CMT earthquake parameters, Phys. Earth planet. Inter., 179,
research was supported by the Southern California Earthquake Cen- 178–190.
ter. SCEC is funded by NSF Cooperative Agreement EAR-0529922 Jaynes, E.T., 1983. Papers on Probability, Statistics, and Statistical Physics,
and USGS Cooperative Agreement 07HQAG0008. The SCEC con- Springer.
tribution number for this paper is 1791. This work was partially Jaynes, E.T., 2003. Probability Theory: The Logic of Science, Cambridge
supported by the National Science Foundation under Grant No.

Univ. Press.
EAR-0941374. Ji, C., Wald, D.J. & Helmberger, D.V., 2002. Source description of the 1999
Hector Mine, California, earthquake, part I: wavelet domain inversion
theory and resolution analysis, Bull. seism. Soc. Am., 92(4), 1192–1207.
REFERENCES Kanamori, H. & Rivera, L., 2006. Energy partitioning during an earthquake,
Aki, K. & Richards, P., 2002. Quantitative Seismology, Freeman. in Radiated Energy and the Physics of Faulting, Geophys. Monograph
Barnhart, W.D. & Lohman, R.B., 2010. Automated fault model discretization Series, pp. 3–13, American Geophysical Union.
for inversions for coseismic slip distributions, J. geophys. Res., 115(B10), King, G., 2007. Earthquake seismology: 4.08—fault interaction, earthquake
B10419, doi:10.1029/2010JB007545. stress changes, and the evolution of seismicity, in Treatise on Geophysics,
Bayes, M., 1763. An essay towards solving a problem in the doctrine of pp. 225–255, Elsevier.
chances by the late Rev. Mr. Bayes, F.R.S. communicated by Mr. Price, in Komatitsch, D. & Vilotte, J.-P., 1998. The spectral element method: an
a letter to John Canton, A.M.F.R.S., Phil. Trans. R. Soc., 53, 370–418. efficient tool to simulate the seismic response of 2D and 3D geological
Beck, J.L. & Au, S.K., 2002. Bayesian updating of structural models and structures, Bull. seism. Soc. Am., 88, 368–392.
reliability using Markov chain Monte Carlo simulation, J. Eng. Mech., Lohman, R.B. & Simons, M., 2005. Some thoughts on the use of In-
128, 380–391. SAR data to constrain models of surface deformation: noise structure
Beck, J.L. & Zuev, K.M., 2013. Asymptotically independent Markov sam- and data downsampling, Geochem. Geophys. Geosyst., 6(1), Q01007,
pling: a new Markov chain Monte Carlo scheme for Bayesian inference, doi:10.1029/2004GC000841.
Int. J. Uncertainty Quantification, 3(5), 445–474. Love, A.H.E., 1906. A Treatise on the Mathematical Theory of Elasticity,
Bishop, C.M., 2006. Pattern Recognition and Machine Learning, Informa- 2nd edn, Cambridge Univ. Press.
tion Science and Statistics, Springer. Mai, M., 2012. Source inversion validation (SIV): current results & devel-
Bonnet, M. & Constantinescu, A., 2005. Inverse problems in elasticity, opments, in Proceedings of the 2012 SCEC Source Inversion Validation
Inverse Probl., 21(2), R1–R50. Workshop, Hilton Palm Springs Resort.
Bouchon, M., Toksoz, N.N., Karabulut, H., Bouin, M.-P., Dietrich, M., Aktar, Marquering, H., Nolet, G. & Dahlen, F.A., 1998. Three-dimensional wave-
M. & Edie, M., 2002. Space and time evolution of rupture and faulting form sensitivity kernels, Geophys. J. R. astr. Soc., 132(3), 521–534.
during the 1999 Izmit (Turkey) earthquake, Bull. seism. Soc. Am., 92(1), Minson, S.E., Simons, M. & Beck, J.L., 2013. Bayesian inversion for finite
256–266. fault earthquake source models I—theory and algorithm, Geophys. J. Int.,
Cervelli, P., Segall, P., Amelung, F., Garbeil, H., Meertens, C., Owen, S., 194(3), 1701–1726.
Miklius, A. & Lisowski, M., 2002. The 12 September 1999 Upper East Minson, S.E., Simons, M., Beck, J.L., Ortega, F., Jiang, J. & Owen, S.E.,
Rift Zone dike intrusion at Kilauea Volcano, Hawaii, J. geophys. Res., 2014. Bayesian inversion for finite fault earthquake source models II—
107(B7), ECV3 1–13. The 2011 great Tohoku-oki, Japan earthquake, Geophys. J. Int., submitted.
Ching, J. & Chen, Y.C., 2007. Transitional Markov chain Monte Carlo Okada, Y., Kasahara, K., Hori, S. & Obara, K., 2004. Recent progress of
method for Bayesian model updating, model class selection, and model seismic observation networks in Japan-Hi-net, F-net, K-NET and KiK-
averaging, J. Eng. Mech., 133, 816–832. net, Earth Planets Space, 56, 15–28.
Dahlen, F.A. & Baig, A.M., 2002. Fréchet kernels for body-wave amplitudes, Pritchard, M.E., Simons, M., Rosen, P.A., Hensley, S. & Webb, F.H., 2002.
Geophys. J. Int., 150(2), 440–466. Co-seismic slip from the 1995 July 30 Mw= 8.1 Antofagasta, Chile,
Dahlen, F.A. & Tromp, J., 1998. Theoretical Global Seismology, Princeton earthquake as constrained by InSAR and GPS observations, Geophys. J.
Univ. Press. Int., 150(2), 362–376.
Delouis, B., Lundgren, P., Salichon, J. & Giardini, D., 2000. Joint inversion Pritchard, M.E., Norabuena, E.O., Ji, C., Boroschek, R., Comte, D.,
of InSAR and teleseismic data for the slip history of the 1999 Izmit Simons, M., Dixon, T.H. & Rosen, P.A., 2007. Geodetic, teleseis-
(Turkey) earthquake, Geophys. Res. Lett., 27(20), 3389–3392. mic, and strong motion constraints on slip from recent southern
Delouis, B., Giardini, D., Lundgren, P. & Salichon, J., 2002. Joint inversion Peru subduction zone earthquakes, J. geophys. Res., 112(B3), B03307,
of InSAR, GPS, teleseismic, and strong-motion data for the spatial and doi:10.1029/2006JB004294.
temporal distribution of earthquake slip: application to the 1999 Izmit Rice, J.R. & Cocco, M., 2007. Seismic fault rheology and earthquake dynam-
mainshock, Bull. seism. Soc. Am., 92, 278–299. ics, in Tectonic Faults: Agents of Change on a Dynamic Earth (Dahlem
Du, Y., Segall, P. & Gao, H., 1994. Dislocations in inhomogeneous media via Workshop Reports), pp. 99–137, MIT Press.
a moduli perturbation approach: general formulation and two-dimensional Savage, J.C., 1987. Effect of crustal layering upon dislocation modeling,
solutions, J. geophys. Res., 99(B7), 13 767–13 779. J. geophys. Res., 92(B10), 10 595–10 600.
Segall, P., 2010. Earthquake and Volcano Deformation, Princeton Univ. where ν(m) is a normalization factor defined as
Press. −1/2
Sekiguchi, H. & Iwata, T., 2002. Rupture process of the 1999 Kocaeli, ν(m) = (2π )−N |Cd |−1/2 Cp (m) . (A2)
Turkey, earthquake estimated from strong-motion waveforms, Bull. seism.
Soc. Am., 92(1), 300–311.
Eq. (A1) can be rearranged by separating quadratic terms from the
Simons, M. & Rosen, P.A., 2007. Interferometric synthetic aperture radar linear terms

geodesy, in Treatise on Geophysics, pp. 391–446, Elsevier. 1
Simons, M. et al., 2011. The 2011 magnitude 9.0 Tohoku-Oki earthquake: p(dobs |m) = ν exp − dT A d − 2 bT d + c dd, (A3)
D pred 2
mosaicking the megathrust from seconds to centuries, Science, 332(6036),
1421–1425. where we define
Talagrand, O. & Courtier, P., 2007. Variational assimilation of meteorolog- −1
ical observations with the adjoint vorticity equation. I: theory, Q. J. R. A = C−1
d + Cp (m)
Meteorol. Soc., 113(478), 1311–1328.
Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic bT = dobs
T
C−1 T
d + g(, m) Cp (m)
−1
approximation, Geophysics, 49(8), 1259–1266.

Tarantola, A., 1988. Theoretical background for the inversion of seismic c = dobs
T
C−1 T −1
d dobs + g(, m) Cp (m) g(, m). (A4)
waveforms including elasticity and attenuation, Pure appl. geophys.,
128(1–2), 365–399.
The matrix A being positive definite, we can then show

Tarantola, A., 2005. Inverse Problem Theory and Methods for Model Pa- 1
rameter Estimation, SIAM. p(dobs |m) = ν exp − c − bT A−1 b
2
Tarantola, A. & Valette, B., 1982. Inverse problems = quest for information,

J. Geophys., 50, 159–170. 1 T
× exp − d − A−1 b A d − A−1 b dd.

Tromp, J., Tape, C. & Liu, Q., 2005. Seismic tomography, adjoint methods, 2
D pred
time reversal and banana-doughnut kernels, Geophys. J. Int., 160(1), 195–
216. (A5)
Virieux, J. & Operto, S., 2009. An overview of full-waveform inversion in
exploration geophysics, Geophysics, 74, WCC1–WCC26. The
integral of the Gaussian in the right and side being equal to
Wald, D.J. & Heaton, T.H., 1994. Spatial and temporal distribution of slip (2π ) N /|A|, it follows that
for the 1992 Landers, California, earthquake, Bull. seism. Soc. Am., 84,
1
668–691. p(dobs |m) = η exp − c − bT A−1 b (A6)
Williams, C.A., Aagaard, B. & Knepley, M.G., 2005. Development of soft- 2
ware for studying earthquakes across multiple spatial and temporal scales where
by coupling quasi-static and dynamic simulations, EOS, Trans. Am. geo- −1/2 −1
phys. Un., 86, AGU abstract #S53A–1072. η(m) = (2π )−N /2 |Cd |−1/2 Cp (m) C + Cp (m)−1 −1/2 .
d
Yabuki, T. & Matsu’ura, M., 1992. Geodetic data inversion using a Bayesian
information criterion for spatial distribution of fault slip, Geophys. J. Int., (A7)
109, 363–375. After substituting eqs (A4) into (A6) and using the following iden-
Yagi, Y. & Fukahata, Y., 2011. Introduction of uncertainty of Green’s func-
tities (Tarantola & Valette 1982)
tion into waveform inversion for seismic source processes, Geophys. J.
−1 −1 −1 −1
Int., 186(2), 711–720. Cd + Cp (m) = C−1
d − Cd
−1
Cd + Cp (m)−1 Cd
Yagi, Y. & Kikuchi, M., 2000. Source rupture process of the Kocaeli, Turkey,
earthquake of August 17, 1999, obtained by joint inversion of near-field −1
−1 −1
= C−1
p − Cp (m) Cd + Cp (m)−1 Cp (m)−1
data and teleseismic data, Geophys. Res. Lett., 27(13), 1969–1972.
Yu, S.-B., Chen, H.-Y. & Kuo, L.-C., 1997. Velocity field of GPS stations in
−1 −1 −1
the Taiwan area, Tectonophysics, 274, 41–59. = Cp (m)−1 C−1
d + Cp (m) Cd (A8)
we then obtain

1
m) T
p(dobs |m) = η(m) exp − dobs − g(,
2
APPENDIX A: POSTERIOR

P RO B A B I L I T Y D E N S I T Y F U N C T I O N −1
× Cd + Cp m)
dobs − g(, , (A9)
(PDF) UNDER THE ASSUMPTION OF
G AU S S I A N U N C E RTA I N T I E S
which demonstrates eq. (8) of Section 2.
Under the Gaussian assumption, we demonstrate that the posterior
PDF is given by eq. (8) of the main text. We essentially follow the
demonstration originally given by Tarantola & Valette (1982) in the A P P E N D I X B : S E N S I T I V I T Y K E R N E L Kμ
more general case when Cp depends on m. If p(d∗ |d) and p(d|m) are I N T H E Q UA S I - S TAT I C C A S E F O R
assumed to be Gaussian PDFs as in eqs (2) and (4), the integrated G E N E R A L 3 - D P E RT U R B AT I O N S
likelihood in eq. (7) is proportional to the convolution between two We describe here the formulation of Kμ for the case of a known
Gaussians: source and a set of distributed receivers. We adopt a perturbation

1 approach which has been extensively used in seismic (e.g. Dahlen
p(dobs |m) = ν(m) exp − (dobs − d)T C−1 d (dobs − d) & Baig 2002; Virieux & Operto 2009) and quasi-static problems
D pred 2
(Du et al. 1994; Cervelli et al. 2002). Assuming that the earth

+ d − g(, m) T Cp (m)−1 d − g(, m) dd, model is characterized by the elastic stiffness tensor cijkl , for a
given source field f(x) and a set of boundary conditions, the pre-
(A1) dicted displacement field is described by the equation (Love 1906;
Segall 2010) shear modulus sensitivity kernel Jμ of the predicted displacement

field u(x) for the rth layer or tectonic region
∂ ∂u k
ci jkl (x) (x) + f i (x) = 0, (B1)
∂xj ∂ xl (Jμ )ri (x) = − Hi j (x, x

) σ jk (x
) n k (x
) d2 x
. (B11)
Sr
where the Einstein summation convention applies to repeated sub-
script indices. Introducing the Green’s function Hij (x, xs ), the solu- Recognizing the fact that geodetic predictions (dpred ) represent only
tion of this equation can be written in the integral form a subset of the total predicted displacement field, we can write

dpred = R · u, (B12)
u i (x) = Hi j (x, xs ) f j (xs ) d3 xs , (B2)
Vs
where R is a sampling operator acting on the complete displacement
where V s is the volume of the source region. Let us now assume field u(x). Similarly, the integral form in eq. (B10) can be discretized
that a perturbation of the earth model to express the model prediction uncertainty () as
ci jkl (x) → ci jkl (x) + δci jkl (x) (B3) = R · δu
leads to a perturbation of the predicted displacement field: = R · Jμ · δ ln μ. (B13)
u i (x) → u i (x) + δu i (x). (B4) From the definition of in eq. (14) of Section 3.1, we can then write
To derive the displacement sensitivity to the medium elastic prop- the sensitivity kernel as
erties, we insert eqs (B3)–(B4) into eq. (B1) and drop second-order Kμ = R · Jμ . (B14)
terms:

∂ ∂δu k

ci jkl (x) (x) + δ f i (x) = 0, (B5) S U P P O RT I N G I N F O R M AT I O N
∂xj ∂ xl
where δf is defined as Additional Supporting Information may be found in the online ver-
sion of this article :
∂ ∂u k
δ f i (x) = δci jkl (x) (x) . (B6)
∂xj ∂ xl Figure S1. Source estimation for a constant unit slip distribution if
Cp is neglected.
Note that this first-order perturbation equation is equivalent to the
Figure S2. One-dimensional marginal posterior PDFs for each
Born approximation in elastodynamics (Dahlen & Tromp 1998; Aki
patch as a function of depth. The marginal probability density his-
& Richards 2002) where δf is the secondary Born source. Since the
tograms are shown in green (a) when the prediction uncertainty
solution of the eq. (B1) is given by eq. (B2), the solution for δu in
is neglected and (b) when the prediction uncertainty is taken into
eq. (B5) can be expressed as

account by including Cp in the inversion problem. The target slip
∂ ∂u k
model is indicated as dashed black lines and mean of the distribution
δu m (x) = Hmi (x, x
) δci jkl (x
) (x ) d3 x
, (B7)
V ∂xj ∂ xl is shown in blue. The results in (b) are obtained using approach (2)
where V denotes the earth model volume. After some manipulations described in Section 4.4 of the main text: Cp is calculated at β = 0
using the divergence theorem, we obtain assuming a uniform unitary slip distribution as a function of depth

(cf. model shown in Fig. S4 for β = 0).
∂ Hmi (x, x
) ∂u k
Figure S3. Evolution of the prediction covariance Cp at each tran-
δu m (x) = − (x ) δci jkl (x
) d3 x
V ∂ x
j ∂ xl
sitional step. The approach (3) described in Section 4.4 of the main

text is used here: Cp is calculated at β = 0 assuming a uniform
∂u k
+ Hmi (x, x
)δci jkl (x
)
n j d2 x
, (B8) unitary slip distribution as a function of depth (cf. model shown in
S ∂ xl Fig. S4 for β = 0).
where nj is the normal to the Earth’s surface. The second integral Figure S4. Slip models used for the calculation of the prediction
vanishes because of homogeneous boundary conditions. If we as- covariance Cp . The approach (3) described in Section 4.4. of the
sume an isotropic medium, assuming only perturbations in the shear main text is used here: the model used for the calculation of Cp at
modulus (μ) while holding the Poisson’s ratio constant, we can then β = 0 assuming a uniform unitary slip distribution as a function of
write depth (cf. model shown for β = 0).

Figure S5. One-dimensional marginal posterior PDFs for each
∂ Hi j (x, x
) σ jk (x
)
δu i (x) = − δμ(x
) d3 x
, (B9) patch as a function of depth. The marginal probability density his-
V ∂ xk
μ(x
)
tograms are shown in green (a) when the prediction uncertainty
where σ ij (x
) = cijkl (x
)uk, l (x
). For practical implementation, we is neglected and (b) when the prediction uncertainty is taken into
can discretize the elastic medium into a limited number of layers account by including Cp in the inversion problem. The target slip
(as done in Section 4) or tectonically parametrized regions (e.g. model is indicated as dashed black lines and mean of the distribution
crust, mantle). If we allow for piecewise variation of μ in such is shown in blue. The results in (b) are obtained using approach (3)
regions, we can then simplify eq. (B9) using again the divergence described in Section 4.4 of the main text: Cp is calculated at β = 0
theorem to obtain the following 2-D surface integral (Du et al. using the posterior mean model for which Cp is neglected (cf. blue
1994) bars in (a) and Fig. S7).

Figure S6. Evolution of the prediction covariance Cp at each tran-
δu i (x) = − δ ln μ r
Hi j (x, x
) σ jk (x
) n rk (x
) d2 x
, (B10) sitional step. The approach (3) described in Section 4.4. of the main
Sr
r
text is used here: Cp is calculated at β = 0 using the posterior mean
where n rk (x
) is the normal to the surface S r delimiting the rth model for which Cp was neglected (cf. Fig. S7 and blue bars in
perturbed region. From this equation, we can directly extract the Fig. S5a).
Figure S7. Slip models used for the calculation of the prediction diagonal elements are proportional to observations (cf. approach
covariance Cp . The approach (3) described in Section 4.4. of the (4) described in Section 4.4 of the main text).
main text is used here: the model used for the calculation of Cp at Figure S10. Slip models used for the calculation of the prediction
β = 0 is the posterior mean model when Cp is neglected (cf. blue covariance Cp . We use a preliminary form of Cp whose diagonal el-
bars in Fig. S5a). ements are proportional to observations (cf. approach (4) described
Figure S8. One-dimensional marginal posterior PDFs for each in Section 4.4 of the main text).
patch as a function of depth. The marginal probability density his- Figure S11. Same as Fig. 7 in the main text but using strictly
tograms are shown in green (a) when the prediction uncertainty positive constraints with a prior p(m) = U(0, 20) M (http://gji.
is neglected and (b) when the prediction uncertainty is taken into oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt517/-/DC1).
account by including Cp in the inversion problem. The target slip
model is indicated as dashed black lines and mean of the distribution
is shown in blue. The results in (b) are obtained using approach (4)
described in Section 4.4 of the main text: We use a preliminary form Please note: Oxford University Press are not responsible for the
of Cp whose diagonal elements are proportional to observations. content or functionality of any supporting materials supplied by
Figure S9. Evolution of the prediction covariance Cp at each tran- the authors. Any queries (other than missing material) should be
sitional step. At β = 0, we use a preliminary form of Cp whose directed to the corresponding author for the article.

Accounting For Prediction Uncertainty When Inferring Subsurface Fault Slip

Uploaded by

Copyright:

Available Formats

Accounting For Prediction Uncertainty When Inferring Subsurface Fault Slip

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Accounting For Prediction Uncertainty When Inferring Subsurface Fault Slip

Uploaded by

Copyright:

Available Formats

Geophysical Journal International

Geophys. J. Int. (2014) doi: 10.1093/gji/ggt517

Accounting for prediction uncertainty when inferring subsurface

Zacharie Duputel,1,∗ Piyush S. Agram,1,† Mark Simons,1 Sarah E. Minson1

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

2007), earthquake energy budget analysis (e.g. Kanamori & Rivera

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

p(dobs |m) = η(m)exp ( − χ (m)). Therefore, whatever the source

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

C = ( − )(   T p() d.

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

This correlation among adjacent patches will not affect inversion

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

inversion in which Cp is neglected and (4) using a preliminary form

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Figure 12. An illustration of the effects of prediction uncertainty on inver-

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

approximation, Geophysics, 49(8), 1259–1266.

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Segall 2010) shear modulus sensitivity kernel Jμ of the predicted displacement

∂xj ∂ xl (Jμ )ri (x) = − Hi j (x, x

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

Downloaded from http://gji.oxfordjournals.org/ by guest on January 30, 2014

You might also like

C = ( − )( T p() d.