1

Fault model of the 2012 doublet earthquake, near the up-2

dip end of the 2011 Tohoku-Oki earthquake, based on a3

near-field tsunami: Implications for intraplate stress state

4

Tatsuya Kubota12

5Corresponding author

6

Email: kubotatsu@bosai.go.jp

78

Ryota Hino2

9

Email: hino@tohoku.ac.jp

1011Daisuke Inazu3

12

Email: inazud@m.kaiyodai.ac.jp

1314

Syuichi Suzuki2

15

Email: syuichi.suzuki.a1@tohoku.ac.jp

1617(Institutional addresses)

18

1

19

Tsukuba, Ibaraki, 305-0006, Japan

National Research Institute for Earth Science and Disaster Resilience, 3-1 Tennodai,

20

2

21

Miyagi, 980-8578, Japan

22

3

23

Technology, 4-5-7 Konan, Minato, Tokyo 108-8477, Japan

Graduate School of Science, Tohoku University, 6-6 Aza-Aoba, Aramaki, Aoba-ku, Sendai,

Department of Marine Resources and Energy, Tokyo University of Marine Science and24
25

Abstract

26

On December 7, 2012, an earthquake occurred within the Pacific Plate near the Japan Trench,

27

which was composed of deep reverse- and shallow normal-faulting subevents (Mw 7.2 and28

7.1, respectively) with a time interval of ~10 s. It had been known that the stress state within29

the plate was characterized by shallow tensile and deep horizontal compressional stresses due

30

to the bending of the plate (bending stress). This study estimates the fault model of the31

doublet earthquake utilizing tsunami, teleseismic, and aftershock data and discusses the stress32

state within the incoming plate and spatiotemporal changes seen in it after the 2011 Tohoku-

33

Oki earthquake. We obtained the vertical extents of the fault planes of deep and shallow

34

subevents as ~45–70 km and ~5 (the seafloor)–35 km, respectively. The down-dip edge of the

35

shallow normal-faulting seismic zone (~30–35 km) deepened significantly compared to what

36

it was in 2007 (~25 km). However, a quantitative comparison of the brittle strength and37

bending stress suggested that the change in stress after the Tohoku-Oki earthquake was too

1

38

small to deepen the down-dip end of the seismicity by ~10 km. To explain the seismicity that39

occurred at a depth of ~30–35 km, the frictional coefficient in the normal-faulting depth range

40

required would have had to be ~0.07 ≤ μ ≤ ~0.2, which is significantly smaller than the typical

41

friction coefficient. This suggests the infiltration of pore fluid along the bending faults, down42

to ~30–35 km. It is considered that the plate had already yielded to a depth of ~35 km before

43

2011 and that the seismicity of the area was reactivated by the increase in stress from the44

Tohoku-Oki earthquake.

45
46

Keywords

47

Ocean bottom pressure gauge, Doublet earthquake, Intraplate earthquake, 2011 Tohoku-Oki

48

earthquake, Bending stress, Fault modeling

49

2

50
51

Introduction
It is well known that the stress state within the incoming Pacific Plate near the52

Japan Trench is characterized by shallow tensile and deep horizontal compressional stresses

53

along a direction perpendicular to the trench axis, separated by a thin aseismic (i.e., stress-

54

neutral) “elastic core”, due to the bending of the plate (bending stress, Figure 1c; e.g.,

55

Chapple and Forsyth 1979). It is also well-known that the number of normal-faulting

56

earthquakes occurring within the plate and near the trench axis increases, following interplate

57

megathrust earthquakes; this is attributed to the increased horizontal tensile stress caused by58

the stress release of the interplate coupling (e.g., Christensen and Ruff 1988; Dmowska and59

Lovison 1988).

60

Recently, Craig et al. (2014) investigated the vertical variation in the centroid depth

61

and fault mechanisms based on global catalogs (M > ~5) and noted that temporal changes in62

the transition depths between normal- and reverse-faulting earthquakes were not detected

63

after the 2011 Tohoku-Oki earthquake. However, Obana et al. (2012; 2014; 2015; 2019)64

studied ocean bottom seismographs and reported that the down-dip limit of shallow normal-

65

faulting earthquakes (M < ~5) within the plate near the trench axis deepened from ~25 to ~3566

km after the same earthquake; they interpreted spatiotemporal change in the intraplate stress67

as having been caused by the release of stress associated with the Tohoku-Oki earthquake.

3

68However, the reason for this inconsistency has not been clarified. Furthermore, the69

relationship between the stress state in the incoming Pacific Plate and the changes in70

coseismic stress due to the Tohoku-Oki earthquake have not yet been quantitatively assessed71

in detail.

72

On December 7, 2012, an Mjma 7.3 earthquake occurred within the Pacific Plate

73

near the Japan Trench, where the extremely large coseismic slip (> ~50 m) was estimated to74

have occurred during the 2011 Tohoku-Oki earthquake (e.g., Iinuma et al. 2012; star in Figure

751b). Detailed teleseismic analyses (Lay et al. 2013; Harada et al. 2013) revealed that this76

earthquake was composed of two M ~7 subevents. According to the global centroid moment77

