# W-phase fast source inversion for moderate to large
earthquakes (1990 - 2010)

Zacharie Duputel, Luis Rivera, Hiroo Kanamori, Gavin Hayes

► To cite this version:

Zacharie Duputel, Luis Rivera, Hiroo Kanamori, Gavin Hayes. W-phase fast source inversion for moderate to large earhquakes (1990 - 2010). Geophysical Journal International, 2012, 189, pp.1125-1147. hal-00839951

HAL Id: hal-00839951

https://hal.science/hal-00839951

Submitted on 4 Jan 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# W phase source inversion for moderate to large earthquakes (1990–2010)

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

Accepted 2012 February 10. Received 2012 February 10; in original form 2011 February 6

## SUMMARY

Rapid characterization of the earthquake source and of its effects is a growing field of interest. Until recently, it still took several hours to determine the first-order attributes of a great earthquake (e.g. $M_w \ge 7.5$), even in a well-instrumented region. The main limiting factors were data saturation, the interference of different phases and the time duration and spatial extent of the source rupture. To accelerate centroid moment tensor (CMT) determinations, we have developed a source inversion algorithm based on modelling of the W phase, a very long period phase (100–1000 s) arriving at the same time as the P wave. The purpose of this work is to finely tune and validate the algorithm for large-to-moderate-sized earthquakes using three components of W phase ground motion at teleseismic distances. To that end, the point source parameters of all $M_w \ge 6.5$ earthquakes that occurred between 1990 and 2010 (815 events) are determined using Federation of Digital Seismograph Networks, Global Seismographic Network broad-band stations and STS1 global virtual networks of the Incorporated Research Institutions for Seismology Data Management Center. For each event, a preliminary magnitude obtained from W phase amplitudes is used to estimate the initial moment rate function half duration and to define the corner frequencies of the passband filter that will be applied to the waveforms. Starting from these initial parameters, the seismic moment tensor is calculated using a preliminary location as a first approximation of the centroid. A full CMT inversion is then conducted for centroid timing and location determination. Comparisons with Harvard and Global CMT solutions highlight the robustness of W phase CMT solutions at teleseismic distances. The differences in $M_w$ rarely exceed 0.2 and the source mechanisms are very similar to one another. Difficulties arise when a target earthquake is shortly (e.g. within 10 hr) preceded by another large earthquake, which disturbs the waveforms of the target event. To deal with such difficult situations, we remove the perturbation caused by earlier disturbing events by subtracting the corresponding synthetics from the data. The CMT parameters for the disturbed event can then be retrieved using the residual seismograms. We also explore the feasibility of obtaining source parameters of smaller earthquakes in the range $6.0 \le M_w < 6.5$. Results suggest that the W phase inversion can be implemented reliably for the majority of earthquakes of $M_w = 6$ or larger.

**Key words:** Tsunamis; Earthquake source observations; Surface waves and free oscillations; Wave propagation; Early warning.

## 1 INTRODUCTION

Considerable effort has been made over the last two decades regarding the design and implementation of tools aimed at fast characterization of earthquake sources. Magnitudes, moment tensors, rupture patterns, shake maps, tsunami excitation and propagation scenarios are now routinely calculated and disseminated by several agencies whenever a significant earthquake occurs. The interest in making

these estimations quickly available is twofold. Authorities and relief agencies can use them to plan and perform rescue and aid operations. Earth scientists, on the other hand, rely upon this information to make critical decisions on re-programming a satellite orbit, or designing a field experiment, etc. The delay in availability of such results is highly variable and depends on the size of the event, its location on the globe and on the type of required result itself. For example, for an $M_w = 5$ earthquake occurring today in Japan or in

<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">³U.S. Geological Survey, National Earthquake Information Center (contracted by Synergetics, Inc.), Golden, CO, USA</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** W phase source inversion results for the 2004 great Sumatra–Andaman Islands earthquake. W phase CMT (WCMT) and Global CMT (GCMT) solutions are shown in the topleft corner. Examples of observed waveforms (black lines) and the corresponding synthetics (red lines) computed from the WCMT solution are presented. The station azimuth ($\phi$) and epicentral distance ($\Delta$) are indicated as well as W phase time window, which are bounded by red dots. The WCMT inversion is based on the ground motion of stations within $\Delta \le 90^\circ$ after applying a bandpass filter in the 1–5 mHz passband. W phase and later arrivals are often very well predicted by the WCMT solution although many channels have instrument problems during or after the surface waves arrivals (most of the broad-band data within $\Delta \le 40^\circ$ are saturated).</figcaption></figure>

California, a fairly reliable magnitude estimation is available within a few minutes after the origin time. As an opposite end-member, it usually takes at least several hours before a reliable slip distribution is determined for a large earthquake ($M_w > 7.0$) even in a well-instrumented region (e.g. Hayes et al. 2011).

Although it is desirable to reduce such a delay, there are severe limitations in this practice. The rupture process of large events can last several minutes, and any analyses based on the beginning of the records can only provide a lower bound to the magnitude estimate or to the extent of the associated rupture. Although it is always possible to use some hypothetical, simplified model to predict the final magnitude from signals generated during the early stages of the rupture, there is increasing evidence for a wide diversity in the nature of seismic sources (Kanamori 2004), which translates into significant variations in narrow-band magnitude estimations. In addition to the rupture duration, it is also necessary to account for wave propagation time, which can be substantial (e.g. 9 min for a direct *P* wave from a shallow event to be recorded at $\Delta = 50^\circ$). The use of nearby or regional records would then appear as desirable. Unfortunately, we do not always have good regional coverage with the existing global networks. Even in cases where regional data are available, the use of seismological records close to the source (say $\Delta < 5^\circ$) can be problematic because such data are more sensitive to variations

in the shallow structure of the Earth. Their use therefore requires an accurate earth model which incorporates 3-D regional heterogeneity (Tsuruoka et al. 2009). Furthermore, if the rupture length is large, a point source approximation can be problematic for the closest stations even when analysing long periods. Finally, regional and even teleseismic records of large events are often saturated across the frequency range. For example, most of the broad-band data available at the Federation of Digital Seismograph Network (FDSN) stations within $40^\circ$ of the 2004 Sumatra–Andaman Islands earthquake rupture are saturated at the arrival time of the surface waves. (Fig. 1).

To overcome these limitations, we have developed a centroid moment tensor (CMT) inversion algorithm using the W phase, a very long period phase (100–1000 s) identified after the 1992 Nicaragua tsunami earthquake by Kanamori (1993). Use of the W phase for fast and accurate quantification of the source properties of great earthquakes is advantageous for several reasons. First, the W phase has a fast group velocity, which facilitates rapid inversion after an earthquake occurs. Second, the major part of the propagating energy associated with the W phase is confined to the mantle, which is less heterogeneous than the crust, thus resulting in a relatively simple waveform. Finally, the W phase has a very long period character, which is essential for the source inversion of large earthquakes.
<table><caption>**Table 1.** Number of events as a function of date and magnitude.</caption><thead><tr><th></th><th>6.50–6.74</th><th>6.75–6.99</th><th>7.00–7.24</th><th>7.25–7.49</th><th>7.50–7.99</th><th>8.00–9.00</th><th>Total</th></tr></thead><tbody><tr><td>1990–1993</td><td>65</td><td>39</td><td>25</td><td>7</td><td>11</td><td>0</td><td>147</td></tr><tr><td>1994–1997</td><td>72</td><td>40</td><td>32</td><td>5</td><td>17</td><td>3</td><td>169</td></tr><tr><td>1998–2001</td><td>53</td><td>35</td><td>23</td><td>10</td><td>13</td><td>3</td><td>137</td></tr><tr><td>2002–2005</td><td>72</td><td>26</td><td>19</td><td>14</td><td>11</td><td>4</td><td>146</td></tr><tr><td>2006–2009</td><td>71</td><td>40</td><td>15</td><td>17</td><td>17</td><td>5</td><td>165</td></tr><tr><td>2010–2010</td><td>18</td><td>14</td><td>6</td><td>9</td><td>3</td><td>1</td><td>51</td></tr><tr><td>Total</td><td>351</td><td>194</td><td>120</td><td>62</td><td>72</td><td>16</td><td>815</td></tr></tbody></table>

In general, long periods increasingly dominate the wavefield as moment grows because the duration of rupture gets longer as the event gets larger. Furthermore, the size of the ruptured fault and the amount of final slip control the tsunami wave height. The long-period wavefield is also useful for the source analysis of tsunami earthquakes, defined by Kanamori (1972). These events are often characterized by an anomalously long-period spectrum at the source, which produces unusually large tsunami height relative to short-period magnitude estimations (≤100 s, e.g. $M_s$). In the specific case of CMT determinations, it is also fundamental to consider very long period waves since the point source approximation is used even for earthquakes rupturing large faults (≥100 km).

The W phase CMT algorithm was initially developed by Kanamori & Rivera (2008) using the inversion of vertical components of ground motion. For brevity, we will refer to this paper as KR. A real-time application at the National Earthquake Information Center (NEIC) using this version of the algorithm was set up by Hayes et al. (2009) for $M_w \ge 5.8$ earthquakes on a global scale. The algorithm has now been extended to all three components of the ground motion and has been deployed in real time at the Institut de Physique du Globe de Strasbourg (IPGS) for testing purposes, at NEIC and at the Pacific Tsunami Warning Center (PTWC) of the National Oceanic and Atmospheric Administration (Duputel et al. 2011). We have also adapted the inversion for use on the regional scale, an application which is currently being tested in California, Japan, Mexico and Taiwan (Rivera & Kanamori 2009; Rivera et al. 2010).

The purpose of this work is to fine-tune and validate the W phase algorithm for large to moderate earthquakes using three-component ground motions at regional and teleseismic distances. To this end, the point source parameters of all $M_w \ge 6.5$ earthquakes that occurred between 1990 and 2010 (815 events) are systematically determined using FDSN, Global Seismographic Network broad-band stations (GSN_BROADBAND) and STS1 global virtual networks. Although our new W phase catalogue is complete to $M_w = 6.5$, we also explore the use of the W phase inversion for all earthquakes between $6.0 \le M_w < 6.5$, to assess whether reliable source parameters can also be obtained for smaller events.

In subsequent sections, we focus on the application of the W phase inversion to produce a complete catalogue of events for all $M_w \ge 6.5$ earthquakes since 1990. Although the approach used for smaller events is the same, we discuss this specific application in a separate section, as the resulting W phase catalogue is no longer complete and as such should be distinguished from our main focus.

## 2 DATA AND PRELIMINARY TREATMENTS

We use three-component broad-band data of moderate to large earthquakes that occurred in the period 1990–2010 and were recorded at regional and teleseismic distances. To have a homogeneous reference

