# Waveform inversion for localized three-dimensional seismic velocity structure in the lowermost mantle beneath the Western Pacific

Kensuke Konishi, Kenji Kawai, Robert J. Geller, Nobuaki Fuji

► To cite this version:

Kensuke Konishi, Kenji Kawai, Robert J. Geller, Nobuaki Fuji. Waveform inversion for localized three-dimensional seismic velocity structure in the lowermost mantle beneath the Western Pacific. Geophysical Journal International, 2014, 199, pp.1245-1267. 10.1093/gji/ggu288 . insu-03581063

HAL Id: insu-03581063

https://insu.hal.science/insu-03581063

Submitted on 19 Feb 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
<figure><img src="image_2.png" /></figure>

# Waveform inversion for localized three-dimensional seismic velocity structure in the lowermost mantle beneath the Western Pacific

Kensuke Konishi,<fnref n="1" /> Kenji Kawai,<fnref n="2" />,<fnref n="3" />,* Robert J. Geller<fnref n="4" /> and Nobuaki Fuji<fnref n="5" />

Accepted 2014 July 21. Received 2014 July 16; in original form 2014 January 1

## SUMMARY

We infer 3-D localized shear velocity structure in the lowermost 400 km of the mantle at the western edge of the Pacific large low shear velocity province (LLSVP) by applying waveform inversion to transverse component body-wave waveforms from the F-net seismic array in Japan. Our data set consists of relatively long period (12.5–200 s) broad-band seismic waveforms of Tonga-Fiji deep focus and intermediate deep earthquakes. We conduct several tests to confirm the robustness of the inversion results. We find two low-velocity zones at the bottom of the target region, with a high-velocity zone in the middle, and a low-velocity zone above the high-velocity zone and contiguous with the two deeper low-velocity zones at a depth of 200–300 km above the core-mantle boundary (CMB). This supports the idea that the Pacific LLSVP may be an aggregation of small upwelling plumes rather than a single large thermochemical pile.

**Key words:** Mantle processes; Composition of the mantle; Body waves; Seismic tomography.

## 1 INTRODUCTION

Global-scale seismic velocity structure in the lowermost mantle, which is both the thermal boundary layer for mantle convection and the chemical boundary between the silicate mantle and the liquid iron alloy outer core (e.g. Poirier 2000), has been investigated mainly using shear waves rather than compressional waves, due to their more nearly homogeneous sampling (see review by Romanowicz 2003). Such studies have found large-scale heterogeneity in the lowermost mantle, called the degree-2 pattern, where two large low shear velocity provinces (LLSVPs) beneath the Pacific and Africa are surrounded by high-velocity regions (e.g. Grand 2002). It has been controversial whether this is due to thermal effects, chemical effects or a combination of both.

Modelling of the mantle convection pattern in the lowermost mantle under isochemical conditions tends to produce smaller thermal plumes which cluster beneath the Pacific and Africa (e.g. Bunge et al. 1998; Schubert et al. 2004). On the other hand, thermochemical effects such as depletion or enrichment of components in incompatible elements could produce a large chemically distinct pile in the lowermost mantle, which would be the root of a superplume

(see Bull et al. 2009). The former and latter are called ‘plume clusters’ and ‘thermochemical piles’, respectively. Data which could elucidate the mode of mantle convection can thus contribute to understanding whether the driving forces of mantle convection are thermal, compositional or a combination of both.

The lowermost mantle beneath the western Pacific is located at the western edge of the Pacific LLSVP (Takeuchi 2007). There have been many studies of the velocity structure beneath the western Pacific, which is well sampled due to the favourable distribution of seismic sources and receivers. A forward modelling study by Garnero et al. (1993) found that no single laterally homogeneous shear velocity model beneath the western Pacific could satisfactorily match all of the observations, suggesting the existence of laterally heterogeneous structure. Wysession et al. (1994) analysed *ScS-S* and *sScS-sS* differential traveltimes using long-period data and obtained laterally heterogenous shear wave velocity structure in the lowermost 300 km of the mantle beneath the western Pacific, which suggested a 3 per cent low-velocity region beneath northeastern Indonesia surrounded by high-velocity regions. Kito et al. (2004) analysed short-period *S* waves recorded by the J-Array and Hi-net array using double array methods and found two *S* velocity jumps (one negative and one positive at depths of 340 and 40 km above the CMB, respectively). They inferred 1-D velocity profiles based on the amplitude ratio between the reference and reflected phases and differential traveltimes.

<footnote n="1">¹School of Earth and Environmental Sciences, Seoul National University, Gwanak, Seoul 151-742, Republic of Korea</footnote>

<footnote n="2">²Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Ookayama 2-12-1, Meguro-ku, Tokyo 152-8551, Japan</footnote>

<footnote n="3">³Earth-Life Science Institute, Tokyo Institute of Technology, Ookayama 2-12-1, Meguro-ku, Tokyo 152-8551, Japan</footnote>

<footnote n="4">⁴Department of Earth and Planetary Science, Graduate School of Science, University of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo 113-0033 Japan.</footnote>

<footnote>E-mail: bob@eps.s.u-tokyo.ac.jp</footnote>

<footnote n="5">⁵Institut de Physique du Globe de Paris, 1 rue Jussieu F-75238 Paris Cedex 05, France</footnote>

<footnote>* Now at: Department of Earth Science and Astronomy, Graduate School of Arts and Sciences, University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Event-receiver geometry, with great circle ray paths. The portions of the great circles which sample D'' are shown in red; plus signs indicate the turning points within D''. Blue triangles show the stations, and red stars show the earthquakes. (b) Event-receiver geometry. Red stars show the locations of earthquakes. Blue triangles show the stations. The voxels for the inversion are shown by red stars with black circles. The yellow star shows the location for the partial derivatives of shallow structure used in Fig. 19. (c) Cross section schematically showing the distribution of events (red star), stations (blue triangle), target region (orange region) and shallow perturbation points used in Fig. 19 (yellow star and underlying points). (d) The green lines show the locations of the six cross sections in Figs 4, 5 and 13. The red dots indicate the locations of voxels for the inversion.</figcaption></figure>

He *et al.* (2006) and Takeuchi *et al.* (2008) suggested a sharp-sided boundary of the Pacific LLSVP based on an analysis of S-ScS and sS-sScS differential traveltimes. He *et al.* (2006) conducted forward modelling and reported negative shear velocity anomalies ranging from 0 to -1 to -13 per cent, suggesting the existence of an ultra-low velocity zone (ULVZ) in the lowermost 220 km of the mantle. He & Wen (2009) reported velocity anomalies varying from -3.0 per cent at the top (at 740 km above the CMB) to -3.5 per cent at 100 km above the CMB, with an average shear velocity reduction of -5 per cent in the bottom 100 km of the mantle in the NW-SE direction. Idehara (2011) suggested that there is a ULVZ beneath only parts of the western Pacific. Takeuchi & Obara (2010) studied the topography of the D'' discontinuity by stacking and migration of shear waveform data recorded by the Hi-net tiltmeter array, suggesting laterally heterogeneous fine structure with a length of about 700 km in the NE-SW direction. Takeuchi (2012) conducted waveform inversion for laterally heterogeneous structure relative to an average shear wave velocity structure in the depth range from the CMB to 400 km above the CMB and found a ridge-like structure at the edge of the Pacific LLSVP.