tensor (GCMT) solution (http://www.globalcmt.org/; Ekström et al. 2012), the first subevent78

had a reverse-faulting mechanism with a depth of ~60 km (Mw 7.2) and the second had a79

normal-faulting mechanism (~20 km, Mw 7.2) with a time interval of ~12 s (red CMT

80solutions in Figures 1b and 1c). Hereafter, this earthquake is referred to as the doublet

81

earthquake, and the first and the second subevents are referred to as subevent 1 and subevent82

2, respectively. Since the fault mechanisms of the two subevents are consistent with bending83

stress, the source process of the doublet earthquake should reflect the intraplate stress state

84

after the Tohoku-Oki earthquake. The vertical extents of each fault will be key to discussing85

the temporal change in the vertical variations of the stress state after the 2011 Tohoku-Oki

4

86

earthquake.

87

Rupture processes related to the doublet earthquake have been investigated88

previously. Lay et al. (2013) and Harada et al. (2013) investigated this earthquake using

89

teleseismic data to estimate the CMT solution and the finite fault model. Teleseismic data is

90

generally a powerful dataset for resolving rupture processes of global earthquakes. However,91

since the teleseismic signals from each subevent overlapped, it is difficult to decompose the92

rupture process of the doublet earthquake precisely, especially for the latter, shallower,

93

subevent.

94

Inazu and Saito (2014) estimated the spatial distribution of the initial sea-surface

95

height change (tsunami source) using far-field tsunami data from ~200–2000 km away from96

the focal area (Figure S1). In contrast to teleseismic data, tsunami data is useful for

97

constraining the rupture process of the subevent 2, since shallow earthquakes generally excite

98

tsunamis or cause seafloor vertical deformation more effectively than deep earthquakes.

99

When the 2012 doublet earthquake occurred, off-line autonomous absolute ocean

100

bottom pressure gauges (PGs) installed near the focal area (< ~ 200 km from the source,

101

Figure 1a) recorded clear tsunami signals. This dataset is useful to constrain the fault model of102

subevent 2, which was difficult to constrain with teleseismic data. In the present study, we103

utilize the tsunami and aftershock data and the results of the teleseismic analysis to estimate

5

104

the finite fault model of the 2012 doublet earthquake, focusing particularly on the vertical105

extent of the fault planes of each subevent. We also discuss the relationship between the106

vertical profile of the intraplate seismicity and its spatiotemporal changes associated with the107

2011 Tohoku-Oki earthquake.

108
109

Data and Methods

110

Tsunami Data

111

We use the near-field PGs installed by Tohoku University (hereafter TPG; e.g., Hino

112

et al. 2014; Kubota et al. 2017a; 2017b) (green inverted triangles in Figure 1a). We also use

113

tsunami data obtained by the off-Kushiro online cabled PGs installed by Japan Agency for

114

Marine-Earth Science and Technology (JAMSTEC) (Hirata et al. 2002) (KPGs, pink

115

triangles), by GPS buoys installed by the Port and Airport Research Institute (PARI) (Kato et116

al. 2005) (yellow squares), and by the Deep-ocean Assessment and Reporting of Tsunamis

117(DART) system (Bernard et al. 2014) (blue diamonds). Detailed information is given in Table

118

1.

119

To retrieve the tsunami waveforms, we remove the ocean-tide component using the120

theoretical tide model (Matsumoto et al. 2000) and apply the filter from Saito (1978). The

121

lowpass filter we applied has a cutoff of 3 min to the TPG records, and the bandpass filter has

6

122

a passband of 3–60 min to the KPG, GPS buoy, and DART records. We apply the lowpass

123

filter to the TPG records to preserve the offset in the pressure change caused by the vertical124

deformation of the seafloor.

125

The filtered records are shown in Figure 2. TPGs first capture down-motion

126

tsunamis with amplitudes of ~−5 cm and then larger up-motion tsunamis with amplitudes of

127~+10 cm. The durations of both the down- and up-motion tsunamis are ~5 min (black dashed

128

lines in Figure 2a). Small fluctuations and changes in the pressure offset are observed at the129

stations near the focal area (Figure 2b). Tsunami amplitudes at the DART and the KPG

130stations are very small (~1 cm, black dashed lines in Figures 2c and 2d). At the former,

131

dynamic pressure changes caused by seismic waves (e.g., Kubota et al. 2017b) are also

132

observed. Tsunami signals are also detected by some GPS buoys (e.g., ~15 cm at station 801,

133

Figure 2e).

134
135

Step-by-Step Approach for Fault Modeling136

To decompose the complex rupture process of the 2012 doublet earthquake, we apply137

a step-by-step procedure to tsunami, teleseismic, and aftershock data. We first estimate an initial138

sea-surface height distribution of the tsunami (hereafter, the tsunami source model) by inverting139

tsunami records. Since seafloor crustal deformations, or tsunamis, are very sensitive to shallow

7

140

earthquakes, we then estimate the fault model of shallow subevent 2 based on the tsunami

141

source model. We also use aftershocks detected around subevent 2 (Obana et al. 2015) to obtain

142

information on the fault geometry. We then calculate the change in residual sea-surface height143

between the tsunami source model and the vertical displacement from subevent 2. Given that144

this residual distribution is caused by subevent 1, we then estimate the fault model of subevent

145

1. Furthermore, because there are large trade-offs between fault size and the focal depth and the146amount of slip in a deeper earthquake, we also use the results of teleseismic analyses (Lay et al.

147

2013; Harada et al. 2013) to obtain prior information on the fault parameters such as fault depth

148

and size for subevent 1.

149
150
151

Tsunami Source Modeling using Near-field Tsunami Records
We estimate the tsunami source model by inverting tsunami records via the152conventional inversion analysis method (e.g., Tsushima et al. 2012; Inazu and Saito 2014;

153

Kubota et al. 2018a). The details of the procedure are identical to those described in Kubota et154

al. (2018a). Before estimating the tsunami source model, however, we conduct a preparatory

155

tsunami simulation using the tsunami source model in Inazu and Saito (2014) (Figure S1).

156

Although the pressure fluctuation at the TPGs were explained roughly, the pressure offset

157changes at the stations closest to the source (G09, TJT1, JFAST, and GJT3) were not8

158

explained at all (Figure S1). This indicates that the latter are not due to coseismic seafloor

159

deformations. Thus, they are probably due to tilts or rotations in the sensors related to the160

seafloor strong ground motion (e.g., Wallace et al. 2016; Kubota et al. 2018a).

161

In calculating the Green’s functions for the tsunami, the unit source elements of the162

seafloor displacements (Kubota et al. 2015; 2018a) are distributed around the focal area. The

163

horizontal dimension of the unit source elements is 20 km × 20 km at a spacing of 10 km

164(overlapping with the adjacent elements), distributed along 260 km (in the EW direction) ×

165

240 km (NS) area. To calculate the sea-surface displacement from the unit source elements of166

the seafloor displacement, we consider the spatial filtering effect due to water depth (Saito

1672019). In the depth filtering process, we assume a sea water depth of 6 km. To simulate the168tsunami, we solve the linear dispersive wave equation in the local Cartesian coordinates (e.g.,

169

Saito 2019). The grid spacing is 2 km with a 1 s time step according to interpolation done via

170

ETOPO1 bathymetry data (Amante and Eakins 2009). We assume that the displacements of171

all unit sources occurred instantaneously and simultaneously. The tsunami propagation

172

velocity expected by the linear long-wave theory is expressed as v = (g0H)1/2 (g0: gravity

173

acceleration constant, H0: water depth). Given the assumed water depth of 6 km and average174

depth of the focal area, the propagation velocity is approximated as ~240 m/s (~15 km/min).

175

Thus, the tsunami propagation distance during the duration of the M~7 earthquake (~10 s) and

9

176

the time interval between two subevents (~10 s) is about 3 km, which is sufficiently small

177

compared to the extent of the tsunami source model (~100 km, Inazu and Saito (2014)). We

178

consider the static pressure offsets related to the calculation of Green’s function of the PGs for

179

permanent seafloor deformation (Tsushima et al. 2012; Kubota et al. (2018a)). The same filter

180

used for the observed records is applied to the simulated waveforms.

181

All data is resampled to 15 s intervals for the inversion. We use different time

182

windows for each station, including tsunami main phase (Table 1, thick black lines in Figure

1833). The smoothing constraint is imposed and its weight is determined based on the trade-off

184

between the weight and reduction in variance between the observed and simulated tsunami

185

waveforms (Figure S2). Since the G09, TJT1, JFAST, and GJT3 stations are located near the186source and are probably affected by seafloor ground shaking (Figure S1), we exclude these187

records from the inversion analysis. Since all TPG stations are located landward of the source

188

region, the constraint for the eastern edge of the tsunami source is likely not very good.

189

Therefore, to improve the source constraint further, we also use the DART and KPG records.

190

Since the amplitudes of TPG records are approximately ten times larger than the DART and191

KPG records, we weight the KPG and DART data at values ten times that of the TPG data.

192

10

193

Fault Modeling of Subevent 2

194

Because we found that the subsidence of the tsunami source was generated by the195

shallow subevent 2 (see Figure 4 and Results), we first estimate a fault model for subevent 2,

196

which best explains the subsidence region of the tsunami source model. We use the grid-

197

search approach proposed by Kubota et al. (2015; 2018b), which estimates an optimum

198

rectangular planar fault model with uniform slip. Because the short-wavelength component

199

disappeared in the tsunami source model due to the smoothing constraint imposed in the200

inversion and the spatial smoothing effect used in the deep-sea region during the tsunami

201

generation (see Figure 4 and Results), we consider the smoothing effect in fault modeling by

202

the following procedure. First, we calculate the seafloor deformation using a fault model203

candidate (a set of unknown parameters) (Okada 1992), and then simulate a tsunami. The

204

simulated waveforms are inverted to obtain the initial sea-surface height distribution, under205

the same conditions used in the inversion for the tsunami source model. Finally, we evaluate

206

the goodness of the tsunami waveform fitting from the fault model candidate, by comparing207

the subsided area of the tsunami source model and the inverted sea-surface height.

208

The geometry of the fault plane is assumed to be on the west-dipping plane of

209

GCMT solution (strike = 189° and dip = 50°), which is consistent with the planar structure of210

the aftershock (Obana et al. 2015). We also assume the rake angle from the GCMT solution (=

11

211

−90°). Because the aftershock alignment is located ~2 km west from the GCMT centroid, the212

fault plane is constrained to pass through the point that is 2 km west from the GCMT centroid

213(hereafter, referred to as the reference point). The unknown parameters are the distance from214

the reference point to both ends of the fault, along the strike (i.e., length) and dip (width)215

direction (L1, L2, W1, and W2; see Figures 5c and 5d). The search range for these parameters is216

determined based on the aftershocks and the horizontal extent of the subsidence area of the217

tsunami source. The fault length (L = L1 + L2) is assumed to be greater than the fault width (W218= W1 + W2), as L > W. The top of the fault plane (defined by parameter W2) is constrained as219

to not extend above the seafloor. The fault model candidate is assessed through variance220

reduction (VR) of the subsided areas between the tsunami source model and the fault model

221

candidate:

222
223

𝑉𝑅 = $1 −+,-./0∑:12)*/34565370 8*;<()*
9+,-./0 >∑:?;<=)*

9@ × 100(%),

(2)

224
225

where 𝑥GHIJKLM and 𝑥GLNOPQPNRM are the displacements of the sea-surface at the ith grid point,

226

from the tsunami source model and the fault model candidate assuming the unit slip,

227

respectively. N is the total number of grid points. The slip amount on fault D is determined so

228

that the VR takes the maximal value. We use the grid points within the subsided area of the

12

229

tsunami source model (blue dashed line in Figure 5a) to calculate the VR. The search range of230

the unknown parameters is listed in Table 2.

231
232

Fault Modeling of Subevent 1

233

In order to estimate the fault model of subevent 1, we use the residual sea-surface

234

height between the tsunami source model and the sea-surface displacement expected from the235

fault model of subevent 2 (hereafter, referred to as the residual height distribution). This is236

because the residual height distribution is expected to correspond to the sea-surface

237

displacement due to subevent 1. Because subevent 1 occurred at the deeper part of the incoming

238

plate, it appears to be difficult to constrain the fault parameters such as fault geometry, depth,

239

size, and slip amount only from the residual height distribution. Meanwhile, the residual height240

will contribute to constrain the horizontal location of the fault. Thus, we use the results of the241

teleseismic analyses in previous studies (GCMT; Lay et al. 2013; Harada et al. 2013) to obtain242

prior information on the fault parameters. We fix the centroid depth, fault geometry (strike, dip,243

and rake angles), and fault dimension (length and width) based on previous teleseismic analyses.

244

We estimate the optimum horizontal location (longitude and latitude) of the fault and the245amount of uniform slip on the rectangular fault using the grid-search approach, as in the fault

246

modeling of subevent 2.

13

247

Because Lay et al. (2013) investigated the CMT solution of the 2012 doublet

248

earthquake using the teleseismic W-phase waveforms and showed and discussed the uncertainty

249

of their estimations in detail, we use the CMT solution proposed by them as prior information

250

for our fault modeling. We assume the centroid depth as 60.5 km and the west-dipping nodal

251

plane with geometry of strike =163° dip = 51°, and rake = 57°. The slip distribution of subevent

252

1 obtained by the teleseismic analysis (Lay et al. 2013; Harada et al. 2013) had a main rupture

253

area with a dimension of L ~ 30 km and W ~ 20 km; therefore, we fix the fault length and width254

as 30 and 20 km, respectively (L1 = L2 = 15 km and W1 = W2 = 10 km). The search range of the255

horizontal location of the centroid is determined based on the evaluation of the uncertainty of256

the horizontal location of the W-phase analysis of Lay et al. (2013) (see fig. S2 in Lay et al.

2572013). The search range is listed in Table 2.

258
259

Results

260

Tsunami Source Model261

We obtained a tsunami source distribution that had a pair of large uplifts and262

subsidences (Figure 3a). The observed waveforms (red lines in Figure 3b) were reproduced

263