catalogue, we start from the Global CMT (GCMT; Dziewonski 1982; Ekström et al. 2005; Ekström & Nettles 2006) database and select all the events with $M_w \ge 6.5$. We use the moment tensor elements provided by GCMT to compute scalar moment $M_0$ using Silver & Jordan (1982)'s and Dahlen & Tromp (1998)'s definition of $M_0 = \sqrt{\sum_{ij} M_{ij} M_{ij}/2}$ and $M_w$ as

$$M_w = \frac{2}{3}(\log_{10}(M_0) - 16.10) \quad (1)$$

with $M_0$ in dyne-cm (Kanamori 1977; Hanks & Kanamori 1979). Table 1 lists the magnitude distribution of these events as a function of time.

Our main criterion for data selection is the availability of a broad-band sensor. Four types of sensor dominate our data set: STS-1, STS-2, KS-5400 and CMG-3T. Data are obtained through NetDC from the data holdings at the Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC) and Geoscope. We use 1 s-sampled data (LHZ) mostly from II, IU, G, GE and MN networks. Some additional stations from other FDSN-affiliated networks are also included to improve spatial coverage. For some stations, several streams are available (different location-IDs). In such cases, we give priority to the longer period sensor. Coverage is quite variable, and depends not only on the event size and its epicentral location, but also on time as a result of the improvement of the worldwide broad-band station distribution since 1990 (Fig. 2).

To design the W phase source inversion, we have to choose the proper time window and passband corner frequencies. Following KR, at a given epicentral distance $\Delta$, the W phase time window extends from the theoretical $P$ arrival time $t_P(\Delta)$ to $t_P(\Delta) + 15 \text{ s deg}^{-1} \times \Delta$. The end time is chosen here to ensure that the time window ends before the arrival of large surface wave trains. The original standard frequency band used in KR for large events ($M_w \ge 8.0$) was 1–5 mHz. However, to have a sufficiently high signal-to-noise ratio for smaller events, it is necessary to gradually shift the passband towards higher frequencies (Hayes et al. 2009, see Table 2). This is related to the well-known behaviour of the background noise steadily growing at periods longer than 200 s. The long-period edge of the bandpass filter is also dictated by the fact that some of the seismometers used become noisy at very long periods. To choose the appropriate bandpass corner frequencies (Table 2) prior to the inversion of each event, we perform a preliminary magnitude estimation by measuring the overall vertical W phase amplitude in the 1–5 mHz bandpass as detailed in Section 3.3.

Following KR, both instrumental deconvolution (to ground displacement) and bandpass filtering (fourth order, causal, Butterworth) are implemented in the time domain as infinite impulse response (IIR) filters. Working in the time domain is very useful for real-time operations since the data can be processed sample by sample as they become available. Moreover, in contrast to traditional frequency-domain deconvolution, it allows the retrieval of the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Data coverage for LH channels extracted from SEED volumes of $M_w \ge 6.5$ events for virtual networks: FDSN, GSN_BROADBAND and STS1 between 1990 and 2010. The number of LH channels is shown in (a) whereas the azimuthal coverage is presented in (b). The total number of available LH channels is shown in black. For stations within epicentral distance $\Delta \le 90^{\circ}$, the selected W phase traces before and after the data screening are presented, respectively, in red and blue.</figcaption></figure>

<table><caption>**Table 2.** Corner frequencies used for butterworth bandpass filtering (fourth order, causal) in W phase inversion when using three components. The frequency passbands used by Hayes *et al.* (2009) were defined for W phase inversion using only vertical components.</caption><thead><tr><th>Magnitude range</th><th>Passband filter, mHz (s)</th></tr></thead><tbody><tr><td>$M_{w-wprel} \ge 8.0$</td><td>1.0–5.0 (200–1000)</td></tr><tr><td>$8.0 > M_{w-wprel} \ge 7.5$</td><td>1.7–6.7 (150–500)</td></tr><tr><td>$7.5 > M_{w-wprel} \ge 7.0$</td><td>2.0–8.3 (120–500)</td></tr><tr><td>$7.0 > M_{w-wprel} \ge 6.5$</td><td>4.0–10.0 (100–250)</td></tr></tbody></table>

W phase on records clipped at the arrival of large-amplitude surface waves. Once cut to length, traces are concatenated to build the data set to be used in the inversion.

Finally, two data sets are defined for each event according to the maximal epicentral distance: $5^{\circ} < \Delta < 50^{\circ}$ and $5^{\circ} < \Delta < 90^{\circ}$. The first data set is available 22 min after the origin time; the second requires an additional 13 min. The reason for using these two distinct data sets will be made clear in the next section.

# 3 METHODOLOGY

## 3.1 Overall real-time operation protocol

The W phase centroid moment tensor (WCMT) inversion is in a way similar to the approach of Dziewonski *et al.* (1981) and Dziewonski & Woodhouse (1983). The three main differences are (1) the time window, (2) the longer periods used in the WCMT inversion and (3) the algorithm employed to determine the best point source location (centroid). The source parameters to be determined are the elements of the seismic moment tensor $\mathbf{f} = [M_{rr}, M_{\theta\theta}, M_{\phi\phi}, M_{r\theta}, M_{r\phi}, M_{\theta\phi}]'$ and the four space-time coordinates of the centroid $\boldsymbol{\eta}_c = [\theta_c, \phi_c, r_c, \tau_c]'$ where $\theta_c$ is the colatitude, $\phi_c$ is the longitude, $r_c$ is the radius and $\tau_c$ is the centroid time. The full WCMT solution

vector can thus be defined as

$$ \mathbf{m} = \begin{pmatrix} \mathbf{f} \\ \boldsymbol{\eta}_c \end{pmatrix}. \qquad (2) $$

Here, we use the term ‘centroid’ following the common practice in source inversion studies, but what we actually determine is the best point source location and the mechanism. Thus, the centroid location and CMT here should be interpreted as the best point source location and the moment tensor, respectively. The centroid $\boldsymbol{\eta}_c$ can be estimated by seeking the point source location that minimizes a quadratic misfit function between the W phase data vector ($\mathbf{d}_w$) and the corresponding synthetic vector ($\mathbf{s}_w$)

$$ \chi(\mathbf{m}) = \frac{1}{2}(\mathbf{s}_w(\mathbf{m}) - \mathbf{d}_w) \cdot (\mathbf{s}_w(\mathbf{m}) - \mathbf{d}_w). \qquad (3) $$

The synthetic seismograms $\mathbf{s}_w$ are obtained from pre-computed Green’s functions calculated using normal mode summation for an epicentral distance range of $0^{\circ} \le \Delta \le 90^{\circ}$ with an interval of $0.1^{\circ}$ and for a depth range of 0–760 km. These Green’s functions are computed using the Preliminary Reference Earth Model (PREM) from Dziewonski & Anderson (1981). The effect of finite-source duration on W phase traces is accounted for by assuming the moment rate function (MRF) to be an isosceles triangle of half duration $h_c$ centred at $\tau_c$. There are two main reasons why the location and origin time estimated from body wave arrivals cannot be assumed as the centroid. First, the errors in hypocentral parameters can be substantial when they are obtained within minutes of the origin time. Secondly, for large earthquakes the hypocentre can be significantly different from the centroid.

Fig. 3 presents the overall algorithm we follow in this study. The horizontal axis represents increasing time and does not include effects of data latency. Let us then suppose that somewhere on the globe an event occurs at $t = t_0$ and at $t_a \sim t_0 + 10$ min or so, we receive from some agency the preliminary epicentral coordinates, depth and origin time. For brevity, we will call this information the preliminary determination of epicentre (PDE). In the context of this
<figure><img src="image_4.png" /><figcaption>**Figure 3.** W phase CMT inversion timeline after an earthquake occuring at $t_0 = 0$. The output level 1 (OL1) is presented in dark grey, the output level 2 (OL2) in blue and the output level 3 (OL3) in red. $N$ is the number of channels after data screening, $\gamma$ is the azimuthal gap and $\Delta$ is the epicentral distance. The horizontal axis represents increasing time together with the data flow (red line) and the GF calculation for the centroid optimization (blue lines). The solutions obtained for the $M_w = 8.3$ 2006 Kuril earthquake corresponding to OL2 and OL3 are presented for $\Delta \le 50^\circ$ and $\Delta \le 90^\circ$. The centroid grid-searches are also displayed on the top of the figure for $\Delta \le 50^\circ$ and for $\Delta \le 90^\circ$. The black cross represents the PDE location whereas the red dot corresponds to the WCMT centroid and the green triangle indicates the GCMT centroid. The isolines correspond to the normalized rms misfit computed during the grid-search. The locations of $M_w \ge 5.0$ aftershocks recorded during the 24 hr after the main shock are indicated by white dots (from the USGS PDE catalogue).</figcaption></figure>

work, the PDE means the preliminary hypocentre location and origin time specified in Global CMT solution files (e.g. the first line in the CMTSOLUTION format). Most of these correspond to the best available USGS solution at the time of GCMT computation. The arrival of this information triggers the operation of our application, which will produce three levels of outputs. Output level 1 (OL1) consists of a preliminary W phase magnitude estimation. Output level 2 (OL2) provides a first solution, $\mathbf{m}_{\text{PDE}}$, by using a preliminary fixed location $\eta_{\text{PDE}}$ as the centroid (i.e. the PDE). A grid-search of the centroid location $\eta_c$ is then performed to obtain Output level 3 (OL3), which includes the full optimal WCMT solution $\mathbf{m}_{\text{OPT}}$.

Chronologically, OL1 is first computed when W phase data for stations within $5^\circ \le \Delta \le 50^\circ$ are available at $t_b \equiv t_0 + 22 \text{ min} (\sim t_0 + t_P(\Delta = 50^\circ) + 15 \text{ s deg}^{-1} \times 50^\circ)$. As described in Section 3.2, the preliminary magnitude delivered in OL1 is estimated from vertical W phase amplitudes in the frequency band 1–5 mHz. It is used to determine the most appropriate initial frequency bandpass prior to the first moment tensor inversion. The frequency passband is further updated using the scalar moment estimate performed in OL2. As detailed in Section 3.4, it is also used to determine the a

priori time-shift and half duration of the MRF. Quickly after OL1, our application delivers OL2, consisting of a first moment tensor solution $\mathbf{m}_{\text{PDE}}$ and the optimum time-shift and half duration of the MRF. If there are enough channels (i.e. number of station $N \ge 30$ and azimuthal gap $\gamma \le 270^\circ$), OL3 is also computed to obtain a preliminary optimized CMT solution. At $t_c \equiv t_0 + 35 \text{ min} (\sim t_0 + t_P(\Delta = 90^\circ) + 15 \text{ s deg}^{-1} \times 90^\circ)$, the second data set for the distance range $5^\circ \le \Delta \le 90^\circ$ is available. Since the surface area covered by this data set is nearly three times larger than that with $\Delta \le 50^\circ$, and is thus far more complete, we perform another inversion including stages OL2 and OL3 (using the final magnitude from the preliminary inversion at $\Delta = 50^\circ$ as a substitution for OL1) which we consider as the final optimum WCMT solution. The solution obtained for $\Delta < 50^\circ$ uses less data but is available 13 min before the solution obtained for $\Delta < 90^\circ$ and as such can be very useful for warning purposes.

