Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29
DOI 10.1186/s40645-017-0143-8

Progress in Earth andPlanetary Science

RESEARCH ARTICLE

Open Access

Adjoint tomography of the crust and upper
mantle structure beneath the Kanto region
using broadband seismograms
Takayuki Miyoshi1,2*, Masayuki Obayashi3, Daniel Peter4, Yoko Tono5 and Seiji Tsuboi6

Abstract
A three-dimensional seismic wave speed model in the Kanto region of Japan was developed using adjoint tomography
for application in the effective reproduction of observed waveforms. Starting with a model based on previous travel time
tomographic results, we inverted the waveforms obtained at seismic broadband stations from 140 local earthquakes in
the Kanto region to obtain the P- and S-wave speeds Vp and Vs. Additionally, all centroid times of the source solutions
were determined before the structural inversion. The synthetic displacements were calculated using the spectral-element
method (SEM) in which the Kanto region was parameterized using 16 million grid points. The model parameters Vp andVs were updated iteratively by Newton’s method using the misfit and Hessian kernels until the misfit between theobserved and synthetic waveforms was minimized. Computations of the forward and adjoint simulations wereconducted on the K computer in Japan. The optimized SEM code required a total of 6720 simulations using
approximately 62,000 node hours to obtain the final model after 16 iterations. The proposed model reveals several
anomalous areas with extremely low-Vs values in comparison with those of the initial model. These anomalies were foundto correspond to geological features, earthquake sources, and volcanic regions with good data coverage and resolution.
The synthetic waveforms obtained using the newly proposed model for the selected earthquakes showed better fit than
the initial model to the observed waveforms in different period ranges within 5–30 s. This result indicates that the model
can accurately predict actual waveforms.
Keywords: Seismic wave speed model, Adjoint tomography, Waveform inversion, Broadband seismogram, Kanto region

Introduction
The Japanese metropolitan area in the Kanto plain,
which contains 40 million inhabitants, has a high disaster risk derived from the earthquake and volcanic activities related to the subduction of two oceanic plates, thePacific (PAC) and Philippine Sea (PHS) plates, beneaththe overriding tectonic plates (Fig. 1). Various types ofearthquakes are associated with the convergence of theseplates in the region (e.g., Hori 2006), such as interplate
earthquakes on the upper surface of the PHS or PAC
slab, intraslab earthquakes within the PHS or PAC slab,
and inland earthquakes within the upper crust. The* Correspondence: miyoshi@bosai.go.jp
1
Earthquake Research Institute, University of Tokyo, 1-1-1 Yayoi, Bunkyo-ku,
Tokyo 113-0032, Japan
2
Present address: Earthquake and Tsunami Research Division, National
Research Institute for Earth Science and Disaster Resilience, 3-1 Tennodai,
Tsukuba, Ibaraki 305-0006, Japan
Full list of author information is available at the end of the article

1703 and 1923 great Kanto earthquakes (M8-class) havebeen considered typical examples of interplate earthquakes on the upper surface of the PHS slab and caused
severe damage as a result of their strong shaking andmajor tsunamis. The 1923 event is the most recent destructive earthquake occurring beneath this area andtook the lives of more than 100,000 people in the Kanto
plain. Modern Japanese society in this metropolitan area
has not experienced large earthquakes occurring in thisarea since the 1923 Kanto earthquake. Active faults havebeen identified in the basin despite the thick basin layer,
the thickness of which reaches a maximum of approximately 4 km (e.g., Suzuki 1996). The active faults are
candidate sources of future inland great earthquakes;
however, it is not difficult to imagine that unknown
faults are buried beneath the metropolitan area. Active
volcanoes to the west of the Kanto plain may produce
volcanic disasters on the same scale as the large eruption© The Author(s). 2017 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0
International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to
the Creative Commons license, and indicate if changes were made.

Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 2 of 20

Fig. 1 Distribution of tectonic plates and topography of the study area. Open squares indicate the NIED F-net stations used in the present study.
Red triangles denote active volcanoes. Thin black lines and the broken line denote active faults and the Median Tectonic Line (MTL), respectively.
AMR, NAM, PAC, and PHS in the left figure indicate the Amur, North America, Pacific, and Philippine Sea plates, respectively. BP and IP in the rightfigure indicate the Boso and Izu Peninsulasof Mt. Fuji in 1707. A detailed model of the threedimensional (3D) wave speed structure would greatly
improve the present understanding of seismotectonics
and volcanotectonics and enable the prediction of strong
motion for future earthquakes occurring beneath themetropolitan area.
The crust and upper mantle structure, seismotectonics, subduction process, and other characteristics
based on local tomography using travel time data of Pand S-waves beneath the Kanto region have been reported (e.g., Kamiya and Kobayashi 2000, 2007; Matsubara et al. 2005; Nakajima et al. 2009; Nakajima andHasegawa 2010). These results have contributed to thepresent geophysical understanding of the Kanto region.
However, difficulties arise when observed seismograms
are modeled using numerical simulations based on thesestructure models because it is unknown how the raybased tomographic results quantitatively reproduce observed waveforms. Therefore, seismic tomography basedon wave theory would be a useful approach to address
such difficulties.
From the viewpoint of strong motion prediction,
Koketsu et al. (2008, 2009) proposed a modeling procedure of a 3D wave speed structure to evaluate theground shaking produced by an earthquake. They
used various types of databases, such as geophysicalexperiments, gravity surveys, geological investigations,and seismic observations of natural events, to develop
a layered model of the Kanto basin. Yoshimoto andTakemura (2014) applied a realistic sedimentary structure of the Kanto basin to the Japan integrated
velocity structure model (Koketsu et al. 2008) toreproduce the excitation and dispersion of the observed surface waves. These models reproduced theobserved waveforms but were not suitable to assess
the tectonic processes associated with seismic andvolcanic activities, because they assumed the layered
structure of a homogeneous isotropic media.

Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Recent progress in supercomputer technology has enabled the calculation of accurate seismograms based onwave theory, which requires a large amount of computing resources. In addition to forward modeling, waveform inversion can be conducted by solving an inverse
problem using efficient algorithms (e.g., Tarantola 1984,
1986; Geller and Hara 1993; Tromp et al. 2005). Waveform tomography has been developed as a powerful approach to infer the earth’s structure. In waveform
tomography, the adjoint method is commonly used. In
this method, the gradients of misfit functions with respect to density, shear wave speed, and compressional
wave speed are simply given by the interaction between
the forward wavefields, which propagate from the source
to the receivers, and the adjoint wavefields, which propagate from the receivers back to the source. For each
earthquake, numerical simulations need to be conducted
only twice: once for the construction of the forward andadjoint wavefields and once to obtain the gradients ofthe misfit function. Based on a combination of thespectral-element and adjoint methods, various wave
speed models have been obtained in previous studies, including models of the Southern California crust (Tape etal. 2007, 2009, 2010), the western Pyrenees (Wang et al.
2016), the Australasian upper mantle (Fichtner et al.
2009), the European upper mantle (Zhu et al. 2012,2015; Zhu and Tromp 2013), the East Asian crust andupper mantle (Chen et al. 2015), and the crust andupper mantle around the Japanese archipelago (Simutė
et al. 2016). Bozdag et al. (2016) conducted the firststudy in which a model of the entire crust and mantle ofthe earth was constructed using the iterative fullwaveform inversion technique. Waveform tomography
has become a feasible method of revealing the heterogeneity of the earth’s interior at various spatial scales.
We used waveform tomography with iterative waveform
inversion to obtain a fine structure model through theadjoint method, referred to as “adjoint tomography”(e.g., Tape et al. 2007, 2010).
In this study, a 3D seismic wave speed model of the
Kanto region was constructed to overcome the problems
mentioned above and enable an analysis of both tectonic
features and waveform reproduction using seismogramsof the period band between 5 and 30 s by adjoint tomography. We demonstrated how the initial and final wave
speed structures could qualitatively reproduce seismograms in each frequency band and described the relation
between the wave speed structure and tectonic phenomena. Although it is beyond the scope of the present work
to predict short-period strong motions for the mitigationof earthquake disasters, the basin structure related tothe excitation of surface waves is briefly discussed in thecontext of achieving strong motion prediction in the
Kanto region.