well. The GPS buoy waveforms, which were not used for the inversion, were also explained.

264

Furthermore, although the offset changes at the stations near the source (G09, TJT1, JFAST,

14

265

and GJT3) were not reproduced, the fluctuations in the calculated waveforms were similar to266

the observations. This is consistent with the idea that these changes were due to the tilting or

267

rotation of the sensors (e.g., Wallace et al. 2016; Kubota et al. 2018a).

268

To investigate the contribution from each subevent on the tsunami source, we

269

compute the sea-surface vertical displacement from the GCMT solution and compared the270

results with the tsunami source model. We use the equation in Okada (1992), which assumes

271

that the rectangular planar fault on the west-dipping nodal plane has a uniform slip. For

272

simplicity, values of length, width, and slip L = 70 km, W = 35 km, and D = 0.7 m,

273

respectively, are used. The spatial pattern of the tsunami source model (Figure 4a) is similar to274that of the combined deformation of subevents 1 and 2 (Figure 4b) and the subsidence area is275

similar to that expected from the CMT solution of subevent 2 (Figure 4d). The uplift of the276

tsunami source model is not consistent with either subevent 1 nor 2 alone (Figures 4c and 4d).

277

Based on this comparison, we conclude that the subsidence is generated by subevent 2 alone278

and that both subevents contribute to the uplift. The deformation expected from subevent 2

279(Figure 4d) has a sharp displacement peak that was not estimated in the tsunami source

280

model. This is probably because the short-wavelength components disappeared due to the281smoothing constraint imposed in the inversion and the spatial smoothing effect used in the282

deep-sea region during the tsunami generation (Saito 2019).

15

283
284

Fault Models of Two Subevents of the 2012 Doublet Earthquake285

We obtained an optimum fault model of subevent 2, which had the highest VR value

286

among all the fault model candidates, with a length of 35 km (L1 = 5 km, L2 = 30 km), width of287

25 km (W1 = 15 km, W2 = 10 km), and fault slip of 1.1 m (VR = 96 %, Figure 5). The seismic

288

moment was 5.5 × 1019 Nm (Mw 7.1, assuming the rigidity of 60 GPa). The VR values for all

289

fault model candidates (in the descending order) are shown in Figure 5b. The VR values for the290

top ten candidates are almost flat and relatively high (Figure 5d). This indicates that the top ten

291

candidates reasonably reconstruct the subsidence of the tsunami source model. Hence, we

292

inspect the model parameters for these candidates to evaluate the estimation error of the293

optimum fault model. They are projected onto the vertical cross section in Figures 5c and 5d

294(thin black lines), and the histograms of the model parameters are shown in Figure 6. Among295

the top ten candidates, most models has the down-dip limit of the fault plane of W1 = 15 km,

296

and all models has fault bottoms shallower than the depth of ~35 km (W1 ≤ 20 km, Figure 6).

297This indicates that the lower end of the fault of subevent 2 should be less than ~35 km.

298

We also calculated the sea-surface height assuming the faults with smaller dimension299

with larger slip (Figure S3). In this calculation, the down-dip end of the fault (W1) and slip

300

amount were changed, and the length (L1, L2) and the up-dip end (W2) of the fault, and seismic

16

301

moment Mo were fixed to those of the optimum model. In the fault models with W < 15 km

302(down-dip depth of fault is shallower than ~25 km), the locations of the western edge of the303

subsided area and of the peak displacement are inconsistent with those of the tsunami source

304

model. This indicates the small fault models are implausible. The subsidence of the tsunami

305

source model could be explained when the small faults are located slightly west of the optimum

306

fault location. However, such faults can be rejected because we used the aftershock distribution

307

of Obana et al. (2014; 2015) to constrain the horizontal location of the fault. It is important to308

use the aftershock distribution for prior information on the fault horizontal location, in order to309

accurately constrain the down-dip depth of the subevent 2 fault.

310

The subsided area calculated from the optimum fault model was consistent with that of311

the tsunami source model (Figures 7a and 7b). We calculate the residual sea-surface height312

between the tsunami source model (Figure 7a) and the sea-surface displacement expected from313

the fault model of subevent 2 (Figure 7b; the residual height distribution).

314

From the residual height distribution utilizing the results of the teleseismic analyses,

315

we constrained the fault model of subevent 1. The result is shown in Figure 8. The slip amount316

on the fault was 2.2 m, and the seismic moment was 7.8 × 1019 Nm (Mw 7.2, rigidity of 60

317

GPa). We obtained the optimum VR of 80 %. The VR values are relatively high and flat for the318

top ten fault model candidates (Figure 8b). The uncertainty of the horizontal location of the

17

319

fault center for these candidates is likely to be ± ~10 km (Figure 8a). This horizontal uncertainty

320

is almost consistent with that estimated by Lay et al. (2013).

321
322

Tsunami and Teleseismic Waveform Simulation from Optimum Fault Model

323

The results of the fault modeling of the 2012 doublet earthquake are summarized in

324

Figure 9. The sea-surface height displacement expected from the optimum rectangular fault

325

models of subevents 1 and 2 (Figure 9a) is calculated by the superposition of the displacements

326

from each fault model (Figures 5a and 8a). The distribution is very similar to that of the tsunami

327

source model (Figure 4). The optimum models has vertical ranges of ~6 (seafloor)–30 km for

328

shallow subevent 2 and ~50–70 km for deep subevent 1 (Figures 9d and 9e).

329

From the combined displacement, we numerically simulate tsunami waveforms (Figure

33010). The simulated waveforms reasonably explained the observed tsunami waveforms well, not331

only the near-field TPGs but also the far-field DART, KPG, and GPS buoys (red lines in Figure332

10b). We calculate the VR between the observed and calculated waveforms from equation (1),333

by using the same time window as in the inversion of tsunami records. We finally obtained VR

334= 79 %.

335

We also simulate the teleseismic P-waves using the fault model parameters for336

comparison with the observed teleseismic waveforms (Figure 11). We assume pure-double-

18

337

couple point sources at the centers of the optimum faults of each subevent. We use the338

calculation programs of Kikuchi and Kanamori (2003). A triangular-shaped source time339

function with the rise time of 6 s is assumed, considering the typical rupture duration of M ~7

340

earthquakes (Figure 9b). After simulating the waveforms of each subevent, we stack the341

simulated waveforms. We assume that the difference of focal times between subevents 1 and 2

342

was 10 s, which is determined by inspecting the waveform similarity of the observed and343

stacked waveforms.

344

We use a 1D multi-layered velocity structure model without the water layer, assuming

345

that the source structure was identical to the receiver structure in Table 3 (blue traces in Figure34611b). The simulated waveforms for each subevent are also shown in Figure S4. The peak timing347

and amplitudes of the first up-motion and the subsequent down-motion waves reasonably fitted348

the observation, although the subsequent phases during 50 – 80 s did not perfectly match. This349

is probably because of the assumptions of the velocity structure and simple source time function.

350

We then simulate the teleseismic waveforms incorporating the water layer and oceanic structure,351

shown in Table 3 (red traces in Figure 11b). The agreement of the subsequent phases improved352

compared with the simulation without the water layer. According to the teleseismic analysis by353

Lay et al. (2013), another smaller normal-faulting subevent was estimated at 40 s after subevent354

1. Hence, it is possible that the third smaller subevent also contributed to the generation of the

19

355

later arrival. It is worth pointing out that it is important to use the teleseismic records to resolve356

the temporal complexity of the doublet earthquake in detail.

357

Our fault model explains both tsunami and teleseismic observations. In addition, our358

fault model of subevent 2 is consistent with the aftershock distribution determined by the ocean359

bottom seismographs. The rupture area of subevent 2 estimated by Harada et al. (2013) was

360

located at the outer-trench region and concentrated in the shallower portion of the plate (z <361

~20 km). The horizontal location of subevent 2 centroid by Lay et al. (2013) (Figure 1b) was

362

also inconsistent with our fault model and with the aftershock locations. The consistency of our363

fault model with the tsunami, teleseismic waveforms, and aftershocks indicates that the step-364

by-step procedure used in this analysis can decompose the complex rupture process of the 2012365

doublet earthquake. We conclude that we can obtain a more comprehensive fault model of the3662012 doublet earthquake, than the one estimated from the teleseismic data alone.

367
368

Discussion

369

Importance of Near-field Tsunami Data for the Fault Modeling

370

Inazu and Saito (2014) showed a tsunami source model of the 2012 doublet earthquake

371

using offshore tsunami stations located more than 200 km from the source area (Figure S1). We

372

compare the tsunami source models of this study with that of Inazu and Saito (2014) (Figure

20

373

12a). The horizontal location of the tsunami source was in agreement with the model of Inazu

374

and Saito (2014), although the amplitudes were lower. The simulated tsunami at the TPG

375

stations using the tsunami source model of Inazu and Saito (2014) are similar to the tsunami

376

peak timing of the observation (Figure 12b). This indicates that the horizontal location of the377

tsunami source is reasonably constrained even when using tsunami stations located far from the378

source (> ~200 km).

379

However, the dominant period and amplitudes of the simulated waveforms are longer

380

and smaller than the observation. This indicates that the spatial resolution of the far-field381

tsunami data is not sufficient for the finite fault model, and thus, to constrain the down-dip limit382

of the fault plane of subevent 2. Using the near-field tsunami data, the constraint of the down-383

dip limit of the fault plane of shallow subevent 2 is improved. This enables us to discuss the384

intraplate stress regime after the 2011 Tohoku-Oki earthquake.

385

The fault dimension and depth of subevent 1 estimated in this study are consistent with386

previous results of the teleseismic waveforms (Lay et al. 2013; Harada et al. 2013). This is387