At a given centroid location and origin time, we invert for the moment tensor. This is a linear, least-square inversion in which we can optionally impose the trace of the seismic moment tensor to be zero ($\text{trace}(\mathbf{M}) = 0$). This first moment tensor estimation is
quasi-instantaneous. Only when thousands of synthetic seismograms are calculated does the computation time become significant. This is the case for OL3, where we want to perform an inversion for each potential location $\eta_c$ on a 3-D grid surrounding the epicentre $\eta_{PDE}$. Despite the computation-time cost, we prefer this approach to those based on local derivatives of eq. (3) (as in Dziewonski et al. 1981; Dziewonski & Woodhouse 1983), because of its robustness in a radially varying earth. In real-time operations, the computation delay can be reduced by computing Green's functions over a grid geometry defined around the initial PDE received at time $t_a$. Since this computation takes place as the W phase is travelling to the stations, by the time the W phase data are available at $t_b$ and $t_c$, all the Green's functions are ready to be used in the inversions.

In practice, for this study where we are dealing with past events, the data are made available to the application instantaneously as event SEED volumes containing both the seismic traces and the instrumental responses.

## 3.2 Data screening

An important task throughout the inversion process is data screening. Both in real time and for further analysis, we need to handle situations in which some bad traces are included in the data set (e.g. noisy or dead channels, bad instrument responses, incomplete traces, etc.). As described in the next paragraphs, screening filters are set up at different stages of the WCMT algorithm to identify and remove such records.

Following KR, we fit the instrumental response of each sensor within a pre-defined frequency band (0.001–100 Hz) to that of a simple electromagnetic velocimeter with three free parameters: natural period, viscous damping and gain factor. If the result of this fit is good enough (misfit smaller than 3 per cent) the response is stored in a look-up table in the form of the coefficients of the recursive filter to be used for deconvolution; otherwise, the corresponding traces are discarded. In general, the volume of rejected data at this stage is extremely small.

The first screening performed is a ‘noise screening’, which is used to reject the noisiest traces. The pre-event displacement power spectral densities (PSDs) are computed for the whole data set using a duration of 3 hr preceding the origin time. We reject any channel for which the average difference between its noise curve and the New High Noise Model (NHNM; Peterson 1993) in the frequency band 1–10 mHz is positive (i.e. very noisy traces).

Any incomplete trace over the interval $[t_P, t_P + 15 \text{ s deg}^{-1} \times \Delta]$ is also removed from the data set. Although in some cases it could be technically feasible to use such a trace (thanks to the time-domain analysis we use), we prefer not to use them to maintain a simple algorithm and to avoid potential data artefacts.

Next, a ‘median screening’ is used to reject traces associated with incorrect instrument responses and to remove glitchy or dead channels. It is applied after performing the time-domain deconvolution and bandpass filtering, according to the following procedure. For each trace $j$, we compute its peak-to-peak value $p_j$ in the W phase time window. From the complete set of $p_j$, we compute the event’s median value $m$, and reject traces with $p_j$ significantly different from $m$ ($p_j < 0.1 \times m$ or $p_j > 3 \times m$). Although this screening method can accidentally reject some good data (e.g. a nodal station), it has the advantage of being completely independent of the misfit in eq. (3) and does not require any forward modelling.

Finally, we apply a ‘misfit screening’ based on the similarity between observed and synthetic W phase traces. After performing

a WCMT inversion, we can compute the rms misfit according to

$$ \rho_i = \frac{\|\mathbf{s}_w^i - \mathbf{d}_w^i\|^2}{\|\mathbf{s}_w\|^2}, \quad (4) $$

where $\mathbf{s}_w^i$ and $\mathbf{d}_w^i$ are, respectively, the synthetic and data traces corresponding to the $i$th channel. The normalization is used to dampen the effect of the event’s magnitude. $\rho_i$ is then compared with a given threshold $\rho_{max}$. Those stations for which $\rho_i > \rho_{max}$ are removed before restarting a new inversion with the reduced data set. Several thresholds corresponding to increasingly more stringent criteria are successively applied. In the present application, we use three consecutive thresholds: $\rho_{max} = 3.0$, $\rho_{max} = 2.0$ and $\rho_{max} = 1.0$.

Fig. 2 presents the fraction of LH channels that remain after applying these data screening filters. The initial number of files extracted from SEED volumes is presented in black, the number of channels selected for $\Delta \le 90^\circ$ is shown in red and the final number of W traces after the screening processes is presented in blue. On average for $M_w \ge 6.5$ earthquakes occuring between 1990 and 2010, 50 per cent of channels are rejected during the data screening process.

In this work, we define ‘disturbed events’ as any earthquake whose signal is contaminated by the large amplitude waveforms of a preceding event. More precisely, they are defined as events occurring within 1 hr of $M_w \ge 6.5$ events, or less than 10 hr after $M_w \ge 7.0$ earthquakes, and which demonstrate a poor station distribution after the data screening process for $\Delta \le 50^\circ$ (i.e. $N < 30$ or $\gamma > 270^\circ$). The standard W phase algorithm is not well suited to model such events because the assumption of an isolated source in time and space is no longer valid. Using this approach, 44 disturbed earthquakes have been recognized and rejected from our catalogue. The list of ‘disturbed event’ over the period 1990–2010 is detailed in Table S1. In Section 4.5, we explore a possible scheme to handle them in real time.

## 3.3 Preliminary W phase magnitude estimation (OL1)

At $t_b = t_0 + t_P(\Delta = 50^\circ) + 15 \text{ s deg}^{-1} \times 50^\circ \sim t_0 + 22 \text{ min}$, W phase traces for all stations within $\Delta < 50^\circ$ are available and the first data set can be built. Before trying a formal inversion for the moment tensor, we perform a first-order fit of the W phase amplitudes as a function of distance and azimuth. Following KR, the idea here is to capture the information carried by the overall vertical amplitude of W phase and to translate it into magnitude.

After instrument correction of vertical component data and bandpass filtering in the 1-5 mHz range, we remove incomplete traces and apply a ‘median screening’ to remove conspicuous outliers. We then measure the peak-to-peak value $p_j$ on each W trace $j$. These amplitudes are then reduced to a common distance ($\Delta = 40^\circ$). This procedure is similar to the Richter Magnitude original definition (Richter 1935). To capture the overall amplitude level while allowing some azimuthal variations due to the mechanism, these reduced amplitudes are matched to a two-lobe azimuthal pattern corresponding to a pure thrust or normal-fault earthquake

$$ p_j = q(\Delta_j)[a - b \cos^2(\Phi_j - \Phi_0)], \quad (5) $$

where $q(\Delta_j)$ is the W phase amplitude decay (see Table 2 of KR) and $a, b$ and $\Phi_0$ the parameters to be determined. Eq. (5) can be solved as a linear least-square problem by inverting for $a - b/2$ (the average amplitude), $b \cos(2\Phi_0)/2$ and $b \sin(2\Phi_0)/2$. The resulting average amplitude $a - b/2$ can then be used as a direct measure of the seismic moment. It is also useful to solve for $\Phi_0$ if we want to
obtain a rough estimate of the fault strike $\phi$. This choice of a two-lobe azimuthal pattern associated with a pure thrust mechanism is motivated by the fact that the W phase algorithm is primarily targeted at the inversion of large tsunamigenic earthquakes.

The purpose of this preliminary estimation is twofold. First, it provides a quick, simple and robust magnitude estimation which is independent of any additional hypothesis and modelling details; this is the first output, OL1, of our algorithm. Secondly, the magnitude so obtained can be used as a proxy for the initial estimate of the duration of the moment rate function to be used in subsequent stages. Note that neither the focal mechanism nor the centroid depth are needed for this preliminary magnitude estimation.

### 3.4 MT inversion at PDE (OL2)

After performing the ‘noise screening’, rejecting channels showing a bad instrumental response fit or with truncated records and applying a ‘median screening’, we perform a first moment tensor inversion. As is typical in moment tensor inversion algorithms, we impose a zero trace to the moment tensor to cope with the poor resolution of the isotropic components for shallow earthquakes (Mendiguren 1977). It constrains the seismic source to have no net volume change.

Besides the waveforms, we need a centroid location, centroid time and an MRF duration. At $t_b \sim t_0 + 22$ min, we use the PDE location as our best guess for the centroid location, whereas the duration is estimated by a scaling law from the seismic moment obtained in OL1

$$h_c = 1.2 \times 10^{-8} \times M_0^{1/3}, \qquad (6)$$

where $h_c$ is in seconds and $M_0$ in dyne-cm. This relation is obtained from the constant stress drop scaling relation $h_c \propto M_0^{1/3}$ (Kanamori & Anderson 1975). The constant of proportionality is set so that an earthquake with $M_0 = 10^{27}$dyne-cm ($M_w = 7.3$) has a half duration of 12s. At $t_c \sim t_0 + 35$ min, we use the duration obtained at OL2 for $\Delta \le 50^\circ$.

We can then select the corresponding Green’s functions from the database, convolve them with the MRF shape and apply the same bandpass filter as applied to the data.

One more parameter is necessary to compute synthetic traces directly comparable to the waveform data: the delay $\tau_c$ between origin time (i.e. PDE) and the centroid time. This is determined with a grid-search by performing several moment tensor inversions for a range of trial values of $\tau_c$. This is an inexpensive operation since changing this delay requires simply a time-shift of the Green’s functions or of the data traces in the opposite direction. As a result of this grid-search, we obtain an optimal delay value. In contrast to $\tau_c$, the MRF duration $h_c$ is poorly constrained by the waveforms, being generally significantly smaller than the long periods of the W phase. We thus use the optimal delay value as a new proxy for $h_c$ (i.e. we assume that $h_c = \tau_c$). With these parameters, we compute three successive moment tensor inversions using an increasingly stringent ‘misfit screening’ with thresholds $\rho_{max} = 3.0$, $\rho_{max} = 2.0$ and $\rho_{max} = 1.0$ for the channel rms misfit in eq. (4). The resulting moment tensor solution is our second-level output, OL2.

### 3.5 Optimized CMT inversion (OL3)