<table><caption>**Table 1.** Earthquakes used in this study.</caption><thead><tr><th>Event #</th><th>Date (Y/M/D)</th><th>Latitude</th><th>Longitude</th><th>Depth</th><th>M<sub>w</sub></th></tr></thead><tbody><tr><td>1</td><td>2000 May 4</td><td>-17.72</td><td>-178.31</td><td>539.8</td><td>6.4</td></tr><tr><td>2</td><td>2000 June 14</td><td>-25.45</td><td>178.38</td><td>615.4</td><td>6.4</td></tr><tr><td>3</td><td>2000 December 18</td><td>-21.11</td><td>-178.98</td><td>655.7</td><td>6.5</td></tr><tr><td>4</td><td>2001 April 28</td><td>-18.07</td><td>-176.68</td><td>367.4</td><td>6.8</td></tr><tr><td>5</td><td>2001 May 26</td><td>-20.25</td><td>-177.65</td><td>413.9</td><td>6.3</td></tr><tr><td>6</td><td>2001 September 12</td><td>-20.84</td><td>-178.90</td><td>634.1</td><td>6.4</td></tr><tr><td>7</td><td>2001 November 5</td><td>-17.12</td><td>-178.96</td><td>579.7</td><td>6.2</td></tr><tr><td>8</td><td>2002 January 2</td><td>-17.63</td><td>178.84</td><td>680.8</td><td>6.1</td></tr><tr><td>9</td><td>2002 June 30</td><td>-22.13</td><td>179.43</td><td>631.6</td><td>6.4</td></tr><tr><td>10</td><td>2002 August 9</td><td>-16.25</td><td>-175.85</td><td>381.3</td><td>6.1</td></tr><tr><td>11</td><td>2002 October 4</td><td>-20.86</td><td>-178.74</td><td>650.8</td><td>6.3</td></tr><tr><td>12</td><td>2002 October 17</td><td>-19.80</td><td>-178.23</td><td>621.9</td><td>6.1</td></tr><tr><td>13</td><td>2002 December 10</td><td>-24.02</td><td>179.28</td><td>538.8</td><td>6.0</td></tr><tr><td>14</td><td>2003 January 4</td><td>-20.72</td><td>-177.32</td><td>394.7</td><td>6.5</td></tr><tr><td>15</td><td>2003 March 14</td><td>-17.50</td><td>-174.77</td><td>282.4</td><td>6.3</td></tr><tr><td>16</td><td>2003 July 27</td><td>-21.09</td><td>-176.12</td><td>215.6</td><td>6.6</td></tr><tr><td>17</td><td>2004 April 9</td><td>-13.22</td><td>167.16</td><td>226.3</td><td>6.4</td></tr><tr><td>18</td><td>2005 March 19</td><td>-21.88</td><td>-179.27</td><td>609.2</td><td>6.3</td></tr><tr><td>19</td><td>2005 August 6</td><td>-19.60</td><td>-175.35</td><td>217.7</td><td>6.0</td></tr><tr><td>20</td><td>2006 January 2</td><td>-19.80</td><td>-177.72</td><td>589.5</td><td>7.1</td></tr><tr><td>21</td><td>2006 February 2</td><td>-17.77</td><td>-178.13</td><td>611.6</td><td>6.7</td></tr><tr><td>22</td><td>2006 February 24</td><td>-17.94</td><td>-179.42</td><td>640.9</td><td>6.1</td></tr><tr><td>23</td><td>2006 February 26</td><td>-23.59</td><td>-179.82</td><td>553.9</td><td>6.4</td></tr><tr><td>24</td><td>2006 June 2</td><td>-20.77</td><td>-178.54</td><td>584.6</td><td>6.0</td></tr><tr><td>25</td><td>2006 June 9</td><td>-17.36</td><td>-178.62</td><td>585.9</td><td>6.1</td></tr><tr><td>26</td><td>2007 May 6</td><td>-19.44</td><td>-179.04</td><td>690.8</td><td>6.5</td></tr><tr><td>27</td><td>2007 October 5</td><td>-25.27</td><td>179.50</td><td>540.8</td><td>6.5</td></tr><tr><td>28</td><td>2007 November 19</td><td>-21.05</td><td>-178.63</td><td>562.5</td><td>6.3</td></tr><tr><td>29</td><td>2008 January 15</td><td>-22.05</td><td>-179.34</td><td>603.3</td><td>6.5</td></tr><tr><td>30</td><td>2008 April 18</td><td>-17.26</td><td>-178.98</td><td>577.8</td><td>6.3</td></tr><tr><td>31</td><td>2008 July 3</td><td>-23.41</td><td>-179.69</td><td>589.9</td><td>6.2</td></tr><tr><td>32</td><td>2009 November 22</td><td>-17.72</td><td>-178.36</td><td>546.4</td><td>5.7</td></tr><tr><td>33</td><td>2010 June 30</td><td>-23.19</td><td>179.26</td><td>581.6</td><td>6.4</td></tr><tr><td>34</td><td>2010 August 16</td><td>-20.74</td><td>-178.67</td><td>604.0</td><td>6.2</td></tr><tr><td>35</td><td>2010 December 28</td><td>-23.48</td><td>-179.72</td><td>571.3</td><td>6.3</td></tr></tbody></table>

Our group conducted waveform inversion for 1-D seismic structure in the lowermost mantle beneath the western Pacific and found an 'S-shaped' pattern of the vertical dependence of the shear velocity model (Konishi *et al.* 2009). We recently extended our methods to allow waveform inversion for 3-D seismic structure, and applied these methods to perform waveform inversion for localized 3-D seismic structure in the lowermost mantle beneath central America using a data set from the US-Array (Kawai *et al.* 2014). Kawai *et al.* (2014) present a full explanation of our methods. In this study, we apply these methods to invert for the 3-D S-velocity structure in the lowermost mantle at the western edge of the Pacific LLSVP.

## 2 DATA AND METHODS

In this section we discuss the data and methods used in this study.

### 2.1 Data