Page 3 of 20

Methods/Experimental
Earthquake and broadband waveform data

After the disastrous 1995 Hyogoken-nanbu (Kobe)
earthquake, numerous seismographs were installed inJapan by the National Research Institute for Earth Science and Disaster Prevention (NIED; renamed the National Research Institute for Earth Science and Disaster
Resilience in 2016). The installation of the broadband
seismograph networks was started as the Fundamental
Research on Earthquakes and Earth’s Interior Anomaly(FREESIA) project, and the current full-range seismograph network of Japan (F-net) provides nationwide
coverage with approximately 70 stations (e.g., Okada etal. 2004). These seismographs are mostly located in tunnels to reduce noise from various sources and observe
the ground velocities of three components. The waveforms are recorded continuously and distributed to the
public through the Internet. The instrument conditions
are systematically monitored to maintain data quality(Kimura et al. 2015). The NIED provides an earthquake
catalog of moment tensor solutions based on F-net seismograms for regional earthquakes in and around theJapanese archipelago (e.g., Kubo et al. 2002). In the F-net
earthquake catalog, the centroid depth and moment tensor are determined; however, other parameters are fixed
to rapid information reported by the Japan Meteorological Agency. The source parameters used in thepresent study are discussed in a later subsection.
The earthquakes were selected from the F-net
earthquake catalog according to the following selection criteria. The events must have (1) occurred between 1998 and 2014, (2) a focal depth shallower
than 80 km, (3) a magnitude in the range of
4.0 ≤ MW ≤ 5.5, and (4) an average variance reduction of waveform fitting between the observed andsynthetic waveforms larger than 85% once the source
parameters have been determined. Based on these criteria, 140 earthquakes were selected for inversion
data, with 18 events used to evaluate the tomographic
results obtained in the present study. Figure 2 shows
the selected events. The selected events were classified into three types from a tectonic viewpoint:interplate earthquakes with the P-axis in the E–W
direction on the PAC slab, interplate earthquakes with
the P-axis in the NW–SE direction on the PHS slab,
and inland earthquakes of the upper crust. However,
events that occurred in the upper crust beneath the
Kanto basin could not be detected.
Broadband waveforms recorded at the F-net stationsin and around the Kanto basin were used. A total of
approximately 4400 seismograms with good quality wereobtained at 16 stations for the 140 considered events.
The velocity seismograms were integrated after removing the response of the seismometer from originalMiyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 4 of 20

Fig. 2 Distribution of the CMT solutions used in the present study for (a) adjoint inversion and (b) the validation of the inversion results. The beachball color corresponds to the centroid depth

records to obtain the displacement time series in threecomponents.
Spectral-element method and initial model

Both synthetic and observed waveforms had to be obtained to conduct adjoint tomography. Synthetic displacements were obtained using the spectral-element
method (SEM), a numerical method to solve the equations of motion in the weak form precisely. In previous
studies, the SEM was first applied to solve a fluid dynamics problem (Patera 1984) numerically and was later
used in seismic wave simulations for a complex 3D
structure at various scales (e.g., Komatitsch and Tromp
2002a, 2002b; Peter et al. 2011). We used the SPECFEM3D Cartesian program package (e.g., Peter et al.
2011), which has been developed as a parallel computing
code to simulate seismic waves on a local scale.
A structure model in the Universal Transverse Mercator coordinate system with horizontal dimensions of
459 km × 492 km and a depth of 120 km was consideredto span the entire Kanto basin area and to include all selected earthquakes and F-net stations. The initial and iterative models assumed only a solid earth without thefluid parts. Some bias may be produced by neglecting

the effect of the solid–fluid interface and the ocean on
the waveforms for offshore events. However, this bias
was expected to be negligibly small for the minimum
period (5 s) considered in this study because of the shortdistances the waves traveled. The presence of the oceans
was mostly relevant for the dispersion of Rayleigh waves
in the case of shallow events, as well as for free-surface
reflected phases, such as PP, SS, and SP (e.g., Komatitsch
and Tromp 2002b).
Absorbing boundary conditions were applied to theartificial boundaries of the horizontal and bottom surfaces (Clayton and Engquist 1977), and the topography
was represented in the uppermost edges of the model.
The initial P- and S-wave speed, Vp and Vs, respectively,
were based on the 3D travel time tomography model by
Matsubara and Obara (2011), which resolved structuresin the Japanese archipelago and off the coast by using
the travel times of offshore events. The density ρ wascomputed from Vp for the initial model using an empirical relationship (e.g., Gardner et al. 1974; Brocher
2005). The relation ρ = Vp/6 + 5/3 was used in this study.
The density in the iterative model was updated using ascaling factor for the Vs–density relation (δ ln ρ)/(δ lnVs) = 0.33 (e.g., Montagner and Anderson 1989; PanningMiyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

and Romanowicz 2006). The attenuation structure wasestimated from the relation with Vs (Olsen et al. 2003);
that is, Qμ was set to scale as 0.05 times Vs (in meters
per second). The attenuation structure Qκ related to thecompressional wave was not taken into account in the
inversion.
The mesh structure of the initial model in the present
study is shown in Fig. 3. The mesh used in the SEM
model contained approximately 230,000 elements with
1.6 × 107 grid points. Each element had 125 (5 × 5 × 5)
Gauss–Lobatto–Legendre points that share grid points
with neighboring elements. The mesh had one doublinglayer at a depth of approximately 40 km (Fig. 3). To accurately map the topography and bathymetry of the region to the mesh, a minimum element size of
approximately 400 m at the upper boundary was chosen,
and the mesh was coarsened to a maximum size of8000 m at the lower boundary. The model area was extended from the area of interest shown in Figs. 1 and 2
as a safeguard against the reflection of waves from theartificial boundaries in addition to the absorbing

Page 5 of 20