After determining the optimum centroid time and MRF half duration, we attempt to find a centroid location which is better than the preliminary location estimate. For this purpose, we set up a 3-D grid-search (latitude–longitude-depth), where each grid-node is used as a potential centroid location and a complete WCMT inversion is made. The rms misfit in eq. (3) is used as an objective function to choose the optimal centroid location. To make certain that the rms values for different centroids are comparable, we must use the same data set. For this reason, we do not apply any additional screening at this level. The typical dimension of the grid is $2.4^\circ \times 2.4^\circ \times 100$ km, centred on the PDE location, and the minimum allowed centroid depth is 12 km. The depth step ($\Delta h$) is variable with the centroid depth ($h$)

<figure><img src="image_5.png" /><figcaption>**Figure 4.** WCMT solutions obtained for $M_w \ge 6.50$ earthquakes from 1990 to 2010. The stations within $\Delta \le 90^\circ$ epicentral distance are used. The complete collection of solutions is presented in Fig. 5 and in the Supporting Information (Figs S2–S11).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Solutions obtained for events with $M_w \ge 7.59$. The events are listed in order of decreasing GCMT magnitude. *N* is the total number of channels after inversion and *γ* is the azimuthal gap (in degrees). GCMT solutions are presented in green, W phase solutions corresponding to OL2 (i.e. for *N* < 30 or *γ* > 270°) are shown in blue and WCMT solutions corresponding to OL3 (i.e. for *N* ≥ 30 and *γ* ≤ 270°) are presented in red. The complete collection of solutions obtained for $M_w \ge 6.5$ earthquakes are listed in the Supporting Information (Figs S1–S11).</figcaption></figure>
<figure><img src="image_7.png" /></figure>

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

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

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

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

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

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

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

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

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

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

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

**Figure 6.** Examples of W phase preliminary magnitude ($M_{w-wprel}$) amplitude fits. The peak-to-peak amplitudes are reduced to a common distance $\Delta = 40^\circ$. A polar representation is used: the angle and radius correspond, respectively, to the station azimuth and amplitude value. The continuous line represents the result of the regression determined from W phase amplitudes within $\Delta \le 50^\circ$. The 12 examples shown are ordered by decreasing GCMT magnitude $M_{w-gcmt}$. GCMT solutions are displayed in light grey for comparison with peak-to-peak amplitude azimuthal distributions. In addition to the good match between $M_{w-wprel}$ and $M_{w-gcmt}$, there is a good agreement between the focal mechanism and the orientation of the two lobes pattern.
(1) $\Delta h = 2$ km for $h \le 25.5$,

(2) $\Delta h = 5$ km for $25.5 \le h \le 50.5$,

(3) $\Delta h = 10$ km for $h \ge 50.5$.

A multiple-scale grid-search is performed for each depth: First, a global exploration of the latitude–longitude space is conducted using a large sampling step (40 km). We then select several locations which represent the best least-squares misfits between observed and calculated waveforms. Another exploration is then performed around these optimal points by increasing the horizontal sampling resolution (10 km). The initial grid size is increased if one of the chosen locations is within one cell of the grid edge. Finally, we choose the centroid depth, latitude and longitude which minimize the rms misfit in eq. (3) and take them as the optimum WCMT centroid (OL3).

## 4 RESULTS

In this section, we present the results of applying the protocol defined earlier to earthquakes with $M_w \ge 6.50$ since 1990 (i.e. 815 events), and systematically compare them with the GCMT solutions. Fig. 4 shows the global distribution of the WCMT mechanisms. In Fig. 5, we present detailed solutions for $M_w \ge 7.59$ events together with GCMT solutions for comparison. In the online Supporting Information, we provide WCMT and GCMT solutions for $M_w \ge 6.5$ earthquakes (Figs S1–S11) and the solutions resulting from the extension to $6.0 \le M_w < 6.5$ events (Figs S13–S36). Solutions for $M_w \ge 6.5$ earthquakes using data within $\Delta < 90^\circ$ are also available at: http://eost.u-strasbg.fr/wphase/MGE65.

### 4.1 Preliminary magnitude estimation

In Fig. 6, 12 examples of distance-corrected amplitude–azimuth fits are presented. In each panel, the continuous line represents the regression whereas the coloured bars indicate the corrected peak-to-peak values at different epicentral distances. The fits are generally good, even for stations at large epicentral distances ($\Delta > 50^\circ$),

<figure><img src="image_19.png" /><figcaption>**Figure 7.** Amplitude ratios (observed/predicted) as a function of azimuth from the epicentre of the 2004 Sumatra–Andaman Islands earthquake. The observed amplitudes are compared with predictions from the GCMT solution indicated in Fig. 1 for the W phase in red and for R1 and R2 Rayleigh waves in blue. Large amplitude ratios observed towards the north-northwest are interpreted in terms of directivity along the rupture direction of ~330°N.</figcaption></figure>

which are not included in the inversion. In addition to the good match between these W phase preliminary magnitudes ($M_{w-wprel}$) and the GCMT magnitudes ($M_{w-gcmt}$), we note that the orientation of the lobes specified by $\Phi_0$ in eq. (5) is often a good estimate of the fault strike.

The assumed two-lobe azimuthal pattern corresponding to a thrust source is clearly less convenient for strike-slip events, since the W phase radiation pattern has four lobes for such earthquakes. However, even if detailed azimuthal variations of the W phase amplitude are not properly reproduced, we note that the average amplitude is sufficiently well modelled for the purpose of preliminary moment magnitude determination, as we can see in Figs 6(h) and (k) for the 2001 Kokoxili and 2010 Haiti earthquakes, respectively.

In Fig. 6(a), we note that the peak-to-peak values observed for the 2004 Sumatra–Andaman Islands earthquake are larger at stations to the north of the epicentre. This amplification cannot be related to structural heterogeneities since W phase mainly propagates through the upper mantle and thus should not be affected by shallow crustal structures. To estimate directivity effects on W phase, we computed synthetics for point source models and simply measure the amplitude ratios (observed/predicted) at each station. For comparison, we also conduct this experiment using R1 and R2 Rayleigh waves after bandpass filtering in the 2.4–10 mHz passband. If there were no effects of directivity, the resulting amplitude ratios shown in Fig. 7 would be unity for all azimuths. We note that both W phase and R1–R2 amplitudes are enhanced towards the north-northwest, which is compatible with the rupture propagation of ~330°N, as discussed for Rayleigh waves by Ammon et al. (2005). However,

<figure><img src="image_20.png" /><figcaption>**Figure 8.** W phase preliminary magnitudes for $M_w \ge 6.5$ earthquakes in the period 1990–2010. The distribution of magnitude difference ($\Delta M_w = M_{w-wprel} - M_{w-gcmt}$) is presented in (a). The comparison of the preliminary magnitude estimates ($M_{w-wprel}$) versus the GCMT magnitude ($M_{w-gcmt}$) is shown in (b). The circles are coloured according to the GCMT centroid depth. Dashed lines indicate $\pm 0.1$ and dot-dashed lines $\pm 0.2$ magnitude units.</figcaption></figure>
<table><caption>**Table 3.** Statistical comparison between WCMT and GCMT solutions for the preliminary W phase magnitude estimation, first and second WCMT inversions. The standard deviation ($\sigma_{M_w}$), the mean and median of the magnitude difference ($\Delta M_w = M_{w-wprel} - M_{w-gcmt}$), as well as the proportion of events showing $\Delta M_w < 0.2$ and $\Delta M_w < 0.1$ are presented.</caption><thead><tr><th></th><th>σ<sub>M<sub>w</sub></sub></th><th>Mean (ΔM<sub>w</sub>)</th><th>Median (ΔM<sub>w</sub>)</th><th>ΔM<sub>w</sub> &lt; 0.2</th><th>ΔM<sub>w</sub> &lt; 0.1</th></tr></thead><tbody><tr><td>Preliminary M<sub>w</sub> estimation</td><td>0.15</td><td>0.09</td><td>0.08</td><td>73 per cent</td><td>43 per cent</td></tr><tr><td>First WCMT solution</td><td>0.05</td><td>0.01</td><td>0.00</td><td>98 per cent</td><td>87 per cent</td></tr><tr><td>Second WCMT solution</td><td>0.04</td><td>0.00</td><td>0.00</td><td>99 per cent</td><td>93 per cent</td></tr></tbody></table>

the azimuthal variation of amplitude ratio for the W phase is not as large as for Rayleigh waves, which indicates that the W phase is much less sensitive to source directivity as a consequence of its high group velocity and long period.

A systematic comparison between $M_{w-wprel}$ and $M_{w-gcmt}$ is shown in Fig. 8; an overall statistical comparison is presented in Table 3. Because of its long-period character, $M_{w-wprel}$ is quite accurate for large earthquakes ($M_w \ge 7.5$). However, as a consequence of the background seismic noise in the 1–5 mHz passband, $M_{w-wprel}$ generally overestimates $M_w$ for smaller events.

## 4.2 First WCMT inversion: $t_b \sim t_0 + 22 \text{ min}$

WCMT solutions obtained at $t_b \sim t_0 + 22 \text{ min}$ (i.e. using stations at $\Delta \le 50^\circ$) are shown in Fig. 5 and from Fig. S1 to S11.

Fig. 9 compares $M_{w-gcmt}$ with $M_{w-gcmt}$ for the complete set of events, using stations within $\Delta \le 90^\circ$. Overall, the magnitudes are estimated extremely well: 98 per cent of inversions yield a $M_{w-gcmt}$ within $M_{w-gcmt} \pm 0.2$, and the standard deviation between the two moment magnitude estimates is $\sigma_{M_w} = 0.05$ as indicated in Table 3.

The reliability of the solution obtained at $\Delta \le 50^\circ$ is particularly important for tsunami warning purposes since it is available only 22 min after origin time. Tsunami earthquakes are of particular interest because they are generally associated with an anomalous moment rate spectrum, which complicates the rapid characterization of their source. The determination of an accurate time-shift is important for these events because of their anomalously long rupture duration. The results obtained for two examples are detailed hereafter.

### 4.2.1 Tsunami earthquake example: 1992 Nicaragua earthquake ($M_w = 7.6$, id = 090292A)

Kanamori (1993) gave the name W phase to a prominent long-period wave observed on several records of this earthquake. The WCMT solutions obtained for this event are included in Fig. 5 (first column of penultimate row) and in Fig. 10. The preliminary W phase magnitude is $M_{w-wprel} = 7.5$. Based on this estimate, ground motions are bandpass filtered in the 2.0–6.7 mHz passband. Fig. 10 shows the W phase solutions obtained at different stages of the algorithm. The dark green focal mechanism in Fig. 10(b) corresponds to the solution obtained with the PDE location using the initial estimate of the centroid timing from eq. (6): $\tau_c = h_c = f(M_{w-wprel}) \sim 16 \text{ s}$. Although $M_{w-gcmt}$ matches $M_{w-gcmt}$ well, the strike is very different from the GCMT solution. The blue mechanism in Fig. 10(c) is the solution obtained after the $\tau_c$ optimization (shown in Fig. 10e). This operation significantly improves the solution and the resulting $\tau_c = 48 \text{ s}$ is comparable to the GCMT estimate. The red mechanism in Fig. 10(d) corresponds to the solution obtained after the centroid location optimization.