because we used the teleseismic analyses as the prior information for fault dimension. In

388

contrast, the horizontal location and fault dimension of subevent 2 differ from the teleseismic389

analyses. The down-dip limit of subevent 2 constrained from the teleseismic analysis (Harada390

et al. 2012) is considerably shallower (< ~20 km) than that of the constrained from tsunami data

21

391

(~40 km). Our fault model has better consistency with the aftershock distribution, for both392

horizontal location and depth range. It is likely that the rupture process of subevent 2 was not393

resolved well from the teleseismic data, because the teleseismic signal from subevent 2 overlaps

394with that of subevent 1. By using the near-field tsunami records, the resolution of subevent 2

395

rupture was considerably improved.

396
397

Temporal Change in Down-dip Limit of Normal-faulting Earthquake

398

The normal-faulting aftershocks around subevent 2 (Obana et al. 2014; 2015) mainly

399

occurred at depths of ~ 30 km, which corresponds the down-dip depth of the optimum fault of400

subevent 2 (Figure 9b). Some normal-faulting seismicity also occurred at depths of ~ 35 km,

401

along the down-dip direction (Figure 9b). According to Obana et al. (2014), the estimation error

402

of the aftershocks is less than 5 km (Figure 2 in Obana et al., 2014). Considering the estimation

403

error of the fault model and aftershock distribution, the down-dip limit where the shallow

404

normal-faulting seismicity can occur around this region is ~ 30–35 km (Figure 13).

405

We also investigate the temporal change of the intraplate seismicity before and after

406

the Tohoku-Oki earthquake (Figure 13). The down-dip limit of subevent 2 fault (~35 km) is

407

clearly ~10 km deeper, compared with the down-dip limit of the normal-faulting seismicity

408

observed in 2007 (<~25 km, Hino et al. 2009). Because both researches use the arrays of the

22

409

ocean bottom seismometers installed just above the focal area, which have identical sensitivity

410

and were distributed with almost identical spatial intervals (~10 km), the detectability in both411

observation periods is expected to be identical and thus the difference between the seismicity

412

depths of the lower limit of the shallow normal-faulting seismicity was confidentially413

significant. In addition, a few deeper (> ~40 km) events were detected in both observations.

414

This also suggests the misdetection of deeper shallow normal-faulting event (~30–40 km)

415

before the Tohoku-Oki earthquake is unlikely to occur.

416

In contrast, the up-dip limits of the subevent 1 fault (~50 km) and the deep reverse-

417

faulting seismicity (~50 km, Obana et al. 2015) are almost equivalent to the deep reverse-

418

faulting seismicity before the Tohoku-Oki earthquake (Seno and Gonzalez 1987; Hino et al.

4192009), although it is difficult to discuss this in detail because of the very low seismicity. By

420

focusing on the shallow normal-faulting seismicity, we discuss the cause of the deepening of421

the down-dip limit of the normal-faulting seismicity.

422

The yield strength of the plate is characterized by the brittle rupture at the shallow

423

portion and ductile failure laws at the deeper parts of the plate (Figure 14a, e.g., Scholtz 1998;

424

Turcotte and Schubert 2002; Hunter and Watts 2016). Based on the Anderson theory of faulting,425

the brittle strength along the horizontal direction normal to the trench axis τxx(z) is expressed as426(e.g., Turcotte and Schubert 2002):

23

427
428

UV(W Y Z1[ )X X\𝜏)) (𝑧) = HQO U]^V(_1LIH
,
U])

(3)

429
430

where ρ0 is the crust density, pw is the pore pressure, z is depth (downward is positive), δ is the431

fault dip angle, and μ is the frictional coefficient. This equation implies that the rock strength is432

proportional to depth z (green line in Figure 14a). Further, by assuming the plate as a rigid two-433

dimensional elastic plane (x- and z-axes are the subducting direction and vertical direction,434

respectively), the vertical distribution of the bending stress along the dip direction σxx(z) is

435

approximated as (e.g., Turcotte and Schubert 2002; Craig et al. 2014; Hunter and Watts 2016)

436
437

ab

𝜎)) (𝑧) = − _1c9 (𝑧 − 𝑧d ),

(4)

438
439

where E is the Young's modulus, C is the curvature of plate bending, ν is the Poisson's ratio,

440

and z0 is the stress-neutral depth (tensile stress is positive, blue line in Figure 14a).

441

At a shallower portion of the plate, where the bending stress exceeds the brittle strength442(the blue background area in Figure 14a), the rock cannot remain elastic and the stress is443

released, or the rock yields, leading to shallow normal-faulting earthquakes. In contrast, elastic

444

behavior is expected in the depth range where the bending stress does not exceed strength; the

24

445

range is termed the elastic core (e.g., Craig et al. 2014; Hunter and Watts 2016). The actual446

deviatoric stress profile within the plate is represented by the red solid line in Figure 14a. The

447

top of the aseismic elastic core, or the bottom of the vertical range of the normal-faulting

448

seismicity, can be defined as a depth where the bending stress and the frictional strength are449

equal. As the top of the elastic core is present at a depth of ~30–35 km and the top of the reverse-

450

faulting seismicity (the bottom of the elastic core) at ~45 km, the stress-neutral plane is expected451

to be located in the depth range between 30–35 and 45 km. This depth range is almost consistent452

with the depths where fault mechanisms flip from the shallow normal-faulting to the deep453

thrust-faulting mechanisms, near the trench axis off NE Japan (~25 – 40 km, Gamage et al.

4542009; Koga et al. 2012).

455

We compare the vertical profiles of the brittle strength and bending stress. Assuming E

456= 80 GPa and ν = 0.25, and C = 2 × 10-7 m-1 (McNutt and Menard 1982), dσxx/dz = EC/(1−ν2)457

is ~15 MPa/km. Assuming the hydrostatic pressure condition pw = ρwg0z (ρw = 1030 kg/m3,

458

seawater density) and ρ0 = 2700 kg/m3 , δ = 50°, and μ = 0.6 (e.g. Byerlee 1978), dτxx/dz is ~ 11459

MPa/km as per equation (3). In this situation, assuming z0 = 40 km, τxx(z) and σxx(z) are equal460

at z ~25 km (Figure 14b). This is consistent with the down-dip limit of normal-faulting

461

seismicity observed before the Tohoku-Oki earthquake (Hino et al. 2009). However,462

considering dσxx/dz = ~15 MPa/km, the stress increment of ~300 MPa is needed to deepen the

25

463

top of the elastic core from 25 to 35 km (orange arrow in Figure 14b). This value is too large

464

compared to the coseismic stress change around subevent 2 expected from the fault model of465

the Tohoku-Oki earthquake by Iinuma et al. (2012), Δσxx ~20 MPa. The expected depth change466

of the down-dip limit of the normal faulting seismicity is only a few km (Figure 14c).

467

It must be considered that τxx(z) and σxx(z) are equal at ~30–35 km depth before the468

Tohoku-Oki earthquake (Figure 14c). In order that the normal-faulting earthquakes occur at a469

depth of 35 km, the brittle strength must be reduced compared to the typical frictional condition

470(Figure 14b). If we assume z0 = 40 ± 5 km (Figure 13), the frictional coefficient of the fault

471

around the lower end of the fault of subevent 2 is μ ~0.07 ± 0.06 so that τxx(z) and σxx(z) are472

equal at z = 35 km (Figure 14c). Even when the top of the elastic core is assumed to be at z =473

30 km, the expected frictional coefficient μ is ~ 0.2 (Figure S5). These values are much smaller474

than that of the typical rocks, but comparable to that estimated for the other incoming plate

475(Craig et al. 2014).

476

Reduction of friction has often been reported in studies on inland earthquakes (e.g.,

477

Yoshida et al. 2018); this reduction can be attributed to the existence of the pore fluid (e.g., Bell478

and Nur 1978). Based on the active seismic survey in the outer-rise region of the Japan trench479(e.g., Fujie et al. 2018), the significant seismic wave velocity reduction and high Vp/Vs area

480

were detected at the shallow part of the subducting plate (<5 km), which are interpreted as being

26

481

the results of pore fluid penetration through the pre-existing bending faults in the shallower part482

of the plate (e.g., Peacock 2001). Considering these studies, it is suggested that that the strength483

reduction within the plate might be related to the pore fluid. Numerical modeling by Faccenda484

et al. (2009) demonstrated that the pore fluid can infiltrate the plate as deep as the lower limit

485of the normal faulting seismicity observed in this study. Cai et al. (2018) also reported the486

serpentinized mantle wedge associated with water infiltration into the subducting plate (down487

to ~35 km) at the Mariana subduction zone.

488

However, majority of the seismicity at depths of 30–35 km is located around the489

subevent 2’s fault (Figure 9) and less in the other portions of the plate. This localization of490

seismicity suggests that the pore fluid, or the strength reduction, is localized within the plate,491

as also suggested by Faccenda et al. (2009) and Obana et al. (2019).

492

The activation of the normal-faulting seismicity at depths of 25–35 km after the493

Tohoku-Oki earthquake can be interpreted as follows: the plate at 25–35 km depths had yielded494

before the Tohoku-Oki earthquake, leading to intrinsically aseismic region; a stress increment495

by the Tohoku-Oki earthquake enhanced the horizontal tensile stress in a broad depth range

496

near the top of the elastic core, which activates seismicity. Less normal-faulting seismicity at

497

depths of 25–35 km during observation from April to June 2007 (Hino et al. 2009) may be498

representative of the long-term-averaged deformation. It is expected that the stressing rate due

27

499

to the bending deformation is lower near the stress-neutral depth, where null-deformation is500

