# Uncertainty estimations for seismic source inversions

Zacharie Duputel, Luis Rivera, Yukitoshi Fukahata, Hiroo Kanamori

► To cite this version:

Zacharie Duputel, Luis Rivera, Yukitoshi Fukahata, Hiroo Kanamori. Uncertainty estimations for seismic source inversions. Geophysical Journal International, 2012, 190, pp.1243-1256. hal-00839956

HAL Id: hal-00839956

https://hal.science/hal-00839956

Submitted on 9 Nov 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
# Uncertainty estimations for seismic source inversions

Zacharie Duputel,<fnref n="1" /> Luis Rivera,<fnref n="2" /> Yukitoshi Fukahata<fnref n="3" /> and Hiroo Kanamori<fnref n="1" />

Accepted 2012 May 21. Received 2012 May 21; in original form 2011 December 1

## SUMMARY

Source inversion is a widely used practice in seismology. Magnitudes, moment tensors, slip distributions are now routinely calculated and disseminated whenever an earthquake occurs. The accuracy of such models depends on many aspects like the event magnitude, the data coverage and the data quality (instrument response, isolation, timing, etc.). Here, like in any observational problem, the error estimation should be part of the solution. It is however very rare to find a source inversion algorithm which includes realistic error analyses, and the solutions are often given without any estimates of uncertainties. Our goal here is to stress the importance of such estimation and to explore different techniques aimed at achieving such analyses. In this perspective, we use the W phase source inversion algorithm recently developed to provide fast CMT estimations for large earthquakes. We focus in particular on the linear-inverse problem of estimating the moment tensor components at a given source location. We assume that the initial probability densities can be modelled by Gaussian distributions. Formally, we can separate two sources of error which generally contribute to the model parameter uncertainties. The first source of uncertainty is the error introduced by the more or less imperfect data. This is carried by the covariance matrix for the data (C<sub>d</sub>). The second source of uncertainty, often overlooked, is associated with modelling error or mismodelling. This is represented by the covariance matrix on the theory, C<sub>T</sub>. Among the different sources of mismodelling, we focus here on the modelling error associated with the mislocation of the centroid position. Both C<sub>d</sub> and C<sub>T</sub> describe probability densities in the data space and it is well known that it is in fact C<sub>D</sub> = C<sub>d</sub> + C<sub>T</sub> that should be included into the error propagation process. In source inversion problems, like in many other fields of geophysics, the data covariance (C<sub>D</sub>) is often considered as diagonal or even proportional to the identity matrix. In this work, we demonstrate the importance of using a more realistic form for C<sub>D</sub>. If we incorporate accurate covariance components during the inversion process, it refines the posterior error estimates but also improves the solution itself. We discuss these issues using several synthetic tests and by applying the W phase source inversion algorithm to several large earthquakes such as the recent 2011 Tohoku-oki earthquake.

**Key words:** Time-series analysis; Inverse theory; Earthquake source observations; Surface waves and free oscillations.

## 1 INTRODUCTION

The estimation of the source parameters is a first step to understand the rupture process of large earthquakes. It is also of great interest to study the relationship between the earthquake and its tectonic and geodynamic environment. The inverse problem can be formulated and solved in various ways depending on the nature of data (e.g. seismological, geodetic), the observation scale (e.g. regional, teleseismic) and the time at which it is performed after the event origin time (i.e. ranging from real time to the study of historical earthquakes). The estimated source models can then be used as in-

puts of various algorithms such as Shakemap computation (Wald et al. 2005), tsunami modelling (Satake 2007) or Coulomb stress transfer calculation (King 2007).

