Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

1

Using tsunami waves reflected at the coast to improve offshore earthquake source2parameters: Application to the 2016 Mw 7.1 Te Araroa earthquake, New Zealand

3

Tatsuya Kubota1, Tatsuhiko Saito1, Yoshihiro Ito2, Yoshihiro Kaneko3,

4

Laura M. Wallace3,4, Syuichi Suzuki5, Ryota Hino5, and Stuart Henrys3

51

National Research Institute for Earth Science and Disaster Resilience, Tsukuba, Japan.

62Disaster Prevention Research Institute, Kyoto Univesity, Uji, Japan.

7

3

GNS Science, Lower Hutt, New Zealand.

8

4

University of Texas Institute for Geophysics, Austin, Texas, USA.

95

Graduate School of Science, Tohoku University, Sendai, Japan.

10
11Corresponding author: Tatsuya Kubota (kubotatsu@bosai.go.jp)

12
13
14

Key Points:
•

15
16

observed by offshore ocean bottom pressure gauges
•

17
18
19

Tsunami from coastal reflected waves associated with the 2016 Te Araroa EQ are clearly

We greatly reduce uncertainties in the centroid location and fault dimensions of the Te
Araroa EQ using tsunami reflected from the coast•

Later tsunami arrivals should be used more widely to extract earthquake source
information in the future1Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

20
21

Abstract
The Te Araroa earthquake occurred on September 2, 2016 (local time) offshore of the22northeastern coast of New Zealand’s North Island (Mw 7.1). When this event occurred, ocean

23

bottom pressure gauges (OBPs), installed ~170 km south of the source area, clearly recorded24

direct tsunami from the source to OBPs (~−1.5 cm), and tsunami from coastal-reflections (~225cm). We estimate the centroid location that best reproduces the OBP waveforms. When using the26

direct wave alone, the centroid location is poorly constrained, with a horizontal uncertainty of

27~100 km. By combining both direct and reflected tsunami waveforms, we obtain a centroid28

location near the Global CMT centroid (~80 km northeast from the coast) with smaller

29

uncertainty (~40 km). We also estimate the earthquake source dimension (length and width) and30

found that the models using coastal reflections require a source dimension larger than ~30 km

31long. Based on the slip distribution obtained by the finite fault inversion, we obtain an energy-

32

based stress drop ΔσE of 1.0 MPa, consistent with typical earthquake stress drop values. This

33

study shows that the information added by coastal reflected tsunami provides much tighter

34

constraints on the centroid location, source dimension, and stress drop of offshore earthquakes,

35

which is difficult to obtain from the onshore seismic data alone. Future studies should utilize the36

information provided by coastal reflected waves to improve earthquake source modeling using

37

ocean bottom pressure data.

2

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

38
39

1. Introduction
On 2 September 2016 at 04:37:55 local time (1 September, 16:37:55 UTC), an Mw 7.1

40earthquake occurred at 36.98°S, 179.52°E, with a depth of 22 km