expected, than the shallowest part of the incoming oceanic plate.

501

On the other hand, if EC/(1−ν2) is small, the static stress change by the Tohoku-Oki

502

earthquake can contribute to deepening the elastic core or the lower limit of the shallow normal

503

faulting seismicity. Supposing that the elastic core is deepened by 10 km due to the coseismic504

static stress change Δσxx of ~20 MPa by the Tohoku-Oki earthquake, EC/(1−ν2) must be ~0.5505MPa/km, which is smaller than the typical elastic condition described above by an order of506

magnitude. However, it is unlikely that the Young’s modulus or the plate curvature are507

significantly reduced by an order of magnitude even if supposing the existence of the pore fluid508

or estimation error of the curvature. Thus, this hypothesis seems unlikely. Although it might be509

possible that the EC/(1−ν2) is small compared to that assumed in this study, its contribution for510

deepening the elastic core associated with the Tohoku-Oki earthquake is not highly significant.

511
512

Conclusions

513

In this study, we estimated the fault model of the intraplate doublet earthquake that514

occurred on December 7, 2012 (subevent 1: a deep reverse-faulting earthquake; subevent 2: a

515

shallow normal-faulting earthquake) strategically utilizing offshore tsunami, aftershocks, and516

the teleseismic records based on the step-by-step analysis procedure. First, the initial sea-

28

517

surface height distribution was estimated by inverting the offshore tsunami records and518

comparing it with the seafloor deformation from the CMT mechanism. It was found that the519

subsidence and uplift areas were generated by subevent 2 and both subevents, respectively.

520

Then, the fault model of each subevent was estimated based on the initial sea-surface height521

model, using information from previous studies. As a result, the vertical extent of the fault plane522

of subevent 2 was obtained as ~5 km (i.e., the seafloor) to 35 km. Finally, we simulated the523

tsunami and teleseismic waveforms from the fault model, which explained the observation well.

524

We compared the tsunami source model obtained from the near-field tsunami data525

acquired at less than 200 km from the epicenter and that from the far-field (>200 km) data. We

526

found that the horizontal location of tsunami source was reasonably constrained, even from the527

far-field tsunami data alone. However, to constrain the finite fault model in more detail, it is528

necessary to use the near-field tsunami records. We also discussed the stress state within the529

plate and its spatiotemporal change after the 2011 Tohoku-Oki earthquake. We found that the530

down-dip limit of the shallow normal-faulting earthquakes was obviously deepened compared531

with that observed in 2007, from 25 to 35 km. However, comparing the coseismic stress change532

by the Tohoku-Oki earthquake and the amount of the bending stress within the plate, the plate533

down to ~35 km in depth should have already yielded before the Tohoku-Oki earthquake, and534

the top of the elastic core been located at ~35 km. Furthermore, as the bending stress around

29

535

the top of the elastic core was much smaller than the rock strength expected from the empirical

536

relationship, the frictional strength in the range of the normal-faulting earthquakes is expected537

to be significantly reduced. The significant strength reduction of the plate suggests pore fluid

538

infiltration down to ~35 km, along the bending faults.

539
540

Abbreviations

541

CMT: Centroid Moment Tensor; DART: Deep-ocean Assessment and Reporting of Tsunamis;

542

GCMT: Global Centroid Moment Tensor; JAMSTEC: Japan Agency for Marine-Earth543

Science and Technology; KPG: Pressure Gauge Installed off Kushiro; PARI: Port and Airport

544

Research Institute; PG: Pressure Gauge; TPG: Pressure Gauge Installed by Tohoku

545

University; VR: Variance Reduction

546
547

Declarations

548

Availability of data and material

549

The TPG data is available in additional file. KPG data were obtained from the Submarine

550