condition. The simulation was accurate up to a minimum period of 2.6 s. Because the mesh structure of theSEM model was maintained throughout the inversion,
the minimum period could be longer than 2.6 s when
the seismic wave speed in the iterative model was slower
than that of the initial model. The initial and ith iterativemodels are hereafter referred to m00 and mi,
respectively.
The following earthquake source parameters were obtained from the F-net moment tensor catalog: the origin
time, epicentral coordinates, centroid depth, moment
tensor, and seismic moment. No information concerning
the rupture duration was included in the F-net catalog,
and the centroid time was fixed at the same value as theorigin time. According to the analysis of the centroid
moment tensor by Matsumura et al. (2006), the time difference between the origin and centroid times was on
the order of a few seconds for local earthquakes inJapan. Obtaining the centroid moment solution using
the 3D model had a high computational cost so only thecentroid time was determined without estimating otherFig. 3 Mesh and shear wave speed structures of the initial model. The color indicates the S-wave speed Vs

Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 6 of 20

source parameters in this study. The difference between
the origin and centroid times was assumed to be theaverage of the time shifts, which were defined as thetime deviation between the observed and synthetic
waveforms estimated using phase correlation for apacket of P-waves. The synthetic displacements calculated using the SEM for the initial model were comparedwith the observed waveforms in the period band of
5–30 s to estimate the lag times. Figure 4 shows a histogram of the time shifts for all events. Almost all time
shifts were positive, indicating that the earthquakes haverupture times of a couple of seconds. Some events had
negative time shifts, and these negative time shifts may
have been caused by the poor estimation of the source
parameters.

Adjoint tomography procedure

This subsection describes the adjoint tomography procedure used in the present study. Figure 5 shows a summary of the present inversion procedure, which was
similar to that used in previous studies (Tape et al. 2010;
Zhu et al. 2015). The adjoint tomography method isbriefly discussed; more detail on the basic theory of adjoint tomography can be found in previous works (e.g.,
Tromp et al. 2005; Tape et al. 2007, 2010; Zhu et al.
2015). The wave speed model was iteratively updated byminimizing the misfit between the observed and synthetic waveforms. The model parameters were the P-and
S-wave speeds Vp and Vs of the grid points, and thenumber of parameters in this model was approximately
1.6 × 107 × 2. The density and attenuation were updated

Fig. 5 Research flow of adjoint tomography in the present study.
The computational cost estimated for the K computer is also shownusing the new Vs value in each iteration according to thevelocity–density and velocity–attenuation scaling relationships, respectively.
The waveforms for such a 3D SEM mesh structure
were synthesized over a duration of 130 s with a time
step of 0.005 s for each iteration to calculate the misfit
function. The misfit function was calculated for the
waveforms in the time windows of the seismic phases.
The FLEXWIN algorithm developed by Maggi et al.
(2009) was used to determine the time windows, inwhich multiple quantitative criteria, such as the crosscorrelation coefficient between the observed and synthetic waveforms, the time lag, and the signal-to-noise(S/N) ratio of the observation, were defined. We selected2630 windows for the period between 20 and 30 s in thefirst iteration and over 4000 windows for the period between 5 and 30 s in the last iteration. The misfit function for the ith window in a single trace was definedusing the multitaper method (Zhou et al. 2004; Tape etal. 2010) as
F i ðmÞ ¼

Fig. 4 Histogram of the time shift for source parameters used in thepresent study. The positive and negative values indicate positive andnegative delay times from each origin time, respectivelyhi ðωÞ ΔT i ðω; mÞ 2dω;σi−∞ H iZ ∞ð1Þwhere σi is the estimated uncertainty associated withthe travel time measurement, ΔT i ðω; mÞ ¼ T obs
i ðωÞ−T iMiyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 7 of 20

ðω; mÞ indicates the frequency-dependent travel timedifference between the observed and synthetic waveforms, hi(ω) is a window in the frequency domain,and Hi is a normalization constant. To minimize themisfit function using the adjoint method, the timereversed technique was applied to this function, andits product was used as the adjoint source in the
backward simulation from each station to the source(e.g., Tromp et al. 2005). After the adjoint wavefield
was calculated, the event kernel was obtained as theinteraction between the forward and adjoint wavefields in each event. This yielded sensitivity kernels
for compressional and shear wave speeds for all
events, as well as an approximate Hessian kernel,
which was regarded as the second derivative of themisfit function.
All event kernels were summed, and smoothing wasapplied to obtain the misfit kernels for Vp and Vs. The
approximate Hessian kernels were used as preconditioners in Newton’s method (e.g., Tape et al. 2007; Fichtner 2010), and the ith model was updated using thegradient g(m) of the misfit function and the step length
α asmiþ1 ¼ mi –αH −1 ðmÞg ðmÞ;ð2Þ

where H(m) is a Hessian kernel. The gradient and theHessian of the misfit function F are given respectively as
follows:
∂F 
g ð mÞ ¼
;
ð3Þ
∂mm

∂2 F 
:ð4Þ
H ð mÞ ¼∂m∂mmHere, g(m) and H(m) were smoothed before updating
the model. Prominent values of the kernels were prevented from appearing around sources and receivers by
the smoothing effect. Because these Hessian kernels represented the approximation of the diagonal elements in
the Hessian matrix, the Hessian kernel acted like anormalization factor within the sensitive area. The adjoint state method (Tarantola 1984; Plessix 2006) was
employed to obtain the gradient, and the approximate
Hessian was calculated as the vector inner product andconvolution of the forward and adjoint accelerations asfollows:Z~ ðxÞ ¼ ∂2 sðx; t Þ⋅∂2 s† ðx; T −t Þdt;
H
ð5Þ
t
t
where s and s† are the forward and adjoint displacement
wavefields, respectively (Zhu et al. 2015). The model was∼
updated using Eq. (2) assuming HðxÞ≈ HðxÞ . The step
length α was set to values ranging from 1 to 5% inincrements of 1% based on the waveform misfits. The
seismograms of a representative subset of 10 earthquakes were calculated for five candidate models, andthe misfit was evaluated for each model to obtain the
optimal step length. The model with the minimum misfit value was selected as the updated model. For example, α = 4% was selected in the first iteration.
Although the principal goal of this study was to obtain
synthetic waveforms that agree well with the observed
waveforms in the period range of 5–30 s, the inversion
was initiated with waveforms in the period band of 20–
30 s to prevent cycle skipping. Period bands, including
shorter periods, were then used in the following stages:
10–30, 8–30, and 5–30 s. The stage was updated after
the convergence of the iterative inversion was achieved
at each stage. The achieved convergence was determinedusing the variance reduction (VR) (e.g., Kubo et al. 2002)
when (i) the difference between the summations of the
VR values from all events in the current and former iterations was less than 1% or (ii) the current VR value waslarger than that of the previous iteration. The variancereduction is defined as:R∞
wi ðtÞ½dðtÞ−sðt; mÞ2 dt
for the ith window wi ðtÞ;V Ri ðmÞ ¼ 1− −∞ R ∞2−∞ wi ðtÞ½dðtÞ dt