<figure><img src="image_21.png" /><figcaption>**Figure 9.** W phase moment magnitudes ($M_{w-wcmt}$) obtained at $t_b \sim t_0 + 22 \text{ min}$ (i.e. using stations within $\Delta \le 50^\circ$). The distribution of magnitude difference ($\Delta M_w = M_{w-wcmt} - M_{w-gcmt}$) is presented in (a). The comparison of $M_{w-wcmt}$ versus the GCMT magnitude ($M_{w-gcmt}$) is shown in (b). W phase solutions corresponding to OL2 (i.e. for $N<30$ or $\gamma > 270^\circ$) are indicated by triangles and WCMT solutions corresponding to OL3 (i.e. $N \ge 30$ and $\gamma \le 270^\circ$) are indicated by circles. The symbols are coloured according to the GCMT centroid depth. Dashed lines indicate $\pm 0.1$ and dot-dashed lines $\pm 0.2$ magnitude units. At $t_b \sim t_0 + 22 \text{ min}$, 87 per cent of $M_{w-gcmt}$ values lie within $M_{w-wcmt} \pm 0.1$ and the standard deviation between the two moment magnitude estimates is $\sigma_{M_w} = 0.05$.</figcaption></figure>

### 4.2.2 Tsunami earthquake example: 2006 Java earthquake ($M_w = 7.7$, id = 200607170819A)

The WCMT solutions obtained for this event are presented in column 4 and row 10 of Fig. 5. The W phase preliminary magnitude for this event is $M_{w-wprel} = 7.4$ (cf. Fig. 6i). This magnitude estimation underestimates $M_w$ but remains acceptable when compared with the $M_{wprel} = 7.2$ estimated in early warning reports (PTWC 2006a,b). The dark green mechanism in Fig. 11(b) shows the first WCMT solution obtained with the PDE location and assuming the initial timing estimates $\tau_c = h_c \sim 14 \text{ s}$ derived from eq. (6). Because of the anomalously long-rupture duration for this earthquake, the
<figure><img src="image_22.png" /><figcaption>**Figure 10.** WCMT solutions calculated at time $t_b \sim t_0 + 22$ min for the 1992 Nicaragua tsunami earthquake. The solutions obtained at different stages of the WCMT algorithm are presented in (b), (c) and (d) and the GCMT solution is shown for comparison in (a). The time-shift ($\tau_c$) grid-search is detailed in (e), where the black cross corresponds to the initial estimate of $\tau_c$ from eq. (6), the red dot corresponds to the optimum $\tau_c$ and the green line indicates the GCMT time delay. Although the strike of solution (b) does not match the WCMT solution well, the solution (c) corresponding to an optimum $\tau_c$ and (d) after centroid position ($\eta_c$) grid-search are very similar to the GCMT mechanism.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 11.** WCMT solutions calculated at time $t_b \sim t_0 + 22$ min for the 2006 Java tsunami earthquake. The solutions obtained at different stages of the WCMT algorithm are presented in (b), (c), (d) versus the GCMT solution is shown for comparison in (a). The time-shift ($\tau_c$) grid-search is detailed in (e), where the black cross corresponds to the initial estimate of $\tau_c$ from eq. (6), the red dot corresponds to the optimum $\tau_c$ and the green line indicates the GCMT time delay. The WCMT solution obtained after centroid location ($\eta_c$) optimization (d) is very similar to the GCMT solution (a).</figcaption></figure>

empirical relation (6) does not hold and $\tau_c$ is clearly underestimated. As a consequence, this initial WCMT solution differs substantially from the GCMT solution. After the time delay optimization, the WCMT solution shown in Fig. 11(c) (blue mechanism) is signif-

icantly changed and the estimated $\tau_c \sim 74$ s agrees well with the GCMT value (cf. Fig. 11e). After a centroid position grid-search, the estimated magnitude $M_{w-wcmt} = 7.76$ is remarkably similar to $M_{w-gcmt}$, as shown in Fig. 11(d).
## 4.3 Second WCMT inversion: $t_c \sim t_0 + 35$ min

The WCMT solutions obtained at $t_c \sim t_0 + 35$ min (i.e. using stations up to $\Delta \le 90^\circ$) correspond to the mechanisms shown in the middle of each frame in Fig. 5, and from Figs S1–S11. The complete collection of solutions is also available at: http://eost.ustrasbg.fr/wphase/MGE65.

Fig. 12 compares $M_{w-wcmt}$ with $M_{w-gcmt}$ for the total set of events using stations within $\Delta \le 90^\circ$. Moment magnitudes are remarkably well determined: 99 per cent of the inversions yield a $M_{w-wcmt}$ within $M_{w-gcmt} \pm 0.2$ and the moment magnitude standard deviation is $\sigma_{M_w} = 0.04$, as presented in Table 3. Fig. 13 presents a comparison between WCMT and GCMT moment tensor elements. Although the dispersion increases for the smaller components, the components $M_{rr}$, $M_{\theta\theta}$, $M_{\phi\phi}$ and $M_{\theta\phi}$ estimated from WCMT and GCMT are very similar. It is well known that the elements $M_{r\theta}$ and $M_{r\phi}$ are difficult to resolve at long period since the relative amplitudes of the associated excitation kernels are very small for shallow earthquakes (Kanamori & Given 1981). A well-known consequence of this difficulty is a tradeoff at long periods between the dip $\delta$ and the scalar moment $M_0$ for shallow earthquakes, such that for pure thrust or pure normal-fault events the product $M_0 \sin \delta$ is well determined, but $M_0$ and $\delta$ are difficult to resolve separately. For $M_w < 7.5$ events, the GCMT estimates for $M_{r\theta}$ and $M_{r\phi}$ are expected to be

<figure><img src="image_24.png" /><figcaption>**Figure 12.** W phase moment magnitudes ($M_{w-wcmt}$) obtained at $t_c \sim t_0 + 35$ min (i.e. using stations within $\Delta \le 90^\circ$). The distribution of magnitude difference ($\Delta M_w = M_{w-wcmt} - M_{w-gcmt}$) is presented in (a). The comparison of $M_{w-wcmt}$ and the GCMT magnitude ($M_{w-gcmt}$) is shown in (b). W phase solutions corresponding to OL2 (i.e. for $N < 30$ or $\gamma > 270^\circ$) are indicated by triangles and WCMT solutions corresponding to OL3 (i.e. $N \ge 30$ and $\gamma \le 270^\circ$) are indicated by circles. The symbols are coloured according to the GCMT centroid depth. Dashed lines indicate $\pm 0.1$ and dot-dashed lines $\pm 0.2$ magnitude units. At $t_c \sim t_0 + 35$ min, 92 per cent of $M_{w-gcmt}$ values lie within $M_{w-wcmt} \pm 0.1$ and the standard deviation between the two moment magnitude estimates is $\sigma_{M_w} = 0.04$.</figcaption></figure>

more robust since they are obtained by combining seismic phases in different frequency bands (body waves, surface waves and mantle waves). However, this trade-off affects the GCMT solutions for large events ($M_w \ge 7.5$) since, in this case, the algorithm is based only on mantle waves filtered at long period. Even though $M_{r\theta}$ and $M_{r\phi}$ are not as well constrained as the other moment tensor components, we finally find the WCMT estimates to be comparable to GCMT.

To quantitatively compare WCMT solutions with GCMT, we use the magnitude difference $\Delta M_w = M_{w-wcmt} - M_{w-gcmt}$ and the angular parameter $\Phi(\mathbf{M}_{wcmt}, \mathbf{M}_{gcmt})$. $\Phi(\mathbf{M}_A, \mathbf{M}_B)$ is the angle of the smallest rotation linking the two sets of $\mathbf{M}_A$ and $\mathbf{M}_B$ principal axes (cf. Appendix). In the case of double-couples, $\Phi(\mathbf{M}_A, \mathbf{M}_B)$ is the angle of the smallest rotation that, when applied to $\mathbf{M}_B$, produces the mechanism of $\mathbf{M}_A$. Fig. 14 shows the variation of $\Delta M_w$ and $\Phi$ as a function of time. Between 1990 and 2010, we find that $\Phi < 20^\circ$ for 85 per cent of the events. A small number of events before 1992 show relatively large $\Phi$ values ($\Phi \ge 50^\circ$). All of these earthquakes are represented by triangles in Fig. 14(b), indicating poor data coverage. There are four events of $M_w > 7.5$ with $\Phi \ge 50^\circ$ or $|\Delta M_w| \ge 0.2$ since 1990:

(1) The first is the 2000 June 4 Enggano earthquake ($M_w = 7.9$, id=060400D, Fig. 5 column 2 and row 6), which is associated with $\Delta M_w = 0.27$ and $\Phi = 82^\circ$. The unusual character of this event is its strike-slip mechanism with a large non-double-couple component having an opposite sign for WCMT and GCMT solutions. The source complexity of this earthquake has been studied by Abercrombie (2003), which suggests a combination of two separate subevents having strike-slip and thrust mechanisms.

(2) The second and third outliers are part of the 2000 Papua New Guinea earthquake sequence. The largest event of this sequence is the 2000 November 16 event ($M_w = 8.0$, id=111600B, Fig. 5 column 1 and row 4). Although $M_{w-wcmt}$ matches $M_{w-gcmt}$ very well ($\Delta M_w = 0.03$), this earthquake shows significant differences between GCMT and WCMT mechanisms ($\Phi = 69^\circ$). The strike-slip mechanism has a northeast dipping plane according to GCMT and a south dipping plane for WCMT. Our results are somewhat more consistent with the body wave solution provided by the Japanese Earthquake Research Institute (ERI) which proposes a strike-slip mechanism on a nearly vertical fault (ERI 2000). The ERI analysis also suggests some source complexity, with a dip-slip subevent occurring 2 min after the main rupture begins.

(3) The last large earthquake of the 2000 Papua New Guinea sequence had a thrust mechanism and occurred on 2000 November 17 ($M_w = 7.5$, id=111700Q, Fig. 5 column 5 and row 7). The match between GCMT and WCMT mechanisms is quite good for this event ($\Phi = 11^\circ$) even though there is a small non-double-couple component in the GCMT solution. The estimated magnitudes are however quite different: $M_{w-wcmt} = 7.52$ and $M_{w-gcmt} = 7.77$ (i.e. $\Delta M_w = -0.25$). Our solution is similar to USGS CMT ($M_w = 7.4$), ERI ($M_w = 7.3$) and Park & Mori (2007) ($M_w = 7.4$). The observed magnitude difference cannot be fully explained by the dip difference between the GCMT and WCMT solutions ($\delta_{wcmt} = 26^\circ$ and $\delta_{gcmt} = 24^\circ$). In fact, running the WCMT inversion by fixing the centroid location and dip values to the GCMT solution gives $M_w = 7.55$ (i.e. $\Delta M_w = -0.22$). As suggested, for example, by Geist (2005), the discrepancy among seismic moment estimates using different passband and time windows can possibly relate to the complexity of this event.

