Kubota et al. Progress in Earth and Planetary Science
https://doi.org/10.1186/s40645-019-0313-y(2019) 6:67

Progress in Earth andPlanetary Science

RESEARCH ARTICLE

Open Access

Fault model of the 2012 doublet
earthquake, near the up-dip end of the2011 Tohoku-Oki earthquake, based on anear-field tsunami: implications for
intraplate stress state
Tatsuya Kubota1,2* , Ryota Hino2, Daisuke Inazu3 and Syuichi Suzuki2

Abstract
On December 7, 2012, an earthquake occurred within the Pacific Plate near the Japan Trench, which wascomposed of deep reverse- and shallow normal-faulting subevents (Mw 7.2 and 7.1, respectively) with a timeinterval of ~10 s. It had been known that the stress state within the plate was characterized by shallow tensile anddeep horizontal compressional stresses due to the bending of the plate (bending stress). This study estimates thefault model of the doublet earthquake utilizing tsunami, teleseismic, and aftershock data and discusses the stressstate within the incoming plate and spatiotemporal changes seen in it after the 2011 Tohoku-Oki earthquake. We
obtained the vertical extents of the fault planes of deep and shallow subevents as ~45–70 km and ~5 (theseafloor)–35 km, respectively. The down-dip edge of the shallow normal-faulting seismic zone (~30–35 km)
deepened significantly compared to what it was in 2007 (~25 km). However, a quantitative comparison of thebrittle strength and bending stress suggested that the change in stress after the Tohoku-Oki earthquake was too
small to deepen the down-dip end of the seismicity by ~10 km. To explain the seismicity that occurred at a depthof ~30–35 km, the frictional coefficient in the normal-faulting depth range required would have had to be ~0.07 ≤μ ≤ ~0.2, which is significantly smaller than the typical friction coefficient. This suggests the infiltration of pore fluid
along the bending faults, down to ~30–35 km. It is considered that the plate had already yielded to a depth of ~35
km before 2011 and that the seismicity of the area was reactivated by the increase in stress from the Tohoku-Oki
earthquake.
Keywords: Ocean bottom pressure gauge, Doublet earthquake, Intraplate earthquake, 2011 Tohoku-Oki earthquake,
Bending stress, Fault modeling

Introduction
It is well known that the stress state within the incoming
Pacific Plate near the Japan Trench is characterized byshallow tensile and deep horizontal compressional stresses
along a direction perpendicular to the trench axis, separated by a thin aseismic (i.e., stress-neutral) “elastic core”,* Correspondence: kubotatsu@bosai.go.jp
1
National Research Institute for Earth Science and Disaster Resilience, 3-1
Tennodai, Tsukuba, Ibaraki 305-0006, Japan
2
Graduate School of Science, Tohoku University, 6-6 Aza-Aoba, Aramaki,
Aoba-ku, Sendai, Miyagi 980-8578, Japan
Full list of author information is available at the end of the article

due to the bending of the plate (bending stress, Fig. 1c;e.g., Chapple and Forsyth 1979). It is also well known thatthe number of normal-faulting earthquakes occurring
within the plate and near the trench axis increases, following interplate megathrust earthquakes; this is attributed tothe increased horizontal tensile stress caused by the stressrelease of the interplate coupling (e.g., Christensen and
Ruff 1988; Dmowska and Lovison 1988).
Recently, Craig et al. (2014) investigated the verticalvariation in the centroid depth and fault mechanisms
based on global catalogs (M > ~5) and noted that temporal changes in the transition depths between normal-© The Author(s). 2019 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.

(2019) 6:67

141˚

142˚

Page 2 of 20

143˚

144˚

145˚

(a)

TPG
KPG
DART
GPS Buoy

39˚

KAMN
P06
P02

KAMS

G09

GJT3
P01 P04 P08
P05
TJT1
P10
P03 MYGI
KPG2
JFAST
A
P07
KCTD KPG1
P12
JP1
805

38˚

807
804
802
803

JP2

0

(b)

-2

GFK

-4

B

801

-6

806
37˚

0 km

50 km

-8

(b)
GCMT #1

Bathymetry [km]Kubota et al. Progress in Earth and Planetary Science

Harada et al. #1

38˚

GCMT #2(c) 108oE from North

0

Harada et al. #2

Tensile37.6˚Lay et al. #1

Stress-neutral

Lay et al. #2

Compressional

1947
0143.8˚

144˚

10

20

144.2˚

30

A(WNW)

144.4˚ 0

20
40

Depth [km]

37.8˚

60

B(ESE)
50

100

80
150

Horizontal Distance [km]Fig. 1 Location map of this study. a Green inverted triangles, yellow squares, pink triangles, and blue diamonds represent TPG, GPS buoy, KPG,
and DART stations, respectively. Small dots denote the aftershocks deduced from the ocean bottom seismographs (Obana et al. 2015). Gray
contour lines indicate the coseismic slip distribution of the Tohoku-Oki earthquake (Iinuma et al. 2012) with 10 m intervals. Black dashed line
represents the trench axis. b Enlarged map around the focal area. Black star shows the epicenter by Japan Meteorological Agency. Red, blue, andgreen CMT solutions are from GCMT, Lay et al. (2013) and Harada et al. (2013), respectively. Gray CMT solution denotes the 1967 mb 4.7
earthquake (Seno and Gonzalez 1987). c The vertical cross-section along the A-B line in a. Aftershocks within the dashed rectangle in a areshown. Black curved line is the plate boundary (Ito et al. 2005). Schematic image of the intraplate bending stress state is also shown (red: downdip compressional stress, blue: extensional stress). The thick line denotes the approximate location of the stress neutral planeand reverse-faulting earthquakes were not detected after
the 2011 Tohoku-Oki earthquake. However, Obana et al.
(2012, 2014, 2015, 2019) studied ocean bottom seismographs and reported that the down-dip limit of shallow
normal-faulting earthquakes (M < ~5) within the platenear the trench axis deepened from ~25 to ~35 km after
the same earthquake; they interpreted spatiotemporalchange in the intraplate stress as having been caused by
the release of stress associated with the Tohoku-Oki
earthquake. However, the reason for this inconsistency
has not been clarified. Furthermore, the relationship between the stress state in the incoming Pacific Plate andthe changes in coseismic stress due to the Tohoku-Oki