Cable Data Center (SCDC), JAMSTEC (http://www.jamstec.go.jp/scdc/top_e.html).

551

Teleseismic data were downloaded from the Data Management Center (DMC) of Incorporated

552

Research Institutions for Seismology (IRIS)

30

553

(http://ds.iris.edu/ds/nodes/dmc/tools/event/3650366). GPS buoy data were obtained by a554

request to PARI. The data from DART, which is jointly operated by NOAA and JMA, was555

provided upon request by Eddie Bernard and Yong Wei from Science Applications

556

International Corporation (SAIC). Aftershock data was provided by Koichiro Obana.

557
558

Competing interests

559

The authors declare that they have no competing interest.

560
561

Funding562

This work was supported by the research project “Research concerning Interaction between

563

the Tokai, Tonankai, and Nankai earthquakes” of the Ministry of Education, Culture, Sports,

564

Science and Technology, Japan, by JSPS KAKENHI [Grant Numbers JP20244070,

565

JP26000002, and JP 19K14818], and by the Sasakawa Scientific Research Grant [2019-2037]

566

of the Japan Science Society.

567
568

Authors' contributions

569

TK designed the study, analyzed the data, interpreted the results, and drafted the manuscript.

570

RH contributed in fault modeling, data interpretation, and revision of the manuscript. DI

31

571

contributed in tsunami modeling and revising the manuscript. SS took part in installation and572

retrieval of the pressure data. All authors approved the final manuscript.

573
574

Acknowledgements

575

We thank the reviewer Thorne Lay and another anonymous reviewer, as well as the editor,

576

Frantisek Gallovic, for their reviews and insightful comments. We thank Koichiro Obana for

577

providing aftershock data and Keisuke Yoshida for providing stress change data due to the578

Tohoku-Oki earthquake. We also thank Eddie Bernard, Yong Wei, Science Applications

579

International Corporation (SAIC), and JMA for providing DART data. We also thank Editage580(www.editage.jp) for English language editing. TK thanks Koichiro Obana and Tatsuhiko

581

Saito for insightful discussions. Without the constructive comments and enormous

582

encouragement from Yusaku Ohta, this paper would not have been published.

583

32

584

Endnotes

585

References

586

Amante C, Eakins BW (2009) ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data

587

Sources and Analysis. National Geophysical Data Center, NOAA, Boulder, Colorado,588

doi:doi.org/10.7289/V5C8276M

589

Bell ML, Nur A (1978) Strength changes due to reservoir-induced pore pressure and stresses

590

and application to Lake Oroville. J Geophys Res Solid Earth 83:4469–4483.

591

doi:10.1029/JB083iB09p04469

592

Bernard E, Wei Y, Tang L, Titov V (2014) Impact of near-field, deep-ocean tsunami593

observations on forecasting the 7 December 2012 Japanese tsunami. Pure Appl Geophys594

171:3483–3491. doi:10.1007/s00024-013-0720-8

595
596

ByerleeJ(1978)

Frictionofrocks.

Pure

Appl

Geophys

116:615–626.

doi:

10.1007/BF00876528

597

Cai C, Wiens DA, Shen W, Eimer M (2018) Water input into the Mariana subduction zone598

estimated from ocean-bottom seismic data. Nature 563:389–392. doi: 10.1038/s41586-

599

018-0655-4

600
601

Chapple WM, Forsyth DW (1979) Earthquakes and bending of plates at trenches. J Geophys
Res 84(B12):6729–6749. doi: 10.1029/JB084iB12p06729

33

602
603

Christensen DH, Ruff LJ (1988) Seismic coupling and outer rise earthquakes. J Geophys Res
Solid Earth 93:13421–13444. doi: 10.1029/JB093iB11p13421

604

Craig TJ, Copley A, Jackson J (2014) A reassessment of outer-rise seismicity and its

605

implications for the mechanics of oceanic lithosphere. Geophys J Int 197:63–89.

606

doi:10.1093/gji/ggu013

607
608

Dmowska R, Lovison LC (1988) Intermediate-term seismic precursors for some coupled
subduction zones. Pure Appl Geophys 126:643–664. doi:10.1007/BF00879013

609

Ekström G, Nettles M, Dziewoński AM (2012) The global CMT project 2004-2010: Centroid-

610

moment tensors for 13,017 earthquakes. Phys Earth Planet Inter 200–201:1–9.

611

doi:10.1016/j.pepi.2012.04.002

612
613

Faccenda M, Gerya TV, Burlini L (2009) Deep slab hydration induced by bending-related
variations in tectonic pressure. Nat Geosci 2:790–793. doi: 10.1038/ngeo656

614

Fujie G, Kodaira S, Kaiho Y, Yamamoto Y, Takahashi T, Miura S, Yamada T (2018) Controlling615

factor of incoming plate hydration at the north-western Pacific margin. Nat Commun

616

9:3844. doi:10.1038/s41467-018-06320-z

617

Gamage SSN, Umino N, Hasegawa A, Kirby SH (2009) Offshore double-planed shallow

618

seismic zone in the NE Japan forearc region revealed by sP depth phases recorded by619

regional networks. Geophys J Int 178:195–214. doi:10.1111/j.1365-246X.2009.04048.x

34

620

Harada T, Murotani S, Satake K (2013) A deep outer-rise reverse-fault earthquake immediately

621

triggered a shallow normal-fault earthquake: The 7 December 2012 off-Sanriku earthquake

622(MW 7.3). Geophys Res Lett 40:4214–4219. doi: 10.1002/grl.50808

623

Hino R, Azuma R, Ito Y, Yamamoto Y, Suzuki K, Tsushima H, Suzuki S, Miyashita M, Tomori

624

T, Arizono M, Tange G (2009) Insight into complex rupturing of the immature bending

625

normal fault in the outer slope of the Japan Trench from aftershocks of the 2005 Sanriku

626

earthquake (Mw = 7.0) located by ocean bottom seismometry. Geochem, Geophys Geosyst

627

10:Q07O18. doi: 10.1029/2009GC002415

628

Hino R, Inazu D, Ohta Y, Ito Y, Suzuki S, Iinuma T, Osada Y, Kido M, Fujimoto H, Kaneda Y

629(2014) Was the 2011 Tohoku-Oki earthquake preceded by aseismic preslip? Examination

630

of seafloor vertical deformation data near the epicenter. Mar Geophys Res 35:181–190.

631

doi: 10.1007/s11001-013-9208-2

632

Hirata K, Aoyagi M, Mikada H, Kawaguchi K, Kaiho Y, Iwase R, Morita S, Fujisawa I, Sugioka

633

H, Mitsuzawa K, Suyehiro K, Kinoshita H, Fujiwara N (2002) Real-time geophysical634

measurements on the deep seafloor using submarine cable in the Southern Kurile

635

subduction zone. IEEE J Ocean Eng 27:170–181. doi: 10.1109/JOE.2002.1002471

636

Hjörleifsdóttir V, Ekström G (2010) Effects of three-dimensional Earth structure on CMT

637

earthquake

parameters.

Phys

Earth

35

PlanetInter

179:178–190.

doi:

638
639
640

10.1016/j.pepi.2009.11.003
Hunter J, Watts AB (2016) Gravity anomalies, flexure and mantle rheology seaward of circumPacific trenches. Geophys J Int 207:288–316. doi: 10.1093/gji/ggw275

641

Iinuma T, Hino R, Kido M, Inazu D, Osada Y, Ito Y, Ohzono M, Tsushima H, Suzuki S,

642

Fujimoto H, Miura S (2012) Coseismic slip distribution of the 2011 off the Pacific Coast643

of Tohoku Earthquake (M9.0) refined by means of seafloor geodetic data. J Geophys Res

644

117:B07409. doi: 10.1029/2012JB009186

645

Inazu D, Saito T (2014) Two subevents across the Japan Trench during the 7 December 2012646

off Tohoku earthquake (Mw 7.3) inferred from offshore tsunami records. J Geophys Res

647

Solid Earth 119:5800–5813. doi: 10.1002/2013JB010892

648

Ito A, Fujie G, Miura S, Kodaira S, Kaneda Y, Hino R (2005) Bending of the subducting oceanic

649

plate and its implication for rupture propagation of large interplate earthquakes off Miyagi,

650

Japan, in the Japan Trench subduction zone. Geophys Res Lett 32:L05310. doi:

651

10.1029/2004GL022307

652

Kato T, Terada Y, Ito K, Hattori R, Abe T, Miyake T, Koshimura S, Nagai T (2005) Tsunami

653

due to the 2004 September 5th off the Kii peninsula earthquake, Japan, recorded by a new654

GPS buoy. Earth Planets Space 57:297–301. doi: 10.1186/BF03352566

655

Kikuchi M, Kanamori H (1991) Inversion of complex body waves—III. Bull Seismol Soc Am

36

656
657
658

81:2335–2350
Kikuchi M, Kanamori H (2003) Note on teleseismic body-wave inversion program.
http://www.eri.u-tokyo.ac.jp/ETAL/KIKUCHI/. Accessed 20 May 2019

659

Koga S, Ito Y, Hino R, Shinohara M, Umino N (2012) Focal mechanisms of small earthquakes

660

within the Pacific Plate near the Japan Trench. Zisin 2 64:75–90 (in Japanese with English

661

abstract). doi:10.4294/zisin.64.75

662

Kubota T, Hino R, Inazu D, Ito Y, Iinuma T (2015) Complicated rupture process of the Mw 7.0

663

intraslab strike-slip earthquake in the Tohoku region on 10 July 2011 revealed by near-field664

pressure records. Geophys Res Lett 42:9733–9739. doi: 10.1002/2015GL066101

665

Kubota T, Hino R, Inazu D, Ito Y, Iinuma T, Ohta Y, Suzuki S, Suzuki K (2017a) Coseismic

666

slip model of offshore moderate interplate earthquakes on March 9, 2011 in Tohoku using

667

tsunami waveforms. Earth Planet Sci Lett 458:241–251. doi: 10.1016/j.epsl.2016.10.047

668

Kubota T, Saito T, Suzuki W, Hino R (2017b) Estimation of seismic centroid moment tensor

669

using ocean bottom pressure gauges as seismometers. Geophys Res Lett 44:10907–10915.

670

doi: 10.1002/2017GL075386

671

Kubota T, Suzuki W, Nakamura T, Chikasada NY, Aoi S, Takahashi N, Hino R (2018a) Tsunami

672

source inversion using time-derivative waveform of offshore pressure records to reduce673

effects

ofnon-tsunami

components.

37

GeophysJInt

215:1200–1214.

doi:

674

10.1093/gji/ggy345

675

Kubota T, Saito T, Ito Y, Kaneko Y, Wallace LM, Suzuki S, Hino R, Henrys S (2018b) Using

676

tsunami waves reflected at the coast to improve offshore earthquake source parameters:

677

Application to the 2016 Mw 7.1 Te Araroa earthquake, New Zealand. J Geophys Res Solid

678

Earth 123:8767–8779. doi: 10.1029/2018JB015832

679

Lay T, Duputel Z, Ye L, Kanamori H (2013) The December 7, 2012 Japan Trench intraplate680

doublet (Mw 7.2, 7.1) and interactions between near-trench intraplate thrust and normal

681

faulting. Phys Earth Planet Inter 220:73–78. doi: 10.1016/j.pepi.2013.04.009

682

Matsumoto K, Takanezawa T, Ooe M (2000) Ocean tide models developed by assimilating

683

TOPEX/POSEIDON altimeter data into hydrodynamical model: A global model and a684

regional model around Japan. J Oceanogr 56:567–581. doi: 10.1023/A:1011157212596

685

McKenzie D, Jackson J, Priestley K (2005) Thermal structure of oceanic and continental

686

lithosphere. Earth Planet Sci Lett 233:337–349. doi: 10.1016/j.epsl.2005.02.005

687

McNutt MK, Menard HW (1982) Constraints on yield strength in the oceanic lithosphere688

derived from observations of flexure. Geophys J Int 71:363–394. doi: 10.1111/j.1365-

689

246X.1982.tb05994.x

690

Mueller S, Choy GL, Spence W (1996) Inelastic models of lithospheric stress - I. Theory and691

application to outer-rise plate deformation. Geophys J Int 125:39–53. doi: 10.1111/j.1365-

38

692

246X.1996.tb06533.x

693

Obana K, Fujie G, Takahashi T, Yamamoto Y, Nakamura Y, Kodaira S, Takahashi N, Kaneda Y,

694

Shinohara M (2012) Normal-faulting earthquakes beneath the outer slope of the Japan

695

Trench after the 2011 Tohoku earthquake: Implications for the stress regime in the696

incoming Pacific plate. Geophys Res Lett 39:L00G24. doi: 10.1029/2011GL050399

697

Obana K, Kodaira S, Nakamura Y, Sato T, Fujie G, Takahashi T, Yamamoto Y (2014)

698

Aftershocks of the December 7, 2012 intraplate doublet near the Japan Trench axis. Earth,699

Planets Sp 66:24. doi: 10.1186/1880-5981-66-24

700

Obana K, Fujie G, Kodaira S, Takahashi T, Yamamoto Y, Sato T, Yamashita M, Nakamura Y,

701

Miura S (2015). Seismicity and structure of the incoming Pacific plate subducting into the702Japan Trench off Miyagi. Abstract T43G-01 presented at 2015 AGU Fall Meeting, San

703

Francisco, CA, 14–18 Dec 2015.

704

Obana K, Fujie G, Takahashi T, Yamamoto Y, Tonegawa T, Miura S, Kodaira S (2019) Seismic705

velocity structure and its implications for oceanic mantle hydration in the trench–outer rise706

of the Japan Trench. Geophys J Int 1629–1642. doi: 10.1093/gji/ggz099

707
708
709

Okada Y (1992) Internal deformation due to shear and tensile faults in a half-space. Bull
Seismol Soc Am 82:1018–1040
Peacock SM (2001) Are the lower planes of double seismic zones caused by serpentine

39

710

dehydration in subducting oceanic mantle? Geology 29:299–302. doi: 10.1130/0091-

711

7613(2001)029<0299:ATLPOD>2.0.CO;2

712
713
714
715
716
717
718
719
720
721

Saito M (1978) An automatic design algorithm for band selective recursive digital filters.
Butsuri Tanko 31:112–135 (in Japanese).
Saito T (2019) Tsunami Generation and Propagation. Springer Japan, Tokyo. doi:10.1007/9784-431-56850-6
Scholz CH (1988) The brittle-plastic transition and the depth of seismic faulting. Geol
Rundschau 77:319–328. doi: 10.1007/BF01848693
Seno T, Gonzalez DG (1987) Faulting caused by earthquakes beneath the outer slope of theJapan Trench. J Phys Earth 35:381–407. doi: 10.4294/jpe1952.35.381
Stein CA, Stein S (1992) A model for the global variation in oceanic depth and heat flow withlithospheric age. Nature 359:123–129. doi: 10.1038/359123a0

722

Tsushima H, Hino R, Tanioka Y, et al (2012) Tsunami waveform inversion incorporating723

permanent seafloor deformation and its application to tsunami forecasting. J Geophys Res

724

117:B03311. doi: 10.1029/2011JB008877

725
726
727

Turcotte DL, Schubert G (2012) Geodynamics, 2nd ed. Cambridge University Pressdoi:10.1017/CBO9780511807442
Wallace LM, Araki E, Saffer D, et al (2016) Near-field observations of an offshore Mw 6.0

40

728

earthquake from an integrated seafloor and subseafloor monitoring network at the Nankai

729

Trough,730

doi:10.1002/2016JB013417

731
732

southwest

Japan.

J

Geophys

Res

Solid

Earth

121:8338–8351.

Wessel P, Smith WHF (1998) New, improved version of generic mapping tools released. Eos,Trans Am Geophys Union 79:579. doi:10.1029/98EO00426

733

Yoshida K, Hasegawa A, Yoshida T, Matsuzawa T (2018) Heterogeneities in stress and734

strength in Tohoku and its relationship with earthquake sequences triggered by the 2011 M9735

Tohoku-Oki earthquake. Pure Appl Geophys. doi:10.1007/s00024-018-2073-9

736

41

Figure legends
141˚

142˚

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]737

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

40
60

B(ESE)
50738

20Depth [km]

37.8˚

100

80
150

Horizontal Distance [km]739

Figure 1. Location map of this study. (a) Green inverted triangles, yellow squares, pink740

triangles, and blue diamonds represent TPG, GPS buoy, KPG, and DART stations,741

respectively. Small dots denote the aftershocks deduced from the ocean bottom seismographs

42

742

(Obana et al. 2015). Gray contour lines indicate the coseismic slip distribution of the Tohoku-

743

Oki earthquake (Iinuma et al. 2012) with 10 m intervals. Black dashed line represents the744

trench axis. (b) Enlarged map around the focal area. Black star shows the epicenter by Japan

745

Meteorological Agency. Red, blue, and green CMT solutions are from GCMT, Lay et al.

746(2013) and Harada et al. (2013), respectively. Gray CMT solution denotes the 1967 mb 4.7

747

earthquake (Seno and Gonzalez 1987). (c) The vertical cross section along the A-B line in748

Figure 1a. Aftershocks within the dashed rectangle in Figure 1a are shown. Black curved line749

is the plate boundary (Ito et al. 2005). Schematic image of the intraplate bending stress state is

750

also shown (red: down-dip compressional stress, blue: extensional stress). The thick line

751

denotes the approximate location of the stress neutral plane.

752

43

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

0

10hPa ~ 10cm 6H2O

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 6H2O

20

801
806

0

0
-30

753

0

30

60

90

120

Time from EQ [min]

-30

0

30

60

90

Time from EQ [min]754

Figure 2. Filtered tsunami records. Light gray and dark gray waveforms are the de-tided

755

waveforms and moving-averaged records, respectively. Red waveforms are the filtered756

records.

757

44

120

(a)
39˚

KAMN
P06
P02

KAMS

G09

[cm]

P04P08GJT3
P05 P10
TJT1
P03 MYGI
JFAST
P07
P12
P01

38˚

15
10
5
0

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

GFK

142˚

143˚

-5
-10

0 km

120

Height [cm]

100
80
60
40
20

Height [cm]

0
60
40
20

144˚

TPG

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

5
4
3
2

758

Used for Inversion

DART
JP1
JP2

1
0
4
3
2
1

KPG

KPG2
KPG1
KCTD

0
80
60
40
20
0

0
0

145˚

Observed
Calculated6

TPG near the epicenter
G09
TJT1
JFAST
GJT3
-10

-15

Contour: 5 cm

(b)
140

50 km

10

20

30

Time from EQ [min]40

GPS Buoy

805
807
804
802
803
801
806
45
-10

0

10

20

30

Time from EQ [min]

40

50

759

Figure 3. Results of tsunami source inversion. (a) Distribution of the tsunami source. Contour

760

line interval is 5 cm. (b) Comparison of the observed (gray) and synthesized (red) waveforms.

761

Thick black lines denote the time window used for the inversion.

762(a) Tsunami source model(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

50 km

0 km

50 km

-15

763

37˚
143˚

764

Figure 4. Sea-surface displacement expected from the GCMT solution. Displacements from765(a) the tsunami source model, (b) the combination of both subevents, (c) subevent 1, and (d)766

subevent 2. The faults are represented by black rectangles.

144˚

145˚ 143˚

46

144˚

145˚

767

(b)

(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

100

10

100

1000

Order
0

(d)

(c)
W2

th
id
ta

lw

20
30

W1

Fault strike

40
50

A(West)
-40

768

L1

To

Depth [km]10

Total length
L2A’(East)-20

0

B(North)20-40

Horizontal Distance [km]-20

B’(South)
0

20

40

Horizontal Distance [km]769

Figure 5. Results for the grid-search of subevent 2. (a) The optimum fault model is shown by

770

the red rectangle. The top 10 fault model candidates are shown by thin black rectangles. (b)771

VR (black) and Mw(red) of the fault model candidates searched in the grid-search, arranged772

in the descending order in terms of VR. The rectangle with blue dashed lines denotes the area773

used for calculating VR. (c, d) Vertical profile along A-A’ and B-B’ lines in Figure 5a.

774

Aftershocks are taken from the rectangular area with the green lines in Figure 5a. The

775

configuration of the fault parameters is also shown.
47

776

Count

10

Top 10 models
W2 (Up-dip)

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

0

Width

10 20 30 40 50

Length

5
0
10

777
778

20

30

40

10

[km]

20

30

40

50

[km]

Figure 6. Histograms of the top 10 model parameters in the grid-search of subevent 2.

779

48

(a) Tsunami source (=#1+#2)(b) Subevent 2 model

38˚

0 km

50 km

0 km

50 km

37˚

(c) Residual distribution
Harada et al.

[cm](d) Subevent 1 model

15

GCMT

10

38˚

5
0
-5
-10
Lay et al.

0 km

50 km

0 km

50 km

-15

780

37˚
143˚

781

Figure 7. Sea-surface height changes from the models in the present study. (a) The tsunami

782

source model with contour intervals of 5 cm. The rectangular area enclosed by the blue

783

dashed lines denotes the area used for the fault modeling of subevent 2. (b) Subevent 2 fault

784

model; the optimum fault is also shown. (c) the residual between (a) and (b). CMT solutions

785

for subevent 1 by GCMT, Lay et al. (2013) and Harada et al. (2013) are also shown. The

786

rectangular area enclosed by the blue dashed lines is the area used for the fault modeling of787

subevent 1. (d) the subevent 1 fault model; the optimum fault is also shown.

144˚

145˚ 143˚

49

144˚

145˚

788

(b)

(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˚

1

145˚

789

Mw

100

10

100

Order

790

Figure 8. Results for the grid-search of subevent 1. (a). The optimum fault model is shown by791

the red rectangle. The top 10 fault candidates are shown by thin black rectangles. A rectangle

792with dot-and-dashed lines denotes the search range for the fault centroid location. (b) VR

793(black) and Mw (red) of the fault model candidates searched in the grid-search, arranged in the794

descending order in terms of VR.

795

50

(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]

2

(b)

34.8

10 20 30 40

Along strike [km]
0

(d)

(e)

Depth [km]

10
20
30
40
50
60
70

A
-60

796

A’
-40

-20

0

20

40

60

Horizontal Distance [km]

B
-60

B’
-40

-20

0

20

40

Horizontal Distance [km]797

Figure 9. Summary of the fault modeling. (a) Sea-surface height distribution calculated from798

the superposition of the those calculate from the fault models of subevents 1 and 2. Locations

799

of the optimum faults are also shown. (b) The moment rate function used for the teleseismic800

waveform calculation. (c) Aftershock distribution on the fault plane of subevent 2.

801

Aftershocks with distances from the fault plane smaller than 5 km are shown. (d, e) Vertical802

profile along A-A’ and B-B’ lines in Figure 9a. Orange and blue rectangles are the optimum

803

fault models of subevents 1 and 2, respectively. Aftershocks are taken from the rectangular
51

804

areas in Figure 9a.

805

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

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

60
40
20

GPS Buoy

805
807
804
802
803
801
806

0

0
-10

806

Observed

6

0

10

20

30

40

-10

Time from EQ [min]

0

10

20

30

40

50

Time from EQ [min]807

Figure 10. Comparison of the observed (gray) and calculated (red) waveforms based on the808superposition of the optimum fault models of subevents 1 and 2.

809

52

(a)

BORG
KIEV
KIV

SFJD

BILL

ABKT
UOSS
FURI

HRV
SSPA
FFC
WCI

KEV

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
IU.SNZO.00
II.MSEY.00
Az=30o Dist=92o Az=57o Dist=76o Az=157o Dist=83o Az=265o Dist=91o Az=339o Dist=63o
II.WRAB.00
II.FURI.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
IU.NWAO.00
II.UOSS.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

810

Time from P-arrival [s]811

Figure 11. Results of the forward calculation of the teleseismic waveforms based on the812

optimum fault model. (a) Locations of the teleseismic stations. (b) Comparison of the813

teleseismic waveforms. Black lines are the observed waveforms, and the synthetic waveforms
53

814

using the velocity structure with and without the water layer are shown by red and blue lines,

815

respectively. The bandpass filter of 10 – 500 s is applied.

816

54

[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

144.5

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

60
40
20

GPS Buoy

805
807
804
802
803
801
806

0

0
-10

817

144.0

Observed
This study
Inazu and Saito (2014)

5
4

143.5

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

0

10

20

30

Time from EQ [min]

40

55

-10

0

10

20

30

Time from EQ [min]

40

50

818

Figure 12. Comparison with the result of Inazu and Saito (2014). (a) Tsunami source model of819

Inazu and Saito (2014). (b) Cross-section of the tsunami sources from this study (red line) and820

Inazu and Saito (2014) (blue) at 37.8°N. (c) Comparison of the observed waveforms (gray)821

and synthesized waveforms. Blue dashed and red lines are synthesized using the tsunami

822

source model of the present study and Inazu and Saito (2014), respectively.

823

824
825

Figure 13. Schematic representation of the vertical-temporal change of seismicity in the outer-826

rise region. Blue and red stars denote the normal-faulting and reverse-faulting small seismicity

827

shown by previous studies (Seno and Gonzalez 1987; Hino et al 2009; Obana et al. 2012; 2014;

8282015). Note that the aftershocks occurred around the Mw 7.6 outer-rise ~40 min after the
56

829

Tohoku-Oki earthquake (Obana et al. 2012) are located ~60 km east from the focal area of the830

2012 doublet earthquake.

831

re
St

th

ng

µ=0.6

Depth

ng
ndi

ss

stre

Be

6mxx 5 20 MPa

40

Depth

gth

ng
di
n
Be6m 5 300 MPa
xx

30

n
Stre

Aseismic (m < o)

ss

re

st

20

µ=0.6

th

TensileSeismic (m > o)sesstrength 30
40

Depth

str

[km]

ric

Ductile832

Seismic (m > o)

20

10

c
ismi
ly se
ntial
ss
stre

(b,c)

200

0

iatric

10

[MPa]
100

0

Dev

strengssn

Be

500 Tensileiatss
ic streBrittle

0

400

Pote

e
str

100

v
De

tr
Devia

g
din

0

(c) µ=0.07[MPa]
200 300µ=0.07

Aseismic(m < o)(b) µ=0.6

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

833

Figure 14. Vertical profile of the strength and the bending stress within the incoming plate. (a)834

A schematic image of vertical profile. Green, gray dashed, and blue lines represent brittle

835

strength, ductile strength, and bending stress, respectively. After the yielding of rock, deviatoric836

stress occurring within the plate is expressed by a red line. Blue and red background colors

837

show the areas where shallow normal-faulting and deep thrust-faulting earthquakes occur,

838

respectively. (b) Stress profile at the shallow part of the plate assuming the typical frictional839

strength condition. Dark red line is the deviatoric stress after the Tohoku-Oki earthquake. Gray

840

hatched area is the region where seismicity becomes active after adding the stress change. (c)841Stress profile assuming the reduced frictional strength. Dark red line is the deviatoric stress

842

after the Tohoku-Oki earthquake, calculated by adding a constant static stress change of 20 MPa.
57

843
844

Tables845

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
P01
P02
P03
P04
P05
P06
P07
P08
P10
P12

37.581238.333138.5006
38.183438.3163
38.3000
38.6338
38.0000
38.2833
38.2500
37.8206

142.7647
142.4167
142.5035
142.3996
142.5657
142.7004
142.5833
142.4486
142.8329
143.1666
142.8996

2245
1038
1109
1056
1265
1412
1269
1064
1424
2066
1635

140
180
180
170
160
150
180
170
140
110
130

0 – 1200
300 – 1500
300 – 1500
300 – 1500
300 – 1500
300 – 1500
300 – 1500
300 – 1500
150 – 1350
0 – 1200
0 – 1200

TPG
TPG
TPG
TPG
TPG
TPG
TPG
TPG
TPG
TPG
TPG

1
1
1
1
1
1
1
1
1
1
1

KAMN
KAMS
MYGI
G09
GJT3
JFAST
TJT1
801
802
803
804
805
806
807
KPG1
KPG2
KCTD

38.8862
38.6347
38.0832
38.4782
38.294837.9336
38.2079
38.2325
39.2586
38.8578
39.6272
40.6333
36.971440.1167
41.7040
42.2365
41.6675

143.3639
143.2621
142.9166
143.7922
143.4811
143.9154
143.7904
141.6836
142.0969
141.8944
142.1867
141.7500
141.6836
142.0667
144.4375
144.8454
144.3409

2360
2246
1697
5500
3260
6799
5744
144
204
160
200
87
137
125
2218
2210
2540

150
130
130
90
90
40
60
240
250
240
320
380
250
320
450
510
440

0 – 1200
0 – 1200
0 – 1200
Not usedNot usedNot usedNot usedNot usedNot usedNot usedNot usedNot usedNot usedNot used
1200 – 2400
1500 – 2700
1200 – 2400

TPG
TPG
TPG
TPG
TPG
TPG
TPG
GPS Buoy
GPS Buoy
GPS Buoy
GPS Buoy
GPS Buoy
GPS Buoy
GPS Buoy
KPG
KPG
KPG

1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
10

58

JP1
JP2
846

40.3777
39.2849a146.1681
145.7845

5125
5183

330
210

1050 – 2250
750 – 1950

DART
DART

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

847
848

Table 2. Search range of the grid-search for fault modeling
Subevent 2
Reference point latitude [°N]
Reference point longitude [°E]
Reference point depth [km]a
Strike [°]a
Dip [°]a
Rake [°]a
L1 [km]b
L2 [km]b
W1 [km]b
W2 [km]b
Slip amount [m]b

Subevent 1

37.77c37.6 – 38.1d
143.81c143.8 – 144.5d
19.5
57.8
189
158
50
59
−90
48
0, 5, 10, ..., 50e15e0, 5, 10, ..., 50
1510, 15, 20, ..., 40e10e10, 15
10
Adjusted so that the VR values become maximal

849aStrike, dip, rake angles, and reference point depth were fixed to the GCMT value.

850b

Rectangular fault model with uniform slip is assumed. The fault length (L = L1 + L2) was

851

assumed to be greater than the fault width (W = W1 + W2).

852

c

A point 2 km west of the GCMT centroid.

853

dHorizontal location of reference point was searched within the range of ±30 km in the EW

854

and NS direction from the centroid location of Lay et al. (2013) (37.82°N, 144.13°E), with855

increments of 5 km.

856eIncrements for the fault length and width are 5 km.
59

15
15

857
858

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

Vp [km/s]

ρ [g/cm3]

H [km]

Source3

Receiver

3

1.50
6.00
8.10
6.00
6.75
8.10

0.00
3.50
4.70
3.50
3.80
4.70

1.00
2.70
3.30
2.70
2.80
3.30

6.0
6.0
Half Space
18.0
18.0
Half Space

859

a

860

km) is assumed for the source structure.

The structure is based on Kikuchi and Kanamori (1991) but the water layer (thickness of 660

Additional files for:Fault model of the 2012 doublet earthquake, near the up-dip
end of the 2011 Tohoku-Oki earthquake, based on a near-fieldtsunami: Implications for intraplate stress state

Tatsuya Kubota12, Ryota Hino2, Daisuke Inazu3, Syuichi Suzuki2

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, Japan3

Department of Marine Resources and Energy, Tokyo University of Marine Science andTechnology, 4-5-7 Konan, Minato, Tokyo 108-8477, Japan1

Observed
Inazu and Saito (2014)
140
120
38˚

[cm]
-15
37˚
143˚

0

0 km

50 km

144˚

60

145˚

21416
21402
KPG2
KPG1

20

21419

0
60

Height [cm]

801

30˚

80

40

50˚

40˚

15

Height [cm]

100

21413

40
20

TPG

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

TPG near the epicenter
G09
TJT1
JFAST
GJT3

0
-10
140˚

150˚

160˚

0

10

20

30

Time from EQ [min]Figure S1. Tsunami source model of Inazu and Saito (2014) and the stations used for theirinversion analysis. Comparison between the observed and simulated tsunami waveforms at
the stations near the source is also shown.
2

40

100

90

VR [%]

80

70

60

40
10-2

0.3

50

10-1

100

101

Smoothing Weight
Figure S2. Trade-off curve between the smoothing weight and the VR used for the inversion
analysis.

3

143˚

144˚

145˚ 143˚

DOWN=-21.1cm

144˚

145˚ 143˚

Length=35km
W1=0km (Width=10km)
VR=85.41%
DOWN=-21.9cm

144˚

145˚

Length=35km
W1=5km (Width=15km)
VR=91.38%
DOWN=-20.7cm

38˚

0 km

50 km

0 km

50 km

0 km

50 km

37˚

Length=35km
W1=10km (Width=20km)
VR=94.98%
DOWN=-19.3cm

Length=35km
W1=15km (Width=25km)
VR=96.41%
DOWN=-18.0cm

Length=35km
W1=20km (Width=30km)
VR=96.12%
DOWN=-16.7cm

38˚

0 km

50 km

0 km

50 km

0 km

50 km

37˚

Length=35km
W1=25km (Width=35km)
VR=94.60%
DOWN=-15.5cm

Length=35km
W1=30km (Width=40km)
VR=92.32%
DOWN=-14.4cm[cm]

Length=35km
W1=35km (Width=45km)
VR=89.62%
DOWN=-13.5cm

15
10

38˚

5
0
-5
0 km

50 km

0 km

50 km

0 km

50 km

-10
-15

37˚

Figure S3. Sea-surface height changes from subevent 2 with fault dimensions changed. The
left-top panel shows the distribution of the tsunami source model, and others show the seasurface height distribution of assumed fault models.

4

(a) Subevent 1
IU.SFJD.00

Obs.
IU.WCI.00

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

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

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

30

60

90

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

30

60

90

0

30

60

90

0

30

60

90

Time from P-arrival [s](b) Subevent 2
IU.SFJD.00

Obs.
IU.WCI.00

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

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

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

30

60

90

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

30

60

90

0

30

60

90

0

30

60

90

Time from P-arrival [s]Figure S4. Result of the forward calculation of the teleseismic waveforms based on the optimum
fault models of (a) subevent 1 and (b) subevent 2.

5

Figure S5. Vertical profile of the strength and the bending stress within the incoming plate,assuming the top depth of elastic core at 30 km.

6Dataset S1. (the CSV file is uploaded separately) The TPG data used in this study.

7