ð6Þ
where d(t) is the amplitude of the data and s(t, m) is thesynthetic waveform of model m. The VR values ranged
from − ∞ to 1, with the latter representing a perfect fit.
The left panels in Fig. 6b–e show the frequency distributions of the VR values over the full time length based on
m00 for four period bands. The waveform fitting of thelonger period bands was better than that of the shorterbands. There was a difference between the observed andsynthetic waveforms even in the longer band (Fig. 6e), so
it was reasonable to first consider the seismograms of
the longer period band.
Computational costs

The present analysis required only 140 pairs of forward
and backward simulations to obtain the gradient of themisfit function in every iteration. However, it still required a high computational cost and extensive time toobtain the final model despite using the adjoint technique. Sufficient computer resources had to be acquired
to conduct the calculations. In both the forward and adjoint simulations, the model region was divided into 256(16 × 16) slices in the horizontal plane, and one processor performed the calculations from bottom to top ineach slice. Parallel computing was used to conduct thecalculation for the 256 slices. The K computer at the
RIKEN Advanced Institute for Computational Science inKobe, Japan, was used for the large computations. The K

Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 8 of 20

Fig. 6 Histograms of the variance reduction VR for the initial and final models m00 and m16. a All windows for the period range of 5–30 s. b Fullwaveform for the period range of 5–30 s. c Full waveform for the period range of 8–30 s. d Full waveform for the period of 10–30 s. e Full waveform
for the period range of 20–30 s. The left, center, and right panels show all events, events that occurred at depths greater than 30 km, and events that
occurred at depths less than 30 km, respectively. The gray bars and black outlines represent the VR values of m00 and m16, respectively. The VR values
can range from − ∞ to 1. The average VR values are also shown in each panel

computer had a very large number of nodes (82,944),
and each node or central processing unit had eightcores. Forward simulations using the SEM on the K
computer have been previously implemented for global
modeling by Tsuboi et al. (2016). They optimized thesource code for the K computer and successfully modeled waveforms with an accuracy of 2 s. In the present

study, the forward and adjoint wavefields were simulatedusing source code optimized in the same manner as
Tsuboi et al. (2016) on 32 nodes (256 cores) with a flat
message passing interface parallel computing model. It
took approximately 10 min to conduct a single forward
simulation and approximately 40 min to obtain a single
event kernel using 256 cores. When sufficient free nodes

Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

were available, 140 forward or adjoint simulations wereconducted simultaneously using 8192 cores. In each iteration, the computational costs to obtain a forward wavefield and a single event kernel were approximately
780 × 8 and 3100 × 8 core hours, respectively. Zhu et al.
(2015) reported that the computational cost of the kernel calculation was twice that of the forward simulation;
however, the present kernel computational cost was fourtimes that of the forward simulation. This highlighted
the well-optimized forward computation concerning theoriginal kernel computation. In general, the cost proportion may depend on the difference between the
optimization levels of the source codes of the forward
and adjoint parts.

Results
Waveform improvement

The final seismic wave speed model was obtained after
16 iterations. Figure 7 shows the line search procedure
for the step length in the first iteration and the iterativetransition of the VR values and the period bands used ineach iteration. In the first iteration, a step length of 4%
was selected by comparing the VR values of five models
to produce the first iterative model m01 (Fig. 7a). The
VR values increased gradually with each iteration(Fig. 7b). Figure 6 shows the results of fitting the synthetic waveform to the observed waveform for m00 andm16. For the waveforms within the time windows, thefitting was improved in the final model for the period
range of 5–30 s (Fig. 6a). Full seismograms were analyzed for the four considered period bands to quantitatively evaluate the waveform fitting. The waveform
fitting was improved by the iterative process in all period
ranges (Fig. 6b–e). The VR value in the period range of
5–30 s was improved by 24%, corresponding to VR
values of − 0.24 and 0.00 for m00 and m16, respectively.
The depth dependence of the waveform fitting was also

Page 9 of 20

analyzed with a depth threshold of 30 km. The fitting
was improved for events occurring at depths both
greater (center panels in Fig. 6) and less (right panels in
Fig. 6) than 30 km. The fitting of the deeper events wasmuch better than that of the shallower events for m16,
as indicated by the VR values. We concluded from theseresults that the proposed iterative process could moreaccurately reproduce waveforms with longer periods or
waveforms of deeper events than it could reproduce
waveforms with shorter periods or waveforms of shallower events.
The observed waveforms for the considered earthquakes were compared with the theoretical ones to perform a detailed investigation of how the waveforms werechanged by the inversion. Figures 8 and 9 show the threecomponents of the displacements of example waveforms
obtained at the JIZ station. The vertical cross sections of
the shear wave speeds of m00 and m16 and their differences along the source and the receiver are also shownin Figs. 8 and 9. As shown in Fig. 8, the radial and vertical components of the m00 waveform differed from thoseof the observed waveform, such as those at the arrival ofsome of the phases near the S-wave, whereas the m16
waveforms reproduced the observation well. Because theshear wave speed between the source and the receiver
was slower than that of the m00 waveform, the m16
waveform was slightly delayed relative to the m00 waveform. As shown in Fig. 9, the three components of thesynthetic waveform for m00 were not well fitted to theobserved ones, and the waveform components for m16
matched better than those of m00 near the arrival of the
S-wave. A faster shear wave speed model between thesource and the receiver was inferred in m16 in comparison with that in m00, and this contributed to the improvement of the waveform models. As shown in Figs. 8
and 9, the theoretical waveforms accurately reproduced
the observed ones near the arrival of the P- and S-waves.

Fig. 7 Waveform fitting in the inversion process. a Step length search in the first iteration. b Waveform misfit between the observed andsynthetic waveforms in each iteration. The fitting value is normalized by the fitting value of the first iteration of each period band. The period
band used in each iteration is also shownMiyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 10 of 20

Fig. 8 Shear wave speed models and waveforms for the earthquake on the PHS (June 2, 2007; MW 4.5). a Source and receiver. b Vertical crosssection of Vs for m00 and m16. c Vertical cross section of the difference between the Vs values for m00 and m16. d Waveform comparisons. The
gray areas in (b) and (c) indicate areas of low data coverage. The crosses and inverted triangles in (b) and (c) show the location of theearthquake and station in (a), respectively. The black, blue, and red lines in (d) indicate the observed, m00, and m16 waveforms, respectively

However, the waveforms around the later arrivals of the
S-wave were not well fitted to the observed ones, because the time windows were set around the P- and Swaves. Later phases must be added as inversion window
data to modify the seismograms, such as the surface
waves derived from the basin; this remains a task for future studies.
Seismic wave speed modelThis section describes the wave speed models obtainedin the present study. The average lateral speeds of m00
and m16 were examined, the wave speeds of m00 and