41(http://www.geonet.org.nz/earthquake/2016p661332, Figure 1a). This earthquake was located42~80 km northeast offshore Te Araroa in the Raukumara Peninsula, the North Island, New

43

Zealand (widely referred to as the Te Araroa earthquake). The earthquake had a normal-faulting

44

mechanism and occurred within the subducting Pacific Plate at the Hikurangi subduction zone

45

(Figure 1b) (Warren-Smith et al., 2018). The centroid moment tensor (CMT) solution of the Te46

Araroa earthquake was estimated using teleseismic data by Global CMT

47(http://www.globalcmt.org, hereinafter, GCMT) and USGS48

(http://earthquake.usgs.gov/earthquakes/eventpage/us10006jbi), and from regional seismic data49

obtained by GeoNet of GNS Science, New Zealand (Ristau, 2008; http://www.geonet.org.nz).

50

Although their depths, strikes, dips, rakes, and seismic moments are similar, their horizontal

51

locations are quite different (Figure 1a, Table 1). Both the GCMT and USGS centroids are52

located ~80 km northeast from the coast, but they are ~20 km apart from each other in the north-53

south direction. The GeoNet centroid is located ~50 km northeast from the GCMT and USGS

54

centroids (~130 km from the coast). The Te Araroa earthquake source parameters, and how it fits

55

into the regional tectonic framework is very poorly characterized, due to the earthquake’s

56

offshore location and a lack of seismic stations in the offshore region to provide regional57

azimuthal coverage. Here, we show a new approach that combines modeling of direct tsunami

58

waves with tsunami waves reflected from the coastline to yield greatly improved constraints on

59

the location and source parameters of the Te Araroa earthquake. We suggest that the same3

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth60

approach can be used to improve understanding of other, poorly characterized offshore

61

earthquakes.

62

It is important to accurately estimate earthquake source parameters to infer earthquake

63

processes, evaluate the physical properties surrounding the rupture area, and to understand how

64

the earthquake plays a role in seismic hazard and the regional geodynamic framework. Many

65

previous studies have assessed the accuracy of the modern worldwide earthquake catalogs,

66

especially CMT catalogs (e.g., Hjörleifsdóttir & Ekström, 2010; Valentine & Trampert, 2012).

67

Seismic waveforms are one of the key pieces of data needed to determine CMT solutions. The

68

seismic waveforms have an advantage of high spatiotemporal resolution, which enables us to69

extract the source time function, or spatiotemporal evolution of the rupture on the fault plane. In

70general, the robustness of the estimated source parameters such as centroid location and centroid

71

time is restricted by various factors, such as the azimuthal coverage of seismic stations, S/N ratio,72

or uncertainties in velocity structure. Although station coverage is improved by teleseismic data,73

the uncertainty in horizontal centroid locations remains tens of km for global events in many

74

parts of the world that are far from seismic networks

75(https://earthquake.usgs.gov/earthquakes/eventpage/terms.php). When an earthquake occurs far

76

from the coast, like the Te Araroa earthquake, the azimuthal coverage of regional seismic

77

stations is usually poor and the uncertainty of the centroid horizontal location is large in the78

direction away from the coast, due to a tradeoff between the centroid time and the horizontal

79

location.

80

As well as onshore seismic data, tsunami data is often analyzed to constrain earthquake81

source parameters (e.g., Gusman et al., 2017; Kubota et al., 2015; 2017a), although this approach82

typically has lower spatiotemporal resolution compared to the seismic waves. On the other hand,4Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth83

because we can use reliable bathymetry data, we can simulate the propagation process of tsunami

84

more accurately than that of seismic waves. Furthermore, because tsunami are not as sensitive to85

the centroid time difference as seismic wave analyses, the tradeoff between the centroid time and86

the horizontal location will be much smaller than that of regional seismic data. Tsunami have87

another advantage in constraining earthquake source dimensions (i.e., fault length and width) as

88these data contain unique information about the area of seafloor deformation (e.g., the extent of89

the tsunami source), because the tradeoff between the earthquake source dimension and rupture

90

velocity is much less significant for tsunami than for seismic waves.

91

When the Te Araroa earthquake occurred, an offshore observation array consisting of four92

ocean bottom pressure gauges (OBPs) spaced ~ 10 km apart, was located ~ 170 km south of the93source area (green triangles in Figure 1a). The long-period signals (period of T ~10 min) were94

observed by four OBPs, at ~25 min and at ~ 90 min after the mainshock (noted by orange and95

yellow dots in Figure 1c, respectively). Because there were no other significant earthquakes and96

aftershocks in this region around the time of the earthquake, these two signals are most likely

97

due to tsunami generated by the Te Araroa earthquake, in spite of the small magnitude (Mw 7.1)

98

and large distance (~ 170 km) of the earthquake from the OBP network. The first tsunami signals

99

are likely to travel directly from the source region (direct tsunami). Because the travel time of the100

secondary signals at ~90 min is too large compared to the direct tsunami, the secondary signals

101

are probably coastal reflected tsunami waves, which propagate in the shallow water region near102

the coast and are reflected along the eastern coastline of the Raukumara Peninsula. We should103

note that such reflected waves are typically excluded from the source estimation analysis and104

only the direct waves are analyzed because it is often difficult to reproduce the reflected waves in105

the observed records with numerical simulations.

5Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

106

This study uses tsunami waveforms observed by the offshore OBP array (Figure 1a) to107

better constrain the centroid location and source dimension of the 2016 Te Araroa earthquake.

108

However, it should be noted that all the available OBPs were located ~ 170 km south of the109

earthquake location. It is difficult to constrain the centroid location and fault dimensions

110

precisely from the direct waves observed by the OBPs, because the ray paths of the direct waves

111

to the closely-spaced OBP array are similar. To overcome this, we analyze not only the direct

112

tsunami wave arrivals, but also the tsunami waves reflected from the coast. Previous studies have

113

recognized the existence of coastal reflected tsunami waves (e.g., Suppasri et al., 2017; Gusman

114

et al., 2017), but they have not previously been used to constrain the earthquake source

115

parameters. This is likely due to the complexity of the reflected waves, which largely arises from116

nonlinear effects in shallow water.

117
118

2. Ocean Bottom Pressure Data and Data Processing

119

The OBP array that we use consists of four autonomous instruments equipped with

120

Absolute Pressure Gauges, deployed in early June 2016, and recovered with acoustic release121

systems mounted on the instruments in late June of 2017. The instruments were located at the122

stations KU16-2, KU16-3, KU16-4, and KU16-5 (Figure 1a, Table 2). Figure 1c shows the123

processed pressure records during the time surrounding the Te Araora earthquake. For the124

processed data, we first removed the ocean-tide components, from the 1-s sampled raw pressure

125

data, using a theoretical tidal model (Matsumoto et al., 2000) (light gray lines in Figure 1c). Then,126

in order to reduce the high frequency components due to seismic or acoustic waves in the127

seawater and beneath the seafloor, we took a moving average with a 60-sec time window (dark128gray lines), and applied a bandpass filter (Saito, 1978) with a passband of 180 – 3600s (red lines).

6

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth129

Before applying the bandpass filter, the static pressure offset was subtracted using the 20-min

130

average before the focal time.

131

We recognize large impulsive signals with pressure changes of ~10 hPa with a dominant132

period T less than ~200 s at the time soon after the focal time (Figure 1c). This is the pressure

133change caused by seismic waves (e.g., An et al., 2017; Kubota et al., 2017b; Saito, 2013; 2017).

134

There are also coherent tsunami signals between the four OBPs, with maximum amplitude of135~−2 cm at ~25 min from the focal time (orange rectangular area in Figure 2) and another set with

136~+1.5 cm at ~90 min (yellow rectangular area). We note that 1 hPa of pressure change is137

equivalent to 1 cm water height change, assuming a water density of 1.03 g/cm3. Based on linear

138

long wave theory (e.g., Satake, 1995), the tsunami propagation velocity v is ! =139

is the water depth. Assuming the average water depth is 2000 m (e.g., average water depth

140

between the source to the OBPs), we obtain v as ~100 m/s (~6 km/min), and then the travel time141

becomes ~28 min assuming the epicentral distance is 170 km. This is consistent with the timing

142of the peak of the first tsunami signals. The duration of the downward wave of the first direct

143

tsunami arrival and the upward wave of the secondary reflected tsunami are ~ 8 – 10 min, which144

corresponds to a horizontal extent of the subsided area of ~ 48 – 60 km, based on an approximate

145

tsunami propagation velocity around the source region (v ~6 km/min).

#$ % where H

146
147
148

3. Tsunami Forward Simulations
In the first step of the analysis, we simulate tsunami records using the CMT solutions

149

estimated by GCMT, GeoNet, and USGS. In this calculation, we assume a rectangular planar

150

fault with uniform slip and calculate the seafloor displacement. The fault length (L) and width

151(W) are assumed based on fault scaling laws of Wells and Coppersmith (1994) (Table 1). The

7

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

152
153

slip amount on the fault plane (D) is calculated by the following relationship:
&=

'(
)*+

,

(1)

154

where Mo is the seismic moment and µ is the rigidity (40GPa). The center of the fault plane

155(longitude, latitude and depth) is assumed to coincide with the centroid of the CMT solution,

156

since, in the case of a rectangular fault with uniform slip, the centroid is located at the fault

157

center. The strike, dip, and rake angles are fixed to those of each CMT solution. The location and

158

extent of the fault planes are shown by colored rectangles in Figure 2a (red: GCMT, blue:

159

GeoNet, green: USGS). Using these rectangular fault models, permanent seafloor vertical160

deformation is calculated using the equations of Okada (1992). We adopted the northeast-161

dipping nodal plane as the fault geometry consistent with double-difference relative hypocenters

162of the smaller pre-cursory events and the mainshock determined by Warren-Smith et al. (2018).

163

The region of predicted seafloor deformation for each CMT solution is quite different

164

between the three solutions (Figure 2a). Thus, the tsunami arrival times should also be different

165

between the three CMT solutions. The maximum subsidence expected from the GCMT, GeoNet,

166

and USGS solutions is 18.6 cm, 26.5 cm, and 24.5 cm, respectively. The diameter of the167

subsided and uplifted areas is ~ 50 km, consistent with the expected source size from the tsunami

168

propagation velocity and the tsunami duration (~ 48 – 60 km).

169

We then calculate the tsunami by solving the linear long wave (LLW) equation in

170

Cartesian coordinates, using a finite difference scheme on a discretized staggered grid (e.g.,

171

Satake, 1995; Saito et al., 2014). We assume that the initial sea surface displacement is equal to172

the seafloor displacement. As the duration of a typical M~7 earthquake is ~10 – 20 s, which has

173

little effect on the resulting tsunami waves, we assumed instantaneous displacement on the174

seafloor. We used 250 m resolution gridded bathymetry data from National Institute of Water8

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth175

and Atmospheric Research (NIWA), in New Zealand released in 2016176

(https://www.niwa.co.nz/our-science/oceans/bathymetry). The grid spacing used is Δx = Δy =177

500 m, and the time step interval is Δt = 1 s. To compare the calculated and observed tsunami

178

calculation, we apply a bandpass filter similar to the filter applied to the observed records.

179

We compare waveforms from the different source models for the station KU16-2 in Figure

180

2b. The maximum amplitude and duration of the direct waves (~25 min) observed by OBPs are181

reasonably explained by all CMT solutions. The arrival time matches the GCMT and GeoNet

182

solutions reasonably well, while the match to the USGS solution is worse (orange rectangular183area). The GCMT and USGS solutions explain the peak timing of the secondary waves at ~90

184

min reasonably well, while the GeoNet solution does not (yellow rectangular area). This185

tendency is also seen in the other OBP waveforms (Figures 2c to 2e). The GCMT solution also

186

matches the tsunami first arrivals observed at the coastal tide gauge at East Cape (Figure 2f),

187

while the GeoNet solution does not. Therefore, we conclude that the GCMT solution’s centroid

188

location provides the best match among the three CMT solutions. We also show the comparison

189of the model results for coastal tide gauges in Figure S1. The calculated waveforms using the190

GCMT solution do not fit the tide gauges, except for East Cape. This is likely due to the lack of

191

finer grids that accurately capture the bathymetry within bays and ports. This suggests that the192

use of the offshore OBP data, which is less sensitive to the fine-scale bathymetry, can greatly

193

improve the accuracy of the source parameter estimation of off-shore earthquakes, which is also

194

discussed by Gusman et al. (2017).

195

We also conduct a tsunami simulation using the finite fault model obtained from196

teleseismic data, provided by USGS (Figure S2,

197

https://earthquake.usgs.gov/earthquakes/eventpage/us10006jbi#finite-fault). Although the arrival9Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

198

time of the direct and reflected waves is slightly improved compared to that from the rectangular

199

fault based on the USGS solution, the model tsunami wave amplitudes are much smaller than200

those observed. This suggests that the slip distribution is not well constrained by the teleseismic

201

data. We also simulate the tsunami using the conjugate nodal plane of the GCMT solution (strike202

of 220°, Figure S3). The calculated seafloor deformation and waveforms are similar to those203

assuming the northeast-dipping plane (Figure 2), but the amplitudes of the first downward wave

204

are slightly larger than observed (due to the steeper dip and larger normal fault component for205

the conjugate plane).

206

In order to test the idea that the secondary tsunami waveform recorded at ~90 min is

207

caused by tsunami reflected from the coast, we conducted an additional tsunami simulation

208

removing the land area (thus the source of reflections does not exist) (Figures 3 and S4). In the209

calculation, we modified the water depth on the land and near-shore, shallow shelf region (H <210

200 m) to equal a constant depth (H = 200m), so that the tsunami does not reflect at the coast

211(i.e., the coastline does not exist) (Figure 3). We used the GCMT solution to calculate the initial212

sea surface height, and kept the other settings in the calculation the same as the original ones.

213

The first direct tsunami wave arrival is similar to that obtained using the true bathymetry (Figure

214

2, compare red and blue lines in Figure 3). However, the secondary tsunami peak at ~ 90 min

215

does not appear in the simulation without the coastline. Even if we assume that the land and the216bathymetry shallower than 10 m are modified to a constant sea depth of 10 m, the observed

217

reflected waves cannot be reproduced (red lines in Figure S4), suggesting that the reflection at218

the coast is required to generate the secondary tsunami waves. Therefore, we conclude that the219

secondary tsunami represents a coastal reflection from the east coast of the Raukumara Peninsula.

220

Furthermore, the small tsunami wave train prior to the secondary peak (at ~ 60 – 90 min) was not

10

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth221

reproduced by the simulation without the land, whereas the simulation that includes bathymetry222

reproduced this well. This suggests the wave train between 60 – 90 min are also from tsunami

223

reflected at the coast.

224

To investigate where the strong reflected tsunami signals at ~90 min are radiating from, we

225

conduct some additional analyses (Figures 4, S4 to S8; the details are in Text S1). First, we226

investigate the snapshots of the tsunami height distribution from the simulations with and227

without the land (Figures S5a and S5b, respectively), the difference between them (Figure S5c),

228

and their maximum absolute height distribution (Figures S6). We find that the main tsunami

229

energy is directly radiating from the source to the middle part of the East Coast of the230

Raukumara Peninsula (Figures S6b) and is trapped near the middle section of the East Coast of231

the Raukumara Peninsula (Figure S6c), indicating that the tsunami largely radiates from the232

middle part of the East Coast of the Raukumara Peninsula. We also conducted additional tsunami

233

calculations in order to confirm that the coastal reflection occurred in the middle section of the234

East Coast of the Raukumara Peninsula (Figure S7)., We divided the Raukumara Peninsula into

235

three segments (northern, middle, and southern) and we assigned land to each segment but

236

removed the land from the other segments (we modified the water depth shallower than 10 m in237

the other segments to a constant depth of 10 m) in turn; northern (Figure S7a1), middle (Figure

238

S7a2), or southern (Figure S7a3; see the details in Text S1). The simulation results best239

reproduce the observed reflected waves at ~90 min when the land is assumed at the middle

240

segment of the Raukumara Peninsula (red lines in Figure S7b2), whereas the simulations that the241land is assigned to the other segments (blue: northern, and green: southern) did not (although the242tsunami at ~ 70 min are similar to the observation). Therefore we conclude that the coastal

11

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth243

reflections are largely radiating from the central part of the eastern coast of the Raukumara

244

Peninsula, and moderately from the northern and southern segments.

245

We undertake ray tracing of the reflected tsunami using the pseudo-bending method (Um246& Thurber, 1987), assuming the reflection points at the middle part of the East Coast of the247

Raukumara Peninsula and the source at the maximum subsidence of the initial sea surface height248

calculated based on the GCMT solution (Figure 4). We find that an 85–90-minute travel time of249

the reflected waves is consistent with the observed data, whereas the travel times are slightly

250

shorter when the source is assumed at the GCMT centroid and at the maximum uplift based on

251

the GCMT solution (Figure S8), indicating the reflected waves at ~ 90 min are mainly generated

252

by the subsidence of tsunami source. We also find that the travel time of the direct tsunami

253

obtained by the ray tracing (dashed line in Figure 4) is approximately 25 min, also consistent254

with the observations.

255

We note that the timing of the peak tsunami (Table 2, noted by small dots in Figure 1c) is

256

slightly different between the OBP stations. By using the time delay between the two OBPs, we

257

can estimate the approximate incident azimuth of the direct and reflected tsunamis. We estimate258

the time delay of tsunami arrivals between the two OBPs, based on cross correlations (Tables S1

259

and S2, details are described in Text S2). The estimated incident azimuths for the direct and260

reflected tsunamis are 71° and 29°, respectively (Figure S9). This suggests that the incident

261

azimuth and travel times of the reflected waves are substantially different from those of the262

direct wave. By using the additional information from the reflected wave, we are able to better

263

constrain the centroid location and the fault dimension.

264
265

Accurately modeling tsunami propagation near the coast might require the consideration of
nonlinear wave propagation (e.g., Satake, 1995; Saito et al., 2014). In order to assess this, we

12Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth266

compared the simulation results derived by the nonlinear long wave (NLL) equation and the267

LLW equation with the same initial tsunami height distribution and found that the results of the268

LLW and NLL simulations are almost identical, for both the direct and the reflected waves

269(Figure S10, details are described in Text S3). We also compare the simulation results using a270finer grid spacing (Δx = Δy = 250 m) (Figure S10) and confirm that the difference in the271

calculated OBP waveforms for grid spacings of 500 m and 250 m are small. However, there are272

some differences between the two models for the coastal tide gauge (East Cape), particularly for273

the later arrivals (>~ 60 min). The nonlinearity is small in the OBP waveforms of the Te Araroa

274

earthquake due to the small tsunami amplitude at the coast. If we assume an earthquake with275

larger magnitude (Mw 8.0) which generates much larger tsunami in the numerical simulation

276(Figure S11, details are described in Text S3), the nonlinear nature of tsunami would start to277

appear in the reflected waves.

278
279
280

4. Determination of the Centroid Location
In this section, we determine the earthquake centroid location by reproducing the observed

281

tsunami using a grid search approach, as the horizontal locations of the GCMT, GeoNet, and

282

USGS solutions are somewhat different. We fix the strike, dip, rake, seismic moment, and283

centroid depth to the GCMT value (GeoNet and USGS values are similar), and we assume the284

rectangular planar fault dimensions using the GCMT solution discussed earlier (Table 1).

285

Seafloor deformation is calculated using the equations of Okada (1992) and the tsunami is286

calculated based on the LLW equation (details for tsunami calculation is described in Text S4).

287

We search for the best-fitting horizontal centroid locations in a 300 km × 300 km region, with 5288

km intervals. The model waveforms are evaluated based on the variance reduction between the13Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth289

observed and calculated waveforms, VR, expressed as:

290

: 234
5016786
1 01: 0 234 9
11

291

VR = 1 −9×100 (%),(2)

where uiobs and uicalc are the observed and calculated waveforms, respectively.

292

At first, we use only the direct waves (time window of 10 – 40 min, orange rectangular293

area in Figure 5d). An optimum solution is obtained for a centroid located at ~15 km north of the294

GeoNet centroid (~ 130 km northeast from the coast), shown by a large red star in Figure 5a. The

295

direct waves are reasonably reproduced (VR = 76.3 %), although the reflected waves are not296(blue lines in Figure 5d). The centroid locations that gave relatively high VR, are the solutions

297

with a VR of more than 90% of the optimum solution’s VR (an area surrounded by green lines in298

Figure 5a), and extend ~100 km in the WSW-ENE direction. This suggests the centroid location299

is not well-constrained by the direct wave alone.

300

To evaluate why the high-VR area spans a larger region aligned in the WSW-ENE

301

direction, we conduct ray-tracing of tsunami from the OBPs to the solutions with high VR

302(Figure S12). Most of the rays from solutions with high VR arrive at the OBPs with similar

303

travel times of 22–24 min. Moreover, the incident azimuth to the OBPs from the solutions in the304high-VR area are similar (~60 – 65°). These make it difficult to resolve the differences in the305

centroid horizontal location. We note that the incident azimuths are consistent with that expected306

from the time delay of the tsunami arrival at the OBP array.

307

We also use tsunami reflected from the coast (70 – 100 min) in addition to the direct308

tsunami (10 – 40 min) to estimate the centroid location (Figure 5b). We obtain a centroid309

location at ~10 km northwest of the GCMT centroid (178.97°E, 37.10°S, ~ 80 km northeast from310

the coast), which is part of the high-VR area estimated by the direct wave alone (Figure 5a). The

14

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth311

extent of the high-VR area is reduced to ~40 km in the WSW–ENE direction and does not312

coincide with the GeoNet and USGS centroids. The calculated reflected waves reasonably

313

reproduced the observed, as well as the direct waves (red lines in Figure 5d, VR = 73.2 %). Note

314

that the optimum solution estimated from the direct wave alone (Figure 5a) did not explain the315

amplitudes and the arrival times of the reflected waves, as well as the forward-calculated

316

waveforms from the GeoNet solution. When we used the reflected wave alone, we observe a317

similar tendency (Figure 5c). The centroid was located at 178.92°E, 37.15°S, 5 km west and 5318

km south of that from both direct and reflected waves. The incident azimuth and travel times of319

the reflected waves are substantially different from the direct waves (Figure S9). Using the320

additional information from the reflected waves, we can reduce the uncertainty in the centroid321

location.

322
323
324

5. Estimation of the Fault Dimension
We then undertake a search for the earthquake source dimensions (Figure 6). We fix the325

centroid depth of the rectangular fault to the GCMT centroid (27.8 km), the seismic moment Mo

326

and fault strike, dip, and rake to the GCMT value, and we fix the centroid location (longitude and327

latitude) to that obtained from the grid search of the centroid location using both the direct and328

reflected waves (Figure 5b). We assume the ratio of source length L to width W such that L/W =329

2, with a rigidity of µ = 40 GPa (note that the source dimension ratio expected from the GCMT

330

solution and the scaling law of Wells & Coppersmith (1994) is ~2; L = 45.1 km and W = 21.6 km,

331

Table 1), and vary the fault length L at 5 km intervals. We show the result using only the direct332

waves in Figure 6 (blue curve). Considering the range of solutions with VRs larger than 90 % of333

the optimum solution, we estimate the source length uncertainty as L = 40 km ± 20 km. When

15

#2

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth334

using both the direct and reflected waves, we find that the optimum source length is L = 50 ± 10

335

km (red curve in Figure 6). In both analyses, the upper limit of the possible source length is336

similar (~60 km). However, the results using both the direct and reflected waves suggest that a337

smaller source dimension (L < 30 km) is not plausible. Moreover, source dimensions of L < 40338

km also seem improbable when using the reflected waves alone (green curve in Figure 6).

339

To evaluate how sensitive the calculated waveforms are to source dimension, we compare340

the calculated waveforms at the station KU16-2 for a variety of source dimensions (Figure 7)341(see Figure S13 for other stations). The tsunami with larger fault dimensions (L > ~ 60 km) have342

smaller amplitude and longer duration for the direct waves, compared to the observations. When343

the fault dimensions are small (e.g., L < ~ 40 km), the amplitudes of the model reflected waves344

are larger compared to the observations, especially for the earlier part of the reflected wave (70 –345

85 min). Fault dimensions with L ~40 – 60 km provide a reasonable fit to the observed tsunami

346

waveform amplitude and duration.

347

By using the coastal reflected tsunami waveforms, we obtain better constraints on the348source dimension. Our source dimension estimated from the scaling rules of Wells and349

Coppersmith (1994) (Table 1) is similar to the solution obtained from both direct and reflected350

waves. In contrast, Warren-Smith et al. (2018) suggested the possibility that the source

351

dimension is larger than expected from typical fault scaling relationships, based on investigation352

of the GeoNet seismograms (see Section 4.2 in Warren-Smith et al. (2018)). Using the range of353

source dimensions obtained from the analysis of both the direct and reflected waves, we354

calculate the stress drop Δσ, based on the conventional relationship of Kanamori and Anderson355(1975), as:356

AB = C

)D

(3)

*+

16

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth357

where c is constant. Assuming Poisson’s ratio of 0.25, we obtain c = 8/3π. Assuming L = 2W and358

using the seismic moment, we rewrite equation (3) as,
AB = C

359

E E'(
*F

.

(4)

360

Substituting the estimated values of L = 40 – 60 km, we obtain a possible range of stress drops361

Δσ = 0.5 – 3.0 MPa. This value is consistent with typical values for earthquake stress drops (~ 1362

– 10 MPa) (e.g., Abe, 1975; Kanamori & Anderson, 1975).

363
364

6. Finite Fault Inversion

365

In order to investigate the stress drop in more detail, we conduct a finite fault inversion to366

estimate the slip distribution (e.g., Satake, 1989), using the direct and reflected tsunami waves.

367

The details of the inversion are shown in Text S5. We assume a rectangular planar fault 80 km

368long and 70 km wide, based on the optimum fault dimension obtained by the grid search. The

369

center of the assumed fault coincides with the optimum centroid obtained by the grid search

370(Figure 5b). We use the fault geometry (strike, dip, and rake) from the GCMT solution. We

371

divided the planar fault into the subfaults with size 5 km × 5 km, and use these to calculate the372Green’s functions for the inversion. In the inversion, we impose a non-negativity constraint

373(Lawson & Hanson, 1978) and a spatial smoothing constraint. The weighting of the smoothing is374

based on the value of Akaike’s Bayesian Information Criterion (ABIC) (Yabuki & Matsu’ura,3751992) (Figure S14).

376

We show the inversion result using both direct and reflected waves in Figure 8 (the results

377

using the direct waves alone and using the reflected waves alone are shown in Figure S15). The

378

centroid of slip calculated from the finite fault model (white star in Figure 8) is located near the379

center of the assumed fault plane, which coincides with the centroid obtained by the grid search.

17

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth380

We obtain a maximum slip of 0.9 m, and the total seismic moment Mo is 4.34 ×1019 Nm (Mw381

7.03, assuming µ = 40 GPa). The main rupture area, defined as the subfaults with slip larger than382

20 % of the maximum slip (86 % of the total moment is concentrated in this area, marked by383

green lines in Figure 8), has an area of 2000 km2 and the average slip amount within the main384

rupture area is ~ 0.5 m. We calculate the shear stress distribution (Figure S16) from the finite

385

fault model, to estimate an energy-based stress drop (ΔσE) (Noda et al., 2013; Ye et al., 2016)386(the details are described in Text S5). We obtain ΔσE of 1.0 MPa, which is included within the387

uncertainty range obtained by the grid search (Δσ ~ 0.5 – 3.0 MPa), but is at the lower end of388

typical earthquake stress drop values (~ 1 – 10 MPa) (e.g., Abe, 1975; Kanamori & Anderson,389

1975; Ye et al., 2016). The low stress drop of our estimation is consistent with Warren-Smith et390

al. (2018), who suggested that the source dimension is larger than expected from typical fault

391

scaling relationships. However, we note that the estimated stress drop is larger than stress drops392

observed in tsunami earthquakes (Kanamori, 1972; Tanioka & Satake, 1996), which are

393

characterized by extremely low stress drops (< ~ 0.5 MPa) (e.g., Ye et al., 2016). This is394

expected given the intra-slab source of the Te Araora earthquake (tsunami earthquakes are

395

typically generated on the shallow plate interface, not within the slab).

396
397

7. Discussion and Conclusions

398

Ocean bottom pressure gauges (OBPs) deployed in the offshore region clearly showed399

small tsunamis associated with the 2016 Te Araroa earthquake (Mw 7.1). In addition to the400

direct-arrival tsunami with ~−1.5 – 2 cm height, tsunami waveforms with ~ +2 cm height401

reflected from the coast were recorded ~60 min after the direct tsunami arrival. By analyzing the402

reflected tsunamis in addition to the direct tsunami, we improved estimates of the centroid18

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth403

location and fault dimensions of the 2016 Te Araroa earthquake. Although the centroid locations

404

of the GCMT and GeoNet CMT solutions are quite different (~ 80 km apart from each other),405

both the tsunami forward models provided a reasonable fit to the observed direct wave. However,406

we found that the GCMT solution reproduced the reflected waves well, while the GeoNet

407

solution did not. When we used both direct and reflected waves, the horizontal extent of the408

possible centroid location was constrained to within ~40 km (compared to ~100 km for the direct409

wave solution). We also undertook a grid search on the source dimensions (source length and410

width). Using only direct waves, we obtained an optimum source dimension of L = 40 km ± 20

411

km. However, by including reflected waves in the estimate of the source dimension, the source

412

dimension was estimated as L = 50 ± 10 km, which suggests that a source dimension < 40 km

413long is unlikely. Our finite fault inversions suggest that the main slip was concentrated at the414

location of the fault obtained by the grid search for the horizontal fault location. Using the finite

415

fault model, we obtained an energy-based stress drop for this earthquake of 1.0 MPa. This value416

is in the range of typical earthquake stress drops, albeit on the low end of these.

417

Combining direct and coastal reflected waves observed by offshore OBPs, we can estimate418

the centroid horizontal location, fault dimensions, and stress drop of the Te Araroa earthquake

419

more accurately than using direct tsunami wave arrivals and onshore seismic station data alone.

420

The effect of nonlinearity in the OBP waveforms is very small in moderate (M ~ 7) earthquakes

421

because of the small tsunami amplitude at the coast. Coastal reflected tsunami waves have not422

been used previously for estimating earthquake source parameters, although the characteristics of423

coastal reflected waves have been noted by some previous studies (Saito et al., 2013; 2014;424

Gusman et al., 2017; Suppasri et al., 2017). Coastal reflected tsunamis waveforms are very

425

useful to improve source parameter estimation, especially for moderate size offshore earthquakes

19Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth426(M < ~ 7), which are often difficult to constrain because of the small amplitude of the signal at427

onshore seismic stations.

428

Since dense and wide seafloor tsunami observation networks have been constructed in the429

deep ocean (e.g., Kaneda et al., 2015; Kawaguchi et al., 2015; Kanazawa et al., 2016; Uehira et430

al., 2016), it is anticipated that later tsunami arrivals including coastal reflected tsunami signals

431

will be recorded more often from these OBP arrays. In the past, when the details of the432bathymetry were not well-known in many locations, reflected waves could not be simulated433

reliably. For example, Figure S17 demonstrates how low-resolution bathymetry data cannot

434

reproduce the reflected tsunami. We calculate tsunami using the ETOPO1 global bathymetry435

data, with a spatial resolution of 1 arcmin (~ 1.5 km) (Amante & Eakins, 2009) (we use the LLW

436

equation with a grid spacing of Δx = Δy = 2 km and an initial sea surface height distribution

437

obtained from the finite fault inversion; Figure 8). The direct tsunami waves (~ 10 – 40 min)438

agree well with the observed waveforms even with the low-resolution bathymetry data, while the439

reflected tsunami waves (~ 70 – 100 min) do not match the observations, demonstrating the440

dependence of reflected waves on the details of the bathymetry. Tsunami simulations including

441

high-resolution bathymetry data, enable us to simulate the reflected tsunami more reliably than442

the past. This is similar to the advances made in seismology that utilize reflected waves and443

scattered seismic waves (coda waves) to improve understanding of the earthquake source (e.g.,444

Aki, 1969; Umino et al., 1995; Engdahl et al., 1998; Abercrombie, 2013). The analysis of later445

tsunami arrivals holds great promise as a new tool to advance future tsunami research.

446
447

Acknowledgments, Samples, and Data

20

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth448

We thank Editor Yehuda Ben-Zion and two anonymous reviewers for their constructive

449

comments. This study was supported by JSPS KAKENHI 26257206 and partially supported by

450

JP26000002. We greatly appreciate the contribution by the captain and crew of the R/V

451

Tangaroa during, TAN1607, and the ocean bottom observation engineering teams from Kyoto

452

University and Tohoku University to obtain the OBP data. We downloaded GeoNet tide gauge453

data, via IOC (Intergovernmental Oceanographic Commission)'s Sea Level Station Monitoring

454

Facility (http://www.ioc-sealevelmonitoring.org/).We acknowledge support for ship time from455

New Zealand’s Ministry of Business Innovation and Employment (MBIE), and support from456

MBIE’s Endeavour fund program. The figures in this paper were prepared using Generic

457

Mapping Tools (GMT) (Wessel & Smith, 1998). The OBP data used in this study is available in458

the supplementary file (Supplementary Dataset S1).

21

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth459

References

460

Abercrombie, R. E. (2013). Comparison of direct and coda wave stress drop measurements for461

the Wells, Nevada, earthquake sequence. Journal of Geophysical Research: Solid Earth,462

118, 1458–1470. http://doi.org/10.1029/2012JB009638.

463

Aki, K. (1969). Analysis of the seismic coda of local earthquakes as scatted waves. Journal of464

Geophysical Research, 74, 615–631. http://doi.org/10.1029/JB074i002p00615

465

Amante, C., & Eakins, B. W. (2009). ETOPO1 1 arc-minute global relief model: Procedures,

466

data sources and analysis. In NOAA Technical Memorandum NESDIS NGDC-24. Boulder,467

CO: National Geophysical Data Center, NOAA. http://doi.org/10.7289/V5C8276M

468
469

Abe, K. (1975). Reliable estimation of the seismic moment of large earthquakes. Journal of
Physics of the Earth, 23, 381–390. http://doi.org/10.4294/jpe1952.23.381

470

An, C., Cai, C., Zheng, Y., Meng, L., & Liu, P. (2017). Theoretical Solution and Applications of471

Ocean Bottom Pressure Induced by Seismic Seafloor Motion. Geophysical Research

472

Letters, 44, 10272–10281. http://doi.org/10.1002/2017GL075137

473

Engdahl, E.R., van der Hilst, R., & Buland, R. (1998). Global teleseismic earthquake relocation

474

with improved travel times and procedures for depth determination. Bullutin of

475

Seismologcal Society of America, 88, 722–743.

476

Gusman, A. R., Satake, K., Shinohara, M., Sakai, S., & Tanioka, Y. (2017). Fault Slip

477

Distribution of the 2016 Fukushima Earthquake Estimated from Tsunami Waveforms.

478

Pure and Applied Geophysics, 174, 2925–2943. http://doi.org/10.1007/s00024-017-1590-2

479

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

480earthquake parameters. Physics of the Earth and Planetary Interiors, 179, 178–190.

481

http://doi.org/10.1016/j.pepi.2009.11.003

22

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

482
483

Kanamori, H. (1972). Mechanism of tsunami earthquakes. Physics of the Earth and Planetary
Interiors, 6, 346–359. http://doi.org/10.1016/0031-9201(72)90058-1

484

Kanamori, H., & Anderson, D. L. (1975). Theoretical basis of some empirical relations in

485

seismology. Bulletin of the Seismological Society of America, 65(5), 1073–1095.

486

Kanazawa, T., Uehira, K., Mochizuki, M., Shinbo, T., Fujimoto, H., Noguchi, S. … Okada, Y.

487(2016). S-net project, cabled observation network for earthquakes and tsunamis, WE2B–3,

488

Paper presented at SubOptic 2016, SubOptic, Dubai.

489

Kaneda, Y., Kawaguchi, K., Araki, E., Matsumoto, H., Nakamura, T., Kamiya, S … Takahashi,490

N. (2015). Development and application of an advanced ocean floor network system for

491

megathrust earthquakes and tsunamis. In Favali P. et al. (Eds.), Seafloor Observatories (pp.

492643–662). Berlin: Springer

493

Kawaguchi, K., Kaneko, S., Nishida, T. & Komine, T. (2015). Construction of the DONET real-

494

time seafloor observatory for earthquakes and tsunami monitoring. In Favali P. et al. (Eds.),

495

Seafloor Observatories (pp. 211–228). Berlin: Springer.

496

Kubota, T., Hino, R., Inazu, D., Ito, Y., & Iinuma, T. (2015). Complicated rupture process of the497

Mw 7.0 intraslab strike-slip earthquake in the Tohoku region on 10 July 2011 revealed by498

near-field pressure records. Geophysical Research Letters, 42, 9733–9739.

499

http://doi.org/10.1002/2015GL066101

500

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

501

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

tsunami waveforms. Earth and Planetary Science Letters, 458, 241–251.

503

http://doi.org/10.1016/j.epsl.2016.10.047

23

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

504

Kubota, T., Saito, T., Suzuki, W., & Hino, R. (2017b). Estimation of Seismic Centroid Moment505

Tensor Using Ocean Bottom Pressure Gauges as Seismometers. Geophysical Research

506

Letters, 44, 10907–10915. http://doi.org/10.1002/2017GL075386

507
508
509

Lawson, C. L. & Hanson, B.J. (1974). Solving Least Squares Problems. Philadelphia, PA:
Society for Industrial and Applied Mathematics. http://doi.org/10.1137/1.9781611971217
Matsumoto, K., Takanezawa, T., & Ooe, M. (2000). Ocean tide models developed by510

assimilating TOPEX/POSEIDON altimeter data into hydrodynamical model: A global

511

model and a regional model around Japan. Journal of Oceanography, 56, 567–581.

512

http://doi.org/10.1023/A:1011157212596

513

Noda, H., Lapusta, N., & Kanamori, H. (2013). Comparison of average stress drop measures for514

ruptures with heterogeneous stress change and implication for earthquake physics.

515

Geophysical Journal International, 193, 1691–1712. http://doi.org/10.1093/gji/ggt074

516

Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space. Bulletin of

517
518

the Seismological Society of America, 82(2), 1018–1040.
Ristau, J. (2008). Implementation of routine regional moment tensor analysis in New Zealand.

519

Seismological Research Letters, 79(3), 400–415. http://doi.org/10.1785/gssrl.79.3.400

520

Saito, M. (1978). An automatic design algorithm for band selective recursive digital filters.

521
522

Butsuri Tanko, 31, 112–135. (in Japanese).
Saito, T. (2013). Dynamic Tsunami generation due to sea-bottom deformation: Analytical523

representation based on linear potential theory. Earth, Planets and Space, 65, 1411–1423.

524

http://doi.org/10.5047/eps.2013.07.004

525

Saito, T. (2017). Tsunami Generation: Validity and Limitations of Conventional Theories.

526

Geophysical Journal International, 210, 1888–1900. http://doi.org/10.1093/gji/ggx275

24Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth527

Saito, T., Inazu, D., Tanaka, S., & Miyoshi, T. (2013). Tsunami coda across the Pacific Ocean

528

following the 2011 Tohoku-Oki earthquake. Bulletin of the Seismological Society of

529

America. 103(2B), 1429–1443. http://doi.org/10.1785/0120120183

530

Saito, T., Inazu, D., Miyoshi, T., & Hino, R. (2014). Dispersion and nonlinear effects in the 2011

531

Tohoku-Oki earthquake tsunami. Journal of Geophysical Research: Oceans, 119, 5160–

532

5180. http://doi.org/10.1002/2014JC009971

533

Satake, K. (1989). Inversion of tsunami waveforms for the estimation of heterogeneous fault

534

motion of large submarine earthquakes: The 1968 Tokachi-oki and 1983 Japan Sea

535

earthquakes. Journal of Geophysical Research. 94(B5), 5627–5636. http://doi.org/

536

10.1029/JB094iB05p05627

537

Satake, K. (1995). Linear and nonlinear computations of the 1992 Nicaragua earthquake tsunami.

538

Pure and Applied Geophysics, 144, 455–470. http://doi.org/10.1007/BF00874378

539

Suppasri, A., Leelawat, N., Latcharote, P., Roeber, V., Yamashita, K., Hayashi, A., … Imamura,

540

F. (2017). The 2016 Fukushima earthquake and tsunami: Local tsunami behavior and541

recommendations for tsunami disaster risk reduction. International Journal of Disaster

542

Risk Reduction, 21, 323–330. http://doi.org/10.1016/j.ijdrr.2016.12.016

543

Tanioka, Y., & Satake, K. (1996). Fault parameters of the 1896 Sanriku tsunami earthquake

544

estimated from tsunami numerical modeling. Geophysical Research Letters, 23(13), 1549–

545

1552. http://doi.org/10.1029/96GL01479

546

Uehira, K., Kanazawa, T., Mochizuki, M., Fujimoto, H., Noguchi, S., Shinbo, T. … Yamada, T.

547(2016). Outline of Seafloor Observation Network for Earthquakes and Tsunamis along the548

Japan Trench (S-net), EGU2016-13832, Paper presented at EGU General Assembly 2016,

549

European Geosciences Union, Vienna, Austria.

25Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

550
551
552

Um, J., & Thurber, C. (1987). A fast algorithm for two-point seismic ray tracing. Bulletin of the
Seismological Society of America, 77, 972–986.
Umino, N., Hasegawa, A., & Matsuzawa T. (1995). sP depth phase at small epicentral distances553

and estimated subducting plate boundary. Geophysical Journal International, 120, 356–

554

366. doi:10.1111/j.1365-246X.1995.tb01824.x

555

Valentine, A. P., & Trampert, J. (2012). Assessing the uncertainties on seismic source

556

parameters: Towards realistic error estimates for centroid-moment-tensor determinations.

557

Physics of the Earth and Planetary Interiors, 210–211, 36–49.

558

http://doi.org/10.1016/j.pepi.2012.08.003

559

Warren-Smith, E., Fry, B., Kaneko, Y., & Chamberlain, C. J. (2018). Foreshocks and delayed560

triggering of the 2016 MW7.1 Te Araroa earthquake and dynamic reinvigoration of its

561

aftershock sequence by the MW7.8 Kaikōura earthquake, New Zealand. Earth and562

Planetary Science Letters, 482, 265–276. http://doi.org/10.1016/j.epsl.2017.11.020

563

Wells, D. L., & Coppersmith, K. J. (1994). New Empirical Relationships among Magnitude,

564

Rupture Length, Rupture Width, Rupture Area, and Surface Displacement. Bulletin of the565

Seismological Society of America, 84(4), 974–1002.

566

Wessel, P., & Smith, W. H. F. (1998). New, improved version of generic mapping tools released.

567

Eos, Transactions American Geophysical Union, 79, 579.

568

http://doi.org/10.1029/98EO00426

569

Williams, C. A., Eberhart-Phillips, D., Bannister, S., Barker, D. H. N., Henrys, S., Reyners, M.,

570& Sutherland, R. (2013). Revised Interface Geometry for the Hikurangi Subduction Zone,

571

New Zealand. Seismological Research Letters, 84(6), 1066–1073.

572

http://doi.org/10.1785/0220130035

26

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth573

Yabuki, T., & Matsu’ura M. (1992). Geodetic data inversion using a Bayesian information574

criterion for spatial distribution of fault slip. Geophysical Journal International, 109, 363–

575

375. http://doi.org/10.1111/j.1365-246X.1992.tb00102.x

576

Ye, L., Lay, T., Kanamori, H., & Rivera, L. (2016). Rupture characteristics of major and great577(Mw ≥ 7.0) megathrust earthquakes from 1990 to 2015: 1. Source parameter scaling578

relationships. Journal of Geophysical Research: Solid Earth, 121(2), 826–844.

579

http://doi.org/10.1002/2015JB012426

27

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

580

581

Table 1. CMT parameters and fault parameters used for forward calculationa582

Parameter

Global CMT

GeoNetUSGS

Origin time (hh:mm:ss UTC)

16:38:15.8d

16:37:55.440

16:37:57.3

Longitude [°E]179.03179.5200

179.06

Latitude [°S]37.19

36.9795

37.26

Depth [km]

27.8

23

25.5

Seismic moment Mo [Nm]5.05 × 1019

6.35 × 1019

4.410 × 1019

Mw

7.1

7.1

7.0

Strike [°]a354 / 220

354 / 224

351 / 229

Dip [°]a

20 / 70

21 / 76

26 / 75

Rake [°]a−134 / −76−138 / −74−145 / −68

Length L [km]b

45.13

48.71

43.13

Width W [km]b

21.58

22.77

20.91

Slip amount D [m]c

1.30

1.43

1.22

Parameters for two conjugate nodal planes are shown. We used the east-dipping nodal plane(strike ~ 354°) for the calculation.

583

b

584

(1994).

585

c

Slip amount (D) is calculated from D = Mo/(µLW), assuming µ = 40 GPa.

586

dNote that the centroid time is shown.

Length (L) and width (W) were assumed based on the scaling law of Wells and Coppersmith

28Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

587

588

Table 2. Locations of the OBP stations and peak timing of tsunamis.

aStationLongitude[°E]

Latitude
[°S]

Depth
[m]

KU16-2178.8729

38.8465

KU16-3178.7555

KU16-4
KU16-5

Peak timing [min]aDirect wave

Reflected wave

2,138

25.4

93.6

38.8914

1,379

27.0

93.6

178.6609

38.7112

1,047

26.8

91.8

178.8950

38.7208

2,450

24.6

91.9

Peak time is measured from the focal time. These are noted by small dots in Figure 1c.

29

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

177˚
(a)

178˚

180˚

-2000

-2
00

0

GCMT

USGS

-37˚

179˚
A

GeoNet

-10

00
0
00

-2
0

R
Pe auknin umsu ara
la

-1

-38˚

00

B

[km]

-3000

0

10
KU16-5
KU16-4-39˚

20

KU16-3

KU16-20 km
Depth [km]030

100 km

(b)

20

A

B

40
-60

10

-30

0
30
Horizontal Distance [km]

(c)

60

90

BPF (180-3600s)

86KU16-2

4

KU16-3

2

KU16-40

KU16-52hPa~2cm

-2
-20
589

0

20

40

30
60

80

100

Lapse Time from Earthquake [min]120

Aftershock depth

East CapeConfidential manuscript submitted to Journal of Geophysical Research: Solid Earth

590

Figure 1. (a) Location map of this study. The green inverted triangles denote the OBP stations

591

used for the analysis. The blue square is the location of the coastal tide gauge at East Cape. The

592

CMT solutions and the centroid locations from GCMT (red), GeoNet (blue), and USGS (green)593

are shown. The green rectangle denotes the location of the finite fault model provided by USGS.

594

The contour lines of the bathymetry are drawn at 1000 m intervals. Small circles denote the595

aftershock distribution determined by GeoNet, until 14 September 2016. Thin gray contour lines596

show the subducting plate interface from Williams et al. (2013), with 10 km depth intervals. (b)597

Vertical cross section along the A-B line in Figure 1a. The thick black curve is the Hikurangi

598

plate interface from the model of Williams et al. (2013). (c) Pressure time series observed at the599

OBPs. Light and dark gray waveforms denote the raw and moving-averaged pressure time series,

600

respectively. Red waveforms are obtained by applying a bandpass filter to the moving-averaged

601

data. In the application of the bandpass filter, we set the parameters of Saito (1978) as Ap=0.50,602

As=5.00, Fs = 1/120 Hz, Fp = 1/180 Hz, and Fl=1/3600 Hz. Small filled dots denote the peak603

timing of the direct (orange) and reflected tsunamis waveform arrivals (yellow).

31

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

604

32

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth605

Figure 2. (a) Initial tsunami height distributions calculated from the CMT solutions. Background606

colors (red and blue) show the initial tsunami height distribution calculated from the GCMT

607

solution. Colored rectangles and contour lines denote the fault models and tsunami source

608

distribution, calculated from the GCMT (red), GeoNet (blue), and USGS (green) solutions

609(contour intervals are 5 cm). Comparison of the forward-calculated waveforms based on the610

CMT solutions, where trace colors match the corresponding fault models used in (a), at the611

stations (b) KU16-2, (c) KU16-3, (d) KU16-4, (e) KU16-5, and (f) East Cape.

33

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

612
613

Figure 3. Comparison of the observed (black) and calculated waveforms based on the GCMT

614

solutions, with (blue) and without (red) the land. In the calculation that excludes land, we

615

modified the water depth in the shallow region (H < 200 m) to a constant depth (H =200 m). The

616

200 m iso-depth contour from the true bathymetry is noted by the dashed line.

34

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

178˚

179˚

180˚

-37˚

-38˚

KU16-2 0 km

50 km

-39˚
Travel time [min]

617

0

20

40

60

80 100e618

Figure 4. Ray tracing of the coastal reflected tsunami waves, assuming a source at the maximum619

subsidence expected from the GCMT solution. Red star denotes the location of the assumed

620

source. The dashed lines denote the ray paths of the direct tsunami waves. Colors of the rays of621

the reflected tsunami waves denote the travel time.

35

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

622
623
178˚

-36˚

179˚

180˚

-179˚

178˚

179˚

180˚

(a) Direct wave(b) Direct and reflected waves

VRmax=76.3%

VRmax=73.2%

-179˚

VR[%]
80
70

-37˚

60
50
40
30

-38˚

20

0 km

100 km200 km(c) Reflected wave
-36˚

0 km

(d)

VRmax=69.5%

100 km200 km10

Observed
Direct wave
Direct and reflected waves
Reflected wave

KU16-2-37˚

KU16-3

KU16-42hPa~2cm
KU16-5-38˚
0 km

100 km200 km
-20

624

0

20

40

60

80

100

120

Lapse Time from Earthquake [min]625

Figure 5. Results of the grid search for the centroid location, using (a) only the direct waves, (b)626both direct and reflected waves, and (c) only the reflected waves. (a–c) Distribution of the VR

627

between the calculated and observed waveforms. The large red star denotes the location of the628

optimum centroid location. The areas surrounded by the green contour are the solutions with629

relatively high VR (> 90 % of the optimum VR). Locations of the CMT centroids are also shown630

by small stars (GCMT: red, GeoNet: blue, USGS: green). (d) Comparison of the waveforms

36Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth631

between the observed (black) and calculated waveforms from direct waves alone (blue), both632

direct and reflected waves (red), and reflected waves alone (green), for each OBP station. The

633

area shown by the colored rectangles are used for the VR calculation.

37

Confidential manuscript submitted to Journal of Geophysical Research: Solid EarthFault Width [km]80

5

10 15 20 25 30 35 40 45 50

70

VR [%]

60
50
40
30
20

Direct wave
Direct and reflected waves
Reflected wave
10 20 30 40 50 60 70 80 90 100
Fault Length [km]

634
635

Figure 6. Result of the grid search for the fault dimension, using both direct and reflected waves

636(red), direct wave alone (blue), and coastal reflected wave alone (green). Note that tsunami was

637

calculated using the LLW equation. Colored dashed lines denote the 90 % of the optimum VR

638

for each analysis.

38

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

KU16-2
L=100km
L=95km
L=90km
L=85km
L=80km
L=75km
L=70km
L=65km
L=60km
L=55km
L=50km
L=45km
L=40kmL=35km
L=30km
L=25km
L=20km
L=15km
L=10km
L=5km

0.5hPa~0.5cm

Observed

-20

0

20

40

60

80

100

120

639

Lapse Time from Earthquake [min]640

Figure 7. Comparison between the observed (thick black curve) and the calculated waveforms

641

for different earthquake source lengths, at the station KU16-2. The waveform calculated from the642

optimum fault solution using both direct and reflected waves (L = 50 km) is shown in red.

643

39

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

-36˚

178˚

179˚

180˚

-179˚

(a) Direct and reflected waves(b)

Maximum slip = 0.9m
Mo(Slip > 0.2×Maximum slip) = 85.7 %Observed

Calculated

GCMT

-37˚

KU16-2

KU16-3

KU16-4-38˚

2hPa~2cm
KU16-5-39˚

644

KU16-4
KU16-3

KU16-5
KU16-20 km

100 km

0.5

1.0

[m]
0.0

-20

0

20

40

60

80

100

120

Lapse Time from Earthquake [min]645

Figure 8. Result of the finite fault inversion using both direct and reflected waves. (a) Slip

646

distribution on the fault (slip amount is represented by color scale). Contour lines denote vertical

647

displacement calculated from the slip distribution (contour intervals are 5 cm). The star denotes648

the location of the centroid calculated from the fault model. The subfaults with slips larger than649

20 % of the maximum slip are marked by right green lines. (b) Comparison of the waveforms

650

between the observed (black) and calculated waveforms (red).

40

Journal of Geophysical Research: Solid Earth
Supporting Information for
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
Tatsuya Kubota1, Tatsuhiko Saito1, Yoshihiro Ito2, Yoshihiro Kaneko3,
Laura M. Wallace3,4, Syuichi Suzuki5, Ryota Hino5, and Stuart Henrys3
1

National Research Institute for Earth Science and Disaster Resilience, Tsukuba, Japan
2
Disaster Prevention Research Institute, Kyoto Univesity, Uji, Japan3
GNS Science, Lower Hutt, New Zealand4
University of Texas Institute for Geophysics, Austin, Texas, USA
5
Graduate School of Science, Tohoku University, Sendai, Japan

Contents of this file
Text S1 to S5
Figure S1 to S17
Table S1 to S2
Dataset S1 (File is uploaded separately)

Introduction
Text S1 describes the detail of the investigation of the reflected waves. Text S2 describes themethod to estimate the time delay of tsunami arrival and the incident azimuth to the OBP array.
Text S3 describes the evaluation test of the nonlinearity based on tsunami simulation. The method
to calculate tsunami for the estimation of the centroid location and fault dimensions is described in
Text S4. The detail of the finite fault inversion is described in Text S5. Comparison of the coastal tidegauges calculated from the CMT solutions are shown in Figure S1. Figure S2 shows the result of thetsunami simulation based on the USGS finite fault model. The simulation result based on the
conjugate nodal plane from the GCMT solution is given in Figure S3. Figure S4 is the results of thetsunami simulation without the land. Figures S5 and S6 show the snapshots and the maximum
tsunami height based on of the tsunami simulations, respectively. Result of the simulation with themodified bathymetry is in Figure S7. The ray tracing of the reflected waves is shown in Figure S8.
Figure S9 shows the result of the array analysis of OBP data to extract tsunami incident azimuth.
The comparison of the forward calculated waveforms based on LLW and NLL equations is in Figure
S10. The result of tsunami simulation assuming Mw 8 earthquake is in Figure S11. Figure S12 is theresult of ray tracing from the centroid locations obtained by grid search using the direct waves.
Figure S13 shows the comparison between the observed waveform and the calculated waveforms
from the different source length. Figure S14 shows the ABIC values used for the finite fault slip
inversion. Figure S15 shows the results of the finite fault inversion. Distribution of the stress drop on
the fault plane is shown in Figure S16.Tables S1 and S2 gives the time delay of tsunami arrival, forthe direct and reflected waves, respectively. The raw OBP data used in this study is separately
uploaded as Dataset S1.

1

Text S1
In order to investigate the where the coastal reflection radiated from, we show snapshotsof the tsunami height distributions with and without the land (Figures S5a and S5b,
respectively) and the distribution of the maximum absolute height (Figures S6a and S6b). We
find that the main tsunami energy from the source is radiating to the middle part of the East
Coast of the Raukumara Peninsula (Figures S6b). In order to retrieve the tsunami radiation
associated with the coastal reflection, we calculate the difference between the tsunami heightfrom the simulations with and without the land (Figure S5c) and its maximum height
distribution (Figure S6c). We find that tsunami energy is trapped near the middle section of theEast Coast of the Raukumara Peninsula (Figures S5c and S6c). These results suggest that themain energy from the reflected waves concentrates at the middle part of the East Coast of the
Raukumara Peninsula.
We also conducted the additional tsunami calculation by modifying the bathymetry, inorder to confirm that the coastal reflection occurred in the middle section of the East Coast of
the Raukumara Peninsula (Figure S7). In the simulation, we modified the water depth shallower
than 10 m to a constant depth of 10 m, but assigning land to each segment in turn; northern(Figure S7a1), middle (Figure S7a2), or southern (Figure S7a3). The simulation result assumingthat the land exists only at the middle part of Raukumara Peninsula (red lines in Figure S7b2)
reproduced the reflected waves at ~ 90 min well, and the simulations assuming the other partof the Raukumara Peninsula (blue lines in Figure S7b1 and green lines in Figure S7b3) did not
explain the reflected waves at ~ 90 min (although the waves at ~ 70 min are reasonably
explained). This result supports the idea that the coastal reflection mainly occurred at the
middle section of the East Coast of the Raukumara Peninsula.
We then undertook the ray tracing of the coastal reflected tsunami waves (Figures 4 andS8). In the ray tracing, we assumed the reflection points at the middle part of the East Coast of
the Raukumara Peninsula. In order to calculate the ray path, we used the pseudo-bending
method proposed by Um and Thurber (1987). In this approach, the starting point and theending point of the ray path are assumed, and the optimum ray path is searched so that thetravel time between two points is minimized. In the calculation, we separately calculated theray path: one is between the source and the coastal reflection point, and the other is between
the reflection point and the OBP. We assumed three source locations, (1) at the maximum
subsidence of the initial sea surface height distribution calculated based on the GCMT solution(Figures 4 and S7a), (2) at the GCMT centroid (Figure S7b), and (3) at the maximum uplift basedon the GCMT solution (Figure S7c). Most of the ray paths are almost the same regardless of theassumed source location, but the travel times take slightly longer when the source is assumed
at the maximum subsidence (~ 85 – 90 min). Considering that the main phases of the reflected
waves are observed at ~ 90 min, the reflected waves are mainly generated by the subsidence oftsunami source. Furthermore, we also calculated the ray paths directly propagating from thesource to the OBP (dashed lines in Figure S5). The travel time calculated of the direct tsunami is
25.4 min when the source is assumed at the maximum subsidence, which is consistent with theobserved travel time (Table 2).

2

Text S2
The peak timing of direct and reflected tsunami waveform arrivals (Table 2, noted bysmall dots in Figure 1c) are slightly different between the OBP stations. By using the time delay
between the two OBPs, we can estimate the approximate incident azimuth of the direct andreflected tsunamis. We calculated the time delay of tsunami arrivals between the two OBPs by
calculating the cross correlation between ith and jth OBPs (CCij), defined as:
CC"# =- %& '() %* '()+,- %& '().

- %* '()+,

.

,

(S1)

where /" 0 denotes the waveform of the ith OBP at 0 = 120, 20 is the data sampling interval (1s), and 3 is the lag time between the two waveforms. We estimated the time delay between the456
ith and jth stations (3"#) for the direct wave (10 – 40 min, Table S1) and the reflected wave (70– 100 min, Table S2), by calculating the lag time 3 which maximize the CCij.
As a result, the first direct tsunami arrival (~ 25 min) occurs at the northeastern-moststation KU16-5, and then arrives at station KU16-2 < 1 min later. Tsunami arrivals are recorded
last at the western stations KU16-3 and KU16-4 almost at the same time, ~ 2 min after the firstarrival at KU16-5. This suggests that the first tsunami propagated from the northeast/eastnortheast of our OBP array. The secondary tsunami arrives first at KU16-5, and is earlier by ~ 1.5
min at the northern stations (KU16-4 and KU16-5) compared to the southern stations (KU16-2
and KU16-3). This suggests that the secondary tsunamis propagated from the north-northeast.
We also show the time delay estimated from the timing of the secondary tsunami in Table S2.
The lag times estimated from the direct waves are almost identical to those calculated from thepeak timing, but those from the reflected waves are different between the cross correlation andpeak timing. This suggests that the reflected waves are more complex than the direct waves,probably because the complex coastal bathymetry distort the waveforms in addition to thetravel time change.
After estimating the time delay, we then estimate the incident azimuth of the tsunami
arrivals based on array analysis. Based on the plane wave approximation, the time delay
789
between the ith and jth stations (3"#: ) can be expressed as:;<=
3"# : = : ∙ ?# − ? " ,(S2)
where : = AB , AD is the slowness vector and ? " = E" , F" is the location of the ith station. We
search for the optimum slowness vector, which minimizes the root mean square (RMS) misfit of
the time delay (RMS(s)), defined as:
RMS : =

J

LMN
789
3"#
− 3"#
:

"P#

O

,

(S3)

K
456
where 3"#is the time delay obtained from the observed tsunami waveform (Tables S1 and S2).

The estimated incident azimuths for the direct and reflected tsunamis are 71° and 29°,
respectively (Figure S9), which is consistent with those expected from the time delay. We also
calculate the water depth expected from the optimum slowness, based on the linear long wave
theory:Q=
where A = : =

R.
ST

=

J
ST N .

,(S4)

ABO + ADO . Using the direct wave, we obtain H = 1771 m, which is almost

equivalent to the average water depth around the OBP array. We also obtain H = 4360 m fromthe reflected waves, which is much deeper than the average depth. This is probably due to thecomplexity of the reflected waves, including the travel time change caused by bathymetry
changes along the coast.

3

Text S3
In general, nonlinearity should be considered to accurately model tsunami propagationnear the coast (e.g., Satake, 1995; Saito et al., 2014), whereas nonlinear effects are very small
when tsunami propagate in the deep ocean. Therefore, we additionally conducted tsunami
calculation considering the nonlinearity. We used the nonlinear long wave (NLL) equation (e.g.,
Satake, 1995; Saito et al., 2014), imposing a bottom friction with Manning’s coefficient of n =0.03 (m−1/3 s). We also included the tsunami inundation in our calculation, as a moving boundary
condition between the land (dry cell) and sea (wet cell) The other details of the NLL tsunami
calculation scheme is described in Saito et al. (2014). The other settings were all the same as
the simulation using the LLW equation.
We show a comparison of the calculated tsunami waveforms using the NLL equations inFigure S10. The waveforms from the LLW (red) and NLL (blue) equations are very similar, notonly for the direct waves but also for the reflected waves. In order to evaluate the agreement
between the model waveforms, we calculate the variance reduction (VRcal), using the following
equation:
VR ;<= = 1 −[ X
Z .
& %& Y%&[
X
Z& %& × %&×100 (%),(S5)

where /"b and /"c is the i-th data of the calculated waveforms A and B, respectively, and N
denotes the number of data used for the calculation. VRcal for four OBP waveforms between
from LLW equation and from NLL equations using a time window of 0 – 60 min and 60 – 120
min are 100 % and 98.5 %, respectively. This suggests that the nonlinearity is small in thissituation.
We also conducted additional tsunami simulations using a finer grid spacing (Δx = Δy =250 m) compared to the original calculation (Δx = Δy = 500 m) (Figure S10). As a result, even if
we use the finer grid spacing of 250 m, the difference between the calculated waveforms using
the LLW and NLL equations is small (i.e., the nonlinearity is very small). We also compare thesimulation result using a grid spacing of 500 m with that using a finer grid spacing (250 m).
Overall, the calculated OBP waveforms are similar between the simulations using the grid
spacing of 500 m and the finer 250 m spacing, although the waveforms at the coastal station,
East Cape, varied slightly. This suggests the offshore OBP waveforms are less affected by thecomplex coastal bathymetry, than the coastal tide gauges.
We find that the effect from nonlinearity is small in the reflected waves observed by the
OBP array. This is probably because of the moderate magnitude and small tsunami amplitude
at the coast. To illustrate this, we conduct additional tsunami simulations, by assuming anearthquake with larger magnitude (Mw 8.0) which generates a much larger tsunami (FigureS11). In this calculation, we assumed the magnitude of Mw 8.0, and estimated the fault length
and width, and slip amount, D, based on the scaling law (L = 132 km, W = 46km, D = 5.2 m). We
used the strike, dip, rake and centroid location of the GCMT solution (Figure S11a). As a result,the direct waves are the same between those calculated by using the LLW (red) and NLL (bluein Figure S11b) equations, because the direct waves propagate in the deep ocean. But thereflected waves, which propagate near the shallow area, are slightly different. The VRcal
between the waveforms calculated by the LLW and NLL equations, using the time window of 0
– 60 and 60 – 120 min was 100 % and 92.9 %, respectively. VRcal for the reflected waves is lower
than that from the original (Mw 7.1) result, suggesting that effects of nonlinearity at the coastcannot be neglected in the larger earthquake.

4

Text S4
For the estimation of the centroid location and the fault dimensions, the tsunami is
calculated based on a linear superposition of pre-computed tsunami Green’s functions fromsmall source elements for seafloor displacement (unit tsunami sources). This method is
proposed by Kubota et al. (2015), in order to minimize the tsunami calculation time. We
distributed 39× 39 unit sources in the area 400 km × 400 km around the focal area. We use
pyramid-like-shaped unit sources with dimensions of 20 km × 20 km, identical to those used in
Kubota et al. (2015). Each has a horizontal spacing of 10 km, and overlaps with the neighboring
ones. After the calculation of tsunami from each unit source, we apply a bandpass filter similar
to that applied to the observed records, to obtain the tsunami Green’s function. Aftercalculating the seafloor deformation using the equations of Okada (1992), the deformation isexpressed by a linear superposition of the unit sources. Then tsunami was calculated by a linear
superposition of the tsunami Green’s function.

5

Text S5
In the finite fault slip inversion, we assumed the rectangular planar fault with the size of
80 km length and 70 km width, which is determined based on the optimum fault dimension
obtained by the grid search analysis. The center of the rectangular planar fault coincides withthe point 178.97°E, 37.10°S, and 27.8 km depth, which is the optimum fault center location
obtained by the grid-search of the fault horizontal location. We used the fault geometry fromthe GCMT solution (strike = 354°, dip = 20°, and rake = –134°). We calculated the seafloor
vertical deformation from each subfault by using the equation of Okada (1992). We assumed
that the sea surface vertical displacement is equal to the seafloor vertical deformation. Tsunami
Green’s functions for the finite slip inversion from each subfault was calculated by thesuperposition of the tsunami Green’s function of the tsunami source (described in Text S4).
In the finite fault inversion, we assumed the observed tsunami waveforms can be
expressed as a superposition of the tsunami Green’s function from each subfault. The
observation equation is expressed as:

f
d
=
i,
αh
e

(S6)

where d is a vector composed of the tsunami data, G is a matrix composed of the tsunami
Green’s function for the finite slip inversion, and u is a vector representing the slip distributionof the fault, which are we want to obtain. A matrix S denotes the smoothing constraint of thefault slip, and α is the weighing of the smoothing constraint.
The weighting of the smoothing constraint was set as α = 5, which was determinedbased on the value of Akaike’s Bayesian Information Criterion (ABIC) (Yabuki and Matsu’ura,
1992) (Figure S14).
We also calculated the location (longitude and latitude) of the centroid (xg = (xg, yg)) fromthe finite fault model, using the following equation:?S =K j & ×8&
"kJ 8 ,

(S7)&

where N denotes the number of the subfaults included in the main rupture area, and xi = (xi, yi)
is the location of the center of the ith subfault included in the main rupture area and ai is the slip
amount of the ith subfault.
We then calculated the shear stress (i.e., the stress drop) distribution, using the equationof Okada (1992), by computing the shear stress change at the center of each subfault (Figure
S16). We also calculated the energy-based stress drop ΔσE (e.g., Noda et al., 2013; Ye et al.,2016), expressed as:

Δmn = s(op (%p qr
s(%p qr,(S8)

where Δσ1 and Δu1 are the distributions of the stress drop and slip amount on the slip area. Thisequation can be simplified as:

Δmn =

[
&tp (o& %&
[ % ,
&tp &

(S9)

where Δσi and ui are the stress drop and slip amount of ith subfault. Using the stress drop andslip amount of all subfaults, we obtain ΔσE as 0.8 MPa. If we use the subfaults inside the main
rupture area (marked by green rectangles in Figure 8), we obtain ΔσE as 1.0 MPa.

6

174˚

177˚

North Cape

Observed GCMT GeoNet USGS(b)
North Cape
-400

0

0
-10 -2000
00
00
-2

-36˚

-177˚

00

-2000

180˚
-30

171˚
-33˚
(a)

Great Barrier Is

Great Barrier Is
Tauranga

-200
0
00
East
-30Cape

East Cape-1000

Gisborne

Tauranga-39˚

Gisborne
Port Napier

Port of Napier

Castlepoint

20cm

300

0

Castlepoint

-42˚
[cm]
-10

-5

0

5

10

-20

0

20

40

60

80

100

120

Lapse Time from Earthquake [min]Figure S1. (a) Location map of the coastal tide gauges. Background colors (red and blue) show the initial tsunami height distribution from the
GCMT solution. Colored rectangles denote the fault models calculated from the CMT solutions. The black dashed lines denote the area shownin Figure 2a. (b) Comparison of the forward-calculated waveforms of the coastal tide gauges (trace colors match the corresponding fault
models used in (a)).

7

177˚
-36˚
(a)

178˚

179˚

180˚

(b)

Observed
USGS FAULT

00

-40

GCMT CMTUSGS CMT

-20
0

0

0

-20
00

00
-3

KU16-2-2000

KU16-3

-3
0

-1

00
0

East Cape

00

-37˚

-179˚

KU16-4-38˚2hPa ~ 2cm
KU16-5

KU16-2

-39˚
[cm]
-10

-5

0

5

10

-20

0

20

40

60

80

100

120

Lapse Time from Earthquake [min]Figure S2. (a) Initial tsunami height distributions calculated from the finite fault model from the teleseismic analysis by USGS
(https://earthquake.usgs.gov/earthquakes/eventpage/us10006jbi#finite-fault). Back ground colors (red and blue) show the initial tsunami
height distribution calculated from the USGS finite fault model. The black rectangle and contour lines denote the fault models and tsunami
source distribution (contour intervals are 5 cm). (b) Comparison of the forward-calculated waveforms from the rectangular model based on the
GCMT solution (red), the USGS finite fault model (blue), and the rectangular model based on the USGS solution (green), with the observed
waveform (black).

8

177˚
-36˚
(a)

178˚

179˚

180˚

NE-dipping

(b)
-40

NW-dipping

00

-20
0

0

0

-20
00

00
-3

KU16-2-2000

KU16-3

-3
0

-1

00
0

East Cape

00

-37˚

-179˚

KU16-4-38˚2hPa ~ 2cm
KU16-5

KU16-2

-39˚
[cm]
-10

-5

0

5

10

-20

0

20

40

60

80

100

120

Lapse Time from Earthquake [min]Figure S3. (a) Initial tsunami height distributions calculated from the GCMT solution, assuming the northwest-dipping nodal plane (strike =220°). Background colors (red and blue) show the initial tsunami height distribution. The black and red rectangles denote the fault models
calculated from the northwest-dipping (strike = 220°) and north-east dipping (strike = 354°) planes, respectively. (b) Comparison of the forwardcalculated waveforms from the northwest-dipping (blue) and northeast-dipping (red) nodal planes.

9

177˚
-36˚

178˚

179˚

180˚

-179˚

(a)

Observed
10m
200m

(b)

-37˚

With land
50m
500m

100m

KU16-2-50-500-100 0-1-200

KU16-3

KU16-4-38˚

2hPa~2cm
KU16-5-39˚

KU16-4

KU16-5

KU16-3

KU16-2
[cm]

-10

-5

0

5

10

-20

0

20

40

60

80

100

120

Lapse Time from Earthquake [min]Figure S4. Comparison of the observed (gray) and calculated waveforms based on the GCMT solutions, with (black) and without the land. In thecalculation that excludes land, we modified the water depth in the shallow region to a constant depth of 10 m (red), 50 m (light blue), 100 m(pink), 200 m (blue), and 500 m (green). The iso-depth contours from the true bathymetry are noted by the dashed lines.

10(a) With land
t=0.0mint=15.0min

t=30.0min

t=45.0min

t=60.0mint=75.0min

t=90.0min(b) Without land (H0=10m)
t=0.0mint=15.0min

t=30.0min

t=45.0min

t=60.0mint=75.0min

t=90.0min

[cm]

10

5

(c) Differencet=0.0min0
t=15.0min

t=30.0min

t=45.0min

t=60.0mint=75.0min

t=90.0min-5-10

Figure S5. Snapshots for the tsunami height from (a) the simulation with the original bathymetry data, (b) the simulation without the land (<200 m), and (c) the difference between the simulation with and without the land. Contours are drawn at 0.5 cm intervals.

11

(a) With land(b) Without land (H0=10m)(c) Difference
[cm]
10
8
6
4
2
0

Figure S6. Distributions of the absolute maximum value of tsunami height, calculated from (a) the simulation with the original bathymetrydata, (b) the simulation with the modified bathymetry data (water depth shallower than 200 m is modified to a constant depth of 200 m), and(c) the difference between these two simulations. The contours are drawn at 0.5 cm intervals.

12

177˚
-36˚

178˚

179˚

180˚

-179˚

(a1)

(b1)

-37˚

Observed

With land

North

KU16-2

KU16-3

KU16-4-38˚

2hPa~2cm
KU16-5-39˚

KU16-4

KU16-5

KU16-3

KU16-2

-20

[cm]

177˚
-36˚

178˚

-10 179˚
-5
0

5180˚
10

-179˚

(a2)

0

40

60

80

100

120

Lapse Time from Earthquake [min]
(b2)

-37˚

20

Observed

With land

MiddleKU16-2

KU16-3

KU16-4-38˚

2hPa~2cm
KU16-5-39˚

KU16-4

KU16-5

KU16-3

KU16-2

-20

[cm]

177˚
-36˚

178˚

-10 179˚
-5
0

5180˚
10

-179˚

(a3)

0

40

60

80

100

120

Lapse Time from Earthquake [min]
(b3)

-37˚

20

Observed

With landSouth

KU16-2

KU16-3

KU16-4-38˚

2hPa~2cm
KU16-5-39˚

KU16-4

KU16-5

KU16-3

KU16-2
[cm]

-10

-5

0

5

10

-20

0

20

40

60

80

100

120

Lapse Time from Earthquake [min]Figure S7. Result of the tsunami calculation with the modified bathymetry, assigning the land
to the (a1) northern, (a2) middle, and (a3) southern segments of the Raukumara Peninsula. The
10 m iso-depth contour are noted by the dashed line. The black lines denote the area with theland assigned. The other segments are modified to a constant water depth of 10 m. (b1, b2, b3)
Comparison of the observed (gray) and calculated waveforms based on the GCMT solutions,using the original bathymetry (black) and the modified bathymetry (colored).

13

178˚

179˚

180˚178˚

(a)

179˚

180˚178˚

(b)

179˚

180˚

(c)

-37˚

-38˚

KU16-2 0 km

50 km

KU16-2 0 km

50 km

KU16-2 0 km

50 km

-39˚
Travel time [min]
0

20

40

60

80 100Travel time [min]
0

20

40

60

80 100Travel time [min]
0

20

40

60

80 100Figure S8. Ray tracing of the coastal reflected tsunami waves, assuming a source at (a) the maximum subsidence expected from the GCMT
solution, (b) the GCMT centroid, and (c) the maximum uplift expected from the GCMT solution. Red star denotes the location of the assumed
source. The dashed lines denote the ray paths of the direct tsunami waves. Colors of the rays of the reflected tsunami waves denote the travel
time from the source.

14

0.3

y slowness [min/km]

0.2(a) Direct wave

0.3

Optimum RMS = 0.09min
Optimum azimuth = 71°
Optimum slowness = 0.13min/km
Optimum velocity = 7.9km/min
Expected water depth = 1771m

0.2

0.1

0.1

0.0

0.0
0.1

-0.1

(b) Reflected wave
Optimum RMS = 0.15min
Optimum azimuth = 29°
Optimum lowness = 0.08min/km
Optimum velocity = 12.4km/min
Expected water depth = 4360m

0.1

-0.1

0.2
-0.2
-0.3
-0.3

0.2
0.0 0.2 0.4 0.6 0.8 1.0

0.3

RMS [min]
-0.2

-0.1

0.0

0.1

x slowness [min/km]

0.2

-0.2

-0.3
0.3
-0.3

0.0 0.2 0.4 0.6 0.8 1.0

0.3

RMS [min]
-0.2

-0.1

0.0

0.1

0.2

0.3

x slowness [min/km]

Figure S9. Distribution of the RMS(s) for (a) direct wave and (b) reflected waves. Color and solid contour lines denote the RMS(s) value. Thinand thick contour lines are drawn by 0.1 and 1 min intervals, respectively. Dashed circles denote the iso-slowness lines with 0.1 min/km
intervals.

15

177˚
-36˚
(a)

178˚

179˚

180˚

-179˚
(b) KU16-2

Direct wave

00

-20

00

0
00
-3
-20

00

Reflected wave

250m

-2000

East Cape

NLL

-3

-1
0

00
0

-37˚

Observed LLW

00
-40

500m

-38˚
3hPa~3cm
KU16-4

KU16-5

KU16-3

KU16-2

-39˚

-20

[cm]
-10

-5

0

5

10

0

20

Direct wave

Reflected wave

Direct wave
250m

500m

500m
3hPa~3cm

20

40

60

80

100

120 -20120

Reflected wave

0

20

40

60

80

100

120

Lapse Time from Earthquake [min](e) KU16-5(f) East Cape

Direct wave

Reflected wave

250m

250m

500m

500m
3hPa~3cm

0100

3hPa~3cm

Lapse Time from Earthquake [min]

-20

80

(d) KU16-4

250m

0

60

Lapse Time from Earthquake [min](c) KU16-3

-20

40

20

40

60

40cm

80

100

Lapse Time from Earthquake [min]120 -20

0

20

40

60

80

100

120

Lapse Time from Earthquake [min]Figure S10. Comparison of the forward calculated waveforms based on the Global CMT
solution using the LLW (red) and NLL (blue) equations, with the computational grids of 500 and
250m.

16

177˚
-36˚
(a)

178˚

179˚

180˚

(b) Modified-GCMT (Mw 8.0)
LLW

-40

00

NLL

KU16-2

-20
0

0

0

-20
00

00
-3

-2000

-3
0

-1

00
0

East Cape

KU16-3
00

-37˚

-179˚

KU16-4-38˚

15hPa~15cm
KU16-5-39˚

KU16-4

KU16-5

KU16-3

KU16-2[cm]-50 -25

0

25

50

-20

0

20

40

60

80

100

120

Lapse Time from Earthquake [min]Figure S6S11. (a) Tsunami source distribution calculated from the modified GCMT solutions, which has a magnitude comparable to an Mw 8.0.
(b) Comparison of the forward calculated waveforms based on the CMT solutions at the OBPs from LLW (red) and NLL (blue) equations. The
observed pressure waveforms associated with the Te Araroa earthquake is also shown by black lines.

17

-36˚

178˚

179˚

180˚

-179˚

100 km200 km-37˚-38˚

KU16-2

-39˚
20

25

30

0 km

[min]
Figure S12. Result of the ray tracing from the solutions in the high-variance reduction (VR)area, obtained by the direct waves alone (small green circles) to the OBP station KU16-2. The
red star denotes the optimum solution obtained from the direct waves.

18

(a) KU16-3(b) KU16-4

L=100km
L=95km
L=90km
L=85km
L=80km
L=75km
L=70km
L=65km
L=60km
L=55km
L=50km
L=45km
L=40kmL=35km
L=30km
L=25km
L=20km
L=15km
L=10km
L=5km

L=100km
L=95km
L=90km
L=85km
L=80km
L=75km
L=70km
L=65km
L=60km
L=55km
L=50km
L=45km
L=40kmL=35km
L=30km
L=25km
L=20km
L=15km
L=10km
L=5km

0.5hPa~0.5cm

Observed

-20

0

20

40

60

80

0.5hPa~0.5cm

Observed

100

120 -20

Lapse Time from Earthquake [min]020(d) East Cape

L=100km
L=95km
L=90km
L=85km
L=80km
L=75km
L=70km
L=65km
L=60km
L=55km
L=50km
L=45km
L=40kmL=35km
L=30km
L=25km
L=20km
L=15km
L=10km
L=5km

L=100km
L=95km
L=90km
L=85km
L=80km
L=75km
L=70km
L=65km
L=60km
L=55km
L=50km
L=45km
L=40kmL=35km
L=30km
L=25km
L=20km
L=15km
L=10km
L=5km

0.5hPa~0.5cm

-20

0

20

40

60

80

Lapse Time from Earthquake [min]

80

100

120

5cm

Observed

100

60

Lapse Time from Earthquake [min](c) KU16-5

Observed

40

120 -20

0

20

40

60

80

100

120

Lapse Time from Earthquake [min]Figure S13. Comparison between the observed (thick black curve) and the calculated
waveforms from the different source length, at the stations (a) KU16-3, (b) KU16-4, (c) KU16-5,
and (d) East Cape. The waveforms calculated from the optimum fault solution using both direct
and reflected waves (L = 50 km) are shown in red.

19

1500

Direct wave
Direct and reflected waves
Reflected wave

ABIC

1000

500

0
0.1

1

10

Smoothing weight α
Figure S14. Values of ABIC for the finite fault inversion using the direct waves alone (blue), thereflected waves alone (green), and both direct and reflected waves (red), plotted as a functionof α. We adopted the smoothing weight of α = 5, denoted by the dashed line.

20Figure S15. Results of the finite fault inversion, using (a) the direct waves alone, (b) both direct
and reflected waves, and (c) the reflected waves alone. (a–c) Slip distribution on the fault (slipamount is represented by color scale). Contour lines denote vertical displacement calculated
from the slip distribution (contour intervals are 5 cm). The star denotes the location of the
optimum centroid location obtained by the grid-search for the horizontal fault location. The
subfaults with slips larger than 20 % of the maximum slip are noted by right green lines. The
GCMT solution is also shown. (d) Comparison of the waveforms between the observed (black)and calculated waveforms from direct waves alone (blue), both direct and reflected waves (red),and reflected waves alone (green), for each OBP station. The area shown by colored rectangle
is used for the inversion.

21

Figure S16. Distribution of shear stress change calculated from the finite fault model using
both the direct and reflected waves (Figure 8). Positive value means the shear stress was
reduced after the earthquake. The subfaults with slips larger than 20 % of the maximum slip arenoted by right green lines.

22

Figure S17. Comparison of the observed (black) and calculated waveforms based on the result of the finite fault inversion (Figure 8), with the
NIWA bathymetry data (red) and ETOPO1 bathymetry data (blue).

23

$%&
Table S1. Time delay of tsunami arrival (!"#
) for the direct waves
Time delay from crossStation i
Station jcorrelation [min]aTime delay from peaktiming [min]b

KU16-5

KU16-20.72

0.75

KU16-5

KU16-3

2.28

2.35

KU16-5

KU16-42.30

2.13

KU16-2

KU16-3

1.60

1.60

KU16-2

KU16-4

1.60

1.38

KU16-3

KU16-4

0.15

−0.22

a

Time delay is measured by cross correlation (CCij). Positive time delay denotes that tsunami
arrives earlier at station i.
bTime delay is measured by peak timing of tsunami, shown in Table 2. Positive time delay
denotes that tsunami arrives earlier at station i.

$%&
Table S2. Time delay of tsunami arrival (!"#
) for the reflected waves
Time delay from crossTime delay from peak
Station i
Station jcorrelation [min]atiming [min]b

KU16-5

KU16-4

0.78

−0.15

KU16-5

KU16-2

1.18

1.67

KU16-5

KU16-3

1.88

1.65

KU16-4

KU16-2

0.07

−1.82

KU16-4

KU16-3

1.25

1.80

KU16-2

KU16-3

0.65

1.60

a

Time delay is measured by cross correlation (CCij). Positive time delay denotes that tsunami
arrives earlier at station i.
bTime delay is measured by peak timing of tsunami, shown in Table 2. Positive time delay
denotes that tsunami arrives earlier at station i.

24