earthquake have not yet been quantitatively assessed indetail.
On December 7, 2012, an Mjma 7.3 earthquake occurred
within the Pacific Plate near the Japan Trench, where theextremely large coseismic slip (> ~50 m) was estimated tohave occurred during the 2011 Tohoku-Oki earthquake(e.g., Iinuma et al. 2012; star in Fig. 1b). Detailed teleseismic analyses (Lay et al. 2013; Harada et al. 2013) revealed
that this earthquake was composed of two M ~7 subevents.
According to the global centroid moment tensor (GCMT)
solution (http://www.globalcmt.org/; Ekström et al. 2012),
the first subevent had a reverse-faulting mechanism with adepth of ~60 km (Mw 7.2) and the second had a normal-

(2019) 6:67

Page 3 of 20

faulting mechanism (~20 km, Mw 7.2) with a time intervalof ~12 s (red CMT solutions in Fig. 1b and c). Hereafter,
this earthquake is referred to as the doublet earthquake,and the first and the second subevents are referred to as
subevent 1 and subevent 2, respectively. Since the fault
mechanisms of the two subevents are consistent withbending stress, the source process of the doublet earthquake should reflect the intraplate stress state after the
Tohoku-Oki earthquake. The vertical extents of each fault
will be key to discussing the temporal change in the vertical variations of the stress state after the 2011 Tohoku-Oki
earthquake.
Rupture processes related to the doublet earthquake
have been investigated previously. Lay et al. (2013) andHarada et al. (2013) investigated this earthquake using
teleseismic data to estimate the CMT solution and thefinite fault model. Teleseismic data is generally a powerful dataset for resolving rupture processes of global
earthquakes. However, since the teleseismic signals fromeach subevent overlapped, it is difficult to decompose
the rupture process of the doublet earthquake precisely,
especially for the latter, shallower, subevent.
Inazu and Saito (2014) estimated the spatial distribution of the initial sea-surface height change (tsunami
source) using far-field tsunami data from ~200–2000 km
away from the focal area (Additional file 1: Figure S1). In
contrast to teleseismic data, tsunami data is useful for
constraining the rupture process of the subevent 2, sinceshallow earthquakes generally excite tsunamis or cause
seafloor vertical deformation more effectively than deep
earthquakes.
When the 2012 doublet earthquake occurred, off-line
autonomous absolute ocean bottom pressure gauges (PGs)
installed near the focal area (< ~200 km from the source,
Fig. 1a) recorded clear tsunami signals. This dataset is useful to constrain the fault model of subevent 2, which wasdifficult to constrain with teleseismic data. In the present
study, we utilize the tsunami and aftershock data and theresults of the teleseismic analysis to estimate the finite fault
model of the 2012 doublet earthquake, focusing particularlyon the vertical extent of the fault planes of each subevent.
We also discuss the relationship between the vertical profileof the intraplate seismicity and its spatiotemporal changes
associated with the 2011 Tohoku-Oki earthquake.

Airport Research Institute (PARI) (Kato et al. 2005) (yellow
squares), and by the Deep-ocean Assessment and Reportingof Tsunamis (DART) system (Bernard et al. 2014) (blue diamonds). Detailed information is given in Table 1.
To retrieve the tsunami waveforms, we remove theocean-tide component using the theoretical tide model(Matsumoto et al. 2000) and apply the filter from Saito(1978). The lowpass filter we applied has a cutoff of 3
min to the TPG records, and the bandpass filter has apassband of 3–60 min to the KPG, GPS buoy, andDART records. We apply the lowpass filter to the TPG
records to preserve the offset in the pressure change
caused by the vertical deformation of the seafloor.
The filtered records are shown in Fig. 2. TPGs first capture down-motion tsunamis with amplitudes of ~−5 cm
and then larger up-motion tsunamis with amplitudes of~+10 cm. The durations of both the down- and up-motion
tsunamis are ~5 min (black dashed lines in Fig. 2a). Small
fluctuations and changes in the pressure offset are observedat the stations near the focal area (Fig. 2b). Tsunami amplitudes at the DART and the KPG stations are very small (~1cm, black dashed lines in Fig. 2c and d). At the former, dynamic pressure changes caused by seismic waves (e.g.,
Kubota et al. 2017b) are also observed. Tsunami signals arealso detected by some GPS buoys (e.g., ~15 cm at station
801, Fig. 2e).

Kubota et al. Progress in Earth and Planetary Science

Methods/Experimental
Tsunami data

We use the near-field PGs installed by Tohoku University(hereafter TPG; e.g., Hino et al. 2014; Kubota et al. 2017a,
2017b) (green inverted triangles in Fig. 1a). We also use
tsunami data obtained by the off-Kushiro online cabled
PGs installed by Japan Agency for Marine-Earth Scienceand Technology (JAMSTEC) (Hirata et al. 2002) (KPGs,
pink triangles), by GPS buoys installed by the Port andStep-by-step approach for fault modelingTo decompose the complex rupture process of the 2012
doublet earthquake, we apply a step-by-step procedure totsunami, teleseismic, and aftershock data. We first estimate an initial sea-surface height distribution of the tsunami (hereafter, the tsunami source model) by inverting
tsunami records. Since seafloor crustal deformations, or
tsunamis, are very sensitive to shallow earthquakes, wethen estimate the fault model of shallow subevent 2 basedon the tsunami source model. We also use aftershocks
detected around subevent 2 (Obana et al. 2015) to obtain
information on the fault geometry. We then calculate thechange in residual sea-surface height between the tsunami
source model and the vertical displacement from subevent
2. Given that this residual distribution is caused by subevent 1, we then estimate the fault model of subevent 1.
Furthermore, because there are large trade-offs between
fault size and the focal depth and the amount of slip in a
deeper earthquake, we also use the results of teleseismic
analyses (Lay et al. 2013; Harada et al. 2013) to obtainprior information on the fault parameters such as fault
depth and size for subevent 1.
Tsunami source modeling using near-field tsunami
records

We estimate the tsunami source model by inverting tsunami
records via the conventional inversion analysis method (e.g.,Kubota et al. Progress in Earth and Planetary Science(2019) 6:67

Page 4 of 20

Table 1 List of tsunami stations used in this study
StationLatitude[°N]

Longitude[°E]

Depth [m]Epicentral distance [km]Inversion time window [s]Instrument

Sampling rate oforiginal data [s]a

GFK

37.5812142.7647

2245

140

0– 200

TPG

1

P01

38.3331142.4167

1038

180

300–1500

TPG

1

P02

38.5006

142.5035

1109

180

300–1500

TPG

1

P03

38.1834

142.3996

1056

170

300–1500

TPG

1

P04

38.3163

142.5657

1265

160

300–1500

TPG

1

P05

38.3000

142.7004

1412

150

300–1500

TPG

1

P06

38.6338

142.5833

1269

180

300–1500

TPG

1

P07

38.0000

142.4486

1064

170

300–1500

TPG

1

P08

38.2833

142.8329

1424

140

150–1350

TPG

1

P10

38.2500

143.1666

2066

110

0–1200

TPG

1

P12

37.8206

142.8996

1635

130

0–1200

TPG

1

KAMN

38.8862

143.3639

2360

150

0–1200

TPG

1

KAMS

38.6347

143.2621

2246

130

0–1200

TPG

1

MYGI

38.0832

142.9166

1697

130

0–1200

TPG

1

G09

38.4782

143.7922

5500

90

Not used

TPG

1

GJT3

38.2948

143.4811

3260

90

Not used

TPG

1

JFAST

37.9336143.91546799

40

Not used

TPG

1

TJT1

38.2080

143.79045744

60

Not used

TPG

1

801

38.2325141.6836

144

240

Not used

GPS Buoy1

802

39.2586

142.0969

204

250

Not used

GPS Buoy1803

38.8578141.8944

160

240

Not used

GPS Buoy1804

39.6272

142.1867

200

320

Not used

GPS Buoy1805

40.6333

141.7500

87

380

Not used

GPS Buoy1806

36.9714141.1856

137

250

Not used

GPS Buoy1807

40.1167

142.0667

125

320

Not used

GPS Buoy1

KPG1

41.7040

144.4375

2218

450

1200–2400

KPG

1

KPG2

42.2365

144.8454

2210

510

1500–2700

KPG

1

KCTD

41.6675

144.3409

2540

440

1200–2400

KPG

10

JP1

40.3777

146.1681

5125

330

1050–2250

DART

15

JP2

39.2849

145.7845

5183

210

750–1950

DART

15

a

Observed records were resampled to 15 s in the inversion for the tsunami source model

Tsushima et al. 2012; Inazu and Saito 2014; Kubota et al.
2018a). The details of the procedure are identical to thosedescribed in Kubota et al. (2018a). Before estimating the tsunami source model, however, we conduct a preparatory tsunami simulation using the tsunami source model in Inazu
and Saito (2014) (Additional file 1: Figure S1). Although thepressure fluctuation at the TPGs was explained roughly, thepressure offset changes at the stations closest to the source(G09, TJT1, JFAST, and GJT3) were not explained at all(Additional file 1: Figure S1). This indicates that the latter is
not due to coseismic seafloor deformations. Thus, they
are probably due to tilts or rotations in the sensors
related to the seafloor strong ground motion (e.g.,
Wallace et al. 2016; Kubota et al. 2018a).

In calculating the Green’s functions for the tsunami,
the unit source elements of the seafloor displacements(Kubota et al. 2015, 2018a) are distributed around the
focal area. The horizontal dimension of the unit source
elements is 20 km × 20 km at a spacing of 10 km (overlapping with the adjacent elements), distributed along
260 km (in the EW direction) × 240 km (NS) area. To
calculate the sea-surface displacement from the unit
source elements of the seafloor displacement, we consider the spatial filtering effect due to water depth (Saito2019). In the depth filtering process, we assume a seawater depth of 6 km. To simulate the tsunami, we solve
the linear dispersive wave equation in the local Cartesian
coordinates (e.g., Saito 2019). The grid spacing is 2 km

Kubota et al. Progress in Earth and Planetary Science(2019) 6:67

Page 5 of 20

6

140

120

Height [cm]

100

80

60

40

20

(a) TPG
KAMS
KAMN
MYGI
P12
P10
P08
P07
P06
P05
P04
P03
P02
P01
GFK

4

2

6

4

2

10hPa ~ 10cm ΔH2O

0
100

80

(b) TPG near the epicenter
60

Height [cm]

60

20

JP1
JP2
1cm

0

0

40

(c) DART

G09
TJT1
JFAST
GJT3

(d) KPG
KPG2
KPG1
KCTD
1cm
(e) GPS Buoy
805
807
804
802
803

40

10cm
10hPa ~ 10cm ΔH2O

20

801
806

0

0
-30

0

30

60

90

Time from EQ [min]120Time from EQ [min]Fig. 2 Filtered tsunami records. The records for a TPGs, b TPGs near the epicenter, c DART, d KPGs, and e GPS buoys are shown. Light gray anddark gray waveforms are the de-tided waveforms and 60s-moving-averaged records, respectively. Red waveforms are the filtered records

with a 1 s time step according to interpolation done via
ETOPO1 bathymetry data (Amante and Eakins 2009).
We assume that the displacements of all unit sources occurred instantaneously and simultaneously. The tsunami
propagation velocity expected by the linear long-wave
theory is expressed as v = (g0H)1/2 (g0: gravity acceleration constant, H0: water depth). Given the assumed
water depth of 6 km and the average depth of the focalarea, the propagation velocity is approximated as ~240m/s (~15 km/min). Thus, the tsunami propagation distance during the duration of the M~7 earthquake (~10s) and the time interval between two subevents (~10 s)
is about 3 km, which is sufficiently small compared tothe extent of the tsunami source model (~100 km, Inazu
and Saito (2014)). We consider the static pressure offsets
related to the calculation of Green’s function of the PGs
for permanent seafloor deformation (Tsushima et al.
2012; Kubota et al. (2018a)). The same filter used for theobserved records is applied to the simulated waveforms.

All data is resampled to 15 s intervals for the inversion. We use different time windows for each station, including tsunami main phase (Table 1, thick black lines
in Fig. 3). The smoothing constraint is imposed and its
weight is determined based on the trade-off between theweight and reduction in variance between the observed
and simulated tsunami waveforms (Additional file 1: Figure S2). Since the G09, TJT1, JFAST, and GJT3 stations
are located near the source and are probably affected byseafloor ground shaking (Additional file 1: Figure S1),
we exclude these records from the inversion analysis.
Since all TPG stations are located landward of thesource region, the constraint for the eastern edge of thetsunami source is likely not very good. Therefore, toimprove the source constraint further, we also use the
DART and KPG records. Since the amplitudes of TPG
records are approximately ten times larger than the
DART and KPG records, we weight the KPG andDART data at values ten times that of the TPG data.

Kubota et al. Progress in Earth and Planetary Science(2019) 6:67

Page 6 of 20

(a)
39˚

KAMN
P06
P02

KAMS

G09
P04P08GJT3
P01
P05 P10
TJT1
P03 MYGI
JFAST
P07
P12

38˚

[cm]
15
10
5
0

KPG2
KPG1
KCTD
805 JP1807
JP2802 804
803
37˚
801
806

GFK

142˚

143˚

-5
-10

0 km

50 km

-15

Contour: 5 cm
144˚

145˚

(b)
140
120

Height [cm]

100
80
60
40
20

Height [cm]

0
60
40
20

TPG

KAMS
KAMN
MYGI
P12
P10
P08
P07
P06
P05
P04
P03
P02
P01
GFK

Observed
Calculated

Used for Inversion

6
5
4
3
2

DART
JP1
JP2

1
0
4
3
2
1

KPG

KPG2
KPG1
KCTD

0
80

TPG near the epicenter
G09
TJT1
JFAST
GJT3

0

60
40
20

GPS Buoy

805
807
804
802
803
801
806

0
-10

0

10

20

30

40

Time from EQ [min]

-10

0

10

20

30

40

50

Time from EQ [min]Fig. 3 Results of tsunami source inversion. a Distribution of the tsunami source. Contour line interval is 5 cm. b Comparison of the observed(gray) and synthesized (red) waveforms. Thick black lines denote the time window used for the inversion

Fault modeling of subevent 2

Because we found that the subsidence of the tsunami
source was generated by the shallow subevent 2 (see Fig.
4 and Tsunami source model subsection in Results andDiscussion section), we first estimate a fault model for
subevent 2, which best explains the subsidence region of
the tsunami source model. We use the grid-search approach proposed by Kubota et al. (2015, 2018b), which
estimates an optimum rectangular planar fault modelwith uniform slip. Because the short-wavelengthcomponent disappeared in the tsunami source model dueto the smoothing constraint imposed in the inversion andthe spatial smoothing effect used in the deep-sea region
during the tsunami generation (see Fig. 4 and Tsunami
source model subsection in Results and Discussionsection), we consider the smoothing effect in fault modeling by the following procedure. First, we calculate theseafloor deformation using a fault model candidate (a setof unknown parameters) (Okada 1992) and then simulate
a tsunami. The simulated waveforms are inverted toKubota et al. Progress in Earth and Planetary Science(2019) 6:67(a) Tsunami source model

Page 7 of 20(b) GCMT #1+#2

38˚

0 km

50 km

0 km

50 km

37˚

(c) GCMT #1[cm](d) GCMT #2

15
10

38˚

5
0
-5
-10
0 km

37˚
143˚

144˚

50 km
145˚ 143˚

0 km

50 km

144˚

-15

145˚

Fig. 4 Sea-surface displacement expected from the GCMT solution. Displacements from a the tsunami source model, b the combination of both
subevents, c subevent 1, and d subevent 2. The faults are represented by black rectangles

obtain the initial sea-surface height distribution, under thesame conditions used in the inversion for the tsunami
source model. Finally, we evaluate the goodness of the tsunami waveform fitting from the fault model candidate, bycomparing the subsided area of the tsunami source modeland the inverted sea-surface height.
The geometry of the fault plane is assumed to be on
the west-dipping plane of GCMT solution (strike = 189°
and dip = 50°), which is consistent with the planar structure of the aftershock (Obana et al. 2015). We also assume the rake angle from the GCMT solution (= −90°).
Because the aftershock alignment is located ~2 km west
from the GCMT centroid, the fault plane is constrained
to pass through the point that is 2 km west from the
GCMT centroid (hereafter, referred to as the reference
point). The unknown parameters are the distance fromthe reference point to both ends of the fault, along the
strike (i.e., length) and dip (width) direction (L1, L2, W1,
and W2; see Fig. 5c and d). The search range for theseparameters is determined based on the aftershocks andthe horizontal extent of the subsidence area of the tsunami source. The fault length (L = L1 + L2) is assumed
to be greater than the fault width (W = W1 + W2), as L
> W. The top of the fault plane (defined by parameterW2) is constrained as not to extend above the seafloor.

The fault model candidate is assessed through variancereduction (VR) of the subsided areas between the tsunami source model and the fault model candidate:2 !
PN  source−Dxcandidatei
i¼1 xiVR ¼ 1−
PN  source 2
k¼1 xi 100ð%Þ;ð1Þ
and xcandidate
are the displacements of thewhere xsourceiisea-surface at the ith grid point, from the tsunami
source model and the fault model candidate assuming
the unit slip, respectively. N is the total number of grid
points. The slip amount on fault D is determined so that
the VR takes the maximal value. We use the grid points
within the subsided area of the tsunami source model(blue dashed line in Fig. 5a) to calculate the VR. The
search range of the unknown parameters is listed inTable 2.
Fault modeling of subevent 1

In order to estimate the fault model of subevent 1, we
use the residual sea-surface height between the tsunami
source model and the sea-surface displacement expectedfrom the fault model of subevent 2 (hereafter, referred(2019) 6:67

Page 8 of 20

(b)

100

(a)
B
VR [%]

38˚

A
A’

90

7.35

80

7.30

70

7.25

60

7.20

50

7.15

40

7.10

30

7.05

20

B’
0 km

7.00

VR
Mw

10

50 km

7.40

6.95
6.90

0

37˚
143˚

144˚

1

145˚

Mw

Kubota et al. Progress in Earth and Planetary Science

10

100

1000

Order
0

(c)

(d)
W2
th
id
ta

lw

20
30

To

Depth [km]10

Total length
L2

L1

W1

Fault strike

40
50

A(West)-40

B(North)A’(East)
-20

0

-40

20

-20

B’(South)
0

20

40

Horizontal Distance [km]Horizontal Distance [km]Fig. 5 Results for the grid-search of subevent 2. a The optimum fault model is shown by the red rectangle. The top ten fault model candidates
are shown by thin black rectangles. b VR (black) and Mw (red) of the fault model candidates searched in the grid-search, arranged in thedescending order in terms of VR. The rectangle with blue dashed lines denotes the area used for calculating VR. c, d Vertical profile along A-A’
and B-B’ lines in a. Aftershocks are taken from the rectangular area with the green lines in a. The configuration of the fault parameters is
also shown

Table 2 Search range of the grid-search for fault modeling
Reference point latitude [°N]

Subevent 2

Subevent 1

37.77c37.6–38.1d

Reference point longitude [°E]

143.81143.8–144.5d

Reference point depth [km]a

19.5

57.8

Strike [°]

189

158

Dip [°]a

50

59
48

a

c

Rake [°]−90

L1 [km]b0, 5, 10, ..., 50e

15

L2 [km]e0, 5, 10, ..., 50

15

W1 [km]b

10, 15, 20, ..., 40e

10

a

b

W2 [km]10, 15

Slip amount [m]bAdjusted so that the VR values become maximal

b

a

e

10

Strike, dip, rake angles, and reference point depth were fixed to the GCMT value
Rectangular fault model with uniform slip is assumed. The fault length (L = L1 + L2) was assumed to be greater than the fault width (W = W1 + W2)
A point 2 km west of the GCMT centroiddHorizontal location of reference point was searched within the range of ± 30 km in the EW and NS direction from the centroid location of Lay et al. (2013)(37.82°N, 144.13°E), with increments of 5 kmeIncrements for the fault length and width are 5 kmbc(2019) 6:67

Page 9 of 20

to as the residual height distribution). This is because
the residual height distribution is expected to correspond to the sea-surface displacement due to subevent 1.
Because subevent 1 occurred at the deeper part of theincoming plate, it appears to be difficult to constrain thefault parameters such as fault geometry, depth, size, andslip amount only from the residual height distribution.
Meanwhile, the residual height will contribute to constrain the horizontal location of the fault. Thus, we use
the results of the teleseismic analyses in previous studies(GCMT; Lay et al. 2013; Harada et al. 2013) to obtainprior information on the fault parameters. We fix thecentroid depth, fault geometry (strike, dip, and rake angles), and fault dimension (length and width) based on
previous teleseismic analyses. We estimate the optimum
horizontal location (longitude and latitude) of the fault
and the amount of uniform slip on the rectangular fault
using the grid-search approach, as in the fault modelingof subevent 2.
Because Lay et al. (2013) investigated the CMT solution of the 2012 doublet earthquake using the teleseismic W-phase waveforms and showed and discussed theuncertainty of their estimations in detail, we use the
CMT solution proposed by them as prior informationfor our fault modeling. We assume the centroid depth as
60.5 km and the west-dipping nodal plane with geometry of strike = 163° dip = 51°, and rake = 57°. The slip
distribution of subevent 1 obtained by the teleseismic
analysis (Lay et al. 2013; Harada et al. 2013) had a main
rupture area with a dimension of L ~ 30 km and W ~ 20km; therefore, we fix the fault length and width as 30
and 20 km, respectively (L1 = L2 = 15 km and W1 = W2
= 10 km). The search range of the horizontal location of
the centroid is determined based on the evaluation of
the uncertainty of the horizontal location of the Wphase analysis of Lay et al. (2013) (see Figure S2 in Lay
et al. 2013). The search range is listed in Table 2.

displacement from the GCMT solution and compared
the results with the tsunami source model. We use theequation in Okada (1992), which assumes that the rectangular planar fault on the west-dipping nodal plane has auniform slip. For simplicity, values of length, width, andslip L = 70 km, W = 35 km, and D = 0.7 m, respectively,
are used. The spatial pattern of the tsunami source model(Fig. 4a) is similar to that of the combined deformation ofsubevents 1 and 2 (Fig. 4b) and the subsidence area issimilar to that expected from the CMT solution of subevent 2 (Fig. 4d). The uplift of the tsunami source model is
not consistent with either subevent 1 nor 2 alone (Fig. 4cand d). Based on this comparison, we conclude that thesubsidence is generated by subevent 2 alone and that both
subevents contribute to the uplift. The deformation expected from subevent 2 (Fig. 4d) has a sharp displacementpeak that was not estimated in the tsunami source model.
This is probably because the short-wavelength components disappeared due to the smoothing constraint imposed in the inversion and the spatial smoothing effectused in the deep-sea region during the tsunami generation(Saito 2019).

Kubota et al. Progress in Earth and Planetary Science

Results and Discussion
Tsunami source model

We obtained a tsunami source distribution that had apair of large uplifts and subsidences (Fig. 3a). The observed waveforms (red lines in Fig. 3b) were reproduced well. The GPS buoy waveforms, which were
not used for the inversion, were also explained. Furthermore, although the offset changes at the stations
near the source (G09, TJT1, JFAST, and GJT3) were
not reproduced, the fluctuations in the calculated
waveforms were similar to the observations. This isconsistent with the idea that these changes were dueto the tilting or rotation of the sensors (e.g., Wallace
et al. 2016; Kubota et al. 2018a).
To investigate the contribution from each subevent on
the tsunami source, we compute the sea-surface vertical

Fault models of two subevents of the 2012 doublet
earthquake

We obtained an optimum fault model of subevent 2,
which had the highest VR value among all the fault
model candidates, with a length of 35 km (L1 = 5 km, L2= 30 km), width of 25 km (W1 = 15 km, W2 = 10 km),and fault slip of 1.1 m (VR = 96 %, Fig. 5). The seismic
moment was 5.5 × 1019 Nm (Mw 7.1, assuming the rigidity of 60 GPa). The VR values for all fault model candidates (in the descending order) are shown in Fig. 5b.
The VR values for the top ten candidates are almost flatand relatively high (Fig. 5d). This indicates that the top
ten candidates reasonably reconstruct the subsidence ofthe tsunami source model. Hence, we inspect the model
parameters for these candidates to evaluate the estimation error of the optimum fault model. They are projected onto the vertical cross-section in Fig. 5c and d(thin black lines), and the histograms of the model parameters are shown in Fig. 6. Among the top ten candidates,
most models have the down-dip limit of the fault plane of
W1 = 15 km, and all models have fault bottoms shallower
than the depth of ~35 km (W1 ≤ 20 km, Fig. 6). This indicates that the lower end of the fault of subevent 2 shouldbe less than ~35 km.
We also calculated the sea-surface height assuming
the faults with smaller dimension with larger slip
(Additional file 1: Figure S3). In this calculation, thedown-dip end of the fault (W1) and slip amount werechanged, and the length (L 1, L2) and the up-dip end(W 2) of the fault, and seismic moment Mo were fixed
to those of the optimum model. In the fault models

(2019) 6:67

Kubota et al. Progress in Earth and Planetary Science

Count

10

Top 10 models
W2 (Up-dip)

Page 10 of 20

with W < 15 km (down-dip depth of fault is shallower
than ~25 km), the locations of the western edge ofthe subsided area and of the peak displacement areinconsistent with those of the tsunami source model.
This indicates the small fault models are implausible.
The subsidence of the tsunami source model could be
explained when the small faults are located slightly
west of the optimum fault location. However, such
faults can be rejected because we used the aftershock
distribution of Obana et al. (2014, 2015) to constrain
the horizontal location of the fault. It is important to
use the aftershock distribution for prior informationon the fault horizontal location, in order to accurately
constrain the down-dip depth of the subevent 2 fault.
The subsided area calculated from the optimum
fault model was consistent with that of the tsunami
source model (Fig. 7a and b). We calculate the residual
sea-surface height between the tsunami source model(Fig. 7a) and the sea-surface displacement expectedfrom the fault model of subevent 2 (Fig. 7b; the residual height distribution).
From the residual height distribution utilizing the results of the teleseismic analyses, we constrained the fault

L2 (South)

5
0
10

Count

10

15

0

W1 (Down-dip)

10 20 30 40 50

L1 (North)

5
0
10

Count

10

20

30

40

Width

0

10 20 30 40 50

Length

5
0
10

20

30

40

10

20

[km]

30

40

50

[km]

Fig. 6 Histograms of the top ten model parameters in the gridsearch of subevent 2(a) Tsunami source (=#1+#2)(b) Subevent 2 model

38˚

0 km

50 km

0 km

50 km

37˚

[cm]

(d) Subevent 1 model(c) Residual distribution
Harada et al.

15

GCMT

10

38˚

5
0
-5
-10
Lay et al.

37˚
143˚

0 km
144˚

0 km

50 km
145˚ 143˚

144˚

50 km

-15

145˚

Fig. 7 Sea-surface height changes from the models in the present study. a The tsunami source model with contour intervals of 5 cm. The
rectangular area enclosed by the blue dashed lines denotes the area used for the fault modeling of subevent 2. b Subevent 2 fault model; theoptimum fault is also shown. c the residual between a and b. CMT solutions for subevent 1 by GCMT, Lay et al. (2013) and Harada et al. (2013)
are also shown. The rectangular area enclosed by the blue dashed lines is the area used for the fault modeling of subevent 1. d the subevent 1
fault model; the optimum fault is also shown(2019) 6:67

Page 11 of 20

model of subevent 1. The result is shown in Fig. 8. The
slip amount on the fault was 2.2 m, and the seismic moment was 7.8 × 1019 Nm (Mw 7.2, rigidity of 60 GPa). We
obtained the optimum VR of 80%. The VR values are relatively high and flat for the top ten fault model candidates(Fig. 8b). The uncertainty of the horizontal location of thefault center for these candidates is likely to be ± ~10 km(Fig. 8a). This horizontal uncertainty is almost consistent
with that estimated by Lay et al. (2013).

We also simulate the teleseismic P-waves using the fault
model parameters for comparison with the observed teleseismic waveforms (Fig. 11). We assume pure doublecouple point sources at the centers of the optimum faultsof each subevent. We use the calculation programs of
Kikuchi and Kanamori (2003). A triangular-shaped source
time function with the rise time of 6 s is assumed, considering the typical rupture duration of M ~7 earthquakes(Fig. 9b). After simulating the waveforms of each subevent, we stack the simulated waveforms. We assume thatthe difference of focal times between subevents 1 and 2
was 10 s, which is determined by inspecting the waveform
similarity of the observed and stacked waveforms.
We use a 1D multi-layered velocity structure model
without the water layer, assuming that the source structure was identical to the receiver structure in Table 3 (bluetraces in Fig. 11b). The simulated waveforms for eachsubevent are also shown in Additional file 1:Figure S4.
The peak timing and amplitudes of the first up-motion
and the subsequent down-motion waves reasonably fitted
the observation, although the subsequent phases during
50–80 s did not perfectly match. This is probably becauseof the assumptions of the velocity structure and simple
source time function. We then simulate the teleseismic
waveforms incorporating the water layer and oceanic
structure, shown in Table 3 (red traces in Fig. 11b). The
agreement of the subsequent phases improved comparedwith the simulation without the water layer. According to
the teleseismic analysis by Lay et al. (2013), anothersmaller normal-faulting subevent was estimated at 40 s

Tsunami and teleseismic waveform simulation from
optimum fault model

The results of the fault modeling of the 2012 doublet earthquake are summarized in Fig. 9. The sea-surface height displacement expected from the optimum rectangular fault
models of subevents 1 and 2 (Fig. 9a) is calculated by thesuperposition of the displacements from each fault model(Figs. 5a and 8a). The distribution is very similar to that of
the tsunami source model (Fig. 4). The optimum models
have vertical ranges of ~6 (seafloor)–30 km for shallow subevent 2 and ~50–70 km for deep subevent 1 (Fig. 9d and e).
From the combined displacement, we numericallysimulate tsunami waveforms (Fig. 10). The simulated
waveforms reasonably explained the observed tsunami
waveforms well, not only the near-field TPGs but also
the far-field DART, KPG, and GPS buoys (red lines
in Fig. 10b). We calculate the VR between the observed and calculated waveforms from Eq. (1), by
using the same time window as in the inversion oftsunami records. We finally obtained VR = 79 %.

(b)

100

(a)

VR [%]

38˚

Search range
for fault location

90

7.35

80

7.30

70

7.25

60

7.20

50

7.15

40

7.10

30

7.05

20
0 km
37˚
143˚

50 km

7.40

7.00

VR
Mw

10

6.95
6.90

0
144˚

145˚

Mw

Kubota et al. Progress in Earth and Planetary Science

1

10

100

Order
Fig. 8 Results for the grid-search of subevent 1. a The optimum fault model is shown by the red rectangle. The top ten fault candidates areshown by thin black rectangles. A rectangle with dot-and-dashed lines denotes the search range for the fault centroid location. b VR (black) andMw (red) of the fault model candidates searched in the grid-search, arranged in the descending order in terms of VR

(2019) 6:67

Page 12 of 20

2

(b)
(a)

1

B’

0

38˚

-5

0 km
37˚
143˚

50 km

144˚

15

20

25

30

(c)
Along dip [km]

15
10
5
0
-5
-10
-15

A’

10

Time [s]

[cm]

B

5

0

11.8

10

19.5

20

27.2

30

Distance ≤ 5km
-10

145˚

0

Depth [km]A

0Moment rate [×1019 Nm/s]

Kubota et al. Progress in Earth and Planetary Science

34.8

10 20 30 40

Along strike [km]
0

(d)

(e)

Depth [km]

10
20
30
40
50
60
70

A
-60

A’
-40

-20

0

20

40

Horizontal Distance [km]

60

B
-60

B’
-40

-20

0

20

40

Horizontal Distance [km]Fig. 9 Summary of the fault modeling. a Sea-surface height distribution calculated from the superposition of the those from the fault models of
subevents 1 and 2. Locations of the optimum faults are also shown. b The moment rate function used for the teleseismic waveform calculation. c
Aftershock distribution on the fault plane of subevent 2. Aftershocks with distances from the fault plane smaller than 5 km are shown. d, eVertical profile along A-A’ and B-B’ lines in a. Orange and blue rectangles are the optimum fault models of subevents 1 and 2, respectively.
Aftershocks are taken from the rectangular areas in aafter subevent 1. Hence, it is possible that the third smaller
subevent also contributed to the generation of the later arrival. It is worth pointing out that it is important to use
the teleseismic records to resolve the temporal complexityof the doublet earthquake in detail.
Our fault model explains both tsunami and teleseismic
observations. In addition, our fault model of subevent 2
is consistent with the aftershock distribution determinedby the ocean bottom seismographs. The rupture area ofsubevent 2 estimated by Harada et al. (2013) was locatedat the outer-trench region and concentrated in the shallower portion of the plate (z < ~20 km). The horizontal location of subevent 2 centroid by Lay et al. (2013) (Fig. 1b)
was also inconsistent with our fault model and withthe aftershock locations. The consistency of our fault
model with the tsunami, teleseismic waveforms, andaftershocks indicates that the step-by-step procedureused in this analysis can decompose the complex rupture process of the 2012 doublet earthquake. We conclude that we can obtain a more comprehensive fault
model of the 2012 doublet earthquake, than the one
estimated from the teleseismic data alone.

Importance of near-field tsunami data for the fault
modeling

Inazu and Saito (2014) showed a tsunami source modelof the 2012 doublet earthquake using offshore tsunami
stations located more than 200 km from the source area
(Additional file 1: Figure S1). We compare the tsunami
source models of this study with that of Inazu and Saito(2014) (Fig. 12a). The horizontal location of the tsunami
source was in agreement with the model of Inazu andSaito (2014), although the amplitudes were lower. The
simulated tsunami at the TPG stations using the tsunami
source model of Inazu and Saito (2014) is similar to thetsunami peak timing of the observation (Fig. 12b). This
indicates that the horizontal location of the tsunami
source is reasonably constrained even when using tsunami stations located far from the source (> ~200 km).
However, the dominant period and amplitudes of thesimulated waveforms are longer and smaller than theobservation. This indicates that the spatial resolution of
the far-field tsunami data is not sufficient for the finite
fault model, and thus, to constrain the down-dip limit of
the fault plane of subevent 2. Using the near-fieldKubota et al. Progress in Earth and Planetary Science

140
120

Height [cm]

100
80
60
40
20

Height [cm]

0
60
40
20

(2019) 6:67

Page 13 of 20

TPG

KAMS
KAMN
MYGI
P12
P10
P08
P07
P06
P05
P04
P03
P02
P01
GFK

Observed

6
5
4
3
2

Calculated

DART
JP1
JP2

1
0
4
3
2
1

KPG

KPG2
KPG1
KCTD

0
80

TPG near the epicenter
G09
TJT1
JFAST
GJT3

0

60
40
20

GPS Buoy

805
807
804
802
803
801
806

0
-10

0

10

20

30

40

Time from EQ [min]

-10

0

10

20

30

40

50

Time from EQ [min]Fig. 10 Comparison of the observed (gray) and calculated (red) waveforms based on the superposition of the optimum fault models of
subevents 1 and 2

tsunami data, the constraint of the down-dip limit of thefault plane of shallow subevent 2 is improved. This enables us to discuss the intraplate stress regime after the2011 Tohoku-Oki earthquake.
The fault dimension and depth of subevent 1 estimatedin this study are consistent with previous results of the teleseismic waveforms (Lay et al. 2013; Harada et al. 2013).
This is because we used the teleseismic analyses as the priorinformation for fault dimension. In contrast, the horizontal
location and fault dimension of subevent 2 differ from theteleseismic analyses. The down-dip limit of subevent 2 constrained from the teleseismic analysis (Harada et al. 2013) is
considerably shallower (< ~20 km) than that of the constrained from tsunami data (~40 km). Our fault model has
better consistency with the aftershock distribution, for bothhorizontal location and depth range. It is likely that the rupture process of subevent 2 was not resolved well from theteleseismic data, because the teleseismic signal from subevent 2 overlaps with that of subevent 1. By using the nearfield tsunami records, the resolution of subevent 2 rupture
was considerably improved.
Temporal change in down-dip limit of normal-faulting
earthquake

The normal-faulting aftershocks around subevent 2(Obana et al. 2014, 2015) mainly occurred at depths of~30 km, which corresponds to the down-dip depth ofthe optimum fault of subevent 2 (Fig. 9b). Some normalfaulting seismicity also occurred at depths of ~35 km,
along the down-dip direction (Fig. 9b). According to
Obana et al. (2014), the estimation error of the aftershocks is less than 5 km (Figure 2 in Obana et al., 2014).
Considering the estimation error of the fault model andaftershock distribution, the down-dip limit where theshallow normal-faulting seismicity can occur around thisregion is ~30–35 km (Fig. 13).
We also investigate the temporal change of the intraplate seismicity before and after the Tohoku-Oki earthquake (Fig. 13). The down-dip limit of subevent 2 fault(~35 km) is clearly ~10 km deeper, compared with thedown-dip limit of the normal-faulting seismicity observed in 2007 (< ~25 km, Hino et al. 2009). Because
both researches use the arrays of the ocean bottom seismometers installed just above the focal area, which haveidentical sensitivity and were distributed with almost
identical spatial intervals (~10 km), the detectability inboth observation periods is expected to be identical andthus the difference between the seismicity depths of thelower limit of the shallow normal-faulting seismicity wasconfidentially significant. In addition, a few deeper (>~40 km) events were detected in both observations. This
also suggests the misdetection of the deeper shallow
normal-faulting event (~30–40 km) before the TohokuOki earthquake is unlikely to occur.

(2019) 6:67

Kubota et al. Progress in Earth and Planetary Science

(a)

BORG
KIEV
KIV

SFJD

HRV
SSPA
WCI

KEV

FFC
BILL

ABKT
UOSS
FURI

Page 14 of 20

COLA RSSD HKTPFO SLBS
POHA

PALK

XMAS
FUNA

MSEY
DGAR COCO WRAB

RAO

NWAO

(b)

SNZO

Obs.
IU.WCI.00

IU.SFJD.00

Syn. (without water layer)
IU.XMAS.10
II.COCO.00

Syn. (with water layer)
II.ABKT.00

Az=5o Dist=74o Az=37o Dist=91o Az=108o Dist=64o Az=231o Dist=66o Az=299o Dist=65o
IU.BILL.00
IU.FUNA.00
II.KIV.00
IU.RSSD.00
II.DGAR.00
Az=15o Dist=32o Az=43o Dist=77o Az=137o Dist=56o Az=252o Dist=80o Az=311o Dist=71o
IU.HRV.00
IU.HKT.00
II.PALK.10
IU.KIEV.00
IU.RAO.10
Az=25o Dist=93o Az=48o Dist=92o Az=146o Dist=75o Az=259o Dist=64o Az=323o Dist=74o
IU.KEV.00
IU.SSPA.00
II.PFO.00
II.MSEY.00
IU.SNZO.00
Az=30o Dist=92o Az=57o Dist=76o Az=157o Dist=83o Az=265o Dist=91o Az=339o Dist=63o
II.FURI.00
II.WRAB.00
IU.COLA.00
IU.SLBS.00
II.BORG.00
Az=32o Dist=47o Az=61o Dist=87o Az=190o Dist=58o Az=286o Dist=96o Az=353o Dist=77o
II.FFC.00
IU.POHA.00
II.UOSS.00
IU.NWAO.00
Az=33o Dist=71o
0

30

60

90

Az=91o Dist=54o Az=203o Dist=74o Az=288o Dist=73o
0

30

60

90

0

30

60

90

0

30

60

90

Time from P-arrival [s]Fig. 11 Results of the forward calculation of the teleseismic waveforms based on the optimum fault model. a Locations of the teleseismic
stations. b Comparison of the teleseismic waveforms. Black lines are the observed waveforms, and the synthetic waveforms using the velocity
structure with and without the water layer are shown by red and blue lines, respectively. The bandpass filter of 10–500 s is applied

Table 3 Velocity structure used for the teleseismic calculationa
Structure# of layerVp [km/s]

Vp [km/s]

ρ [g/cm3]

H [km]

Source

3

1.50

0.00

1.00

6.0

6.00

3.50

2.70

6.0

8.10

4.70

3.30

Half space

6.00

3.50

2.70

18.0

6.75

3.80

2.80

18.0

8.10

4.70

3.30

Half space

Receivera3

The structure is based on Kikuchi and Kanamori (1991) but the water layer(thickness of 6 km) is assumed for the source structure

In contrast, the up-dip limits of the subevent 1 fault(~50 km) and the deep reverse-faulting seismicity (~50km, Obana et al. 2015) are almost equivalent to the deep
reverse-faulting seismicity before the Tohoku-Oki earthquake (Seno and Gonzalez 1987; Hino et al. 2009), although it is difficult to discuss this in detail because of thevery low seismicity. By focusing on the shallow normalfaulting seismicity, we discuss the cause of the deepeningof the down-dip limit of the normal-faulting seismicity.
The yield strength of the plate is characterized by thebrittle rupture at the shallow portion and ductile failurelaws at the deeper parts of the plate (Fig. 14a, e.g.,
Scholtz 1988; Turcotte and Schubert 2002; Hunter andWatts 2016). Based on the Anderson theory of faulting,the brittle strength along the horizontal direction normal

Kubota et al. Progress in Earth and Planetary Science(2019) 6:67

Page 15 of 20

[cm]

(a)

15
10
5
0
-5
-10
-15

38˚

(b)

(b) 37.8oN
Inazu and Saito (2014)
This study

0 km
37˚
143˚

144˚

(c)

140
120

Height [cm]

100
80
60
40
20

Height [cm]

0
60
40
20

50 km
145˚

TPG

KAMS
KAMN
MYGI
P12
P10
P08
P07
P06
P05
P04
P03
P02
P01
GFK

143.5

144.0

[cm]
30
20
10
0
-10
-20
-30

144.5

Observed
This study
Inazu and Saito (2014)
5
4

DART

JP1

3
2

JP2

1
0
4
3
2
1

KPG

KPG2
KPG1
KCTD

0
80

TPG near the epicenter
G09
TJT1
JFAST
GJT3

0

60
40
20

GPS Buoy

805
807
804
802
803
801
806

0
-10

0

10

20

30

40

Time from EQ [min]

-10

0

10

20

30

40

50

Time from EQ [min]Fig. 12 Comparison with the result of Inazu and Saito (2014). a Tsunami source model of Inazu and Saito (2014). b Cross-section of the tsunami
sources from this study (red line) and Inazu and Saito (2014) (blue) at 37.8°N. c Comparison of the observed waveforms (gray) and synthesized
waveforms. Blue dashed and red lines are synthesized using the tsunami source model of the present study and Inazu and Saito(2014), respectively

to the trench axis τxx(z) is expressed as (e.g., Turcotteand Schubert 2002):
2μ ρ0 g 0 z−pw
τ xx ðzÞ ¼;
sin2δ þ μð1− cos2δ Þð2Þwhere ρ0 is the crust density, pw is the pore pressure, z
is depth (downward is positive), δ is the fault dip angle,
and μ is the frictional coefficient. This equation implies
that the rock strength is proportional to depth z (green
line in Fig. 14a). Further, by assuming the plate as a rigid

two-dimensional elastic plane (x- and z-axes are the subducting direction and vertical direction, respectively),
the vertical distribution of the bending stress along thedip direction σxx(z) is approximated as (e.g., Turcotteand Schubert 2002; Craig et al. 2014; Hunter and Watts
2016)
σ xx ðzÞ ¼ −EC
ðz−z0 Þ;1−ν2

ð3Þwhere E is the Young’s modulus, C is the curvature ofplate bending, ν is the Poisson’s ratio, and z0 is the(2019) 6:67

Kubota et al. Progress in Earth and Planetary Science

Page 16 of 20

Fig. 13 Schematic representation of the vertical-temporal change of seismicity in the outer-rise region. Blue and red stars denote the normalfaulting and reverse-faulting small seismicity shown by previous studies (Seno and Gonzalez 1987; Hino et al. 2009; Obana et al. 2012, 2014, 2015).
Note that the aftershocks occurred around the Mw 7.6 outer-rise ~40 min after the Tohoku-Oki earthquake (Obana et al. 2012) are located ~60
km east from the focal area of the 2012 doublet earthquake

normal-faulting earthquakes. In contrast, elastic behavior
is expected in the depth range where the bending stressdoes not exceed strength; the range is termed the elastic
core (e.g., Craig et al. 2014; Hunter and Watts 2016).
The actual deviatoric stress profile within the plate isrepresented by the red solid line in Fig. 14a. The top of(b) μ=0.6
0400

500 Tensile

Depth

[km]

s

B

Δσxx ≈ 20 MPa

ngth

Aseismic (σ < τ)

40

in
end

Stre

μ=0.6

30

tres

gs

c

th

ng

BΔσ ≈ 300 MPa
xx

ss

re

St

Depth

g

in

d
en

20

μ=0.6

ismi

ss

strength 30
40

Depth

tre

cs

DuctilessrestTensile

ly se

20

10

stre

(b,c)

tri

10

via

Be

ntial
Pote

tre

gs

in
nd

200

0

iatricSeismic (σ > τ)

ss

[MPa]
100

0

Seismic (σ > τ)100(c) μ=0.07

Devstreng[MPa]
200 300

De

Brittle0

ess

tric str

Devia

th

Stress, Strength
Compressional
Tensile[km](a) Schematic

μ=0.07

Aseismic(σ < τ)stress-neutral depth (tensile stress is positive, blue line
in Fig. 14a).
At a shallower portion of the plate, where the bending
stress exceeds the brittle strength (the blue backgroundarea in Fig. 14a), the rock cannot remain elastic and thestress is released, or the rock yields, leading to shallow

Fig. 14 Vertical profile of the strength and the bending stress within the incoming plate. a A schematic image of vertical profile. Green, graydashed, and blue lines represent brittle strength, ductile strength, and bending stress, respectively. After the yielding of rock, deviatoric stressoccurring within the plate is expressed by a red line. Blue and red background colors show the areas where shallow normal-faulting and deep
thrust-faulting earthquakes occur, respectively. b Stress profile at the shallow part of the plate assuming the typical frictional strength condition.
Dark red line is the deviatoric stress after the Tohoku-Oki earthquake. Gray hatched area is the region where seismicity becomes active after
adding the stress change. c Stress profile assuming the reduced frictional strength. Dark red line is the deviatoric stress after the Tohoku-Oki
earthquake, calculated by adding a constant static stress change of 20 MPa

(2019) 6:67

Page 17 of 20

the aseismic elastic core, or the bottom of the verticalrange of the normal-faulting seismicity, can be defined
as a depth where the bending stress and the frictional
strength are equal. As the top of the elastic core ispresent at a depth of ~30–35 km and the top of thereverse-faulting seismicity (the bottom of the elasticcore) at ~45 km, the stress-neutral plane is expected to
be located in the depth range between 30–35 and 45
km. This depth range is almost consistent with thedepths where fault mechanisms flip from the shallow
normal-faulting to the deep thrust-faulting mechanisms,
near the trench axis off NE Japan (~25–40 km, Gamage
et al. 2009; Koga et al. 2012).
We compare the vertical profiles of the brittle strengthand bending stress. Assuming E = 80 GPa and ν = 0.25,
and C = 2 × 10-7 m-1 (McNutt and Menard 1982), dσxx/dz= EC/(1−ν2) is ~15 MPa/km. Assuming the hydrostaticpressure condition pw = ρwg0z (ρw = 1030 kg/m3, seawater
density) and ρ0 = 2700 kg/m3 , δ = 50°, and μ = 0.6 (e.g.,
Byerlee 1978), dτxx/dz is ~11 MPa/km as per Eq. (2). In thissituation, assuming z0 = 40 km, τxx(z) and σxx(z) are equalat z ~25 km (Fig. 14b). This is consistent with the downdip limit of normal-faulting seismicity observed before the
Tohoku-Oki earthquake (Hino et al. 2009). However, considering dσxx/dz = ~15 MPa/km, the stress increment of~300 MPa is needed to deepen the top of the elastic corefrom 25 to 35 km (orange arrow in Fig. 14b). This value is
too large compared to the coseismic stress change aroundsubevent 2 expected from the fault model of the
Tohoku-Oki earthquake by Iinuma et al. (2012), Δσxx
~20 MPa. The expected depth change of the down-dip limitof the normal faulting seismicity is only a few km (Fig. 14c).
It must be considered that τxx(z) and σxx(z) are equalat ~30–35 km depth before the Tohoku-Oki earthquake(Fig. 14c). In order that the normal-faulting earthquakes
occur at a depth of 35 km, the brittle strength must be
reduced compared to the typical frictional condition(Fig. 14b). If we assume z0 = 40 ± 5 km (Fig. 13), thefrictional coefficient of the fault around the lower end of
the fault of subevent 2 is μ ~0.07 ± 0.06 so that τxx(z)and σxx(z) are equal at z = 35 km (Fig. 14c). Even
when the top of the elastic core is assumed to be at z
= 30 km, the expected frictional coefficient μ is ~0.2
(Additional file 1: Figure S5). These values are muchsmaller than that of the typical rocks but comparable to that
estimated for the other incoming plate (Craig et al. 2014).
Reduction of friction has often been reported in studies on inland earthquakes (e.g., Yoshida et al. 2018); thisreduction can be attributed to the existence of the pore
fluid (e.g., Bell and Nur 1978). Based on the active seismic survey in the outer-rise region of the Japan trench(e.g., Fujie et al. 2018), the significant seismic wave velocity reduction and high Vp/Vs area were detected at theshallow part of the subducting plate (< 5 km), which areinterpreted as being the results of pore fluid penetration
through the pre-existing bending faults in the shallower
part of the plate (e.g., Peacock 2001). Considering thesestudies, it is suggested that that the strength reduction
within the plate might be related to the pore fluid.
Numerical modeling by Faccenda et al. (2009) demonstrated that the pore fluid can infiltrate the plate as deepas the lower limit of the normal faulting seismicity observed in this study. Cai et al. (2018) also reported theserpentinized mantle wedge associated with water infiltration into the subducting plate (down to ~35 km) at
the Mariana subduction zone.
However, the majority of the seismicity at depths of 30–
35 km is located around the subevent 2’s fault (Fig. 9) andless in the other portions of the plate. This localization ofseismicity suggests that the pore fluid, or the strength reduction, is localized within the plate, as also suggested byFaccenda et al. (2009) and Obana et al. (2019).
The activation of the normal-faulting seismicity at
depths of 25–35 km after the Tohoku-Oki earthquake
can be interpreted as follows: the plate at 25–35 km
depths had yielded before the Tohoku-Oki earthquake,leading to intrinsically aseismic region; a stress increment by the Tohoku-Oki earthquake enhanced the horizontal tensile stress in a broad depth range near the topof the elastic core, which activates seismicity. Less
normal-faulting seismicity at depths of 25–35 km duringobservation from April to June 2007 (Hino et al. 2009)
may be representative of the long-term-averaged deformation. It is expected that the stressing rate due tothe bending deformation is lower near the stress-neutral
depth, where null-deformation is expected, than theshallowest part of the incoming oceanic plate.
On the other hand, if EC/(1−ν2) is small, the static stresschange by the Tohoku-Oki earthquake can contribute todeepening the elastic core or the lower limit of the shallow
normal faulting seismicity. Supposing that the elastic coreis deepened by 10 km due to the coseismic static stresschange Δσxx of ~20 MPa by the Tohoku-Oki earthquake,EC/(1−ν2) must be ~0.5 MPa/km, which is smaller than
the typical elastic condition described above by an orderof magnitude. However, it is unlikely that the Young’s
modulus or the plate curvature are significantly reduced
by an order of magnitude even if supposing the existenceof the pore fluid or estimation error of the curvature.
Thus, this hypothesis seems unlikely. Although it might
be possible that the EC/(1−ν2) is small compared to thatassumed in this study, its contribution to deepening theelastic core associated with the Tohoku-Oki earthquake is
not highly significant.

Kubota et al. Progress in Earth and Planetary Science

Conclusions
In this study, we estimated the fault model of the intraplate doublet earthquake that occurred on December 7,

Kubota et al. Progress in Earth and Planetary Science(2019) 6:67

2012 (subevent 1: a deep reverse-faulting earthquake;subevent 2: a shallow normal-faulting earthquake) strategically utilizing offshore tsunami, aftershocks, and theteleseismic records based on the step-by-step analysis
procedure. First, the initial sea-surface height distribution was estimated by inverting the offshore tsunami records and comparing it with the seafloor deformationfrom the CMT mechanism. It was found that the subsidence and uplift areas were generated by subevent 2 andboth subevents, respectively. Then, the fault model ofeach subevent was estimated based on the initial seasurface height model, using information from previous
studies. As a result, the vertical extent of the fault planeof subevent 2 was obtained as ~5 km (i.e., the seafloor)
to 35 km. Finally, we simulated the tsunami and teleseismic waveforms from the fault model, which explained
the observation well.
We compared the tsunami source model obtained fromthe near-field tsunami data acquired at less than 200 km
from the epicenter and that from the far-field (> 200 km)
data. We found that the horizontal location of the tsunami
source was reasonably constrained, even from the far-fieldtsunami data alone. However, to constrain the finite fault
model in more detail, it is necessary to use the near-fieldtsunami records. We also discussed the stress state within
the plate and its spatiotemporal change after the 2011
Tohoku-Oki earthquake. We found that the down-dip
limit of the shallow normal-faulting earthquakes was obviously deepened compared with that observed in 2007,
from 25 to 35 km. However, comparing the coseismic
stress change by the Tohoku-Oki earthquake and theamount of the bending stress within the plate, the plate
down to ~35 km in depth should have already yielded before the Tohoku-Oki earthquake, and the top of the elastic
core been located at ~35 km. Furthermore, as the bending
stress around the top of the elastic core was much smaller
than the rock strength expected from the empirical relationship, the frictional strength in the range of thenormal-faulting earthquakes is expected to be significantly
reduced. The significant strength reduction of the platesuggests pore fluid infiltration down to ~35 km, along thebending faults.

Supplementary information
Supplementary information accompanies this paper at https://doi.org/10.
1186/s40645-019-0313-y.
Additional file 1: Figure S1. Tsunami source model of Inazu and Saito(2014) and the stations used for their inversion analysis. Comparison
between the observed and simulated tsunami waveforms at the stations
near the source is also shown. Figure S2. Trade-off curve between thesmoothing weight and the VR used for the inversion analysis. Figure S3.
Sea-surface height changes from subevent 2 with fault dimensionschanged. The left-top panel shows the distribution of the tsunami source
model, and others show the sea-surface height distribution of assumed
fault models. Figure S4. Result of the forward calculation of the

Page 18 of 20

teleseismic waveforms based on the optimum fault models of (a)subevent 1 and (b) subevent 2. Figure S5. Vertical profile of the strengthand the bending stress within the incoming plate, assuming the top
depth of elastic core at 30 km.
Additional file 2: Dataset S1. The TPG data used in this study.

Abbreviations
CMT: Centroid moment tensor; DART: Deep-ocean Assessment andReporting of Tsunamis; GCMT: Global Centroid Moment Tensor;
JAMSTEC: Japan Agency for Marine-Earth Science and Technology;
KPG: Pressure Gauge Installed off Kushiro; PARI: Port and Airport ResearchInstitute; PG: Pressure gauge; TPG: Pressure gauge installed by Tohoku
University; VR: Variance reduction

Acknowledgements
We thank the reviewer Thorne Lay and another anonymous reviewer, as wellas the editor, Frantisek Gallovic, for their reviews and insightful comments.
We thank Koichiro Obana for providing aftershock data and Keisuke Yoshidafor providing stress change data due to the Tohoku-Oki earthquake. We
thank Eddie Bernard, Yong Wei, Science Applications International
Corporation (SAIC), and JMA for providing DART data. We also thank JAMSTEC and PARI for providing the KPG and GPS buoy data, respectively. We
also thank Editage (www.editage.jp) for English language editing. The figuresin this paper were prepared using Generic Mapping Tools (Wessel & Smith
1998). TK thanks Koichiro Obana and Tatsuhiko Saito for insightful discussions. Without the constructive comments and enormous encouragementfrom Yusaku Ohta, this paper would not have been published.

Authors' contributions
TK designed the study, analyzed the data, interpreted the results, and drafted
the manuscript. RH contributed to fault modeling, data interpretation, andrevision of the manuscript. DI contributed to tsunami modeling and revising
the manuscript. SS took part in the installation and retrieval of the pressuredata. All authors approved the final manuscript.

Funding
This work was supported by the research project “Research concerning
Interaction between the Tokai, Tonankai and Nankai earthquakes” of the
Ministry of Education, Culture, Sports, Science and Technology, Japan, by
JSPS KAKENHI [Grant Numbers JP20244070, JP26000002, and JP 19K14818],
and by the Sasakawa Scientific Research Grant [2019-2037] of the Japan
Science Society.

Availability of data and materials
TPG data is available in Additional file. KPG data were obtained from theSubmarine Cable Data Center (SCDC), JAMSTEC (http://www.jamstec.go.jp/
scdc/top_e.html). Teleseismic data were downloaded from the DataManagement Center (DMC) of Incorporated Research Institutions for
Seismology (IRIS) (http://ds.iris.edu/ds/nodes/dmc/tools/event/3650366). GPS
buoy data were obtained by a request to PARI. The data from DART, which
is jointly operated by NOAA and JMA, was provided upon request by Eddie
Bernard and Yong Wei from Science Applications International Corporation(SAIC). Aftershock data was provided by Koichiro Obana.

Competing interests
The authors declare that they have no competing interests.
Author details
National Research Institute for Earth Science and Disaster Resilience, 3-1
Tennodai, Tsukuba, Ibaraki 305-0006, Japan. 2Graduate School of Science,
Tohoku University, 6-6 Aza-Aoba, Aramaki, Aoba-ku, Sendai, Miyagi 980-8578,
Japan. 3Department of Marine Resources and Energy, Tokyo University ofMarine Science and Technology, 4-5-7 Konan, Minato, Tokyo 108-8477,
Japan.
1

(2019) 6:67

Page 19 of 20

References
Amante C, Eakins BW (2009) ETOPO1 1 Arc-Minute Global Relief Model:
Procedures, Data Sources and Analysis. National Geophysical Data Center,
NOAA, Boulder, Colorado, doi:doi.org/10.7289/V5C8276M
Bell ML, Nur A (1978) Strength changes due to reservoir-induced pore pressure
and stresses and application to Lake Oroville. J Geophys Res Solid Earth 83:
4469–4483. https://doi.org/10.1029/JB083iB09p04469
Bernard E, Wei Y, Tang L, Titov V (2014) Impact of near-field, deep-ocean tsunami
observations on forecasting the 7 December 2012 Japanese tsunami. Pure
Appl Geophys 171:3483–3491. https://doi.org/10.1007/s00024-013-0720-8
Byerlee J (1978) Friction of rocks. Pure Appl Geophys 116:615–626. https://doi.
org/10.1007/BF00876528
Cai C, Wiens DA, Shen W, Eimer M (2018) Water input into the Mariana
subduction zone estimated from ocean-bottom seismic data. Nature 563:
389–392. https://doi.org/10.1038/s41586-018-0655-4
Chapple WM, Forsyth DW (1979) Earthquakes and bending of plates at trenches.
J Geophys Res 84(B12):6729–6749. https://doi.org/10.1029/JB084iB12p06729
Christensen DH, Ruff LJ (1988) Seismic coupling and outer rise earthquakes. J
Geophys Res Solid Earth 93:13421–13444. https://doi.org/10.1029/
JB093iB11p13421
Craig TJ, Copley A, Jackson J (2014) A reassessment of outer-rise seismicity andits implications for the mechanics of oceanic lithosphere. Geophys J Int 197:
63–89. https://doi.org/10.1093/gji/ggu013
Dmowska R, Lovison LC (1988) Intermediate-term seismic precursors for somecoupled subduction zones. Pure Appl Geophys 126:643–664. https://doi.org/
10.1007/BF00879013
Ekström G, Nettles M, Dziewoński AM (2012) The global CMT project 2004-2010:
Centroid-moment tensors for 13,017 earthquakes. Phys Earth Planet Inter
200–201:1–9. https://doi.org/10.1016/j.pepi.2012.04.002
Faccenda M, Gerya TV, Burlini L (2009) Deep slab hydration induced by bendingrelated variations in tectonic pressure. Nat Geosci 2:790–793. https://doi.org/
10.1038/ngeo656
Fujie G, Kodaira S, Kaiho Y, Yamamoto Y, Takahashi T, Miura S, Yamada T (2018)
Controlling factor of incoming plate hydration at the north-western Pacificmargin. Nat Commun 9:3844. https://doi.org/10.1038/s41467-018-06320-z
Gamage SSN, Umino N, Hasegawa A, Kirby SH (2009) Offshore double-planed
shallow seismic zone in the NE Japan forearc region revealed by sP depth
phases recorded by regional networks. Geophys J Int 178:195–214. https://doi.org/10.1111/j.1365-246X.2009.04048.x
Harada T, Murotani S, Satake K (2013) A deep outer-rise reverse-fault earthquake
immediately triggered a shallow normal-fault earthquake: The 7 December
2012 off-Sanriku earthquake (MW 7.3). Geophys Res Lett 40:4214–4219.
https://doi.org/10.1002/grl.50808
Hino R, Azuma R, Ito Y, Yamamoto Y, Suzuki K, Tsushima H, Suzuki S, Miyashita M,
Tomori T, Arizono M, Tange G (2009) Insight into complex rupturing of theimmature bending normal fault in the outer slope of the Japan Trench fromaftershocks of the 2005 Sanriku earthquake (Mw = 7.0) located by ocean
bottom seismometry. Geochem, Geophys Geosyst 10:Q07O18. doi: 10.1029/
2009GC002415
Hino R, Inazu D, Ohta Y, Ito Y, Suzuki S, Iinuma T, Osada Y, Kido M, Fujimoto H,
Kaneda Y (2014) Was the 2011 Tohoku-Oki earthquake preceded by aseismicpreslip? Examination of seafloor vertical deformation data near the epicenter.
Mar Geophys Res 35:181–190. https://doi.org/10.1007/s11001-013-9208-2
Hirata K, Aoyagi M, Mikada H, Kawaguchi K, Kaiho Y, Iwase R, Morita S, Fujisawa I,
Sugioka H, Mitsuzawa K, Suyehiro K, Kinoshita H, Fujiwara N (2002) Real-time
geophysical measurements on the deep seafloor using submarine cable in
the Southern Kurile subduction zone. IEEE J Ocean Eng 27:170–181. https://doi.org/10.1109/JOE.2002.1002471
Hunter J, Watts AB (2016) Gravity anomalies, flexure and mantle rheology
seaward of circum-Pacific trenches. Geophys J Int 207:288–316. https://doi.
org/10.1093/gji/ggw275
Iinuma T, Hino R, Kido M, Inazu D, Osada Y, Ito Y, Ohzono M, Tsushima H, Suzuki
S, Fujimoto H, Miura S (2012) Coseismic slip distribution of the 2011 off the
Pacific Coast of Tohoku Earthquake (M9.0) refined by means of seafloor
geodetic data. J Geophys Res 117:B07409. doi: https://doi.org/10.1029/
2012JB009186
Inazu D, Saito T (2014) Two subevents across the Japan Trench during the 7
December 2012 off Tohoku earthquake (Mw 7.3) inferred from offshore

tsunami records. J Geophys Res Solid Earth 119:5800–5813. https://doi.org/10.
1002/2013JB010892
Ito A, Fujie G, Miura S, Kodaira S, Kaneda Y, Hino R (2005) Bending of the
subducting oceanic plate and its implication for rupture propagation of large
interplate earthquakes off Miyagi, Japan, in the Japan Trench subduction
zone. Geophys Res Lett 32:L05310. https://doi.org/10.1029/2004GL022307
Kato T, Terada Y, Ito K, Hattori R, Abe T, Miyake T, Koshimura S, Nagai T (2005)
Tsunami due to the 2004 September 5th off the Kii peninsula earthquake,Japan, recorded by a new GPS buoy. Earth Planets Space 57:297–301. https://doi.org/10.1186/BF03352566
Kikuchi M, Kanamori H (1991) Inversion of complex body waves—III. Bull Seismol
Soc Am 81:2335–2350
Kikuchi M, Kanamori H (2003) Note on teleseismic body-wave inversion program.
http://www.eri.u-tokyo.ac.jp/ETAL/KIKUCHI/.
Koga S, Ito Y, Hino R, Shinohara M, Umino N (2012) Focal mechanisms of small
earthquakes within the Pacific Plate near the Japan Trench. Zisin 2(64):75–90(in Japanese with English abstract). https://doi.org/10.4294/zisin.64.75
Kubota T, Hino R, Inazu D, Ito Y, Iinuma T (2015) Complicated rupture process of
the Mw 7.0 intraslab strike-slip earthquake in the Tohoku region on 10 July
2011 revealed by near-field pressure records. Geophys Res Lett 42:9733–9739.
https://doi.org/10.1002/2015GL066101
Kubota T, Hino R, Inazu D, Ito Y, Iinuma T, Ohta Y, Suzuki S, Suzuki K (2017a)
Coseismic slip model of offshore moderate interplate earthquakes on March
9, 2011 in Tohoku using tsunami waveforms. Earth Planet Sci Lett 458:241–251. https://doi.org/10.1016/j.epsl.2016.10.047
Kubota T, Saito T, Ito Y, Kaneko Y, Wallace LM, Suzuki S, Hino R, Henrys S (2018b)
Using tsunami waves reflected at the coast to improve offshore earthquake
source parameters: Application to the 2016 Mw 7.1 Te Araroa earthquake,
New Zealand. J Geophys Res Solid Earth 123:8767–8779. https://doi.org/10.
1029/2018JB015832
Kubota T, Saito T, Suzuki W, Hino R (2017b) Estimation of seismic centroid
moment tensor using ocean bottom pressure gauges as seismometers.
Geophys Res Lett 44:10907–10915. https://doi.org/10.1002/2017GL075386
Kubota T, Suzuki W, Nakamura T, Chikasada NY, Aoi S, Takahashi N, Hino R(2018a) Tsunami source inversion using time-derivative waveform of offshorepressure records to reduce effects of non-tsunami components. Geophys J
Int 215:1200–1214. https://doi.org/10.1093/gji/ggy345
Lay T, Duputel Z, Ye L, Kanamori H (2013) The December 7, 2012 Japan Trench
intraplate doublet (Mw 7.2, 7.1) and interactions between near-trench
intraplate thrust and normal faulting. Phys Earth Planet Inter 220:73–78.
https://doi.org/10.1016/j.pepi.2013.04.009
Matsumoto K, Takanezawa T, Ooe M (2000) Ocean tide models developed byassimilating TOPEX/POSEIDON altimeter data into hydrodynamical model: A
global model and a regional model around Japan. J Oceanogr 56:567–581.
https://doi.org/10.1023/A:1011157212596
McNutt MK, Menard HW (1982) Constraints on yield strength in the oceanic
lithosphere derived from observations of flexure. Geophys J Int 71:363–394.
https://doi.org/10.1111/j.1365-246X.1982.tb05994.x
Obana K, Fujie G, Kodaira S, Takahashi T, Yamamoto Y, Sato T, Yamashita M,
Nakamura Y, Miura S (2015). SeismSSicity and structure of the incoming
Pacific plate subducting into the Japan Trench off Miyagi. Abstract T43G-01
presented at 2015 AGU Fall Meeting, San Francisco, CA, 14–18 Dec 2015.
Obana K, Fujie G, Takahashi T, Yamamoto Y, Nakamura Y, Kodaira S, Takahashi N,
Kaneda Y, Shinohara M (2012) Normal-faulting earthquakes beneath theouter slope of the Japan Trench after the 2011 Tohoku earthquake:
Implications for the stress regime in the incoming Pacific plate. Geophys Res
Lett 39:L00G24. doi: https://doi.org/10.1029/2011GL050399
Obana K, Fujie G, Takahashi T, Yamamoto Y, Tonegawa T, Miura S, Kodaira S(2019) Seismic velocity structure and its implications for oceanic mantle
hydration in the trench–outer rise of the Japan Trench. Geophys J Int 1629–
1642. doi: https://doi.org/10.1093/gji/ggz099
Obana K, Kodaira S, Nakamura Y, Sato T, Fujie G, Takahashi T, Yamamoto Y (2014)
Aftershocks of the December 7, 2012 intraplate doublet near the Japan
Trench axis. Earth, Planets Sp 66:24. https://doi.org/10.1186/1880-5981-66-24
Okada Y (1992) Internal deformation due to shear and tensile faults in a halfspace. Bull Seismol Soc Am 82:1018–1040
Peacock SM (2001) Are the lower planes of double seismic zones caused byserpentine dehydration in subducting oceanic mantle? Geology 29:299–302.
https://doi.org/10.1130/0091-7613(2001)029<0299:ATLPOD>2.0.CO;2
Saito M (1978) An automatic design algorithm for band selective recursive digital
filters. Butsuri Tanko 31:112–135 (in Japanese)

Kubota et al. Progress in Earth and Planetary Science

Received: 31 August 2019 Accepted: 13 November 2019

Kubota et al. Progress in Earth and Planetary Science(2019) 6:67

Saito T (2019) Tsunami Generation and Propagation. Springer Japan, Tokyo.
https://doi.org/10.1007/978-4-431-56850-6
Scholz CH (1988) The brittle-plastic transition and the depth of seismic faulting.
Geol Rundschau 77:319–328. https://doi.org/10.1007/BF01848693
Seno T, Gonzalez DG (1987) Faulting caused by earthquakes beneath the outerslope of the Japan Trench. J Phys Earth 35:381–407. https://doi.org/10.4294/
jpe1952.35.381
Tsushima H, Hino R, Tanioka Y et al (2012) Tsunami waveform inversion
incorporating permanent seafloor deformation and its application to tsunami
forecasting. J Geophys Res 117:B03311. https://doi.org/10.1029/2011JB008877
Turcotte DL, Schubert G (2002) Geodynamics, 2nd ed. Cambridge University
Press. https://doi.org/10.1017/CBO9780511807442
Wallace LM, Araki E, Saffer D et al (2016) Near-field observations of an offshore
Mw 6.0 earthquake from an integrated seafloor and subseafloor monitoring
network at the Nankai Trough, southwest Japan. J Geophys Res Solid Earth
121:8338–8351. https://doi.org/10.1002/2016JB013417
Wessel P, Smith WHF (1998) New, improved version of generic mapping tools
released. Eos, Trans Am Geophys Union 79:579. https://doi.org/10.1029/
98EO00426
Yoshida K, Hasegawa A, Yoshida T, Matsuzawa T (2018) Heterogeneities in stressand strength in Tohoku and its relationship with earthquake sequences
triggered by the 2011 M9 Tohoku-Oki earthquake. Pure Appl Geophys.
https://doi.org/10.1007/s00024-018-2073-9

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.

Page 20 of 20