The lowermost mantle beneath the western Pacific (Fig. 1) is the target region for the 3-D inversion in this study. The combined availability of sources—intermediate and deep focus earthquakes—and receivers determines which areas can be studied. Since seismicity is unlikely to change significantly, further expansion of the studyable areas thus depends on the installation of new seismic arrays.
We use the transverse components of broad-band waveform data (obtained by rotating the N–S and E–W components) for 35 events (Table 1, Fig. 1) from the Japanese F-net (obtained from the data centre of the National Research Institute for Earth Science and Disaster Prevention of Japan). The data are observed at epicentral distances of 54° < Δ < 76°. We deconvolve the instrument response, apply a bandpass filter to the data and construct data sets for the passband 0.005 to 0.08 Hz (i.e. for the period range, 12.5–200 s). Note that the same filter is also applied to the synthetics. We gather a total of 2240 records which include S, ScS and the other phases which arrive between them. Two selection criteria for data quality are then used. First, we compute the ratio of the maximum amplitude of the data and the corresponding synthetic, and eliminate records for which the ratio is greater than 2 or less than 0.5. Secondly, records are rejected if the correlation coefficient between the data and synthetics is less than 0.5. We use the above criteria to select 1597 records (the ray paths shown in Fig. 1). We use the ‘auto-pick’ method to make static corrections for the effect of structure outside the target region (see Fuji et al. 2010; Kawai & Geller 2010; Konishi et al. 2012; Kawai et al. 2014). We use the reciprocal of the maximum amplitude in each trace as the weighting factor for each trace.

## 2.2 Synthetics and partials

We compute synthetics for each source at each receiver using the direct solution method (DSM, e.g. Geller & Ohminato 1994; Kawai et al. 2006) for the anisotropic PREM model (Dziewonski & Anderson 1981). The horizontal dependence of the wavefield is given by vector spherical harmonics and separation of variables is then

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (Top) AIC values for each model. Red: inversion for ‘NEW 1-D’, Green: 3-D inversion with respect to PREM and Blue: 3-D inversion with respect to ‘NEW 1-D’. (Bottom) Same as top panel, but showing results for only the first 10 CG vectors. The horizontal axis (k) shows the number of CG basis vectors used in each inversion.</figcaption></figure>

<table><caption>**Table 2.** Variance and AIC for each model.</caption><thead><tr><th>Model</th><th>Variance (per cent)</th><th>AIC</th></tr></thead><tbody><tr><td>PREM</td><td>82.4</td><td>—</td></tr><tr><td>PREM with time shift</td><td>72.8</td><td>444</td></tr><tr><td>‘NEW 1-D’ with time shift</td><td>63.6</td><td>434</td></tr><tr><td>CG3 (w.r.t. PREM)</td><td>60.6</td><td>418</td></tr><tr><td>CG3 (w.r.t. ‘NEW 1-D’)</td><td>59.9</td><td>412</td></tr></tbody></table>

used to obtain the following (see Kawai et al. 2006, for details). The DSM obtains the solution of the weak form of the equation of motion by directly solving the Galerkin weak form of the equation of motion:

$$ (\omega^2 \mathbf{T} - \mathbf{H} + \omega \mathbf{R})\mathbf{c} = -\mathbf{g}, \quad (1) $$

where c is the vector of coefficient of the trial functions, T is the mass matrix, H is the stiffness matrix, g is the force vector and R enforces continuity conditions at fluid-solid boundaries, which exist only for the spheroidal case.

When our group first started performing inversion of surface wave waveform data for 3-D Earth structure we used modal superposition to compute the synthetics and computed the partial derivatives of each eigenfrequency and eigenfunction to compute the partial derivatives of the synthetic seismograms (Hara et al. 1991). We then realized that this time-consuming step could be greatly streamlined by directly using the Born approximation to compute the partials (Hara et al. 1993); we called this method the DSM. We began applying this approach to body waves shortly thereafter. Note that the matrix on the left-hand side of eq. (1) is only factorized once for each frequency (see Kawai et al. 2006).

Inversion for the 3-D structure in the target region—limited to only S velocity in this study—is conducted by an inversion using the first order Born approximation as formulated by Geller & Hara (1993). We express the *i*-component of the displacement for the *k*th earthquake at the point *r* as $u_i^{(k)}(r)$, and the *i*-component of the back-propagated displacement at *r* excited by a point force in

<figure><img src="image_5.png" /><figcaption>**Figure 3.** ‘NEW 1-D’ (red) and PREM (black).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Results of the inversion with respect to PREM for the 3-D shear wave velocity structure in the lowermost mantle beneath the western Pacific. The reference model for the velocity perturbation is PREM, and the colour bar has units of percentage. (a) Map view for the depth range (relative to the CMB) shown at the top of each panel. (b) Cross sections. The vertical axes show height in km above the CMB.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** 3-D model inferred using 'NEW 1-D' as the initial model. The reference model for the figure is PREM. Other details same as Fig. 4.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** (Top) Coefficients of correlation between CG basis vectors of inversions with respect to PREM and ‘NEW 1-D’. (Bottom) Same as top panel, but showing results for only the first five CG vectors.</figcaption></figure>

the j-direction at the pth station as $\eta_{ij}^{(p)}(\mathbf{r})$. Note that $\mathbf{r}$ will in general be a point in the Earth’s interior (the ‘target’ of the inversion).

The partial derivatives in the wavefield formulation are expressed as follows:

$$ \left\{ \frac{\partial u_i^{(k)}[\mathbf{r}^{(p)}]}{\partial m_l} \right\}^* = \int_V (\omega^*)^2 [u_j^{(k)}]^* [\rho^{(l)}]^* \eta_{ji}^{(p)} dV - \int_V [u_{j,q}^{(k)}]^* [C_{jqr s}^{(l)}]^* \eta_{ri,s}^{(p)} dV, \quad (2) $$

where $\rho^{(l)} = \delta\rho\delta(r - r_0)$ or $C_{ijpq}^{(l)} = \delta C_{ijpq} \delta(r - r_0)$ and $u_{i,j} = \frac{\partial u_i}{\partial x_j}$ is a locally Cartesian derivative (see Geller & Hara 1993, for details). Since the transverse component of waveforms is used in this study, we take into account toroidal–toroidal coupling and neglect other coupling.

### 2.2.1 3-D synthetics based on the Born approximation

We compute 3-D synthetics for voxel perturbations including anisotropy and anelasticity using the first-order Born approximation as follows:

$$ u_{3D}(r) = u_{1D}(r) + \sum_{l} \frac{\partial u(r)}{\partial m_l} \delta m_l, \quad (3) $$

where $u_{3D}$ is a synthetic seismogram for the 3-D Earth model and $u_{1D}$ is a synthetic for a 1-D Earth model, such as PREM (Dziewonski & Anderson 1981).

Partial derivatives with respect to 3-D structure are computed using the formulation of Geller & Hara (1993). Details of our inversion methods are given by Kawai et al. (2014). There is some error due to the use of the Born approximation in both the individual synthetic seismograms for the perturbed model and in the 3-D model obtained by the inversion. Quantifying the extent of this error is an important topic for future research.

The source parameters (moment tensors, centroids and half durations) are fixed to the Global CMT solution. We convolve boxcar moment rate functions with the appropriate half durations with the synthetic seismograms and their partial derivatives.

The source parameters (moment tensors, centroids and half durations) are fixed to the Global CMT solution. We convolve boxcar moment rate functions with the appropriate half durations with the synthetic seismograms and their partial derivatives.