m16 were compared, and the absolute wave speeds inthe depth slices and their deviations from the average
were assessed. First, some distinct characteristics of the
wave speed models are discussed. Figure 10a shows theaverage wave speeds of m00 and m16 at each depth. The
average shear and compressional wave speeds of m16
were both lower than those of m00 (Fig. 10a). The wave
speeds of m00 obtained at the grid points were comparedwith those of m16 to examine the differences between
m00 and m16 in detail. The Vs and Vp values of m16 were
distributed widely around the Vs and Vp values of m00,
respectively (Fig. 10b, c). Based on a comparison of theMiyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 11 of 20

Fig. 9 Same as Fig. 8 for the wave speed models and waveforms for the earthquake on the PAC (May 13, 2014; MW 5.0). a Source and receiver.
b Vertical cross section of Vs for m00 and m16. c Vertical cross section of the difference between the Vs values for m00 and m16. d Waveformcomparisons. The gray areas in (b) and (c) indicate areas of low data coverage. The crosses and inverted triangles in (b) and (c) show the locationof the earthquake and station in (a), respectively. The black, blue, and red lines in (d) indicate the observed, m00, and m16 waveforms, respectively

shear wave speed (Fig. 10b), inversion tended to reduce
low speeds and increase high speeds.
Figure 11 shows horizontal cross sections of the absolute wave speeds of m00 and m16 and the differences between the wave speeds of the two models at
depths of 5.0, 10.0, and 40.0 km. The S- and P-wave
speeds Vs and Vp of m16 differed by approximately 20
and 10%, respectively, from those of m00. The overall
patterns of the new Vs values were essentially unchanged from m00 at depths of 10.0 and 40.0 km.
However, the absolute values of Vs were higher andlower in the fast and slow areas, respectively. Therefore, the boundaries between the low- and high-speed

areas in m16 were better defined than those in m00.
At a depth of 5.0 km, the total trend of the new Vsvalue was changed in the eastern part of the model.
The change appeared to be independent of the wave
speed in m00. The Vp value of m16 was lower thanthat of m00 throughout the region, and a low-wave
speed area appeared along the western coast of the
Boso Peninsula (Fig. 11b). Figure 12 shows the deviation of the wave speeds of m00 and m16 from theaverage speeds at depths of 5.0, 10.0, and 40.0 km.
Two areas with extremely low-Vs were detected in
m16 at a depth of 5.0 km (Fig. 12a). The first was beneath the Boso Peninsula (A in Fig. 12a), and theMiyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 12 of 20

Fig. 10 Results of seismic wave speed changes. a Average wave speed of Vs and Vp concerning depth. Black and red lines indicate m00 and m16,
respectively. b Comparison of Vs values of m00 and m16. c Comparison of Vp values of m00 and m16

Fig. 11 Depth slices concerning the seismic wave speed model. a Vs and b Vp models. The left, center, and right panels show the wave speed of
m00, the wave speed of m16, and the difference between the wave speeds of m00 and m16. Open squares and red triangles represent F-net
stations and volcanoes, respectively. In (a), active faults and the MTL are shown as black lines and a broken line, respectively

Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 13 of 20

Fig. 12 Depth slices of the perturbation of the seismic wave speed from the average speed. a Vs and b Vp models. The left and right panels
indicate the perturbations of the wave speeds for m00 and m16, respectively. Areas labeled A to E in the right figures indicate low-wave speed
anomalies. Other aspects of the figure are the same as in Fig. 11

second was NW of the first (B in Fig. 12a). The firstarea was not detected clearly in m00 and appeared tobe contiguous to the low-Vs area at a depth of
10 km. At a depth 40 km, low seismic wave speed
patterns were observed more clearly in m16 than in
m00 and were distinct in the E–W (C in Fig. 12a)and N–S (D and E in Fig. 12a) directions. Regarding
the Vp deviation shown in Fig. 12b, low-wave speeds
were obtained in the eastern part of the model region, whereas the western part of the model region
showed high-wave speeds at depths of 5.0 and
10.0 km. The deviation of Vp from the average had asimilar pattern to that of Vs at a depth of 40 km. The

following section discusses the tectonic significance of
these extremely low-wave speed areas.

Discussion
Validation of proposed model

One of the purposes of the model developed in thepresent study is to predict the waveforms of earthquakesthat occur beneath the Kanto region. To verify the proposed model, 18 additional event data that were notconsidered in the inversion were examined. Additional
forward simulations concerning m00 and m16 were conducted for the 18 extra events. The full-waveform fittingof the synthetic waveforms to the observation data isMiyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 14 of 20

Fig. 13 Waveform fitting for extra events. a 5–30 s, b 8–30 s, c 10–30 s, and d 20–30 s. The gray bars and black outlines indicate the VR values of
m00 and m16, respectively. The average VR values are also shown in each panel

shown in Fig. 13. The final model m16 better reproduced
the observation than the initial model m00, and the synthetic waveforms in the longer period range showed better fitting than those in the shorter period range. The
VR value in the period range of 5–30 s was improved by
47%, corresponding to change from − 0.97 to − 0.50 for
the m00 and m16 waveforms, respectively. Figure 14
shows a comparison of the observed and synthetic waveforms for an earthquake occurring in the upper crust
and a vertical cross section of the shear wave speed
structure. The transverse and vertical components of the
S-waves for m16 were better fit than those for m00 to theobserved waveform. The proposed model did not reproduce the observed waveforms after body waves, such as
surface waves derived from the basin (Fig. 14) because
no windows were set around the surface waves, whichindicated that the wave speed model could still be improved with future research.
Figure 15 shows the event profiles of the vertical components of waveforms at different epicentral distances
for the events used in the inversion and an extra eventthat was not used in the inversion. The synthetic waveforms for m16 were better fit than those for m00 to theobserved waveforms except for the waveforms after Swave arrivals (Fig. 15b). The proposed model still overestimates the amplitudes of the local surface waves. However, the inversion appeared to have accurately
reproduced the body waves. These results indicated thatthe proposed model had great potential to model the
waveforms of past and future earthquakes and that theredetermination of source parameters using the time
shift was an adequate technique for correcting the centroid time in the inversion. Furthermore, this improved
model was obtained using the waveforms of 140 sources
recorded at only 16 broadband stations.
Data coverage and resolution

The parameters of the seismic wave speed were revised
at all grid points throughout the model region. It is wellknown that there is only waveform sensitivity aroundthe ray paths between sources and receivers. The misfit
kernels were calculated using the theoretical m16 waveforms and deviation waveforms of m16 to determine thedata coverage in the present study. A 1.0 s phase shiftanomaly concerning the theoretical waveforms of m16
was included, and the misfit kernels were calculated byassuming the observed waveforms were perturbed waveforms. The obtained misfit kernels of Vp and Vs wereregarded to represent the data coverage of Vp and Vs, respectively. The data coverage results are shown in Fig. 16.
The kernels were normalized by a subjective value of
1.5 × 10−11. Excluding the coverage of Vp at a depth of
40 km (Fig. 16), the data used in the present study covered nearly the entire region shown in Figs. 11 and 12.
As discussed in previous waveform inversion studies,
traditional checkerboard tests with a workable resolution
had a high computational cost (e.g., Zhu et al. 2015). In
the present study, a point-spreading function using the
Hessian and model perturbation was considered (e.g.,
Fichtner and Trampert 2011; Zhu et al. 2015). The

Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 15 of 20