(4) The last outlier is the 2009 Samoa event ($M_w = 8.1$, id=200909291748A, Fig. 5 column 4 and row 3). The differences
<figure><img src="image_25.png" /><figcaption>**Figure 13.** W phase moment tensor components obtained at $t_c \sim t_0 + 35$ min (i.e. using stations within $\Delta \le 90^\circ$). The comparisons of the moment tensor elements (in dyne-cm) of WCMT (vertical axis) and GCMT (horizontal axis) are given for (a) $M_{rr}$, (b) $M_{\theta\theta}$, (c) $M_{\phi\phi}$, (d) $M_{r\theta}$, (e) $M_{\theta\phi}$ and (f) $M_{\phi\phi}$. Both WCMT and GCMT solutions are obtained by inverting for the deviatoric moment tensor. The blue crosses indicate that the moment tensor component of WCMT and GCMT are of opposite sign, and green symbols, the same sign. See also Fig. 14(b).</figcaption></figure>

observed between the WCMT and GCMT mechanisms ($\Phi = 51^\circ$, $\Delta M_w = -0.03$) can be related to the source complexity of this earthquake. As noted by Lay *et al.* (2010), this earthquake is a highly complex event starting with an outer rise subevent ($M_w = 8.1$) and is followed within 2 min by two major interplate subevents (both with $M_w = 7.8$).

The great-circle distances between WCMT, GCMT and PDE locations are presented in Fig. 15. The directions between these locations are given in Fig. S12 and the WCMT centroid locations are listed from Table S2 to S3. Compared with the PDE location, the WCMT location obtained after grid-search is closer to the GCMT centroid for most events: the distance between the GCMT and WCMT centroid is less than 50 km for 91 per cent of the events. A significant part of WCMT centroid mislocations (with respect to GCMT centroids) can be related to limited spacial resolution of the very long periods being used. In Fig. 15(b), the largest distances are observed for earthquakes in South America. For events in the western part of the Americas, the GCMT centroids are on average shifted 30 km to the west with respect to WCMT locations (cf. Fig. S12b). As discussed by Hjörleifsdóttir & Ekström (2010), the GCMT locations in this region are biased ~15 km to the west due to unmodelled structural heterogeneity Dziewonski & Anderson (1981). We performed synthetic tests considering a strongly unbalanced network geometry in the east-west direction (poor coverage to the west) and adding long-period noise to PREM synthetics. Our results showed

a centroid shift of about 15 km to the east. This effect can thus contribute to the systematic east-west bias observed between GCMT and WCMT centroid locations in this region. We also tested a data weighting scheme proportional to the inverse of the azimuthal density of stations. This approach strongly dampens the effects of an unbalanced network, reducing the final WCMT-GCMT offsets to less than our Green's function grid spacing.

Focal depth is important for local warning purposes since it has a great influence on near-field tsunami wave height. For a far-field tsunami however, the depth dependence is reduced, as shown by Okal (1988). In fact, to the first order the size of great earthquakes is well correlated to the maximum amplitude of far-field tsunami waves regardless of the event depth (Abe 1979, 1989). In this global application, WCMT solutions are available 22–35 min after the origin time. This delay corresponds to tsunami arrivals at distances greater than 250 km and a very accurate depth estimation is no longer essential at this distance range. However, as the scalar moment and fault dip could vary with depth, the centroid depth determination has been included in the algorithm. The depth differences between WCMT and GCMT locations are presented in Fig. 16. Although the depth resolution is limited because of the very long wavelengths involved, we note in Fig. 16(c) that the depth differences between WCMT and GCMT solutions are clearly reduced after the spatial grid-search. The average difference between WCMT and GCMT centroid depths is +9.6 km and 90 per cent of the inversions yield a depth within the GCMT value ±11 km.
<figure><img src="image_26.png" /><figcaption>**Figure 14.** Difference between WCMT and GCMT solutions. The magnitude difference $\Delta M_w$ is presented in (a), the angle $\Phi(M_{wcmt}, M_{gcmt})$ is shown in (b) whereas (c) displays the azimuthal gap $\gamma$ (blue) and the number of channels (red) after the data screening. In (a) and (b), the W phase solutions corresponding to OL2 (i.e. for $N < 30$ or $\gamma > 270^{\circ}$) are indicated by triangles and WCMT solutions corresponding to OL3 (i.e. $N \ge 30$ and $\gamma \le 270^{\circ}$) are indicated by circles. The symbols are coloured according to the GCMT moment magnitude $M_{w-gcmt}$. Events prior to 1992 generally have poor data coverages.</figcaption></figure>

### 4.3.1 Example: 2004 Sumatra–Andaman Islands earthquake ($M_w = 9.2$, id = 122604A)

The Sumatra–Andaman Islands 2004 megathrust earthquake is the largest event within the time period considered in this study. The solutions obtained at 50° and 90° are presented in the top left-hand corner of Fig. 5. The dip of the final ($\Delta \le 90^{\circ}$) WCMT best double-couple is similar to GCMT ($\delta_{wcmt} = 7^{\circ}$ versus $\delta_{gcmt} = 8^{\circ}$). The centroid grid-search is presented in Fig. 17(a). The depth difference at the optimum centroid is $\Delta r_c = r_{c-gcmt} - r_{c-wcmt} \approx +15$ km, with a horizontal offset of about 60 km, which amounts to 5 per cent of the rupture length. This offset is difficult to interpret since the point source centroid location for such a large and complex rupture is somewhat meaningless. Fig. 1 shows the waveform fit associated with the corresponding optimum WCMT solution. The fits are quite good in the W phase time window, while the later surface wave arrivals are affected by instrument saturation due to

their large amplitudes. The WCMT inversion for this event leads to a magnitude difference of 0.15 units: $M_{w-wcmt} = 9.15$ while $M_{w-gcmt} = 9.00$. As discussed by Ammon *et al.* (2005), Park *et al.* (2005) and Lay *et al.* (2005), the amplitudes observed at very long periods (i.e. >500 s) indicate a moment magnitude larger than 9.0 by 0.1–0.3 units depending on the assumed fault dip. The magnitude difference between GCMT and WCMT can only be partly explained by the dip difference (i.e. $\delta_{gcmt} = 8^{\circ}$ and $\delta_{wcmt} = 7^{\circ}$). By running the W phase inversion with a fixed dip $\delta_{wcmt} = \delta_{gcmt} = 8^{\circ}$ (the dip of the shallow plane of the GCMT best double-couple) and using the GCMT centroid location, $M_{w-wcmt}$ is still greater than $M_{w-gcmt}$ by more than 0.1 magnitude units. It is thus reasonable to explain the greater $M_{w-wcmt}$ estimates at 50° and 90° by the longer period character of the WCMT algorithm. The very long periods considered in our approach (~600 s) are certainly more representative of the entire rupture process than those used for GCMT, which is limited to shorter periods (~300 s for mantle waves).
<figure><img src="image_27.png" /><figcaption>**Figure 15.** Distances between WCMT, GCMT and PDE locations. The great-circle distance $\Delta x_{PDE}$ between the PDE location and the GCMT centroid is shown on the map (a) and in the histogram with blue bars (c). The distance $\Delta x_c$ between GCMT and WCMT centroid locations is presented in (b) and in the histogram with red bars (c).</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Figure 16.** Depth difference between GCMT and WCMT locations before and after grid-search. The difference $\Delta r_{PDE} = r_{c-gcmt} - r_{PDE}$ between the PDE and the GCMT depth (i.e. before grid-search) is shown in (a) and in the histogram with blue bars (c). The depth difference $\Delta r_c = r_{c-gcmt} - r_{c-wcmt}$ between GCMT and WCMT centroids (i.e. after grid-search) is presented in (b) and in the histogram with red bars (c).</figcaption></figure>

### 4.3.2 Example: 2010 Haiti earthquake ($M_w = 7.0$, id = 201001122153A)

The 2010 Haiti earthquake was the deadliest earthquake since the 2004 Sumatra–Andaman Islands event, with more than 300 000

fatalities according to the official estimates (USGS 2010). This highlights the fact that even moderate size earthquakes ($M_w \le 7.0$) can cause major human casualties if they occur near large population centres with poor building construction practices. Although the proximity of the event to populated areas prevents any early
<figure><img src="image_29.png" /><figcaption>**Figure 17.** Centroid location optimizations for (a) the 2004 Sumatra–Andaman Islands earthquake, (b) the 2010 Haiti event and (c) the 2010 Maule–Chile earthquake. We illustrate the latitude/longitude grid-search at the WCMT optimum centroid depth by showing the normalized rms contour-lines. The black cross represents the PDE epicentre whereas the red dot corresponds to the WCMT optimum point source location. The GCMT centroid is indicated by a green triangle. The locations of $M_w \ge 5.0$ aftershocks recorded during the 24 hr after the main shock are indicated by white dots (from the USGS PDE catalogue).</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 18.** W phase source inversion results for the 2010 Haiti earthquake. W phase CMT (WCMT) and Global CMT (GCMT) solutions are shown in the topleft corner. Examples of observed waveforms (black lines) and the corresponding synthetics (red lines) computed from the WCMT solution are presented. The station azimuth ($\phi$) and epicentral distance ($\Delta$) are indicated as well as W phase time windows, bounded by red dots. The WCMT inversion is based on the ground motion recorded at stations within $\Delta \le 90^\circ$ after applying a bandpass filter in the 2–8.3 mHz passband. W phase and later arrivals are very well predicted by the WCMT solution.</figcaption></figure>

warning alerts based on global scale systems, a quick determination of the earthquake source is of major interest for planning rapid rescue operations. The WCMT solutions obtained at 50° and 90° are presented in column 4 and row 2 of Fig. S5. The waveform fits presented in Fig. 18 shows that the W phase—and later arrivals—are

very well predicted by the WCMT solution. The W phase estimate $M_{w- WCMT} = 7.04$ yields a magnitude difference of $\Delta M_w = 0.01$ with respect to GCMT and we note also that the mechanisms are very similar. Fig. 17(b) presents the centroid position grid-search. Both the GCMT and WCMT best point source locations are shifted to
the north of the aftershock cloud. The great-circle distance between WCMT and GCMT centroids is about 20 km, with a small depth difference of $\Delta r_c \approx -1.5$ km.