## 2.3 Inverse problem formulation

We formulate the waveform inversion problem in the time domain as follows. We define $\delta d$ to be the residual (the difference of the observed seismogram, $d_{\text{OBS}}$, and the synthetic seismogram for the initial model, $d_{\text{INIT}}$):

$$ \delta d = d_{\text{OBS}} - d_{\text{INIT}}. \quad (4) $$

The inverse problem is usually written as:

$$ A\delta m = \delta d, \quad (5) $$

where $A$ is the $N \times M$ matrix of partial derivatives, $N$ is the number of observations, $M$ is the number of model parameters and $\delta m$ is the perturbation to the initial model. Since the number of unknowns (the number of elements of $\delta m$) is smaller than the number of data points (number of elements of $\delta d$) ($M < N$), it is well known that eq. (5) implies the following minimization rather than strict equality:

$$ |A\delta m - \delta d|^2 = \text{minimum.} \quad (6) $$

Note that weighting can be included in the above formulation, but we omit details here. In this study, we solve the waveform inversion inverse problem using the conjugate gradient (CG) method. See Appendix A of Kawai et al. (2014) for details. We solve the inverse problem by the CG method using Akaike’s Information Criterion (AIC) as the stopping criterion. Note that the partials are not updated during this process.

## 2.4 Variance

When we use a relatively small amount of waveform data, we can directly compute $A\delta m$. However, note that we can easily compute the total variance $|\delta d - A\delta m|^2$, as follows:

$$ \begin{aligned} |\delta d - A\delta m|^2 &= (\delta d - A\delta m)^T (\delta d - A\delta m) \\ &= (\delta d^T - (A\delta m)^T)(\delta d - A\delta m) \\ &= \delta d^T \delta d - \delta d^T A\delta m - (A\delta m)^T \delta d + (A\delta m)^T A\delta m \\ &= \delta d^T \delta d - (A^T \delta d)^T \delta m - \delta m^T (A^T \delta d) \\ &\quad + \delta m^T (A^T A) \delta m. \end{aligned} \quad (7) $$

Since $A^T \delta d$ and $A^T A$ are $M \times 1$ and $M \times M$ matrices, respectively, the above computation will always be feasible.

# 3 INVERSION RESULTS

The S-wave velocities at points 400 km or more above the CMB are fixed to PREM, while those in the target region (within 400 km of the CMB) are the unknown parameters. We assume that the density $\rho$ is not perturbed and consider only perturbations to the shear modulus, $\mu$. The relation between the perturbation to the shear modulus $\Delta\mu$ and the perturbation to the shear velocity $\Delta\beta$ is given by

$$ \Delta\beta = \frac{\Delta\mu}{2\rho\beta}. \quad (8) $$
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Recovered models for checkerboard heterogeneity input model for each depth. The input model is shown at top.</figcaption></figure>

We use two different 1-D initial models for the structure in the target region. One is PREM, and the second is ‘NEW 1-D’, which is discussed below. The study area is at latitudes between -5° and 20° and longitudes between 150° and 175°. We divide the studied volume into 5° × 5° × 50 km voxels. The orientation of the voxelization is chosen on the basis of the ray path directions through the target region (Fig. 1b). The number of unknown parameters

is 144. We conduct inversion using the first *n* basis vectors obtained by CG. We choose the value of *n* which minimizes AIC (Akaike 1977). We show the AIC values for the empirical redundancy parameter $\alpha$ = 1000 (see Kawai et al. 2014). The total number of data points is 175 514. The AIC values in Table 2 are thus calculated using $ND = 175514/1000$. Defining the variance of the data to be 100 per cent, the variance and AIC for each model are
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Recovered model for checkerboard input model used in Fig. 7 with white noise with a variance of 70 per cent added to the synthetics.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Input models for block model test. Each figure shows the input velocity perturbation in the lowermost 400 km of the mantle.</figcaption></figure>

shown in Fig. 2 and Table 2. The variance using PREM as the initial model is 82.4 per cent; thus the variance reduction is 17.6 per cent (i.e. 100–82.4 per cent). A further variance reduction of 9.6 per cent is achieved by making the static corrections. The variance for the 3-D model obtained by the inversion using PREM as the starting model, which minimizes AIC using the first three CG basis vectors, is 60.6 per cent, as shown in Fig. 2 and Table 2. The average and standard deviation of the time shifts (static corrections) was −1.46 ± 2.01 s.

In order to study the extent to which the results of the inversion might depend on the starting model, we invert for an alternate 1-D model, which we call ‘NEW 1-D’, and use that model, in place of PREM, as the starting model in the target region. Waveform inversion for ‘NEW 1-D’ is conducted by the CG method using one CG vector. The obtained 1-D model is shown in Fig. 3. A 0.05 km s⁻¹ low S-velocity zone (relative to PREM) is found in the depth range from about 2500–2750 km. An increase in the S velocity of about

0.05 km s⁻¹ is found, peaking around a depth of 2800–2850 km, which is consistent with the phase transition from pv to ppv expected for average composition models such as pyrolite. The variance for ‘NEW 1-D’ with time shifts (i.e. after static corrections) is 63.6 per cent. Note that ‘NEW 1-D’ is not intended as a realistic model (e.g. there is an unphysical velocity discontinuity at the top). The 3-D model obtained using ‘NEW 1-D’ as the initial model achieves a further variance reduction of 3.7 per cent. The 3-D models in Figs 4 and 5 were both obtained using the first three CG basis vectors ($n = 3$).

The 3-D models obtained by inverting with respect to PREM and ‘NEW 1-D’ are shown in Figs 4 and 5, respectively. The models have the same general pattern, especially in the deeper depth ranges. There are two marked low-velocity zones (indicated by arrows) immediately above the CMB (0–50 km) NW and SE of the high-velocity zone (panel h of Figs 4a and 5a). The low and high-velocity perturbations reach around −9 and 4 per cent, respectively. A similar pattern of velocity anomalies is seen in the depth range from 50 to 100 km above the CMB (panel g of Figs 4a and 5a), although the pattern shifts somewhat to the NW and the magnitudes of the anomalies become smaller. In the next shallower depth range from 100 to 150 km above the CMB (panel f of Figs 4a and 5a), the velocity in the central part (indicated by an arrow) of the target region is 9 per cent slow. In the remaining depth ranges (panels a-d of Figs 4a and 5a), we can see similar patterns to that of 100–150 km above the CMB (panel f of Figs 4a and 5a), but with decreasing magnitudes.

Six cross sections for the inversion with respect to PREM and for the inversion with respect to the ‘NEW 1-D’ model are shown in Figs 4(b) and 5(b), respectively. The slow region (indicated by arrows in panel i of Figs 4b and 5b) immediately above the CMB located at the northern end of the target region corresponds to the
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Recovered model for block input model in Fig. 9a (a) for 3, (b) 5 and (c) 10 CG basis vectors, respectively.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Recovered model for block input model in Fig. 9b (a) for 3, (b) 5 and (c) 10 CG basis vectors, respectively.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 12. Recovered model using the model obtained by the 3-D inversion shown in Fig. 4 as the input model.</figcaption></figure>