Fig. 14 Same as Fig. 8 for the wave speed models and waveforms for the extra shallow event (Feb. 19, 2012; MW 4.9). a Source and receiver.
b Vertical cross section of Vs for m00 and m16. c Vertical cross section of the difference between the Vs values for m00 and m16. d Waveformcomparisons. The gray areas in (b) and (c) indicate areas of low data coverage. The crosses and inverted triangles in (b) and (c) show the locationof the earthquake and station in (a), respectively. The black, blue, and red lines in (d) indicate the observed, m00, and m16 waveforms, respectively

gradient of the misfit for model m + δm perturbed frommodel m was considered to examine the model resolution. The gradient of the misfit is given byg ðm þ δmÞ≈g ðmÞ þ H ðmÞδm:

ð7Þ

If model m is the preferred model, the gradient g(m)nearly equals zero, and Eq. (7) can be rewritten asg ðm þ δmÞ≈H ðmÞδm:

ð8Þ

According to Eq. (8), we must access the full Hessian
H(m) or obtain the gradient of the misfit function on

the model m + δm to know the resolution. In this study,∼HðmÞ was used to avoid additional computational costs.
This approach was expected to provide a first approximate guess of the spreading. Future work is needed to
conduct a more thorough resolution analysis (e.g.,
Fichtner 2010; Zhu et al. 2015). Although the gradientof the misfit generally represented the pattern of the difference between the old and new models, the action∼HðmÞδm of the Hessian with respect to the model perturbation was considered as a substitute for g(m + δm)based on Eq. (8). Two cases of the point-spreading test
were conducted based on Eq. (8) in the area of interest,Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 16 of 20

Fig. 15 Seismogram profiles of the data and extra events. a Vertical waveforms for the earthquake used in the inversion (June 2, 2007; MW 4.5). b
Vertical waveforms for the extra earthquake not used in the inversion (Feb. 19, 2012; MW 4.9). Black traces indicate observed waveforms, and blue
and red traces indicate m00 and m16 waveforms, respectively. The absolute maximum amplitudes (units: 10−6 m) are shown in the right panels.
Black curves indicate the theoretical travel time curves of the P- and S-waves obtained based on a 1D model (Kubo et al. 2002). The CMT solutions and stations are also shownas shown in Fig. 17. The perturbation to the final model
m16 was applied at depths of 40 and 5 km in the lowwave speed area. In both of these test cases, the perturbation patterns remained unchanged and indicated thatthe proposed model was well constrained. However, it
must be considered that the position of the maximum∼HðmÞδm shifted to the surface in Fig. 17 because of thelarge values of the Hessian near the surface. Further∼more, the absolute value of the maximum HðmÞδm israther small (Fig. 17a, b). Although the wide sensitivityof the misfit kernels required numerous iterations to
obtain a sufficiently accurate model, the Hessian wasapplied as a preconditioning operator in the present
study to avoid slow recovery.

Tectonic interpretation of the proposed modelThis subsection focuses on the shear wave structure
shown in Fig. 12a and presents a tectonic interpretationof the proposed model. A low-wave speed zone in theE–W direction was observed at a depth of 40 km on
m16 (C in Fig. 12a) and identified as the Vp structure on
m00 and m16 (Fig. 12b). In the proposed model, thisanomalously low-wave speed was lower and clearer thanthat in m00. As reported in previous studies, the lowwave speed anomaly around a depth of 40 km was considered to be serpentinized peridotite in the mantle
above the subducted PHS plate beneath the Kanto region (e.g., Kamiya and Kobayashi 2000; Matsubara et al.
2005; Nakajima et al. 2009). Additional low-wave speed

Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 17 of 20

Fig. 16 Data coverage concerning seismic wave speed. The left and right figures indicate the coverage of Vs and Vp, respectively. The kernels
were normalized by a subjective value of 1.5 × 10−11. Open squares and red triangles indicate F-net stations and volcanoes, respectively

anomalies were identified in the northern part of the
Boso Peninsula (D in Fig. 12a) and along the volcano
chain at a depth of 40 km (E in Fig. 12a). Nakajima et al.
(2009) interpreted anomaly D as the serpentinizedmantle of the PHS slab above the PAC slab, which was
associated with lack of interplate earthquakes on the
PAC slab. Anomaly E has also been described inprevious studies (e.g., Matsubara et al. 2005; Nakajima etal. 2009) and was interpreted as the result of the dehydration reaction caused by the subduction of the PAC
plate. The anomaly may be associated with volcanic generation and activity.
The low-wave speed anomaly areas at depths of 5.0
and 10.0 km (A and B in Fig. 12a) mentioned in aMiyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Page 18 of 20

Fig. 17 Point-spreading test. a Model perturbation of Gaussian distribution at a depth of 40 km. b Response H(m)δm of (a). c Model perturbation∼
at a depth of 10 km. d Response H(m)δm of (c). HðmÞ≈H ðmÞ was assumed in this test

previous section are discussed here. The shear wave
speed at anomaly A was approximately 10% lower than
the average speed beneath the Boso Peninsula. This
trend was identified in the horizontal slices of the compressional wave speed (Fig. 12b). The area of anomaly A
roughly matched the structure of the Kanto basin (e.g.,
Suzuki 1996) and the upper part of the subducting PHS
plate (e.g., Ishida 1992). Anomaly B extended in the
NW–SE direction along the active faults. It was continuous to the north, and its eastern edge met a high-wave
speed anomaly near the TSK station. These low-wave
speed features were in good agreement with the topography of the lowland (Fig. 1). Anomaly B was also contiguous to the northwestern part, and its strike appeared
to change to the NE–SW direction in the northern areaof the Izu Peninsula. This implied that it was related to
the lateral bending of the Median Tectonic Line andzonal geologic terrains associated with the collision tectonics of the opening of the Japan Sea in the Middle
Miocene and the migration of the PHS plate that occurred until 6 Ma (Takahashi and Saito 1997). Although
the travel time tomography results have not revealed
such an anomaly at a shallow depth, adjoint tomographyusing waveforms demonstrated the characteristics of ananomaly associated with the surface geology. The
present results still have low resolution because theperiod range of waveforms must be longer than thoseused to study the case of California (e.g., Tape et al.
2010) to provide a detailed structure. A detailed modelwith high resolution could be obtained by using moredata and seismograms of shorter period range to contribute to strong motion predictions by reproducing entire waveforms, including surface waves generated fromthe basin structure.

Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