### 4.3.3 Example: 2010 Maule-Chile earthquake ($M_w = 8.8$, id = 201002270634A)

The 2010 Maule-Chile earthquake, being the first $M_w \ge 8.5$ event since the W phase algorithm was established in real-time operation, is particularly interesting for evaluating the WCMT method (Duputel et al. 2010). The solutions obtained at 50° and 90° are presented in the second frame from left on the top row of Fig. 5. The waveform fit is shown in Fig. 19 for a sample of channels within $\Delta \le 90°$. The W phase and later arrivals are very well predicted by the WCMT solution. We also note that some channels are saturated at the surface wave arrivals. The saturation affects most of the broadband data within $\Delta \le 15°$, but also distant stations like CRZF-LHZ, which is located at $\Delta \approx 84°$. Thanks to the small amplitudes of the W phase and since the instrument response deconvolution is performed by using a time-domain IIR filter, the W phases are not affected at all by the saturation of the seismometer. The final $M_{w-wcmt}$ estimation is similar to $M_{w-gcmt}$ with a magnitude difference of $\Delta M_w = 0.06$. This difference can be fully explained by the dip difference between GCMT and WCMT solutions of $\delta_{wcmt} = 14°$

and $\delta_{gcmt} = 18°$. The centroid search is shown in Fig. 17(c) for the optimal centroid depth (30 km). The great-circle distance between WCMT and GCMT centroids is 40 km with a depth difference of $\Delta r_c \approx -7$ km.

## 4.4 The extension of W phase to lower magnitudes

We follow the same approach as outlined earlier to invert for the W Phase CMT solution for all earthquakes with $6.0 \le M_w < 6.5$ (1727 events). Though the application of the inversion here is identical to the analysis of $M_w \ge 6.5$ events, we treat these smaller earthquakes separately because our approach no longer produces a complete catalogue with this magnitude threshold, nor does it perform as consistently well as for larger events—at least not in an automated sense. Here, we summarize results for these inversions and discuss levels of completeness for this analysis. Detailed solutions obtained for $6.0 \le M_w < 6.5$ events together with GCMT solutions for comparison are presented from Fig. S13 to S36.

Fig. 20 shows the variation of $\Delta M_w$ and $\Phi$ as a function of time for $6.0 \le M_w < 6.5$ events. Of the 1727 earthquakes from the GCMT catalogue falling within this magnitude range, 1665 can be considered undisturbed by the criteria outlined in Sections 3.2 and 4.5). Of these, 1511 events produce a WCMT solution whose magnitude is within $\pm 0.2$ units of the GCMT solution—approximately

<figure><img src="image_31.png" /><figcaption>**Figure 19.** W phase source inversion results for the 2010 great Maule-Chile earthquake. W phase CMT (WCMT) and Global CMT (GCMT) solutions are shown in the topleft corner. Examples of observed waveforms (black lines) and the corresponding synthetics (red lines) computed from the WCMT solution are presented. The station azimuth ($\phi$) and epicentral distance ($\Delta$) are indicated as well as W phase time windows, bounded by red dots. The WCMT inversion is based on the ground motions recorded at stations within $\Delta \le 90°$ after applying a bandpass filter in the 1–5 mHz passband. W phase and later arrivals are very well predicted by the WCMT solution. For some channels like EFI-LHZ or CRZF-LHZ, the surface waves are affected by instrument problems, although the W phase signal itself is not affected.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>**Figure 20.** Difference between WCMT and GCMT solutions for $6.0 \le M_w < 6.5$ events. The magnitude difference $\Delta M_w$ is presented in (a) and the angle $\Phi(M_{wcmt}, M_{gcmt})$ is shown in (b). In (a) and (b), the W phase solutions corresponding to OL2 (i.e. for $N < 30$ or $\gamma > 270^\circ$) are indicated by triangles and WCMT solutions corresponding to OL3 (i.e. $N \ge 30$ and $\gamma \le 270^\circ$) are indicated by circles. The symbols are coloured according to the GCMT moment magnitude $M_{w-gcmt}$.</figcaption></figure>

91 per cent of events. Of the remaining 9 per cent, one-third (52 earthquakes) did not converge on a W phase solution because no station passed the inversion signal-to-noise screening filters. At $\Delta = 50^\circ$, these statistics change only slightly; 90 per cent of undisturbed events obtain W phase magnitudes within $\pm 0.2$ units of GCMT, and solutions were absent for just less than 4 per cent of all earthquakes.

As expected, at these smaller magnitudes inverted mechanisms show larger variability between W phase and GCMT solutions than do results for larger events. 65 per cent of earthquakes with $M_w < 6.5$ demonstrate $\Phi < 20^\circ$, while as many as 12 per cent show $\Phi > 50^\circ$ (versus 85 per cent and less than 1 per cent for the $M_w \ge 6.5$ data set). 89 per cent of solutions lie within 50 km of

(a) Standard W phase algorithm

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

(b) Modified W phase algorithm for disturbed events

<figure><img src="image_34.png" /><figcaption>**Figure 21.** WCMT solutions obtained for the 2009 Vanuatu earthquake sequence. The results shown in (a) are computed using the standard W phase algorithm in which we assume isolated sources in time and space. The WCMT solutions shown in (b) are obtained using the modified W phase algorithm for ‘disturbed events’ in which the perturbations caused by the disturbing events are removed from the data inverted for later solutions. The WCMT solutions are displayed in red for (1) 200910072203A, in blue for (2) 200910072218A, in dark green for (3) 200910072250A and in black for (4) 200910072313A. The Global CMT solutions are shown using small light green beach balls for comparison.</figcaption></figure>
GCMT centroid locations—very similar to the results for events with $M_w \ge 6.5$. Interestingly, average depth differences for these smaller events are just 6.5 km when compared to GCMT solutions, and 90 per cent of the solutions obtain depth estimates within 12 km of GCMT. These results suggest closer alignment with the GCMT results than for the $M_w \ge 6.5$ data set.

## 4.5 Disturbed events

As noted in Section 3.2, events that occur soon after another large earthquake are problematic. We define ‘disturbed events’ as earthquakes occurring within 1 hr of $M_w \ge 6.5$ earthquakes or less than 10 hr after $M_w \ge 7.0$ events. These events have poor station azimuthal coverage after performing the W phase data screening for $\Delta \le 50^\circ$ (i.e. $N < 30^\circ$ or $\gamma > 270^\circ$). While the magnitude of a ‘disturbed event’ is usually smaller than a preceding event (since most of them are aftershocks), this is not necessarily a rule, because there is no consideration of the size of earthquakes in the definition of ‘disturbed events’. This makes such a definition particularly adaptable for real-time operations of W phase inversions. To retrieve the WCMT solution of such events in real time, Hayes et al. (2009) proposed to modify the time window and bandpass filter to perform a CMT inversion based on surface wave data. We explore here an alternative approach in which we compute the synthetics for the disturbing (preceding) event and subtract them from the data to produce the residual trace. Then, we run the WCMT algorithm on the residual trace to obtain the source parameters of the disturbed event, as we do for a normal earthquake.

To illustrate this technique, we show a remarkable example for the 2009 October 7 Vanuatu earthquake sequence, in which four $M_w \ge 6.8$ events occurred over a period of less than 2 hr. The four earthquakes of this sequence are 200910072203A: $M_w = 7.6$; 200910072218A: $M_w = 7.8$; 200910072250A: $M_w = 6.8$ and 200910072313A: $M_w = 7.4$ (cf. Table S1). For brevity, we will refer to these earthquakes as (1), (2), (3) and (4), respectively. To estimate the source parameters of the ‘disturbed events’ (2), (3) and (4), we use the following procedure: We consider that the WCMT solution of event (1) obtained for $\Delta \le 50^\circ$ is well constrained and is not perturbed at all by the waveforms associated with the later events. We calculate the full synthetics corresponding to this solution and subtract them from the ground motions at each station. Using the corresponding residuals as the input data, we perform a WCMT inversion for (2). We then subtract the synthetics corresponding to the superposition of (1) and (2) from the data and perform another inversion to obtain the source parameters of (3). We finally use the same procedure to estimate the CMT of (4) by removing the synthetics corresponding to the combination of (1), (2) and (3).

The solutions shown in Fig. 21(a) are obtained using the standard W phase algorithm without removing the perturbations as described earlier. They show significant differences from the focal mechanisms and magnitudes of GCMT, particularly for events (2) and (3). The centroid solutions resulting from the modified W phase algorithm for ‘disturbed events’ are displayed in Fig. 22. After WCMT inversion, the four mechanisms are well aligned along the North New Hebrides Trench and the data are very well fit, as shown in Fig. 23. In Fig. 21(b), the new WCMT solutions are compared to GCMT. The W phase moment magnitudes match GCMT very well. The focal mechanisms are also in good agreement, although the solution obtained for (3), the smallest earthquake of the 2009 Vanuatu sequence, is probably affected by the misfit resulting from inversions of events (1) and (2).

<figure><img src="image_35.png" /><figcaption>**Figure 22.** WCMT solutions obtained for the 2009 October Vanuatu earthquake sequence using the modified W phase algorithm for ‘disturbed events’. The events that make up this sequence are (1) 200910072203A which is shown in red, (2) 200910072218A in blue, (3) 200910072250A in green and (4) 200910072313A in black. (a) shows the normalized moment rate functions and (b) shows the focal mechanisms with the corresponding centroid locations. The station distribution is indicated on the globe at the right-hand corner of (b). To determine the mechanisms of the ‘disturbed events’ (2), (3) and (4), the WCMT algorithm is modified by correcting for the disturbance of preceding event(s).</figcaption></figure>

The simple scheme detailed here can be easily used to retrieve in real time the moment tensor as well as the centroid position and timing of ‘disturbed events’. However, this procedure does not work well if the size of the ‘disturbed event’ is small compared to the preceding earthquake. In such cases, the W phase amplitude of the ‘disturbed event’ can be small with respect to the misfit associated with the preceding solution. Difficulty also arises when the origin time difference between the two events is extremely small (e.g. less than 12 min). In such cases, a simultaneous inversion of the two point source parameters could be a better option.

## 5 CONCLUSION

The W phase source inversion provides a robust estimation of CMT parameters for moderate to large earthquakes (i.e. $M_w \ge 6.5$). An efficient data screening is used to automatically reject bad traces from the data set. The WCMT algorithm provides three levels of output: (1) a simple and robust preliminary W phase estimation of $M_w$, (2) a first moment tensor calculation assuming a PDE location as the centroid position and (3) a full CMT inversion that includes the centroid position determination (latitude, longitude and depth).