low-velocity region shown in panel f of Figs 4(b) and 5(b). There is a high-velocity anomaly immediately above the CMB beneath the low-velocity anomaly in the depth range of 100–300 km around the centre of the target region (indicated by arrows in panel ii of Figs 4b and 5b). The two low-velocity anomalies (indicated by arrows in panel h of Figs 4b and 5b) immediately above the CMB extend up to 100 km above the CMB (panels i, iii, v and vi of Figs 4b and 5b). As the low-velocity anomaly which emerges in the depth range of 100–150 km (indicated by the arrow in panel f of Figs 4b and 5b) above the fast anomaly immediately above the CMB is contiguous with the two low-velocity anomalies immediately above the CMB (panels v and vi of Figs 4b and 5b), the laterally heterogeneous structure can be described as ‘tower-shaped’ low-velocity anomalies surrounded by high-velocity regions.

## 3.1 Comparison of basis vectors

If we use all of the CG vectors, they theoretically form a complete set (ignoring the effect of round-off errors). However, the data set does not have nearly enough resolving power to justify using more than a small number of CG basis vectors. In this study, we terminate the CG expansion using AIC, but in any case some termination is required. In order to further study the dependence of our 3-D models on the initial model, we show the correlation coefficients between the CG basis vectors for the 3-D inversions with respect to PREM and the 3-D inversion with respect to ‘NEW 1-D’ in Fig. 6.

We define the angle (in the model space) between the ith CG vectors for the respective inversions as follows:

$$ \cos \theta_i = \frac{p_i^{\text{PREM}} \cdot p_i^{\text{NEW 1-D}}}{|p_i^{\text{PREM}}||p_i^{\text{NEW 1-D}}|} \quad (9) $$

where $p_i$ is the ith CG vector. As shown in Fig. 6 the angle is a bit less than 45° for i = 1, and is about 60° for i = 2. The angles between the respective third through fifth CG vectors range between these two values. It therefore inevitably follows that the 3-D models shown in Figs 4 and 5 will be somewhat different, as they are linear combinations of different basis vectors.

## 3.2 Resolution tests

We conduct a synthetic resolution (‘checkerboard’) test (Fig. 7) to examine the ability of our methods to resolve a synthetic structure model under ideal conditions. Since the input model is synthetic, there are no static corrections. Synthetic seismograms were calculated using the Born approximation. We confirm that waveform inversion can resolve the lateral heterogeneity well in all depth ranges for this region, although the resolution in the depth range from 100 to 200 km above the CMB is relatively low. Waveform inversion can resolve laterally heterogeneous structure from waveforms propagating only in a relatively limited range of azimuths because of the large amount of information contained in the waveforms, which are linearly independent. As a further test, we add white noise with a variance of 70 per cent to the synthetic data and repeat the checkerboard test. As shown in Fig. 8 this test is basically able to recover the input model, although the resolution is somewhat poor at the edges of the target region.

We also conduct a synthetic test for the block models shown in Fig. 9. We consider two models which have 1 per cent positive and negative velocity perturbations in the lowermost 400 km, respectively. (Strictly speaking, we make positive and negative perturbations of exactly equal amplitude to the shear modulus, so the amplitudes of the positive and negative velocity perturbations differ
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Results of the inversion for the 3-D shear wave velocity structure without static corrections. Details same as Fig. 4.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 14. Results of the first three jack-knife tests.</figcaption></figure>
<table><caption>**Table 3.** Variance for PREM and for the three models shown in Fig. 14.</caption><thead><tr><th>Dataset #</th><th>PREM</th><th>3-D</th></tr></thead><tbody><tr><td>1</td><td>73.7</td><td>60.2</td></tr><tr><td>2</td><td>72.4</td><td>58.8</td></tr><tr><td>3</td><td>73.1</td><td>59.8</td></tr></tbody></table>

slightly.) The input model and the starting model are both fixed to PREM above the lowermost 400 km of the mantle. The obtained models are shown in Figs 10 and 11, respectively. For a small number of CG vectors we are unable to fully recover the block input model, but as the number of CG vectors is increased the recovery is reasonably good. Thus the limitation on recovering a model in the inversion of actual data is not our methods *per se* but rather the limited number of events and receivers, which limits the number of CG basis vectors that can be used in the actual inversions.

<figure><img src="image_19.png" /><figcaption>**Figure 15.** (a) Three groups of events. (b) Three groups of stations. (c) Three groups of turning points. (d–f) The three groups of ray paths for events (d), stations (e) and turning points (f).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 16.** The inversion results for the first, second and third data sets separated by event as shown in Figs 15(a) and (d). The numbers at the left indicate the group of events for each set of panels.</figcaption></figure>

<table><caption>**Table 4.** Variance for PREM and for the three models shown in Fig. 16.</caption><thead><tr><th>Event #</th><th>PREM</th><th>3-D</th></tr></thead><tbody><tr><td>1</td><td>67.9</td><td>57.5</td></tr><tr><td>2</td><td>72.6</td><td>58.1</td></tr><tr><td>3</td><td>81.1</td><td>55.6</td></tr></tbody></table>

As a further test, in order to confirm that data set has certain resolution power for the obtained model, we conduct a synthetic resolution test with the 3-D model shown in Fig. 4 as the input model. The synthetics for this test are noise-free. Fig. 12 shows the result from inverting the synthetic data set, which satisfactorily recovers the input model. We do not present model resolution matrix results here, but they are roughly comparable to those shown by Kawai *et al.* (2014) for the regions beneath central America.
<figure><img src="image_21.png" /><figcaption>**Figure 17.** The inversion results for the first, second and third data sets separated by station as shown in Figs 15(b) and (e). The numbers at the left indicate the group of stations for each view.</figcaption></figure>

<table><caption>**Table 5.** Variance for PREM and for the three models shown in Fig. 17.</caption><thead><tr><th>Station #</th><th>PREM</th><th>3-D</th></tr></thead><tbody><tr><td>1</td><td>68.7</td><td>58.8</td></tr><tr><td>2</td><td>80.5</td><td>59.1</td></tr><tr><td>3</td><td>67.8</td><td>57.1</td></tr></tbody></table>

### 3.3 Influence of time shift