Conclusions
A 3D seismic wave speed model capable of reproducing
observed seismograms using adjoint tomography beneath the Japanese metropolitan area in the Kanto region was successfully developed. The optimized SEM
code on the K computer required a total of 6720 simulations using approximately 62,000 node hours to obtainthe final model through 16 iterations. The fit of the synthetic waveforms to the observed waveforms was improved using the proposed model, by 24 and 47%
relative to the fit of the original model with full waves in
the period range of 5–30 s for inversion data and extra
data, respectively. These results indicated that the proposed model had the potential to evaluate waveforms of
past and future earthquakes reliably and that source relocation using the time shift was an adequate techniquefor correcting centroid time in the inversion. Additionally, this improved model was obtained using the waveforms of 140 sources recorded at only 16 broadband
stations. The wave speeds Vs and Vp of the proposedmodel were changed by approximately 20 and 10%, respectively, from those of the initial model. The proposedmodel showed strong heterogeneity compared with theinitial model at each depth slice. Areas with anomalously
low-wave speeds were detected, and detection was
supported by the analysis of data coverage and point
spreading. In combination with the results of previous
studies, our results indicated that the detected extreme
anomalies were located at the bottom of the Kanto basin
on top of the PHS slab at a depth of 5 km, at theserpentinized mantle above the PHS slab at a depth of
40 km, and at the partial melting beneath the volcanoes
at a depth of 40 km.
Abbreviations
FREESIA: Fundamental Research on Earthquakes and Earth’s Interior Anomaly;
NIED: National Research Institute for Earth Science and Disaster Resilience;
PAC: Pacific; PHS: Philippine Sea; SEM: Spectral-element method
Acknowledgements
We appreciate the valuable discussions we had with Kazuto Ando, Daisuke
Suetsugu, and Nozomu Takeuchi. We used the seismic waveform andearthquake catalog provided by the NIED F-net. Digital map 50 m gridelevation data and the J-EGG500 were provided by the Geospatial Information
Authority of Japan and the Japan Oceanographic Data Center, respectively. The
large computing of the forward and adjoint simulations was conducted on the
K computer at the RIKEN Advanced Institute for Computational Science in Kobe,
Japan. We used the open source program package SPECFEM3D Cartesian and
FLEXWIN from the Computational Infrastructure for Geodynamics (CIG;
geodynamics.org) in this study. GMT software (Wessel and Smith 1998) wasused to make most of the figures in this paper.
Funding
This work was supported by the HPCI System Research project (Proposal
number hp130013), JSPS KAKENHI Grant Number 16K21699, and MEXT
KAKENHI Grant Number 15H05832.
Authors’ contributions
TM designed this study, conducted most of the computing and analysis,
and wrote the manuscript. MO and DP assisted with the data processing
and large computing. YT took part in the analysis of the observed data.

Page 19 of 20

ST assisted with the large computing. All authors read and approved thefinal manuscript.
Competing interests
The authors declare that they have no competing interests.

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.
Author details
1
Earthquake Research Institute, University of Tokyo, 1-1-1 Yayoi, Bunkyo-ku,
Tokyo 113-0032, Japan. 2Present address: Earthquake and Tsunami Research
Division, National Research Institute for Earth Science and Disaster Resilience,
3-1 Tennodai, Tsukuba, Ibaraki 305-0006, Japan. 3Department of Deep Earth
Structure and Dynamics Research, Japan Agency for Marine-Earth Scienceand Technology, 2-15 Natsushima-cho, Yokosuka, Kanagawa 237-0061, Japan.
4
Physical Sciences and Engineering Division, King Abdullah University of
Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia.
5
Research and Development Bureau, Ministry of Education, Culture, Sports,
Science and Technology, 3-2-2 Kasumigaseki, Chiyoda-ku, Tokyo 100-8959,
Japan. 6Center for Earth Information Science and Technology, Japan Agency
for Marine-Earth Science and Technology, 3173-25 Showa-machi,
Kanazawa-ku, Yokohama, Kanagawa 236-0001, Japan.
Received: 21 March 2017 Accepted: 21 September 2017

References
Bozdag E, Peter D, Lefebvre M, Komatitsch D, Tromp J, Hill J, Podhorszki N,
Pugmire D (2016) Global adjoint tomography: first-generation model.
Geophys J Int 207:1739–1766. doi:10.1093/gji/ggw356
Brocher T (2005) Empirical relations between elastic wavespeeds and density inthe Earth’s crust. Bull Seism Soc Am 95:2081–2092
Chen M, Niu F, Liu Q, Tromp J, Zheng X (2015) Multiparameter adjoint
tomography of the crust and upper mantle beneath East Asia: 1. Modelconstruction and comparisons. J Geophys Res Solid Earth 120:1762–1786.
doi:10.1002/2014JB011638
Clayton R, Engquist B (1977) Absorbing boundary conditions for acoustic andelastic wave equations. Bull Seism Soc Am 67:1529–1540
Fichtner A (2010) Full seismic waveform modelling and inversion. SpringerVerlag, Heidelberg
Fichtner A, Kennett B, Igel H, Bunge H (2009) Full seismic waveform tomography for
upper-mantle structure in the Australasian region using adjoint methods.
Geophys J Int 179:1703–1725. doi:10.1111/j.1365-246X.2009.04368.x
Fichtner A, Trampert J (2011) Resolution analysis in full waveform inversion. Geophys
J Int 187:1604–1624. doi:10.1111/j.1365-246X.2011.05218.x
Gardner GHF, Gardner LW, Gregory AR (1974) Formation velocity anddensity—the diagnostic basics for stratigraphic traps. Geophys 39:770–780
Geller R, Hara T (1993) Two efficient algorithms for iterative linearized inversion ofseismic waveform data. Geophys J Int 115:699–710. doi:10.1111/j.1365-246X.
1993.tb01488.x
Hori S (2006) Seismic activity associated with the subducting motion of the Philippine
Sea plate beneath the Kanto district, Japan. Tectonophysics 417:85–100
Ishida M (1992) Geometry and relative motion of the Philippine Sea Plate andPacific Plate beneath the Kanto-Tokai District, Japan. J Geophys Res
97:489–513. doi:10.1029/91JB02567
Kamiya S, Kobayashi Y (2000) Seismological evidence for the existence ofserpentinized wedge mantle. Geophysical Res Lett 27:819–822
Kamiya S, Kobayashi Y (2007) Thickness variation of the descending Philippine
Sea slab and its relationship to volcanism beneath the Kanto-Tokai district,central Japan. J Geophys Res 112:B06302. doi:10.1029/2005JB004219
Kimura T, Murakami H, Matsumoto T (2015) Systematic monitoring of
instrumentation health in high-density broadband seismic networks. Earth
Planet Space 67:55. doi:10.1186/s40623-015-0226-y
Koketsu K, Miyake H, Afnimar TY (2009) A proposal for a standard procedure ofmodeling 3-D velocity structures and its application to the Tokyo
metropolitan area, Japan. Tectonophysics 472:290–300
Koketsu K, Miyake H, Fujiwara H, Hashimoto T (2008) Progress towards a Japan
integrated velocity structure model and long-period ground motion hazard

Miyoshi et al. Progress in Earth and Planetary Science (2017) 4:29