By applying the inversion for all $M_w \ge 6.5$ events between 1990 and 2010, this work highlights the robustness of the WCMT inversion. It is important to point out that all of these results were
<figure><img src="image_36.png" /><figcaption>**Figure 23.** Comparisons between synthetics and data for the 2009 October Vanuatu earthquake sequence. Examples of observed waveforms (black lines) and the corresponding synthetics (red lines) computed from the WCMT solutions of events (1), (2), (3) and (4) (cf. Fig. 22) are presented after bandpass filtering in the 2–8.3 mHz passband. The station azimuth ($\phi$) and epicentral distance ($\Delta$) are indicated relative to the centroid of (4). The W phase time windows are bounded by red dots for event (1), blue dots for (2), green dots for (3) and black dots for (4). W phase and later arrivals associated with the four earthquakes are very well predicted by the WCMT solutions.</figcaption></figure>

obtained without any human interaction, such that only fully automated solutions are presented. The WCMT algorithm can thus operate online in a completely automated manner. The first WCMT solution, obtained 22 min after the earthquake origin time, is generally sufficiently well constrained for tsunami warning purposes. The second inversion, performed 13 min later, provides very robust estimates, which can be used for a variety of purposes including tsunami modelling, finite fault inversions or damage assessments. Completeness of the WCMT catalogue for recent years implies that wherever an $M_w > 6.5$ event occurs today, the present status of the global network is such that the W phase inversion can provide a reliable estimation of $M_w$ and the focal mechanism within 35 min of the earthquake origin time.

We also show that, while not complete to $M_w = 6.0$, the W phase inversion obtains reliable source parameters for the majority of events of this magnitude or larger. In recent years, we

could automatically determine the source parameters of almost all undisturbed events.

In this study, we have defined ‘disturbed events’ as any earthquake contaminated by the large amplitude disturbance caused by a preceding event. The standard W phase algorithm is not suited for such events because the assumption of an isolated source in time and space is no longer valid. In this work, we discuss a possible solution to handle them in real time. However, it is still difficult to establish a general scheme to estimate the source parameters of ‘disturbed events’ quickly and we must further explore methods to cope with such difficult situations.

The recent 2011 Tohoku-oki earthquake ($M_w = 9.0$) is not included in the present catalogue because it occurred after the 1990–2010 time period considered in this study. A detailed report of real-time results and further analyses of this event are presented in Duputel *et al.* (2011).
ACKNOWLEDGMENTS

This work uses Federation of Digital Seismic Networks (FDSN) seismic data and CMT solutions from the Global CMT catalogue. 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, J. Braunmiller, D. Garcia and an anonymous reviewer for their helpful comments on the manuscript.

REFERENCES

Abe, K., 1979. Size of great earthquakes of 1837–1974 inferred from tsunami data, *J. geophys. Res.*, **84**, 1561–1568.

Abe, K., 1989. Quantification of tsunamigenic earthquakes by the Mt scale, *Tectonophysics*, **166**(1–3), 27–34.

Abercrombie, R.E., 2003. The June 2000 $M_w$7.9 earthquakes south of Sumatra: deformation in the India–Australia Plate, *J. geophys. Res.*, **108**, 2018, doi:10.1029/2001JB000674.

Ammon, C.J. et al., 2005. Rupture process of the 2004 Sumatra-Andaman Earthquake, *Science*, **308**, 1133–1139.

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

Duputel, Z., Rivera, L. & Kanamori, H., 2010. W-phase: lessons from the February 27, 2010 Chilean earthquake, *Seismol. Res. Lett.*, **81**, 544.

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.

Dziewonski, A., 1982. Harvard Centroid Moment Tensor Project, Available at: http://www.globalcmt.org.

Dziewonski, A. & Woodhouse, J.H., 1983. Studies of the seismic source using normal-mode theory, in *Earthquakes: Observation Theory and Interpretation: Notes from the International School of Physics “Enrico Fermi”* (1982: Varenna, Italy), pp. 45–137, North-Holland Publ. Co., Amsterdam.

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.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth model, *Phys. Earth planet. Inter.*, **25**(4), 297–356.

Ekström, G. & Nettles, M., 2006. Global Centroid Moment Tensor Project, Available at: http://www.globalcmt.org.

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.

ERI, 2000. EIC Seismological Note 94, Technical report.

Geist, E., 2005. Triggering of tsunamigenic aftershocks from large strike-slip earthquakes: analysis of the November 2000 New Ireland earthquake sequence, *Geochem. Geophys. Geosyst*, **6**, doi:10.1029/2005GC000935.

Hanks, T.C. & Kanamori, H., 1979. A moment magnitude scale, *J. geophys. Res.*, **84**, 2348–2350.

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.

Hayes, G., Earle, P., Benz, H., Wald, D. & Briggs, R., 2011. 88 hours: The US Geological Survey National Earthquake Information Center response to the March 11, 2011 $M_w$ 9.0 Tohoku earthquake, *Seismol. Res. Lett.*, **82**, 481–493.

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

Kagan, Y.Y., 1991. 3-D rotation of double-couple earthquake sources, *Geo-phys. J. Int.*, **106**(3), 709–716.

Kanamori, H., 1972. Mechanism of tsunami earthquakes, *Phys. Earth planet. Inter.*, **6**, 356–359.

Kanamori, H., 1977. The energy release in great earthquakes, *J. geophys. Res.*, **82**(20), 2981–2987.

Kanamori, H., 1993. W phase, *Geophys. Res. Lett.*, **20**(16), 1691–1694.

Kanamori, H., 2004. The diversity of the physics of earthquakes, *Proc. Japan Acad.*, **80**, 297–316.

Kanamori, H. & Anderson, D.L., 1975. Theoretical basis of some empirical relations in seismology, *Bull. seism. Soc. Am.*, **65**, 1073–1095.

Kanamori, H. & Given, J.W., 1981. Use of long-period surface waves for rapid determination of earthquake-source parameters, *Phys. Earth planet. Inter.*, **27**, 8–31.

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

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

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.

Mendiguren, J.A., 1977. Inversion of surface wave data in source mechanism studies, *J. geophys. Res.*, **82**(5), 889–894.

Okal, E.A., 1988. Seismic parameters controlling far-field tsunami amplitudes: a review, *Nat. Hazards*, **1**(1), 67–96.

Park, J. et al., 2005. Earth’s free oscillations excited by the 26 December 2004 Sumatra-Andaman earthquake, *Science*, **308**(5725), 1139–1144.

Park, S.-C. & Mori, J., 2007. Triggering of earthquakes during the 2000 Papua New Guinea earthquake sequence, *J. geophys. Res.*, **112**(B3), B03302, doi:10.1029/2006JB004480.

Peterson, J., 1993. Observations and modelling of background seismic noise, Open-file Report 93-322, U.S. Geological Survey.

PTWC, 2006a. Tsunami Bulletin Number 001, Issued at 0836Z 17 Jul 2006, Technical report.

PTWC, 2006b. Tsunami Bulletin Number 002, Issued at 1108Z 17 July 2006, Technical report.

Richter, C.F., 1935. An instrumental earthquake magnitude scale, *Bull. seism. Soc. Am.*, **25**(1), 1–32.

Rivera, L. & Kanamori, H., 2009. Using W phase for regional tsunami warning and rapid earthquake hazard assessment, in *Proceedings of the 2nd International Workshop on Earthquake Early Warning*, Kyoto, Japan.

Rivera, L., Kanamori, H. & Duputel, Z., 2010. The potential of the W-phase algorithm for regional tsunami warning in Chile, in *Proceedings of AGU Chapman Conference in Chile*.

Silver, P.G. & Jordan, T.H., 1982. Optimal estimation of scalar seismic moment, *Geophys. J. R. astr. Soc.*, **70**, 755–787.

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.

USGS, 2011. *Report on the 2010 Haiti Earthquake*, available at http://earthquake.usgs.gov. (last accessed 2011 January).

APPENDIX: ANGULAR DISTANCE BETWEEN TWO MOMENT TENSORS

To compare the orientations of the moment tensors $\mathbf{M}_A$ and $\mathbf{M}_B$, we calculate the minimum angle necessary to bring the principal axes of $\mathbf{M}_B$ into coincidence with those of $\mathbf{M}_A$. This can be done as follows. Let $\mathbf{P}(\mathbf{M})$ be the proper orthogonal matrix whose columns correspond to the principal axes of the moment tensor $\mathbf{M}$. Because of the ambiguity in the principal axes directions there are four such possible matrices $\mathbf{P}_1(\mathbf{M})$, $\mathbf{P}_2(\mathbf{M})$, $\mathbf{P}_3(\mathbf{M})$ and $\mathbf{P}_4(\mathbf{M})$. We can fix $\mathbf{P}_1(\mathbf{M}_A)$, compute the four angles

$$ \phi_k = \cos^{-1} \frac{\text{trace}[\mathbf{P}_1(\mathbf{M}_A) \mathbf{P}_k^\dagger(\mathbf{M}_B)] - 1}{2}, \quad (A1) $$

and define $\Phi(\mathbf{M}_A, \mathbf{M}_B) = \min[\phi_1, \phi_2, \phi_3, \phi_4]$ as the angular distance between $\mathbf{M}_A$ and $\mathbf{M}_B$. This value ranges from 0° to 120°. This definition is equivalent to the one given by Kagan (1991).
## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article:

**Supplement.** The Table S1 gives the complete list of *disturbed events* considered in this study for $M_w \ge 6.5$ (see Section 4.5 in the main text). The complete collection of solutions obtained for $M_w \ge 6.5$ between 1990 and 2010 after rejecting the *disturbed events* are shown in order of decreasing GCMT magnitude from Fig. S1 to Fig. S11. Each frame represents an earthquake for which we show the GCMT solution in green, and the W phase solutions incorporating stations within $\Delta \le 50^\circ t_b$ and within $\Delta \le 90^\circ$ at $t_c$. WCMT mechanisms corresponding to OL2 (i.e. if $N < 30$ or $\gamma > 270^\circ$) are presented in blue while those corresponding to OL3 (i.e. if $N \ge 30^\circ$ and $\gamma \le 270^\circ$) are shown in red.

$N$ is the total number of channels after inversion and $\gamma$ is the azimuthal gap (in degrees). WCMT solutions for $M_w \ge 6.5$ earthquakes using data within $\Delta < 90^\circ$ are also available on the website http://eost.u-strasbg.fr/wphase/MGE65. The great-circle distances between WCMT, GCMT and PDE locations for $M_w \ge 6.5$ earthquakes are presented in Fig. S12 and the complete list of WCMT centroid locations are listed from Table S2 to Table S15. The solutions resulting from the extension to $6.0 \le M_w < 6.5$ events (see Section 4.4 in the main text) are ordered in decreasing GCMT magnitude $M_{w-gcmt}$ from Fig. S13 to Fig. S36.

Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.