We address the possibility of artefacts due to the time shifts (static corrections). To address this question we conduct an inversion relative to PREM without making any time shifts. A data set of 1597 time windows was used, and a variance reduction of 9.6 per cent (relative to synthetics for PREM when no time shift was made) was achieved. The resulting model is shown in Fig. 13. The model has essentially the same character as the models obtained
<figure><img src="image_22.png" /><figcaption>**Figure 18.** The inversion results for the first, second and third data sets separated by turning point as shown in Figs 15(c) and (f). The numbers at the left indicate the group of turning points for each view.</figcaption></figure>
<table><caption>**Table 6.** Variance for PREM and for the three models shown in Fig. 18.</caption><thead><tr><th>Bottom #</th><th>PREM</th><th>3-D</th></tr></thead><tbody><tr><td>1</td><td>70.5</td><td>60.2</td></tr><tr><td>2</td><td>73.4</td><td>57.2</td></tr><tr><td>3</td><td>72.4</td><td>52.8</td></tr></tbody></table>

by the inversion with time shifts shown in Fig. 4. Thus, we conclude that the pattern of the 3-D velocity models is not an artefact of the time shifting, and that hypothetical errors in the time shifting procedure might cause small errors in the absolute velocities, but will not change the basic spatial pattern of the models.

## 3.4 Jack-knife test

We take random subsets of 50 per cent of our data set (jack-knife) and invert for velocity structure. We conduct 10 such inversions. The results of three such inversions are shown in Fig. 14. The features of all three models are almost identical to those in Fig. 4. The

results for the 0–50 km layer obtained by the remaining seven inversions are also nearly identical (figure not shown). Table 3 shows the variance for each data set for the initial model (PREM) after the time shift, and the variance for the 3-D model obtained by the inversion. Note that the variance reduction for each data set is roughly comparable.

## 3.5 Dependence on data set

We now consider the dependence on events, stations and turning points by dividing the data set as shown in Fig. 15. We divide the data set into three separate subsets by source location (Figs 15a and d), station location (Figs 15b and e) and turning point (Figs 15c and f). Fig. 16 and Table 4 show inversion results for the three groups separated by event location. The models vary but we can see the same patterns, especially in deeper depth ranges. Fig. 17 and Table 5 show inversion results for the three groups separated by station location, and Fig. 18 and Table 6 show results for the three groups separated by turning point. It is true that we can see some dependence on events, stations and turning points, but the models in Figs 16–18 show the same general patterns as Fig 4.

<figure><img src="image_23.png" /><figcaption>**Figure 19.** Panels show the value of the cosine between the partial derivative for the *i*-th point in the target region and the partial derivative for shallow structure. The data set for the respective panels consists of one waveform for one event (a), all waveforms for one event (b), all waveforms for one station (c) and all waveforms for all events (d). The number of waveforms in the data set is 1, 59, 29 and 1597, respectively.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure 20.** Initial 1-D model for inversion to test influence of starting model outside the target region, shown as a perturbation to PREM.</figcaption></figure>

## 3.6 Effect of structure outside target region

In order to examine the possible effects of shallow structure, such as the slab near Japan, on the inversion results, we checked the independence of the partial derivatives for perturbations in such regions and those in the target region in the lowermost mantle. We compute partial derivatives at intervals of 100 km with respect to structure in the depth range 100–500 km beneath the yellow star in Fig. 1(b). The location of the target region and the shallow perturbation points is shown in Fig. 1(c). We compute the cosines between the ith partial derivative in the target region and the jth partial derivative in the shallower part as shown in Fig. 19 (see discussion of fig. 18 of Kawai et al. 2014, for details.) We select one event (2003 July 27) and station (KSK) for use in Figs 19(a)–(c). We can see that when all events and all stations are taken together (Fig. 19d) the partials are close to orthogonal. The data set in this study is thus large enough that neglecting effects from structure outside the target region will not greatly affect the inversion results. On the other hand, as shown by Fig. 19(b), if we were processing data for only one event at a time, as is the case in stacking studies, for example, then there could be significant artefacts due to structure outside the

target region. Note that we have not considered the possibility of trade-off between CMB topography and lower mantle structure (e.g. Takeuchi 2005; Colombi et al. 2012).

We also conduct another experiment to study the effect of 1-D structure outside the target region on the results of the inversion. We define a 1-D model which is roughly based on SH18CEX (Takeuchi 2007, 2012) from the Earth's surface to 400 km above the CMB, but is fixed to PREM in the lowermost 400 km of the mantle (Fig. 20). Using this starting model, we obtain the 3-D model shown in Fig. 21. We can see that the starting model outside the target region has a relatively small effect on the results of the inversion. The models in Figs 4(a) and 21 are broadly similar, although there is some difference in the amplitudes of the lateral variations.

## 3.7 Dependence on grid spacing

As shown in Fig. 14, the jack-knife tests yield inversion results which are highly reproducible regardless of which set of randomly chosen waveforms are used. Thus, the inversions are yielding reliable and repeatable models. On the other hand, there may be questions about how finely the model should be sampled, or whether smoothness constraints should be applied.

As one step towards resolving this issue we conducted the three inversions whose results are shown in Figs 22(a), (b) and (c). The inversion in Fig. 22(a) uses twice as coarse a grid spacing in each horizontal direction, the inversion in Fig. 22(b) uses twice as coarse a grid spacing in each horizontal direction and vertically, and the inversion in Fig. 22(c) uses twice as coarse a grid spacing vertically. The coarsening of the grids was performed using rough approximations. The models in Fig. 22 show the same basic features as those in Fig. 4, but are of course smoother. These features are also generally consistent with other tomographic studies (e.g. Takeuchi 2007).

## 3.8 Example profile

Fig. 23 shows an example of a profile of observed waveforms (red), synthetics for PREM with time shifts (green) and the final model (CG3) obtained by the inversion (green). Improvement in the fit to

<figure><img src="image_25.png" /><figcaption>**Figure 21.** Results of the inversion for the 3-D shear wave velocity structure with respect to the initial model shown in Fig. 20. Details same as Fig. 4(a).</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure 22.** (a) Same as Fig. 4(a), but using twice as large a grid interval horizontally. (b) Same as Fig. 4(a), but using twice as large a grid interval vertically. (c) Same as Fig. 4(a), but using twice as large a grid interval horizontally and vertically.</figcaption></figure>
<figure><img src="image_27.png" /></figure>

**Figure 23.** A record section of the transverse component waveforms for the event of 2006 June 9. The traveltime is reduced at an apparent velocity of 8.3 km s⁻¹. A bandpass filter at periods between 12.5 and 200 s is applied. The red, green and blue traces show observed waveforms, synthetics for the initial model (PREM) with time shifts, and synthetics for the final model (CG3), respectively. ScS and other phases associated with complex structure such as the D'' discontinuity are almost entirely overlapped at the epicentral distances used in this study for this period band.

the core phases can be seen. Note however that this profile should be regarded as strictly ancillary to the waveform inversion using the whole data set.

# 4 DISCUSSION