map. Proceedings of the 14th world conference on earthquake engineering,pp S10–038
Komatitsch D, Tromp J (2002a) Spectral-element simulations of global seismic
wave propagation-I. Validation. Geophys J Int 149:390–412
Komatitsch D, Tromp J (2002b) Spectral-element simulations of global seismic
wave propagation-II. Three-dimensional models, oceans, rotation and selfgravitation. Geophys J Int 150:303–318
Kubo A, Fukuyama E, Kawai H, Nonomura K (2002) NIED seismic moment tensor
catalogue for regional earthquakes around Japan: quality test andapplication. Tectonophysics 356:23–48
Maggi A, Tape C, Chen M, Chao D, Tromp J (2009) An automated time-window
selection algorithm for seismic tomography. Geophys J Int 178:257–281.
doi:10.1111/j.1365-246X.2009.04099.x
Matsubara M, Hayashi H, Obara K, Kasahara K (2005) Low-velocity oceanic crust at
the top of the Philippine Sea and Pacific plates beneath the Kanto region,
central Japan, imaged by seismic tomography. J Geophys Res 110:B12304.
doi:10.1029/2005JB003673
Matsubara M, Obara K (2011) The 2011 off the Pacific coast of Tohoku
Earthquake related to a strong velocity gradient with the Pacific plate. Earth
Planet Space 63:663–667
Matsumura M, Ito Y, Kimura H, Obara K, Sekiguchi S, Hori S, Kasahara K (2006)
Development of accurate and quick analysis system for source parameters(AQUA). Zishin 59:167–184 (in Japanese with English abstract)
Montagner J, Anderson D (1989) Petrological constraints on seismic anisotropy.
Phys Earth Planet Inter 54:82–105
Nakajima J, Hasegawa A (2010) Cause of M ∼ 7 intraslab earthquakes beneath
the Tokyo metropolitan area, Japan: possible evidence for a vertical tear at
the easternmost portion of the Philippine Sea slab. J Geophys Res 115:
B04301. doi:10.1029/2009JB006863
Nakajima J, Hirose F, Hasegawa A (2009) Seismotectonics beneath the Tokyo
metropolitan area, Japan: effect of slab-slab contact and overlap on
seismicity. J Geophys Res 114:B08309. doi:10.1029/2008JB006101
Okada Y, Kasahara K, Hori S, Obara K, Sekiguchi S, Fujiwara H, Yamamoto A (2004)
Recent progress of seismic observation networks in Japan -Hi-net, F-net, KNET and KiK-net-. Earth Planet Space 56:xv-xxviii
Olsen KB, Day SM, Bradley CR (2003) Estimation of Q for long-period (>2 sec)
waves in the Los Angeles basin. Bull Seismol Soc Am 93:627–638
Panning M, Romanowicz B (2006) A three-dimensional radially anisotropic modelof shear velocity in the whole mantle. Geophys J Int 167:361–379
Patera AT (1984) A spectral element method for fluid dynamics: laminar flow in achannel expansion. J Comput Phys 54:468–488
Peter D, Komatitsch D, Luo Y, Martin R, Goff N, Casarotti E, Loher P, Magnoni F, Liu Q,
Blitz C, Nissen-Meyer T, Basini P, Tromp J (2011) Forward and adjoint simulationsof seismic wave propagation on fully unstructured hexahedral meshes. Geophys J
Int 186:721–739. doi:10.1111/j.1365-246X.2011.05044.x
Plessix RE (2006) A review of the adjoint-state method for computing thegradient of a functional with geophysical applications. Geophys J Int 167:
495–503. doi:10.1111/j.1365-246X.2006.02978.x
Simutė S, Steptoe H, Cobden L, Gokhverg A, Fichtner A (2016) Full-waveform
inversion of the Japanese Islands region. J Geophys Res Solid Earth 121:
3722–3741. doi:10.1002/2016JB012802
Suzuki H (1996) Geo1ogy of the Koto deep borehole observatory and geo1ogicalstructure beneath the metropolitan area, Japan. Rep Natl Res Inst Earth Sci
Disaster Prev 56:77–123 (in Japanese with English abstract)
Takahashi M, Saito K (1997) Miocene intra-arc bending at an arc-arc collisionzone, central Japan. Island Arc 6:168–182
Tape C, Liu Q, Maggi A, Tromp J (2009) Adjoint tomography of the southern
California crust. Science 325:988–992. doi:10.1126/science.1175298
Tape C, Liu Q, Maggi A, Tromp J (2010) Seismic tomography of the southern
California crust based on spectral-element and adjoint methods. Geophys J
Int 180:433–462. doi:10.1111/j.1365-246X.2009.04429.x
Tape C, Liu Q, Tromp J (2007) Finite-frequency tomography using adjoint
methods—Methodology and examples using membrane surface waves.
Geophys J Int 168:1105–1129. doi:10.1111/j.1365-246X.2006.03191.x
Tarantola A (1984) Inversion of seismic reflection data in the acoustic
approximation. Geophysics 49:1259–1266
Tarantola A (1986) A strategy for nonlinear elastic inversion of seismic reflection
data. Geophysics 51:1893–1903
Tromp J, Tape C, Liu Q (2005) Seismic tomography, adjoint methods, time
reversal and banana-doughnut kernels. Geophys J Int 160:195–216.
doi:10.1111/j.1365-246X.2004.02453.x

Page 20 of 20

Tsuboi S, Ando K, Miyoshi T, Peter D, Komatitsch D, Tromp J (2016) A 1.8 trillion degreesof-freedom, 1.24 petaflops global seismic wave simulation on the K computer. Int JHigh Perform Comput Appl 30:411–422. doi:10.1177/1094342016632596
Wang Y, Chevrot S, Monteiller V, Komatitsch D, Mouthereau F, Manatschal G,
Sylvander M, Diaz J, Ruiz M, Grimaud F, Benahmed S, Pauchet H, Martin R(2016) The deep roots of the western Pyrenees revealed by full waveform
inversion of teleseismic P waves. Geology 44:475–478. doi:10.1130/G37812.1
Wessel P, Smith WHF (1998) New, improved version of the Generic Mapping
Tools released. EOS Trans AGU 79:579. doi:10.1029/90EO00426
Yoshimoto K, Takemura S (2014) Surface wave excitation at the northern edge ofthe Kanto Basin, Japan. Earth Planet Space 66. doi:10.1186/1880-5981-66-16
Zhou Y, Dahlen F, Nolet G (2004) Three-dimensional sensitivity kernels for surface wave
observables. Geophys J Int 158:142–168. doi:10.1111/j.1365-246X.2004.02324.x
Zhu H, Bozdag E, Peter D, Tromp J (2012) Structure of the European upper
mantle revealed by adjoint tomography. Nat Geosci 5:493–498.
doi:10.1038/NGEO1501
Zhu H, Bozdag E, Tromp J (2015) Seismic structure of the European upper mantle
based on adjoint tomography. Geophys J Int 201:18–52. doi:10.1093/gji/ggu492
Zhu H, Tromp J (2013) Mapping tectonic deformation in the crust and upper
mantle beneath Europe and the north Atlantic Ocean. Science 341:871–875.
doi:10.1126/science.1241335