Despite their importance, these source inversion results are often lacking of uncertainty estimations and the inversion algorithms themselves generally do not include realistic error analyses. The importance of having such estimates has been earlier stressed by Dziewonski et al. (1981) and significant efforts have been made in geodetic finite fault inversion studies (Yabuki & Matsu'ura 1992; Fukahata & Wright 2008; Sudhaus & Jónsson 2009). It is however rare to find a source inversion solution based on seismological

<footnote n="1">¹Seismological Laboratory, California Institute of Technology, Pasadena, CA, USA. E-mail: zacharie@gps.caltech.edu</footnote>

<footnote n="2">²Institut de Physique du Globe de Strasbourg, IPGS - UMR 7516, CNRS and Université de Strasbourg (EOST), France</footnote>

<footnote n="3">³Disaster Prevention Research Institute, Kyoto University, Japan</footnote>
data including such error analysis (e.g. Ide et al. 1996; Yagi & Fukahata 2011). Our goal here is to discuss how to take errors explicitly into account in seismic source inversion problems.

In this perspective, we use the W phase source inversion algorithm developed by Kanamori & Rivera (2008). As discussed by Duputel et al. (2012), this algorithm provides fast and robust centroid moment tensor estimations for moderate to large earthquakes. Besides the standard application, its simplicity and versatility make it suitable for diverse applications. We incorporate a more formal linearized error analysis into the algorithm and discuss the above by applying it to several large earthquakes.

## 2 FORMULATION OF THE SOURCE INVERSION PROBLEM

We consider here the simple linear case of a point source inversion at a given centroid location. The elements available to solve this problem are (1) the observables **d** (the data **d**<sub>obs</sub> corresponding to measurements of **d**) and (2) the theory **d** = **G****m** which relates the observables **d** to a model **m** using a linear operator **G**. In case of W phase inversion, **m** contains the moment tensor elements, **d**<sub>obs</sub> are the displacement traces and **G** are the so called Green's functions.

The source inversion problem is formulated here by using the least-squares criterion using a Bayesian formulation as proposed for example, by Tarantola & Valette (1982) or Yabuki & Matsu'ura (1992). We assume that the initial probability densities can be modelled by Gaussian distributions. As a consequence of the forward problem being linear, the posterior errors on the model parameters are ensured to be Gaussian. Moreover, the covariance matrices play a central role in the formulation of uncertainties and the inversion results will strongly depend on the information they provide.

In our problem, we can separate two sources of error. On one side we have the uncertainty introduced by imperfect data. This information is provided by the probability density $\rho_D(\mathbf{d})$ which is defined by the actual observations $\mathbf{d}_{\text{obs}}$ and the covariance matrix $\mathbf{C}_d$

$$ \rho_D(\mathbf{d}) = ((2\pi)^N \det \mathbf{C}_d)^{-1/2} \exp\left(-\frac{1}{2}(\mathbf{d} - \mathbf{d}_{\text{obs}})' \mathbf{C}_d^{-1} (\mathbf{d} - \mathbf{d}_{\text{obs}})\right), \quad (1) $$

where $N$ is the total number of data samples.

The second source of error, often ignored, is associated with the modelling uncertainties or mismodelling. We assume that for a given source model **m**, instead of predicting exactly the data **d**, we have a fuzzy theory described by a conditional probability density

$$ \rho_T(\mathbf{d}|\mathbf{m}) = ((2\pi)^N \det \mathbf{C}_T)^{-1/2} \exp\left(-\frac{1}{2}(\mathbf{d} - \mathbf{G}\mathbf{m})'\mathbf{C}_T^{-1}(\mathbf{d} - \mathbf{G}\mathbf{m})\right). \quad (2) $$

associated with a covariance matrix $\mathbf{C}_T$. The conditional probability density $\rho_T(\mathbf{d}|\mathbf{m})$ gives us the probability of obtaining the data vector **d** given the model **m**. $\rho_T(\mathbf{d}|\mathbf{m})$ is introduced, for example, by Tarantola & Valette (1982) and Yagi & Fukahata (2008).

Following Tarantola & Valette (1982) and Tarantola (2005), if we take the conjunction of the two states of information described in eqs (1)–(2) and integrate over the data space, we obtain a marginal probability density $\sigma_M(\mathbf{m})$ which is the solution of the inverse problem describing a *posteriori* information in the model space

$$ \sigma_M(\mathbf{m}) = k \exp\left(-\frac{1}{2}(\mathbf{G}\mathbf{m} - \mathbf{d}_{\text{obs}})'\mathbf{C}_D^{-1}(\mathbf{G}\mathbf{m} - \mathbf{d}_{\text{obs}})\right), \quad (3) $$

where $k$ is a normalization factor and

$$ \mathbf{C}_D = \mathbf{C}_T + \mathbf{C}_d. \quad (4) $$

This shows that, under the Gaussian assumption, the mismodelling and observational uncertainties are combined by simply adding the corresponding covariance matrices. $\sigma_M(\mathbf{m})$ in eq. (3) can be rewritten explicitly as a Gaussian distribution (Tarantola 2005)

$$ \sigma_M(\mathbf{m}) = ((2\pi)^N \det \tilde{\mathbf{C}}_M)^{-1/2} \exp\left(-\frac{1}{2}(\mathbf{m} - \tilde{\mathbf{m}})'\tilde{\mathbf{C}}_M^{-1}(\mathbf{m} - \tilde{\mathbf{m}})\right), \quad (5) $$

where $\tilde{\mathbf{m}}$ is the point at the maximum of the posterior Gaussian $\sigma_M(\mathbf{m})$ and $\tilde{\mathbf{C}}_M$ is the posterior covariance matrix

$$ \tilde{\mathbf{m}} = \tilde{\mathbf{C}}_M G' \mathbf{C}_D^{-1} \mathbf{d}_{\text{obs}}, \quad \tilde{\mathbf{C}}_M = (G' \mathbf{C}_D^{-1} G)^{-1}. \quad (6) $$

Eq. (6) is a particular case of the result found by Jackson (1979) and Tarantola (2005) without prior information on the model parameters.

In source inversion problems, like in many other field of geophysics, a common assumption is to consider the data covariance ($\mathbf{C}_D$) as diagonal or even proportional to the identity matrix ($\mathbf{I}$). In this work, we demonstrate the importance of using a more realistic form for $\mathbf{C}_D$. In this paper, we focus on how to estimate the observational covariance ($\mathbf{C}_d$) and the mismodelling covariance ($\mathbf{C}_T$) and discuss its implications for the solution of the inverse problem.

## 3 OBSERVATIONAL ERRORS

In the case of seismological observations at long period, the data uncertainty is mostly related to the background seismic noise which steadily grows at long period. As indicated by Sorrells (1971), the long period noise is mainly related to atmospheric pressure disturbances. The noise level varies depending on many factors such as the station location and the quality of the instrument isolation. It depends also on the orientation of the record since the horizontal components are often noisier than the vertical ones.

To assess the effect of long period noise on W phase solutions, one possibility is to add actual noise to synthetic waveforms and to make an inversion using the noisy synthetic data set. For this purpose, we use the following procedure: the synthetic seismograms are first computed by normal mode summation and are convolved with the instrument response at each station. We then add raw noise data to the resulting signals. We use the continuous noise records of networks II, IC, IU and G which have been extracted using the Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC). Considering the network geometry and focal mechanism depicted in Fig. 1, we computed the noisy synthetic data sets corresponding to different moment magnitudes (i.e. $M_w = \{6.0, 7.0, 8.0\}$). The moment rate function is assumed here to be an isosceles triangle with a half duration given by

$$ h_c = 1.05 \times 10^{-8} \times M_0^{1/3}, \quad (7) $$

where $M_0$ is in dyn-cm and $h_c$ in seconds. This scaling law is obtained empirically by fitting the scalar moment $M_0$ to $h_c$ for all $M_w \ge 6.0$ in the Global CMT (GCMT) catalogue since 2003 (Ekström et al. 2005). The noisy synthetic data set thus obtained is then used as input of the W phase algorithm.

In the standard W phase centroid moment tensor (WCMT) algorithm, we assume that $\mathbf{C}_d = \sigma^2\mathbf{I}$ where $\mathbf{I}$ is the identity matrix and $\sigma$ is an estimate of the data error. As shown in Fig. 2, if we filter the data at very long period using a fixed 1–5 mHz passband, the signal to noise ratio clearly decreases as the event gets smaller and the match between the WCMT solution and the actual source model gets worse accordingly. As illustrated in Fig. 3, to improve the W phase solutions for smaller events, one possibility is to increase
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Source-station geometry assumed for the synthetic experiment. The star indicates the centroid location and the blue triangles indicate the station locations. The focal mechanism used to compute the synthetics is shown in green.</figcaption></figure>

the signal-to-noise ratio by changing the frequency passband and by performing an efficient data screening to reject the excessively noisy stations. Although this approach is robust and provides very good results at different scales (regional and teleseismic) in real-

time or in post-mortem studies (Hayes et al. 2009; Duputel et al. 2012), it has two main disadvantages. The first one is that it needs a preliminary magnitude to choose the frequency passband used to filter the data. The second is that it doesn't provide accurate error estimates on the source model parameters. This is illustrated in Fig. 4 which shows the posterior uncertainty obtained if we assume a large noise level of 70 per cent between 1 and 5 mHz in the case of the $M_w = 7.0$ event presented in Fig. 2(b) (therefore $\sigma = 0.7 \times \|d\| \approx 0.01$ mm). The error estimates of the model parameters and the correlation between them are presented in Fig. 4(b) where the uncertainty is indicated by showing the nodal planes corresponding to a population of 1000 source models drawn from the posterior probability density given in eq. (5). Although $\sigma \approx 0.01$ mm clearly overestimates the noise level on the data traces (see Fig. 2b), we note that the posterior uncertainty is clearly underestimated since the actual model does not fall within the error bars.

We can improve the solution by taking into account the data uncertainty for each station instead of assuming $C_d = \sigma^2 I$. This approach should allow us to improve the solution and posterior error estimates while using the complete data set and a fixed frequency band (no preliminary magnitude needed). As a first guess, we can consider a diagonal $C_d$ whose elements correspond to a measurement of the noise level $\sigma_d^n$ at each station $n$

$$ (C_d^n)^{ij} = (\sigma_d^n)^2 \delta^{ij}. \qquad (8) $$

In this study, $\sigma_d^n$ is estimated by measuring the pre-event noise level. The results obtained by considering this diagonal $C_d$ are shown in Fig. 5 in the case of the $M_w = 7.0$ event presented in Fig. 2(b) after a bandpass filtering with the 1–5 mHz passband. The solution in Fig. 5(b) obtained by considering the diagonal $C_d$ in eq. (8) shows a better match to the actual model than the solution

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Synthetic experiment using a fixed 1–5 mHz passband. We calculated three noisy synthetic data sets corresponding to different event magnitudes: (a) $M_w = 6.0$, (b) $M_w = 7.0$ and (c) $M_w = 8.0$. The noisy traces are obtained by adding ambient seismic noise to synthetic seismograms which are calculated for the focal mechanism indicated in green. The W phase CMT (WCMT) solution obtained for each data set after bandpass filtering in the 1–5 mHz passband is shown in red. Examples of noisy synthetic traces (black lines) and the corresponding noise free synthetics (green lines) are presented. The noise free synthetics have been slightly shifted up to distinguish the two waveforms. For the $M_w = 8.0$ earthquake, noisy traces are not much different from pure synthetics for a majority of stations. For smaller events, the data is clearly contaminated by the long period noise and the WCMT solutions are significantly different from the actual model.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Synthetic experiment using a magnitude dependent frequency passband and a data screening. We calculated three noisy synthetic data sets corresponding to different event magnitudes: (a) $M_w = 6.0$, (b) $M_w = 7.0$ and (c) $M_w = 8.0$. The noisy traces are obtained by adding ambient seismic noise to synthetic seismograms which are calculated for the focal mechanism indicated in green. For each data set, the W phase CMT (WCMT) solution obtained after bandpass filtering and data screening is shown in red. The passband is shifted toward higher frequencies for smaller events to reduce the long period noise contamination. A data screening is performed to reject the noisy stations from the data set. Examples of noisy synthetic traces (black lines) and the corresponding noise free synthetics (green lines) are presented. The noise free synthetics have been slightly shifted up to distinguish the two waveforms. The WCMT solutions are very similar to actual source models. Note however that by using a high frequency bandpass for small events (i.e. a 7–20 mHz passband for $M_w \le 6.0$), we are moving away from the standard W phase inversion since the effect of shallow heterogeneities may become non-negligible when using actual data filtered at such short periods.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Results of the WCMT inversion with the assumption that $C_d = \sigma^2 I$ where $\sigma$ is taken simply as $\sigma = 0.70 * \|d\| \approx 0.1$ mm in the case of the $M_w = 7.0$ event presented in Fig. 2(b). We show in (a) the WCMT solution and the actual source model for comparison. The posterior uncertainty associated with this solution is presented in (b). We superimposed on the WCMT mechanism a population of 1000 nodal planes drawn from the posterior probability density. The WCMT significantly differs from the actual solution and the posterior errors on the moment tensor elements are obviously underestimated.</figcaption></figure>

in Fig. 4 obtained with $C_d = \sigma^2 I$. This is not surprising since using eq. (8) is equivalent to weight the data as a function of noise level at each station. However, as observed in Fig. 4(b), there is an obvious underestimation of uncertainty in Fig. 5(c). Because the background seismic noise is the only source of error introduced in this synthetic experiment, some information is clearly missing in the present formulation of $C_d$ given by eq. (8).

## 4 DATA OVERSAMPLING: THE IMPORTANCE OF COVARIANCE OFF-DIAGONAL TERMS

As discussed in the previous section, the diagonal form of $C_d$ in eq. (8) leads to an obvious under-estimation of posterior uncertainties. To identify what is missing in this formulation, we can take a closer look at the data. Fig. 6 shows the BFO vertical displacement seismograms which have been recorded during the 2004 Sumatra-Andaman Islands earthquake. The frequency range used here is the 1–5 mHz passband which is used for the determination of the source mechanism of large events using W phase waveforms. Since the WCMT algorithm uses ‘LH’ channels sampled at 1 Hz, the data traces are clearly oversampled. This is demonstrated on Fig. 6 by drawing red dots every 10 samples. The data points are strongly correlated and, under such circumstances, we cannot neglect the interdependence of observational errors (Fukahata & Wright 2008). Neglecting these correlations by assuming a diagonal $C_d$ can seriously bias the error estimates of WCMT inversions. Since we are working at teleseismic distances, the interstation distance is large enough to neglect any correlation due to spatial oversampling.

There are basically two different ways to address this issue. The first approach is to significantly decimate the data so that the hypothesis of independent data samples (i.e. diagonal $C_d$) becomes approximately valid. However, since a battery of low pass filters should be applied to reduce the data sampling, this will add a significant delay which is not suited for W phase fast source inversions, even without considering timing complications. Moreover, this first approach results in a reduction of information originally included in
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Results of the WCMT inversion which takes into account the background noise level at each station in the case of the $M_w = 7.0$ event presented in Fig. 2(b). The blue mechanism in (b) corresponds to the WCMT solution obtained by assuming a $C_d$ whose diagonal elements are determined by measuring the pre-event noise levels for each trace [i.e. the diagonal $C_d$ formulated in eq. (8)]. The posterior uncertainty associated with this solution is presented in (c). We superimposed on the WCMT mechanism a population of 1000 nodal planes drawn from the posterior probability density. For comparison, we show in (a) the actual source model and the WCMT solution obtained using the standard W phase algorithm in which we assume $C_d = I$ (cf. Fig. 2b). Although the WCMT solution in (b) is closer to the actual mechanism than the WCMT solution in (a), the posterior errors on the moment tensor elements are obviously underestimated.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** W phase displacement recorded on LHZ BFO channel during the 2004 Sumatra-Andaman Islands earthquake. After removal of instrument response, the signal is bandpass filtered using a 1–5 mHz passband. Red circles are plotted every 10 samples (i.e. every 10 s) to illustrate the data oversampling.</figcaption></figure>

the data and is not desirable in source inversion problem. The other possibility is to drop the hypothesis of a diagonal $C_d$ by including the data correlation in the construction of $C_d$, which leads to non-diagonal terms. This approach is explored here by introducing a block-diagonal $C_d$ which is given by the following expression for a

station n:

$$ (C_d^n)^{ij} = (\sigma_d^n)^2 \exp(-|\Delta t^{ij}|/t_0), \quad (9) $$

where $t_0$ represents a characteristic correlation duration and $\Delta t^{ij}$ equals the time difference between samples i and j. A similar expression is used by Tarantola & Valette (1982) and Fukahata & Wright (2008) to take into account spatial correlation of error. In practice, $t_0$ in eq. (9) is chosen as the shortest period content available after filtering the data (i.e. in our case we assume $t_0 = 200$ s). An alternative approach to this parametric representation of the data correlation is to estimate the covariance matrix directly from background noise data samples (as suggested by Gouveia & Scales 1998; Sambridge 1999; Tarantola 2005). However, for the very long periods considered here, such calculation would necessitate very long time-series to have a statistical significance.

Fig. 7 shows the results of employing eq. (9) for an inversion based on the noisy synthetic data set generated for a $M_w = 7.0$ earthquake filtered with the 1–5 mHz passband (cf. Fig. 2b). As done in Figs 4 and 5, the posterior error are indicated in Fig. 7(c) by showing a random population of probable nodal planes. The posterior uncertainty shown here is much more realistic than those obtained when neglecting temporal correlations (e.g. Figs 5b or c). Interestingly, taking into account the off-diagonal covariance components not only enhances the posterior errors estimates but it also improves the solution itself since we note a pretty good match between the solution and the original model. By using oversampled seismograms at long period, the information from the data is

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Results of the WCMT inversion which takes into account (1) the background noise level at each station and (2) the oversampling of the W phase traces in the case of the $M_w = 7.0$ event presented in Fig. 2(b). The blue mechanism in (b) corresponds to the WCMT solution obtained by assuming a block-diagonal $C_d$ as formulated in eq. (9). The posterior uncertainty associated with this solution is presented in (c). We superimposed on the WCMT mechanism a population of 1000 nodal planes drawn from the posterior probability density. For comparison, we show in (a) the actual source model and the WCMT solution obtained using the standard W phase algorithm in which we assume $C_d = I$ (cf. Fig. 2b).</figcaption></figure>
artificially overweighted. If we perform the source inversion while neglecting the correlation of errors in oversampled traces, a model that excessively fits observed data would be selected. On the contrary, if we consider the effect of data covariance, the correlation of errors are reasonably taken into account which avoids biased inversion results.

## 5 MODELLING UNCERTAINTIES

In source inversion practices, we often know that our modelling is imperfect. Source mislocation, unmodelled finiteness, incorrect prior fault geometry, oversimplified Earth models are among many sources of modelling errors having a potential effect on source inversion results. To introduce this information into the inverse problem, one possibility is to allow for a non-deterministic theory incorporating such fuzziness. As discussed in Section 2, this can be done by defining a probability density $\rho_T(\mathbf{d}|\mathbf{m})$ describing the theoretical relation between model parameters $\mathbf{m}$ and data $\mathbf{d}$ as well as modelling uncertainties. In the literature, the mismodelling is often neglected, which is roughly equivalent to assuming a delta probability distribution

$$ \rho_T(\mathbf{d}|\mathbf{m}) = \delta(\mathbf{d} - \mathbf{Gm}), \quad (10) $$

where $\delta$ is the Dirac delta function. Although this hypothesis can be true in some particular cases, eq. (10) does not generally hold. In most applications, a significant portion of the data misfit can actually be attributed to modelling uncertainties due to inadequate parametrization or to imperfect physical theories with (too) simplistic hypotheses.

In the particular case of the WCMT inversion, the mismodelling can be caused by the centroid mislocation, the source finiteness for very large earthquakes ($M_w \ge 9.0$, e.g. 2004 Sumatra-Andaman Island earthquake), the source complexity (e.g. 2009 Samoa great earthquake sequence, Beavan et al. 2010; Lay et al. 2010) or the large amplitude disturbances caused by a preceding event (e.g. 2009 Vanuatu earthquake sequence). Among these different sources of mismodelling, we focus here on the modelling error associated with the mislocation of the centroid position. Another source of mismodelling which can be taken into account is the possibility of having an incorrect Earth model. However, the W phase mainly propagates through the mantle and thus barely affected by shallow crustal structures. As discussed by Duputel et al. (2012), the effect of directivity due to the source finiteness is visible on the W phase recordings for $M_w > 9.0$ earthquakes (e.g. for the Sumatra-Andaman 2004 earthquake, the amplitudes are enhanced in the direction of the rupture propagation). However, the W phase is not as affected as traditional Rayleigh waves because of its high group velocity and long period.

Let us then assume an erroneous centroid location $\bar{\mathbf{x}}$, knowing, for example, that this location is wrong mainly in some direction. This information is provided by the probability density $\rho_x(\mathbf{x})$ and the covariance $C_x$ describing the uncertainty in the centroid position $\mathbf{x}$

$$ \rho_x(\mathbf{x}) = ((2\pi)^N \det C_x)^{-1/2} \exp\left(-\frac{1}{2}(\mathbf{x} - \bar{\mathbf{x}})^T C_x^{-1}(\mathbf{x} - \bar{\mathbf{x}})\right). \quad (11) $$

In practice, the 3 × 3 matrix $C_x$ carries the uncertainty in the centroid location. Let also $\mathbf{d} = \mathbf{d}(\mathbf{x})$ represent the relation between the centroid position $\mathbf{x}$ and the corresponding predicted data vector $\mathbf{d}$. We assume now that the mismodelling $\rho_T(\mathbf{d}|\mathbf{m})$ in eq. (2) is entirely attributable to the centroid mislocation provided by $\rho_x(\mathbf{x})$. Using the Jacobian rule $\rho_x(\mathbf{x}) = \rho_T(\mathbf{d}) \frac{\partial \mathbf{d}}{\partial \mathbf{x}}$ (where $\frac{\partial \mathbf{d}}{\partial \mathbf{x}}$ is the Jacobian determinant of the transformation $\mathbf{d} = \mathbf{d}(\mathbf{x})$), we can write the

mismodelling covariance as

$$ C_T = \int [\mathbf{d}(\mathbf{x}) - \bar{\mathbf{d}}] [\mathbf{d}(\mathbf{x}) - \bar{\mathbf{d}}]' \rho_x(\mathbf{x}) d\mathbf{x}, \quad (12) $$

where $\bar{\mathbf{d}} = \int \mathbf{d}(\mathbf{x}) \rho_x(\mathbf{x}) d\mathbf{x}$ is the average data vector.

As a first attempt, we can estimate the centroid mislocation covariance $C_T$ by using the following approach. We first calculate a preliminary source model $\bar{\mathbf{m}}$ by performing a WCMT inversion at the centroid location $\bar{\mathbf{x}}$ and generate a random population of locations from $\rho_x(\mathbf{x})C_x$ being given as an input to describe the location uncertainty). For the preliminary model $\bar{\mathbf{m}}$ and each of these locations $\mathbf{x}_k$, we can compute the predicted data $\mathbf{d}(\mathbf{x}_k)$. We then use the second moment of this population as a proxy for the covariance

$$ (C_T)^{ij} = \frac{1}{L} \sum_{k=1}^{L} [d'(\mathbf{x}_k) - \bar{d}'^i][d'(\mathbf{x}_k) - \bar{d}'^i], \quad (13) $$

where $L$ is the total number of locations which have been generated randomly and $\bar{d}'^i$ the empirical mean $\bar{d}'^i = \frac{1}{L} \sum_{k=1}^{L} d'(\mathbf{x}_k)$.

On Fig. 8, this approach is illustrated with a synthetic example where the location error is elongated in the north-south direction. After calculating the synthetic data set using the actual source model shown in green, we filter the waveforms in the 1-5 mHz frequency passband. Since we focus here on the effect of the mismodelling (i.e. considering only $\rho_T(\mathbf{d}|\mathbf{m})$ and neglecting $\rho_D(\mathbf{d})$), we use the synthetic data set without adding background noise to the seismograms. As a result of the centroid position uncertainty, the source is mislocated 0.6° to the south of the actual centroid. The red mechanism is obtained by performing a WCMT inversion which neglects the centroid mislocation information [i.e. ignoring $C_T$ in eq. (4)]. This solution is clearly affected by the mislocation of the centroid. This source model is then used to calculate several data sets corresponding to 1000 random locations drawn from a probability density reflecting the centroid location error (which in this case is larger in the north-south direction). We then use eq. (13) to calculate $C_T$ from this population of data sets. The WCMT solution obtained by including the resulting covariance is shown in orange on Fig. 8. The match between the WCMT solution and the actual source model is indeed clearly improved when the centroid mislocation is taken into account. This approach is time-consuming but it highlights the improvements that can be achieved by incorporating the mismodelling information into the inverse problem. As presented in Section 2, the solution of the inverse problem is the posterior probability density $\sigma_M(\mathbf{m})$ given by eq. (3)-(5). As soon as the covariance $C_D$ does not reflect the true observational errors ($C_D$) and the true modelling uncertainties ($C_T$), the solution $\sigma_M(\mathbf{m})$ is affected and the most probable model $\bar{\mathbf{m}}$ at the maximum of $\sigma_M(\mathbf{m})$ is biased.

Knowing how the moment tensor estimates can be improved by considering the mislocation information, we can now seek a more efficient way to calculate $C_T$. In the vicinity of the centroid location $\bar{\mathbf{x}}$, we can write up to the first order:

$$ \mathbf{d}(\mathbf{x}) = \mathbf{d}(\bar{\mathbf{x}}) + \nabla \mathbf{d}(\bar{\mathbf{x}})(\mathbf{x} - \bar{\mathbf{x}}). \quad (14) $$

Since we are working at very long periods, the centroid mislocation is generally small compared to seismic wavelengths being used and $\mathbf{d}(\mathbf{x})$ is smooth in the neighbourhood of $\bar{\mathbf{x}}$. Thus, we can write $\bar{\mathbf{d}} = \mathbf{d}(\bar{\mathbf{x}})$ and insert eq. (14) in eq. (12) to show

$$ C_T = [\nabla \mathbf{d}(\bar{\mathbf{x}})] C_x [\nabla \mathbf{d}(\bar{\mathbf{x}})]', \quad (15) $$

where $C_x$ is the centroid position covariance in eq. (11) defined as

$$ C_x = \int (\mathbf{x} - \bar{\mathbf{x}})(\mathbf{x} - \bar{\mathbf{x}})' \rho_x(\mathbf{x}) d\mathbf{x}. \quad (16) $$
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Synthetic WCMT inversion incorporating the centroid mislocation information by using the mismodelling covariance $C_T$. The data vector is setup by generating a synthetic data set corresponding to the green focal mechanism (i.e. the actual source model) with a distribution of stations indicated by the blue triangles. The WCMT inversion is performed by assuming a wrong centroid location (red star) which is shifted 0.6° to the south of the actual centroid location (green star). The mismodelling covariance leading to the orange mechanism is computed empirically from eq. (13) using 1000 data sets corresponding to 1000 random centroid locations (black dots). This distribution of location is drawn from a probability density $\rho_X(x)$ of covariance $C_X$ which is used to calculate the mismodelling covariance in eq. (15) associated with the blue mechanism. The WCMT solution resulting from the standard W phase algorithm which neglects the centroid mislocation is shown for comparison in red.</figcaption></figure>

In this section, it is implicitly assumed that the modelling error can be adequately described by a Gaussian probability density (i.e. $\rho_T(d|m)$ in eq. (2) is Gaussian). This hypothesis is not supported in the general case, especially if the data model relationship is non-linear. However, since we can reasonably consider the input location uncertainty $\rho_X(x)$ to be Gaussian, this assumption is valid here as long as eq. (14) stands in the vicinity of the centroid location.

In practice, the data gradient $\nabla d(\bar{x})$ can be estimated by using pre-computed derivatives of $G$ with respect to the source location $\mathbf{x}$

$$[\nabla d(\bar{x})]^{ij} = \frac{\partial G^{ik}}{\partial x_j}(\bar{x}) \bar{m}^k, \qquad (17)$$

where $\bar{m}$ is the moment tensor obtained after performing a preliminary WCMT inversion at the centroid location $\bar{x}$.

The blue mechanism in Fig. 8 shows the results obtained by using eq. (15) for $C_T$ in the previous Californian synthetic application. The resulting WCMT solution is very similar to the actual source model. A good match is also observed using the empirical covariance in eq. (13), but the computation of $C_T$ using eq. (15) is more efficient.

## 6 APPLICATION

To understand the centroid mislocation described by the covariance $C_X$ in the calculation of the mismodelling covariance $C_T$, we compared the location estimates coming from different catalogues. Duputel et al. (2012) recently applied the WCMT algorithm to the set of $M_w \ge 6.5$ earthquakes between 1990 and 2010. Fig. 9(a) shows

the differences between the resulting centroid location and the ones from Harvard/Global CMT (GCMT) catalogue (Dziewonski 1982; Ekström & Nettles 2006). We note that GCMT and WCMT solutions are often separated by several tens of kilometres. As discussed by Duputel et al. (2012), the spatial resolution of the centroid estimated from the W phase algorithm is limited because it involves very long wavelengths. Important location uncertainties can also be observed using shorter periods data. Fig. 9(b) shows the differences between the USGS hypocentre (for brevity we will refer here to the USGS hypocentre as the PDE) and the GCMT centroid location for $5.5 \le M_w \le 6.0$ events. Since these are small earthquakes, the distance between the hypocentre and the centroid location should not exceed 10 km, which is much smaller than what is observed here. This is obviously the consequence of the uncertainty on hypocentre and centroid locations. The PDE location estimates is based on body wave arrivals which are generally affected by unmodelled shallow heterogeneities. On the other hand, as discussed by Hjörleifsdóttir & Ekström (2010), the GCMT centroid location can be biased by the Earth model used to calculate the Green's functions in certain regions of the world (e.g. North and South America).

Given the variability of the centroid location discussed earlier, it is desirable to introduce the location uncertainty in CMT inverse problems. As a first application, we performed WCMT inversions by fixing the centroid to the hypocentre and account for the additional mislocation due to the distance between the rupture initiation and the centroid position. The aim here is to have robust moment tensor estimates regardless of the actual centroid
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Errors in the centroid location. The difference between WCMT and GCMT centroid locations for $M_w \ge 6.5$ earthquakes between 1990 and 2010 (771 events) is shown in (a) and the difference between USGS epicentre (PDE) and GCMT centroids for $5.0 \le M_w \le 6.0$ earthquakes between 2009 and 2010 (705 events) is shown in (b). The arrows point from the WCMT centroid in (a) and from the PDE location in (b) to the GCMT centroid position.</figcaption></figure>

location. We then focus on the 2009 Vanuatu earthquake sequence for which a significant part of the mismodelling can be attributed to the large amplitude disturbances caused by preceding events. We finally discuss the 2011 Tohoku-oki $M_w = 9.0$ earthquake and in particular the uncertainty on the fault dip and the seismic moment associated with the centroid location variability. As is done previously during the synthetic experiments, the data are filtered in the 1–5 mHz passband after proper instrument response correction.

## 6.1 Inversions with the centroid fixed at the preliminary hypocentre

To reduce the perturbations due to the source mislocation, most applications attempt to search for an optimum point source location in addition to moment tensor elements. The centroid position is generally obtained by minimizing a quadratic misfit function which in practice can be performed using various optimization methods such as the steepest descent algorithm (e.g. Dziewonski et al. 1981) or a simple grid search such as in the WCMT algorithm. Performing a spatial grid-search is roughly equivalent to assume a uniform probability density to describe the location uncertainty in the region

where the centroid is explored. As discussed earlier, we know however that the centroid location inverted following such procedures can produce significant errors.

In this section, we follow an alternative strategy: the centroid position is fixed to the PDE location and we introduce a mismodelling covariance $C_T$ to describe the centroid location uncertainty as detailed in Section 5. This approach is of particular interest for warning purposes since our primary goal here is to have better moment tensor estimates regardless of the precise centroid location. As a matter of fact, it has direct implications in grid-based realtime determination of moment tensors (e.g. GRiD MT, Tsuruoka et al. 2009). In this particular context, taking into account the modelling error due to the centroid mislocation should allow a dramatic decrease in the number of point source to be inverted continuously and simultaneously.

In this section, we assume an isotropic horizontal centroid location uncertainty (i.e. $C_x = \sigma_x I$): we use a standard deviation $\sigma_x = 1.0^\circ$ which is appropriate given the distance between the PDE and the centroid location for large earthquakes. The average distance between PDE and WCMT centroid for $M_w \ge 8.0$ earthquakes between 1990 and 2010 is about $0.8^\circ$ (the same value is obtained for the standard deviation between PDE and GCMT locations).

### 6.1.1 2001 Peru earthquake ($M_w = 8.4$)

The GCMT solution for this event is shown in green on Fig. 10. We keep the centroid fixed at the PDE location and perform two WCMT inversions with and without incorporating realistic error analyses. This event is actually quite challenging for this purpose because of a large distance of 150 km between the PDE and the GCMT centroid. This is caused by an unidirectional propagation of the rupture toward the southeast with a significant slip away from the epicentre (Ruegg et al. 2001; Giovanni et al. 2002).

The WCMT solution obtained without any formal error analysis (i.e. neglecting $C_T$ and assuming $C_d = \sigma^2 I$) is shown in red on Fig. 10. Although the magnitude is not much affected, the resulting mechanism significantly differs from GCMT as a consequence of the large distance between the PDE and the centroid location. Moreover, as observed in the synthetic test shown on Fig. 4, the associated posterior uncertainty is clearly under-estimated.

On the contrary, the solution presented in blue on Fig. 10 shows a pretty good match to the GCMT solution and to the WCMT solution proposed by Duputel et al. (2012) with more realistic estimates of the posterior error on the source model parameters. This solution is obtained by taking into account the block diagonal $C_d$ incorporating noise level measurement per-trace and off-diagonal terms as well as the mismodelling covariance $C_T$ associated with the centroid location uncertainty. This shows that even if we assume a point source shifted more than 150 km away from the centroid location, the resulting moment tensor is well resolved with a realistic estimate of its uncertainty.

### 6.1.2 2001 Kokoxili earthquake ($M_w = 7.8$)

The 2001 Kokoxili earthquake is a continental strike-slip earthquake associated with a mostly unidirectional rupture propagation along a long fault of at least 400 km with a significant slip eastward of the hypocentre (Tocheport et al. 2006). This event is interesting in the present context because of the large distance of about 200 km, between the PDE and the centroid location, as shown in Fig. 11. As for the 2001 Peru earthquake, the
<figure><img src="image_11.png" /><figcaption>**Figure 10.** W phase source inversion results for the 2001 Peru earthquake. The GCMT solution for this earthquake is shown in green. The WCMT solutions are shown on the map and detailed on the right with the associated posterior uncertainties. The inversions are performed with a source location fixed at the PDE (red star) which is shifted 150 km to the northeast of the GCMT centroid (green star). The WCMT solution obtained assuming $C_d = \sigma^2 I$ while neglecting $C_T$ is shown in red. The WCMT solution presented in blue is obtained by combining the block diagonal $C_d$ in eq. (9) with the mismodelling covariance $C_T$ of eq. (15). The black dots shown on the map are drawn randomly from the probability density $\rho_X(x)$ describing the centroid location uncertainty associated with the covariance $C_X$ used in the calculation of $C_T$. The distribution of used stations is indicated by the blue triangles on the globe.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** W phase source inversion results for the 2001 Kokoxili earthquake. The GCMT solution for this earthquake is shown in green. The WCMT solutions are shown on the map and detailed on the right with the associated posterior uncertainties. The inversions are performed with a source location fixed at the PDE (red star) which is shifted 214 km to the northeast of the GCMT centroid (green star). The WCMT solution obtained assuming $C_d = \sigma^2 I$ while neglecting $C_T$ is shown in red. The WCMT solution presented in blue is obtained by combining the block diagonal $C_d$ in eq. (9) with the mismodelling covariance $C_T$ of eq. (15). The black dots shown on the map are drawn randomly from the probability density $\rho_X(x)$ describing the centroid location uncertainty associated with the covariance $C_X$ used in the calculation of $C_T$. The distribution of used stations is indicated by the blue triangles on the globe.</figcaption></figure>

WCMT inversions are performed by fixing the centroid to the PDE location with and without incorporating proper data covariance $C_D$.

The solution in red on Fig. 11 is obtained by neglecting $C_T$ and assuming $C_d = \sigma^2 I$. Because of the large difference between the

PDE location used in the inversion and the actual centroid position, the resulting magnitude and mechanism are very different from the GCMT solution with a clear underestimation of the associated posterior uncertainty. The solution obtained automatically by Duputel *et al.* (2012) (yellow) is also different from the GCMT solution as
<figure><img src="image_13.png" /><figcaption>**Figure 12.** W phase source inversion results for the 2009 Vanuatu earthquake sequence. The GCMT solution for (1) the *M*<sub>w</sub> = 7.6 earthquake at 22:03 UTC is indicated in (a), (2) the WCMT solution obtained for the *M*<sub>w</sub> = 7.8 event at 22:18 UTC is shown in (b) and (3) the WCMT solution obtained for the *M*<sub>w</sub> = 7.4 earthquake at 23:13 UTC is presented in (c). In (d) are shown some examples of observed waveforms (black lines) and the corresponding synthetics (red lines) computed from the GCMT solution in (a) and the WCMT solutions shown in (b), and (c). The station azimuth (φ) and epicentral distance (Δ) are indicated relatively to the PDE location of (c). The W phase time windows are bounded by blue circles for the event (2) in (b) and yellow circles for the earthquake (3) in (c).</figcaption></figure>
a consequence of the poor optimum location obtained after grid-search.

If we include a more formal error analysis during the inversion by combining the block diagonal $C_d$ in eq. (9) with the mismodelling covariance $C_T$ of eq. (15), the resulting WCMT solution shown in

blue on Fig. 11 shows a very good match to the GCMT solution. We note also that the posterior uncertainty obtained for this event is larger than the one obtained for the 2001 Peru earthquake. This is consistent with the smaller signal-to-noise ratio for the 2001 Kokoxili event because of its smaller magnitude ($M_w = 7.8$).

<figure><img src="image_14.png" /><figcaption>**Figure 13.** W phase source inversion results for $M_w \ge 7.0$ disturbed events of the 2009 Vanuatu earthquake sequence. The solution obtained for the second event (2) of this sequence at 22:18 UTC is shown in (a), and the one obtained for the third event (3) occurring at 23:13 UTC is presented in (b). The GCMT solution for these earthquakes are shown in green. The WCMT solutions are shown on the map and detailed on the right with the associated posterior uncertainties. The inversions are performed with source locations fixed at the PDE (red star) which are shifted respectively 83 km in (a) and 14 km in (b) to the northeast of the GCMT centroid (green star). The WCMT solutions obtained assuming $C_d = \sigma^2 I$ while neglecting $C_T$ are shown in red. The WCMT solutions presented in blue and yellow are obtained by including a more formal error analysis during the inversion. It includes the covariance $C_T$ of eq. (15) and a block diagonal covariance $C_d$ as described in the main text. The black dots shown on the maps are drawn randomly from the probability density $\rho_X(x)$ describing the centroid location uncertainties associated with the covariance $C_x$ used in the calculation of $C_T$. The distributions of used stations are indicated for both event by the blue triangles on the globes.</figcaption></figure>
## 6.2 Some unusual earthquakes

### 6.2.1 2009 Vanuatu earthquake sequence

The 2009 October 7 Vanuatu earthquake sequence consisted of three $M_w \ge 7.0$ events occurring in less than 5 hr. Due to the short interval between the three events, the W phase waveforms of the second and third earthquakes are contaminated by the large amplitude disturbances caused by earlier events and the standard WCMT inversion does not produce reliable results.

To cope with this situation, we can remove the perturbation caused by earlier events by subtracting the corresponding synthetics from the data in the W phase time window. As discussed by Duputel et al. (2012), this approach is quite straightforward because we run the WCMT algorithm for the disturbed event as an usual earthquake but using the residual traces instead of the original data vector. However, the uncertainty in the resulting solution may be increased since the data misfit associated with the disturbing event source model is added to the background noise level.

We follow here the procedure described earlier by including a formal error analysis during the inversion. As done in Section 6.1, the centroid location is fixed to the PDE location and we use a mismodelling covariance $C_T$ associated with the source mislocation as written in eq. (15). We also include a block diagonal covariance with a form similar to $C_d$ in eq. (9) but with different diagonal elements. We first compute a solution for the disturbed event considering $C_d = \sigma^2 I$ and the residual traces as the data vector. We then use the misfit at each station $n$ as a proxy to determine the diagonal elements $\sigma_d^n$. The new form of $C_D$ used here therefore incorporates (i) the mismodelling due to the centroid mislocation (by using $C_T$), (ii) the background noise level and uncorrected perturbations of earlier events (by including the diagonal elements $\sigma_d^n$) and (iii) the data oversampling (by considering the decaying exponential in $C_d$).

This approach is illustrated in Fig. 12 for the 2009 October 7 Vanuatu earthquakes. This sequence began with (1) a $M_w = 7.6$ earthquake occurring at 22:03 UTC followed 15 min later by (2) a $M_w = 7.8$ earthquake and (3) a $M_w = 7.4$ event at 23:13 UTC. The solution for the $M_w = 7.8$ earthquake (2) is shown in Fig. 12(b) assuming the GCMT model of the $M_w = 7.6$ event (1) shown in Fig. 12(a) to compute the synthetics which are subtracted from the data. We then assess the moment tensor of the $M_w = 7.4$ earthquake (3) shown in Fig. 12(c) after removing the perturbations of the former (1) and (2) events.

The WCMT solutions shown in Figs 12(b) and (c) incorporates the full covariance $C_D$. As shown in Fig. 12(d), there is a very good match between data and waveforms predicted from these point sources. In Fig. 13, they are compared with GCMT and the WCMT solutions obtained by neglecting $C_T$ and assuming $C_d = \sigma^2 I$. On top of having more robust estimates of the posterior uncertainty, taking a more formal error analysis during the inversion allows us to improve the solutions themselves, in particular, by taking into account the additional data error due to the misfit associated with former disturbing events.

### 6.2.2 2011 Tohoku-oki earthquake ($M_w = 9.0$)

The Tohoku-oki earthquake ($M_w = 9.0$) is the largest event during the three years since the W phase algorithm was developed. Thus, it is of particular interest to evaluate the benefit of incorporating realistic error analysis in WCMT inversions. As discussed in Duputel et al. (2011), although the different techniques (e.g. WCMT, GCMT, USGS CMT) provided quite similar solutions in real-time or in post-mortem mode, $M_w$ estimates range from 8.8 to 9.1, and fault dip values are distributed between 10° and 20°, depending on the assumed initial depth. We explore here the possibility

<figure><img src="image_15.png" /><figcaption>**Figure 14.** W phase source inversion results for the 2011 Tohoku-oki earthquake. The GCMT solution for this earthquake is shown in green and the WCMT solution obtained in Duputel et al. (2011) after spacial grid-search using an optimum low-noise data set is presented in yellow. The WCMT solutions obtained in this study are shown in red and blue on the map and detailed on the right with the associated posterior uncertainties. The inversions are performed with a source location fixed at the GCMT centroid (green star). The WCMT solution obtained assuming $C_d = \sigma^2 I$ while neglecting $C_T$ is shown in red. The WCMT solution presented in blue is obtained by combining the block diagonal $C_d$ in eq. (9) with the mismodelling covariance $C_T$ of eq. (15). The black dots shown on the map are drawn randomly from the probability density $\rho_x(x)$ describing the centroid location uncertainty associated with the covariance $C_x$ used in the calculation of $C_T$. The distribution of used stations is indicated by the blue triangles on the globe.</figcaption></figure>
of narrowing and assessing this uncertainty by including the full covariance matrix $C_D$ during the inversion.

The solution obtained in Duputel et al. (2011) for an optimum low-noise data set is shown in yellow on Fig. 14 and the GCMT solution is shown for comparison in green. Although the two mechanisms are similar we note that the magnitudes differ and that the GCMT centroid is shifted to the south of the WCMT optimum location. The solution shown in red on Fig. 14 is obtained if we perform a WCMT inversion with a source location fixed at the GCMT centroid. Although it does not diverge from the GCMT solution and WCMT solution of Duputel et al. (2011), we note that the mechanisms have a different orientation and that the magnitude obtained here is slightly smaller than the GCMT one. We interpret this variability in point source solutions as the consequence of the large fault dimension and of the associated uncertainty in the centroid location. Moreover, as shown in Fig. 14, this uncertainty is clearly underestimated if we assume $C_d = \sigma^2 I$ and neglect $C_T$.

To improve these estimates, we introduce a covariance $C_x$ reflecting the centroid location uncertainty along the fault with standard deviations of 140 km along strike, 70 km along dip and 5 km in

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Posterior uncertainty on the W phase source inversion results for the 2011 Tohoku-oki earthquake. 3000 source models are randomly generated from the posterior Gaussian. The moment magnitude distribution of this population is shown at the top of the figure and the shallow dip histogram is shown below.</figcaption></figure>

depth (cf. black dots in Fig. 14). The associated mismodelling covariance $C_T$ is then calculated using eq. (15). Fig. 14 shows the results obtained if we still assume the GCMT centroid location but incorporate $C_T$ and $C_d$ of eq. (9). The resulting mechanism is similar to GCMT and to the WCMT solution from Duputel et al. (2011). In terms of magnitude, the value $M_w = 9.05$ obtained in this study is between $M_w = 9.02$ from Duputel et al. (2011) and the GCMT estimate $M_w = 9.08$. The histograms on Fig. 15 gives the uncertainty we have on $M_w$ and on the shallow fault dip. Reasonable estimates of $M_w$ range from 9.0 to 9.1 with dip variations between 10° and 12°.

## 7 CONCLUSION

We advocate performing a more formal error analysis in seismic source inversion problems. In this work, we assume Gaussian initial probability densities such that the data uncertainty can be carried by the covariance matrix $C_D$ which should be included into the error propagation analysis. We identify three ingredients to be incorporated in $C_D$ to have an accurate description of input uncertainties. First, there is the noise level which we can measure for each data trace. It forms the non-uniform diagonal of the observational error covariance $C_d$. Secondly, we must include the non-diagonal terms in $C_d$ due to the interdependence of observational errors. Finally we must account for the modelling error which can significantly contribute to errors in the solution.

Although the approach presented here can be applied to much more general inverse problems, we focus on CMT inversion studies based on long period seismological observations. In this perspective, we incorporate a more formal error analysis into the W-phase source inversion algorithm. The forward problem being linear at a given centroid location, the posterior uncertainty on CMT parameters can be modelled by Gaussian distributions (since we assume Gaussian prior probability densities). Given the long period character of the W phase, the observational error is mostly related to the background seismic noise which steadily grows with period. Moreover, we cannot neglect the interdependence of observational errors since the data traces are heavily oversampled. To account for the resulting strong correlation between neighbour data samples, we incorporate non-diagonal terms in $C_d$ by considering a decaying exponential characterized by a correlation duration chosen as the shortest period content after filtering (e.g. ≈200 s). Concerning the modelling error, we consider the case of the mismodelling associated with the mislocation of the centroid position. We also considered the amplitude disturbances caused by preceding events during the 2009 Vanuatu earthquake sequence.

The results obtained by performing W phase inversions on synthetic and actual data sets show the importance of incorporating realistic covariance components during CMT inversions. First, it improves the error estimates on source model parameters. Second, it improves the solution itself notably because using an accurate covariance $C_D$ prevents over-fitting of the data traces. Other sources of mismodelling which could be taken into account are the source finiteness, the source complexity for very large earthquakes or the possibility of having an incorrect Earth model.

## ACKNOWLEDGMENTS

This work uses Federation of Digital Seismic Networks (FDSN) seismic data and CMT solutions from the Global CMT catalog. The Incorporated Research Institutions for Seismology (IRIS) Data Management System (DMS) was used to access the data. This work
made use of the Matplotlib python library. We thank J.-J. Lévêque, M. Sambridge and an anonymous reviewer for their helpful comments on the manuscript.

## REFERENCES

Beavan, J., Wang, X., Holden, C., Wilson, K., Power, W., Prasetya, G., Bevis, M. & Kautoke, R., 2010. Near-simultaneous great earthquakes at Tongan megathrust and outer rise in September 2009, *Nature*, **466**, 959–963.

Duputel, Z., Rivera, L., Kanamori, H., Hayes, G.P., Hirshorn, B. & Weinstein, S., 2011. Real-time W Phase inversion during the 2011 Off the Pacific Coast of Tohoku Earthquake, *Earth Planets Space*, **63**, 535–539.

Duputel, Z., Rivera, L., Kanamori, H. & Hayes, G.P., 2012. W phase source inversion for moderate to large earthquakes, *Geophys. J. Int.*, **189**(2), 1125–1147.

Dziewonski, A., 1982. Harvard Centroid Moment Tensor Project, http://www.globalcmt.org (last accessed 2011 December 1).

Dziewonski, A., Chou, T.A. & Woodhouse, J.H., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, *J. geophys. Res.*, **86**, 2825–2852.

Ekström, G. & Nettles, M., 2006. Global Centroid Moment Tensor Project, http://www.globalcmt.org (last accessed 2011 December 1).

Ekström, G., Dziewonski, A.M., Maternovskaya, N.N. & Nettles, M., 2005. Global seismicity of 2003: centroid–moment-tensor solutions for 1087 earthquakes, *Phys. Earth planet. Inter.*, **148**(2-4), 327–351.

Fukahata, Y. & Wright, T.J., 2008. A non-linear geodetic data inversion using ABIC for slip distribution on a fault with an unknown dip angle, *Geophys. J. Int.*, **173**(2), 353–364.

Giovanni, M.K., Beck, S.L. & Wagner, L., 2002. The June 23, 2001 Peru earthquake and the southern Peru subduction zone, *Geophys. Res. Lett.*, **29**, 2018, doi:10.1029/2002GL015774.

Gouveia, W.P. & Scales, J.A., 1998. Bayesian seismic waveform inversion: parameter estimation and uncertainty analysis, *J. geophys. Res.*, **103**(B2), 2759–2779.

Hayes, G.P., Rivera, L. & Kanamori, H., 2009. Source inversion of the W-Phase: real-time implementation and extension to low magnitudes, *Seismol. Res. Lett.*, **80**, 817–822.

Hjörleifsdóttir, V. & Ekström, G., 2010. Three-dimensional Earth structure on CMT earthquake parameters, *Phys. Earth planet. Inter.*, **179**, 178–190.

Ide, S., Takeo, M. & Yoshida, Y., 1996. Source process of the 1995 Kobe earthquake: determination of spatio-temporal slip distribution by Bayesian modeling, *Bull. seism. Soc. Am.*, **86**(3), 547–566.

Jackson, D.D., 1979. The use of a priori data to resolve non-uniqueness in linear inversion, *Geophys. J. R. astr. Soc.*, **57**(1), 137–157.

Kanamori, H. & Rivera, L., 2008. Source inversion of W phase: speeding up seismic tsunami warning, *Geophys. J. Int.*, **175**(1), 222–238.

King, G., 2007. Fault interaction, earthquake stress changes, and the evolution of seismicity, in *Treatise on Geophysics*, pp. 225–255, Elsevier, Amsterdam.

Lay, T., Ammon, C.J., Kanamori, H., Rivera, L., Koper, K. & Hutko, A.R., 2010. The 2009 Samoa-Tonga great earthquake triggered doublet, *Nature*, **466**(7309), 964–968.

Ruegg, J.-C., Olcay, M. & Lazo, D., 2001. Co-, post- and pre(?)-seismic displacements associated with the Mw 8.4 southern Peru earthquake of 23 June 2001 from continuous GPS measurements, *Seismol. Res. Lett.*, **72**(6), 673–678.

Sambridge, M., 1999. Geophysical inversion with a neighbourhood algorithm II. Appraising the ensemble, *Geophys. J. Int.*, **138**(3), 727–746.

Satake, K., 2007. Tsunamis, in *Treatise on Geophysics*, pp. 483–511, Elsevier, Amsterdam.

Sorrells, G.G., 1971. A preliminary investigation into the relationship between long-period seismic noise and local fluctuations in the atmospheric pressure field, *Geophys. J. R. astr. Soc.*, **26**(1-4), 71–82.

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**, 389–404.

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

Tarantola, A. & Valette, B., 1982. Inverse problems = quest for information, *J. Geophys.*, **50**, 159–170.

Tocheport, A., Rivera, L. & Van der Woerd, J., 2006. A Study of the 14 November 2001 Kokoxili Earthquake: history and geometry of the rupture from teleseismic data and field observations, *Bull. seism. Soc. Am.*, **96**(5), 1729–1741.

Tsuruoka, H., Kawakatsu, H. & Urabe, T., 2009. GRiD MT (grid-based real-time determination of moment tensors) monitoring the long-period seismic wavefield, *Phys. Earth planet. Inter.*, **175**(1-2), 8–16.

Wald, D.J., Worden, B.C., Quitoriano, V. & Pankow, K.L., 2005. *ShakeMap Manual: Technical Manual, User's Guide, and Software Guide*, available at: http://pubs. usgs. gov/tm/2005/12A01 (last accessed 2011 December 1).

Yabuki, T. & Matsu'ura, M., 1992. Geodetic data inversion using a Bayesian information criterion for spatial distribution of fault slip, *Geophys. J. Int.*, **109**, 363–375.

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.

Yagi, Y. & Fukahata, Y., 2011. Introduction of uncertainty of Green's function into waveform inversion for seismic source processes, *Geophys. J. Int.*, **186**(2), 711–720.