We found two low-velocity anomalies immediately above the CMB located near the NW and SE edges of our study area (panel h of Figs 4a and 5a), and a low-velocity anomaly located at the centre of the target region at depths of 100–150 km above the CMB (panel f of Figs 4a and 5a). These anomalies are contiguous, producing a ‘tower-shaped’ low-velocity anomaly with a height of 300 km (panels v and vi of Figs 4b and 5b). There is a high-velocity anomaly between the slow regions in panel h of Figs 4(a) and 5(a), while the regions in panels a to e of Figs 4(a) and 5(a) above the two slow regions in panel h of Figs 4(a) and 5(a) have high velocity. The low-velocity anomaly located at the horizontal and vertical centre of the target region (indicated by the arrow in panel f of Figs 4a and 5a) is also surrounded by high-velocity anomalies.

He *et al.* (2006) studied the lowermost mantle structure at the boundary of the Pacific LLSVP using forward modelling for waveforms whose ray paths are in the NW-SE direction, and reported a 220 km thick 1 per cent low-velocity anomaly with a 30 km thick 13 per cent ULVZ surrounded by a high-velocity region, with the boundary between the high and low-velocity regions in the NE-SW direction (i.e. perpendicular to the direction of the ray paths). On the other hand, our model shows a 300 km high tower-shaped

low-velocity structure with two roots located in the NW-SE direction. The height of the low-velocity region is consistent with the structure obtained by forward modelling by He *et al.* (2006). He & Wen (2009) conducted forward modelling using waveforms whose ray paths are in the NW-SE direction and found a 740 km thick 3 per cent low-velocity anomaly with a few hundred km thick 5 per cent ULVZ. Our results show low-velocity anomalies of about 8 per cent concentrated immediately above the CMB. The thick 3 per cent low-velocity anomaly which He & Wen (2009) obtained seems to be generally consistent with the depth integrated anomaly for our model.

Takeuchi & Obara (2010) conducted differential traveltime analysis with respect to the ScS bounce point at the CMB using Hi-net waveforms whose ray paths are in the NW-SE direction and found a low-velocity barrel roof structure whose ridge is in the same direction as the ray paths. Although their obtained low-velocity structure is qualitatively consistent with our models when integrated in the NW-SE direction, the direction of their ridge is different from ours. This may be due to the fact that their methods implicitly assume laterally homogeneous structure along the ray path (in the NW-SE direction).

Takeuchi (2012) inverted waveforms for laterally heterogeneous global structure and found about 2.5 per cent average low-velocity structure in the 400 km layer immediately above the CMB in this region. Since our model has 0.5 and 1.5 per cent low-velocity anomalies for the lowermost 400 and 150 km of the mantle, respectively, the global tomographic model of Takeuchi (2012) is generally consistent with our results. The slight difference between the two models could be due to a difference in sensitivity, data sets used, model parametrization, etc. As the ScS phases which were used to resolve the lowermost 400 km mantle by Takeuchi (2012) are most sensitive to structure near the turning depth, the apparent average velocity obtained in the lowermost 400 km of the mantle can be lower than the actual average velocity when low-velocity regions exist immediately above the CMB, as was inferred by the present study.

The target region in this study is located at the edge of the LLSVP beneath the western Pacific that was suggested by previous tomographic studies (e.g. Grand 2002; Takeuchi 2012). The 300 km ‘tower-shaped’ low-velocity anomaly surrounded by high-velocity anomalies found by this study is consistent with the evidence for a sharp-sided boundary of the Pacific LLSVP reported by He *et al.* (2006) and Takeuchi *et al.* (2008), but note that we found structure which is heterogeneous inside the Pacific LLSVP.

If the two low-velocity anomalies immediately above the CMB found by this study are interpreted as due to temperature anomalies, the region up to 100 km above the CMB would reach the CMB temperature of 3800 K (Kawai & Tsuchiya 2009), but the temperature could be even higher than 3800 K in the lowermost 50 km of the mantle due to internal heating, because a temperature deviation of about 1500 K from the mantle adiabat at the CMB of 2500 K (Brown & Shankland 1981) would produce only around a 5 per cent velocity reduction for the post-perovskite phase (Wentzcovitch *et al.* 2006). It is important to note that the two low-velocity zones found by this study immediately above the CMB are laterally separated. This suggests that the low-velocity region at the western edge of the Pacific LLSVP might be an aggregate of smaller thermal plumes rather than a single large chemically distinct pile. The ‘tower-shaped’ low-velocity anomalies suggest that upwelling flows from two spatially separated sources of hot material immediately above the CMB merge into a single zone of hot material at a height of about 100–150 km above the CMB, which might be related to an
upwelling flow. The location of the hypothetical upwelling flow is consistent with geological studies which suggest that the location of the source of hotspots is at the edge of LLSVPs (e.g. Burke *et al.* 2008). Also, the location of the prominent low-velocity anomaly in panel f of Fig. 4(a) coincides with the location of the Caroline hotspot shown in fig. 10(a) of Wessel & Kroenke (1998).

Since the absolute values of the velocity anomalies inferred by this study are less certain than the relative values, some caution should probably be exercised in discussing the geodynamic implications of the absolute values. However, taking the 9 per cent low-velocity anomalies found by this study at face value, it is difficult to explain this only in terms of simple thermal effects. An anomaly of about 4 per cent could be due to chemical anomalies such as iron-rich material (Tsuchiya & Tsuchiya 2006). Such chemical anomalies could be produced as drift by mantle flow in the LLSVP. Then, isolated chemically distinct regions could be intensively heated due to heat flux from the core and radioactive heating.

## 5 CONCLUSION

We conducted waveform inversion for 3-D structure in the lowermost mantle beneath the western Pacific. We analysed a data set obtained by a dense array network of seismometers, the Japanese F-net, for deep and intermediate earthquakes which occurred near Tonga and Fiji. The 3-D model obtained for D'' beneath the Western Pacific is 'tower-shaped'. The tower shape indicates upwelling flows at the CMB, suggesting that the western edge of the LLSVP consists of an aggregation of small thermal plumes.

We conducted an analysis of the dependence of the inversion results on various factors. We found that the results varied somewhat for data sets consisting of waveforms for different events, stations or turning points, but are basically robust.

## ACKNOWLEDGEMENTS

F-net data were obtained from servers at NIED. This research was supported by grants from the JSPS (Nos 24840020 and 25400442) and the Korean government (MEST, No. 2009-0092790).

## REFERENCES

* Akaike, H., 1977. An extension of the method of maximum likelihood and the Stein's problem, *Ann. Inst. Statist. Math.*, **29**, 153-164.

* Brown, J.M. & Shankland, T.J., 1981. Thermodynamic parameters in the Earth as determined from seismic profiles, *Geophys. J. R. astr. Soc.*, **66**, 579-596.

* Bull, A.L., McNamara, A.K. & Ritsema, J., 2009. Synthetic tomography of plume clusters and thermochemical piles, *Earth planet. Sci. Lett.*, **278**, 152-162.

* Bunge, H.-P., Richards, M., Lithgow-Bertelloni, C., Baumgardner, J., Grand, S.P. & Romanowicz, B., 1998. Time scales and heterogeneous structure in geodynamic Earth models, *Science*, **280**, 91-95.

* Burke, K., Steinberger, B., Torsvik, T.H. & Smethurst, M.A., 2008. Plume generation zones at the margins of large low shear velocity provinces on the core-mantle boundary, *Earth planet. Sci. Lett.*, **265**, 49-60.

* Colombi, A., Nissen-Meyer, T., Boschi, L. & Giardini, D., 2012. Seismic waveform sensitivity to global boundary topography, *Geophys. J. Int.*, **191**, 832-848.

* Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference earth model, *Phys. Earth planet. Inter.*, **25**, 297-356.

* Fuji, N., Kawai, K. & Geller, R.J., 2010. A methodology for inversion of broadband seismic waveforms for elastic and anelastic structure and its

application to the mantle transition zone beneath the Northwestern Pacific, *Phys. Earth planet. Inter.*, **180**, 118-137.

* Garnero, E., Grand, S.P. & Helmberger, D.V., 1993. Low P wave velocity at the base of the mantle, *Geophys. Res. Lett.*, **20**, 1843-1846.

* Geller, R.J. & Hara, T., 1993. Two efficient algorithms for iterative linearized inversion of seismic waveform data, *Geophys. J. Int.*, **115**, 699-710.

* Geller, R.J. & Ohminato, T., 1994. Computation of synthetic seismograms and their partial derivatives for heterogeneous media with arbitrary natural boundary conditions using the direct solution method, *Geophys. J. Int.*, **116**, 421-446.

* Grand, S.P., 2002. Mantle shear-wave tomography and the fate of subducted slabs, *Phil. Trans. R. Soc. Lond., A.*, **360**, 2475-2491.

* Hara, T., Tsuboi, S. & Geller, R.J., 1991. Inversion for laterally heterogeneous earth structure using a laterally heterogeneous starting model: preliminary results, *Geophys. J. Int.*, **104**, 523-540.

* Hara, T., Tsuboi, S. & Geller, R.J., 1993. Inversion for laterally heterogeneous upper mantle Swave velocity structure using iterative waveform inversion, *Geophys. J. Int.*, **115**, 667-698.

* He, Y. & Wen, L., 2009. Structural features and shear-velocity structure of the "Pacific Anomaly", *J. Geophys. Res.*, **114**, B02309, doi:10.1029/2008JB005814.

* He, Y., Wen, L. & Zheng, T., 2006. Geographic boundary and shear wave velocity structure of the "Pacific anomaly" near the core-mantle boundary beneath western Pacific, *Earth planet. Sci. Lett.*, **244**, 302-314.

* Idehara, K., 2011. Structural heterogeneity of an ultra-low-velocity zone beneath the Philippine Islands: implications for core-mantle chemical interactions induced by massive partial melting at the bottom of the mantle, *Phys. Earth planet. Inter.*, **184**, 80-90.

* Kawai, K. & Geller, R.J., 2010. Waveform inversion for localized seismic structure and an application to D'' structure beneath the Pacific, *J. geophys. Res.*, **115**, B01305, doi:10.1029/2009JB006503.

* Kawai, K. & Tsuchiya, T., 2009. Temperature profile in the lowermost mantle from seismological and mineral physics joint modeling, *Proc. Natl. Acad. Sci. USA*, **106**, 22 119-22 123.

* Kawai, K., Takeuchi, N. & Geller, R.J., 2006. Complete synthetic seismograms up to 2 Hz for transversely isotropic spherically symmetric media, *Geophys. J. Int.*, **164**, 411-424.

* Kawai, K., Konishi, K., Geller, R.J. & Fuji, N., 2014. Methods for inversion of body-wave waveforms for localized three-dimensional seismic structure and an application to D'' structure beneath Central America, *Geophys. J. Int.*, **197**, 495-524.

* Kito, T., Krüger, F. & Negishi, H., 2004. Seismic heterogeneous structure in the lowermost mantle beneath the southwestern Pacific, *J. geophys. Res.*, **109**, B09304, doi:10.1029/2003JB002677.

* Konishi, K., Kawai, K., Geller, R.J. & Fuji, N., 2009. MORB in the lowermost mantle beneath the western Pacific: evidence from waveform inversion, *Earth planet. Sci. Lett.*, **278**, 219-225.

* Konishi, K., Kawai, K., Geller, R.J. & Fuji, N., 2012. Waveform inversion of broad-band body wave data for the S-velocity structure in the lowermost mantle beneath the Indian subcontinent and Tibetan plateau, *Geophys. J. Int.*, **191**, 305-316.

* Poirier, J.-P., 2000. *Introduction to the Physics of the Earth's Interior*, Cambridge Univ. Press.

* Romanowicz, B., 2003. Global mantle tomography: progress status in the past 10 years, *Ann. Rev. Earth planet. Sci.*, **31**, 303-328.

* Schubert, G., Masters, G., Olson, P. & Tackley, P.J., 2004. Superplumes or plume clusters?, *Phys. Earth planet. Inter.*, **146**, 147-162.

* Takeuchi, N., 2005. Finite boundary perturbation theory for the elastic equation of motion, *Geophys. J. Int.*, **160**, 1044-1058.

* Takeuchi, N., 2007. Whole mantle SH velocity model constrained by waveform inversion based on three dimensional Born kernels, *Geophys. J. Int.*, **169**, 1153-1163.

* Takeuchi, N., 2012. Detection of ridge-like structures in the Pacific large low-shear-velocity province, *Earth planet. Sci. Lett.*, **319**, 55-64.

* Takeuchi, N. & Obara, K., 2010. Fine-scale topography of the D'' discontinuity and its correlation to volumetric velocity fluctuations, *Phys. Earth planet. Inter.*, **183**, 126-135.
Tsuchiya, T. & Tsuchiya, J., 2006. Effect of impurity on the elasticity of perovskite and postperovskite: velocity contrast across the postperovskite transition in (Mg, Fe, Al)(Si, Al) O₃, *Geophys. Res. Lett.*, **33**, L12S04, doi:10.1029/2006GL025706.

Takeuchi, N., Morita, Y., Xuyen, N.D. & Zung, N.Q., 2008. Extent of the low-velocity region in the lowermost mantle beneath the western Pacific detected by the Vietnamese Broadband Seismograph Array, *Geophys. Res. Lett.*, **35**, L05307, doi:10.1029/2008GL033197.

Wentzcovitch, R.M., Tsuchiya, T. & Tsuchiya, J., 2006. MgSiO₃ postperovskite at D'' conditions, *Proc. Natl. Acad. Sci. U.S.A.*, **103**, 543–546.

Wessel, P. & Kroenke, L.W., 1998. The geometric relationship between hot spots and seamounts: implications for Pacific hot spots, *Earth planet. Sci. Lett.*, **158**, 1–18.

Wysession, M.E., Bartkó, L. & Wilson, J.B., 1994. Mapping the lowermost mantle using core-reflected shear waves, *J. Geophys. Res.*, **99**, 13 667–13 684.