Ohara et al.
Progress in Earth and Planetary Science(2023) 10:35
https://doi.org/10.1186/s40645-023-00565-z

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open AccessComplex evolution of the 2016 Kaikoura
earthquake revealed by teleseismic body waves
Kenta Ohara1, Yuji Yagi2* , Shinji Yamashita1, Ryo Okuwaki2,3, Shiro Hirano4 and Yukitoshi Fukahata5

Abstract
The 2016 Kaikoura earthquake, New Zealand, ruptured more than a dozen faults, making it difficult to prescribe
a model fault for analyzing the event by inversion. To model this earthquake from teleseismic records, we used
a potency density tensor inversion, which projects multiple fault slips onto a single model fault plane, reducing
the non-uniqueness due to the uncertainty in selecting the faults’ orientations. The resulting distribution of potencyrate density tensors is consistent with observed surface ruptures. In its initial stage, the rupture propagated northeastward primarily at shallow depths. Later, the rupture propagated northeastward at greater depths beneath a gap
in reported surface ruptures. The main rupture phase started in the northeastern part of the Kekerengu fault after 50 s
and propagated bilaterally to the northeast and southwest. The non-double-couple component grew to a large
fraction of the source elements as the rupture went through the junction of the Jordan Thrust and the Papatea fault,
which suggests that the rupture branched into both faults as it back-propagated toward the southwest. The potency
density tensor inversion sheds new light on the irregular evolution of this earthquake, which produced a fault rupture
pattern of unprecedented complexity. Our source model of the 2016 Kaikoura earthquake (e.g., back-rupture propagation) could prompt research to determine a more realistic model with segmented faults using near-field data.
Keywords Earthquake dynamics, Waveform inversion, Body waves, Earthquake source observation
1 Introduction
On 13 November 2016, the moment magnitude (Mw)
7.8 Kaikoura earthquake struck in the South Island of
New Zealand near the boundary between the Pacificand Australia plates (Fig. 1a) (Dziewonski et al. 1981;
Ekström et al. 2012). Field studies reported that theearthquake produced a complex set of surface ruptures*Correspondence:
Yuji Yagi
yagi-y@geol.tsukuba.ac.jp
1
Graduate School of Science and Technology, University of Tsukuba,
Tsukuba, Ibaraki 305‑8572, Japan
2
Faculty of Life and Environmental Sciences, University of Tsukuba,
Tsukuba, Ibaraki 305‑8572, Japan3
Mountain Science Center, University of Tsukuba, Tsukuba, Ibaraki
305‑8572, Japan
4
College of Science and Engineering, Ritsumeikan University, Kusatsu,
Shiga 525‑8577, Japan
5
Disaster Prevention Research Institute, Kyoto University, Uji, Kyoto
611‑0011, Japanof more than 12 faults (Hamling et al. 2017; Stirling
et al. 2017; Litchfield et al. 2018). The surface rupture area, extending a total length of ∼ 165 km, canbe divided into south and north sections separated
by a gap of about 30 km with no mapped surface ruptures between the northeast end of the Conway-Charwell fault and the southwest end of the Manakau fault(Litchfield et al. 2018) (Fig. 1b). The south section
involved the Humps fault and the Conway-Charwell
fault with mixed dextral and reverse faulting (Litchfield et al. 2018) (Fig. 1b). The north section displayed
a linear set of surface ruptures with mixed vertical anddextral displacements on the Manakau fault, the Upper
Kowhai fault, the Jordan Thrust, the Kekerengu fault,
and the Needles fault (Litchfield et al. 2018). In addition, surface rupture with mixed sinistral and reverse
offsets occurred on the west-dipping Papatea fault,
which extends southward nearly orthogonal to the linear rupture set near the junction of the Kekerengu fault© The Author(s) 2023. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Ohara et al. Progress in Earth and Planetary Science(2023) 10:35

Fig. 1 a Tectonic setting of the study region. The dashed lines
represent the plate boundary (Bird 2003). The arrow denotes
the plate motion of the Pacific plate relative to the fixed Australia
plate in NUVEL 1A (DeMets et al. 1994). The star marks the mainshock
epicenter (Lanza et al. 2019). b Seismotectonic summary of the study
region of the 2016 Kaikoura earthquake. The left and right beachball shows the obtained total moment tensor and the Global
Centroid Moment Tensor (Dziewonski et al. 1981; Ekström et al. 2012)
for the mainshock, respectively. Black dots represent aftershocks
during the week after the mainshock (Lanza et al. 2019). Gray, orange,
blue, and green lines indicate surface ruptures of the 2016 Kaikoura
earthquake from the New Zealand Active Faults Database (Langridge
et al. 2016). The black line represents the assumed model plane.
Background contours display topography/bathymetry (Mitchell
et al. 2012). HmF–Humps fault zone, CCF–Conway-Charwell fault,
MF–Manakau fault, UKF–Upper Kowhai fault, JT–Jordan Thrust, PF–
Papatea fault, KF–Kekerengu fault, NF–Needles fault

and the Jordan Thrust (Litchfield et al. 2018) (Fig. 1b).
Aftershocks were distributed throughout the zone ofsurface ruptures (Lanza et al. 2019) (Fig. 1b).
The Global Centroid Moment Tensor (GCMT) solution
for the mainshock (Dziewonski et al. 1981; Ekström et al.
2012) indicates oblique reverse faulting (Fig. 1). A multiple-point-source inversion using the records of longrange seismographs (teleseismic waveforms) detects four
subevents, consisting of three oblique strike-slip subevents and one thrust subevent (Duputel and Rivera 2017),
indicating that the earthquake ruptured multiple faults
with different faulting mechanisms. Finite-fault inversions using seismic data alone (Bai et al. 2017; Hollingsworth et al. 2017; Zhang et al. 2017; Zheng et al. 2018) or
using both seismic and geodetic data (Cesca et al. 2017;
Holden et al. 2017; Wang et al. 2018b) commonly find theinitial rupture episode during the first ~ 60 s, followedby the main rupture episode, and discuss the rupture
propagated toward northeast from the epicenter in both
episodes. Notably, the field surveys identify large displacements with sinistral and reverse offsets at the Papatea fault (Clark et al. 2017; Hamling et al. 2017; Stirling
et al. 2017; Litchfield et al. 2018).

Page 2 of 12

However, finite-fault inversions using only teleseismic
body waves, which can estimate the overall rupture propagation process during an earthquake, have not identified
subevents with focal mechanism corresponding to that
Papatea fault rupture (Bai et al. 2017; Hollingsworth et al.
2017; Zhang et al. 2017).
Finite-fault inversions in previous studies estimatedthe rupture process under the assumption that the rupture unilaterally propagates northeastward (Bai et al.
2017; Cesca et al. 2017; Holden et al. 2017; Hollingsworth et al. 2017; Zhang et al. 2017; Wang et al. 2018b;
Zheng et al. 2018). Such small constraints on the unilateral rupture scenario may not always be appropriate
for the earthquake modeling in a complex fault zone,
which sometimes involves irregularity, e.g., the rupturein the backward direction as a part of its bilateral propagation that is initiated as a secondary rupture episode(e.g., Gallovič et al. 2020; Yamashita et al. 2022a; Yagi
et al. 2023). The assumption of unilateral northeastward
rupture propagation can make the interpretation of theinversion results more difficult. Indeed, it is difficult to
explain how a right-lateral strike-slip rupture propagating in a northeast direction along the Jordan Thrust
could have triggered a reverse fault rupture on the Papatea fault, situated in the extensional quadrant of the focal
mechanism of the Jordan Thrust. Therefore, there shouldstill be a room to investigate whether alternatives to theunilateral rupture scenario proposed by finite-fault inversions exist.
As the 2016 Kaikoura earthquake includes multiple
faults and complex fault geometries, finite-fault inversion
assuming one or a few simplified model fault planes may
produce erroneous inversion results due to modelingerrors caused by the inappropriate assumed fault geometries (Shimizu et al. 2020). This motivates using a potency
density tensor inversion (PDTI) (Shimizu et al. 2020),
that is not requiring assumption of fault geometries, toestimate the rupture process of the 2016 Kaikoura earthquake along an assumed model plane. The PDTI incorporates the uncertainty of the Green’s function in the data
covariance matrix (Yagi and Fukahata 2011) and introduces the Akaike’s Bayesian information criterion (ABIC)(e.g., Akaike 1980, Yabuki and Matsu’ura 1992, Satoet al. 2022), making it possible to perform stable inversion using seismic source model with a high degree of
freedom in the rupture direction (e.g., Hicks et al. 2020;
Yamashita et al. 2022a; Yagi et al. 2023).
In this study, we applied the PDTI to teleseismic
P-waveforms of the 2016 Kaikoura earthquake to simultaneously estimate the rupture propagation and the focal
mechanism variation. It revealed a source process consisting of an initial and a main rupture episode. The initial rupture propagates northeast from the hypocenterOhara et al. Progress in Earth and Planetary Science(2023) 10:35

and breaks shallow and deep parts of the source area;
deep rupture occurs where there is no surface rupture
reported. Then, the main rupture begins at the northeast
end of the Kekerengu fault and propagates bilaterally tothe northeast and southwest. The southwestward rupture(backward rupture) branches out and propagates into the
Jordan Thrust and the Papatea fault.

Page 3 of 12

into the data covariance matrix (Yagi and Fukahata 2011)
and then evaluates the relative weights of informationfrom observed data and prior constraints using Akaike’s
Bayesian information criterion (ABIC) (Akaike 1980,
Yabuki and Matsu’ura 1992, Sato et al. 2022). This inversion formulation reduces the effect of modeling errors
caused by uncertainties in fault geometry and Green’s
function and allows stable estimates of the seismic source
process even when the predefined model plane deviates
2 Methods, data, and modelingfrom the true fault plane (Shimizu et al. 2020). The PDTIThe PDTI of teleseismic P-waveforms has been devel- has been effectively applied to earthquakes for which it
oped to mitigate the effect of modeling errors due to an is difficult to assume a reasonable fault model (Okuwaki
inaccurate model fault geometries (Shimizu et al. 2020). et al. 2020, 2021; Tadapansawut et al. 2021; Yamashita
Teleseismic P-waveforms are sensitive to perturbations et al. 2021, 2022a; Okuwaki and Fan 2022). The PDTI isin the focal mechanism but insensitive to errors in the thus suitable for analyzing the 2016 Kaikoura earthquake.
source location, which is confirmed by both the synThe data covariance matrix, which introduces the error
thetic tests and real applications (e.g., Shimizu et al. 2020; term of the data covariance matrix, is a function of theTadapansawut et al. 2022; Yamashita et al. 2022b). There- model parameter (Yagi and Fukahata 2011). This meansfore, it is important to incorporate the focal mechanism that the inverse problem to be solved in PDTI becomeschange during the rupture propagation when building nonlinear. Following Yagi and Fukahata (2011), we set inia seismic source model to robustly estimate the rupture tial model parameters and solve the nonlinear problem
process (Shimizu et al. 2020). Whereas the conventional iteratively by improving the model parameters (Shimizu
finite-fault inversion method estimates the slip distribu- et al. 2020).
tion along the fault plane as a potency density distribuFor the PDTI, we used the teleseismic P-waveforms
tion, the PDTI method adopted in this study describes (vertical component) from 48 stations at epicentral disthe fault slip along the assumed model plane as a super- tances of 30°–100° downloaded from the Data Manageposition of five basis double-couple (Kikuchi and Kan- ment Center of the Incorporated Research Institutions
amori 1991) and then estimates the rupture evolution for Seismology (IRIS-DMC) (Fig. 2a). We converted(including perturbations in the focal mechanism) as a the raw waveform data to velocity waveforms at a samspatiotemporal distribution of the potency-rate density pling interval of 0.8 s. We calculated Green’s functions
tensor. Thus, the seismic waveform uj observed at a sta- at a sampling interval of 0.1 s by the method of Kikuchition j is given by
and Kanamori (Kikuchi and Kanamori 1991). We used
CRUST2.0 (Bassin et al. 2000) as a 1-D structure model5
uj (t) =Gqj (t, ξ ) + δGqj (t, ξ ) ∗ Ḋq (t, ξ )dξ + ebj (t),around the source (see Additional file 1: Table S1) and
q=1 S
set the value of t*, which controls the inelastic attenuawhere Gqj is the Green’s function of the q th basis double- tion of P-waves, to 1 s. The ray parameters and geometricouple moment tensor, δG qj is the modeling error on Gqj cal spreading factors were calculated based on the ak135(Yagi and Fukahata 2011), Ḋq is the potency-rate density reference velocity model (Kennett et al. 1995). We manufunction for the q th component of the basis double-cou- ally picked and aligned the P-wave first motion to corple moment tensor, ebj is a background and instrumental rect the travel-time deviations due to 3-D earth structure
Gaussian noise, ξ represents a position on the assumed (e.g., Fan and Shearer 2015). The effect of uncertainty ofmodel plane ( S ), and ∗ denotes the convolution operator underground structure was mitigated by introducing thein the time domain. Following Yagi and Fukahata (2011), error term of the Green’s function into the data covarithe background noise level is assessed by referring to the ance matrix (Yagi and Fukahata 2011).
Because the high-frequency component of the teleseispre-signal data, and the variance of the Green’s function
mic
body waveforms is effectively suppressed owing toerror is set to be proportional to the amplitude of thethenaturallow-pass filtering caused by inelastic attenuaGreen’s function.
tion,
given
resamplinginterval of 0.8 s, the waveforms are
Because this inversion allows any type of a deviatoric
only
littleaffectedbyaliasing(see Additional file 1: Fig.
focal mechanism on the assumed model plane, inforS1).
Conversely,
applyingalow-pass
filter that includes
mation about the fault geometry can be extracted fromananti-aliasing
filterincreases
theoff-diagonalcompothe observed data (Shimizu et al. 2020). To stably invert
nent
of
the
data
covariancematrix(Yagi
andFukahatasuch a high degree-of-freedom seismic source model, the
PDTI introduces the error term of the Green’s function

Ohara et al. Progress in Earth and Planetary Science(2023) 10:35

Page 4 of 12

Fig. 2 Summary of inversion results. a Azimuthal equidistant projection of the station distribution used in the inversion. The star denotes
the epicenter (Lanza et al. 2019). Triangles denote station locations; the waveforms for the four stations indicated with red triangles are shown in b.
The circles represent epicentral distances of 30° and 100°. b Observed (upper black trace) and synthetic (lower red trace) waveforms at the stations
marked in red in a. Station codes and maximum amplitudes are shown at the top. c Potency density tensors on the assumed model plane. The
map view in the top panel shows the top row of tensors on the assumed model plane, represented by the black line, and gray lines indicate
surface ruptures (Langridge et al. 2016). The profile in the bottom panel shows the tensors on the assumed model plane. Note that the beach ballsin the map are shown as a lower-hemisphere projection in the map and as a cross-section view from the southeast side in the bottom panel. Beachballs in the bottom panel are colored based on a Frohlich diagram (Frohlich 2001), in which blue is reverse faulting (T), green is strike-slip faulting(SS), red is normal faulting (N), and gray is other. The star denotes the hypocenter (Lanza et al. 2019). d Moment-rate function2011), making it difficult to stably invert the data covariance matrix. Therefore, we do not make any low-passfiltering of real and synthetic data and are able to reproduce observed waveforms well (Fig. 2b, Additional file 1:
Fig. S2).
We adopted a hypocenter location at 172.95° E, 42.62°
S, and 15 km depth (Lanza et al. 2019). We established
a 200 km × 35 km vertical model plane striking NE–SW(230°) to represent surface ruptures (Langridge et al.
2016; Hamling et al. 2017; Stirling et al. 2017; Litchfield et al. 2018) and aftershock activity (Lanza et al.
2019) (Fig. 1b). We set a maximum rupture velocity of
2.6 km/s to allow for the northeastward migration of thehigh-frequency source at about 2.0 km/s indicated by
P-waveform back-projection (Xu et al. 2018). The slip on
the model plane was expanded by linear B-spline functions in space with an interval of 10 km and 5 km in
the strike and dip directions, respectively, and by linear
B-spline functions in time with an interval of 0.8 s witha maximum duration of 60 s for each source element,
which is long enough to detect possible re-rupture and/or back-rupture propagation (Holden et al. 2017; Hicks
et al. 2020). The total duration of the event was set to95 s. The ABIC can prevent overfitting, even using large
number of model parameters (Sato et al. 2022).
We applied a time-adaptive smoothing constraint that
adjusts the smoothing strength in inverse proportion
to the changing amplitude of the potency-rate function(Yamashita et al. 2022b). This constraint can mitigate
the problem of oversmoothing during the main rupture,
which obscures the results (Yamashita et al. 2022b).

3 Results
We estimated the distribution of potency-rate density tensors on the assumed vertical model plane andthen time-integrated them at each source element toyield the spatial distribution of potency density tensors shown in Fig. 2c. This figure shows an area of high

Ohara et al. Progress in Earth and Planetary Science(2023) 10:35

potency density centered around 110 km northeast of
the epicenter on the shallow part of the model plane.
The dominant focal mechanism (with relatively large
potency density) differs along the length of the fault
plane, being oblique reverse slip for 50–120 km andstrike-slip for 130–150 km northeast of the epicenter(Fig. 2c). The total seismic moment is 1.1 × ­1021 Nm(Mw 8.0). The moment-rate function of the wholesource process, obtained by calculating the seismic
moment rate of the best-fitted double-couple source atPage 5 of 12

each sampling time, is around 1.0 × ­1019 Nm/s until 55 sfrom the origin time and then increases rapidly, reaching 6.0 × ­1019 Nm/s at 66 s (Fig. 2d).
Figure 3 shows selected snapshots of the potency-ratedensity tensors on a cross section of the model fault
plane; see Additional file 1: Fig. S3 for the full set of snapshots. Figure 4 shows a map view of the strike angles of
the nodal planes of these tensors along the top of themodel plane. During the first 10 s from the origin time, astrike-slip rupture striking about 25° clockwise from theFig. 3 Selected snapshots of potency-rate density tensors a before 55 s and b after 60 s. Beach balls are shown in cross-section view fromthe southeast side of the assumed model plane. The background color is scaled with the maximum potency-rate density during 0–55 s for (a)and 60–74 s for (b); note that the scales differ for the two plots. The star denotes the hypocenter (Lanza et al. 2019). Black bars are the locationsof the surface faults (Langridge et al. 2016) projected onto the model planeOhara et al. Progress in Earth and Planetary Science(2023) 10:35

Page 6 of 12

increased with time after 64 s. The rupture gradually
weakened after 70 s and ceased at 95 s. The inverted solution well explains the teleseismic P-waveforms (Fig. 2b,
Additional file 1: Fig. S2).

Fig. 4 Map views showing selected snapshots of strikes ofthe potency-rate density tensors (cross-marks) in the top rowof the model plane. Right-lateral or northwest-dipping nodal planesof tensors with relatively large potency-rate density are emphasized.
Note that the color scale changes after 55 smodel plane propagated to the northeast of the epicenter(Figs. 3a and 4). Hereafter, we refer to the origin time as
0 s. The rupture then propagated further northeastward
on the shallow part of the model plane, changing to an
oblique reverse focal mechanism. This shallow rupture
stagnated at about 40 km northeast of the epicenter after
20 s; however, a deeper rupture continued on the model
plane between 20 and 30 s, reaching 70 km northeast of
the epicenter. An isolated reverse rupture occurred at
25–30 s near the ground surface around the epicenter.
During 30–45 s, an oblique reverse rupture appearednear the ground surface about 70 km northeast of theepicenter and propagated northeast; between 45 and
50 s, the rupture propagation pattern was obscure (seeAdditional file 1: Fig. S3a).
After 50 s, the main rupture emerged near the groundsurface about 110 km northeast of the epicenter andpropagated bilaterally to the northeast and southwest(Fig. 3 and Additional file 1: Fig. S4). Between 50 and
55 s, the dominant focal mechanisms were mixed reverseand strike-slip with the right-lateral nodal plane oriented
about 40° clockwise from the model plane (Figs. 3a and
4). The northeastward rupture, a strike-slip rupture striking about 10° counterclockwise from the model plane,
propagated through the shallow part of the model planeand reached the edge of the model plane at about 68 s
(Figs. 3b and 4). The southwestward rupture reachedabout 70 km northeast of the epicenter by 70 s (Fig. 3b).
Between 60 and 64 s, it was dominantly strike-slip near
the ground surface and reverse in the deep part of themodel plane (Fig. 3b). The reverse slip component

4 Reproducibility and sensitivity tests
We performed a numerical experiment to test the stability and reproducibility of our potency-rate density tensor
distribution. We used the obtained source model as an
input model and generated input synthetic waveforms forthe 48 stations used in the analysis. As an error of Green’s
function, we assumed random Gaussian noise with zeromean and a standard deviation of 3% of maximum amplitude of each calculated Green’s function. As the background noise, we assumed random Gaussian noise withzero mean and a standard deviation of 1 μm. The input
synthetic waveform was produced by adding backgroundnoise to the convolved Green’s function with error andthe input slip-rate function (see Additional file 1: Fig. S5).
The input synthetic waveforms were inverted with thesame settings used with the real waveforms.
We performed a structure sensitivity test using the 1-D
structure model CRUST1.0 (Laske et al. 2013) for thesource region instead of CRUST2.0 (Bassin et al. 2000)(see Additional file 1: Table S2). We estimated the rupture evolution using the same observed dataset and thesame inversion settings as for our preferred modeling.
We also performed another sensitivity test projecting
rupture process onto a horizontal model plane. We established a 200 km × 70 km horizontal model plane striking
NE–SW (230°) to represent surface ruptures (Langridge
et al. 2016; Hamling et al. 2017; Stirling et al. 2017; Litchfield et al. 2018) and aftershock activity (Lanza et al.
2019) (see Additional file 1: Figs. S8 and S9). The slip on
the model plane was expanded by bilinear B-spline functions in space with an interval of 10 km. The hypocentral
depth was 10 km, where rupture mainly detected in theanalysis using vertical model plane (Figs. 2 and 3). We
used the same observed dataset and the same inversion
settings as for our preferred modeling using the verticalplane.
Both the reproducibility and structure sensitivity tests
successfully reproduced the robust features in our preferred model: These included the initial strike-slip rupture during the first 10 s, the northeast-propagating
oblique reverse rupture at varying depths between 10
and 30 s and re-appearing near the ground surface about
70 km northeast of the epicenter, and the main bilateral
rupture starting about 110 km northeast of the epicenteraround 50 s with a strike-slip rupture propagating northeast and an oblique-slip rupture propagating southwest (see Additional file 1: Figs. S6, S7, and S8). In theOhara et al. Progress in Earth and Planetary Science(2023) 10:35

reproducibility test, however, potency-rate density tensors between 25 and 35 s at 50–100 km northeast of theepicenter show reverse faulting in the input model but
strike-slip in the output model, and potency-rate density
distribution between 40 and 50 s is different between the
input and output models (Additional file 1: Fig. S6). As
these features are not robustly reproduced in this exercise, we do not discuss them in the following section.
The sensitivity test using the horizontal model plane
reproduced the lateral variation of rupture evolution
in the preferred source model using the vertical model
plane: the initial strike-slip rupture during the first
10 s, the northeast-propagating oblique reverse rupture between 10 and 30 s, and the main bilateral rupture
from about 110 km northeast of the epicenter (see Additional file 1: Figs. S9 and S10). However, the sensitivity
test shows different potency-rate density tensors fromthose in the preferred source model: Potency-rate density tensors between 25 and 30 s at about 70 km northeast of the epicenter had P-axes rotated counterclockwise
by about 90° from those in the preferred source model,
and potency-rate density tensors between 60 and 66 s
at 70–120 km northeast of the epicenter show reverse
faulting in the sensitivity test, whereas those in the preferred source model show strike-slip and oblique slip. In
general, the teleseismic waveforms contain depth phases,
but the horizontal model plane cannot reproduce thesevariations, which may have contributed to the discrepancy between the vertical and horizontal model planes(Shimizu et al. 2020; Yamashita et al. 2022b).

5 Discussion
Our result shows that the rupture process of the 2016
Kaikoura earthquake can be divided into initial andmain rupture episodes: the initial rupture propagatednortheastward; the main rupture propagated bilaterally from 110 km northeast of the epicenter, involving
backward rupture propagation toward the epicenter.
The total moment tensor, which is obtained by integrating the potency-rate density tensors in space and time,shows oblique reverse faulting, which is consistent with
the GCMT solution (Fig. 1). Our estimated total seismic
moment is 1.1 × ­1021 Nm (Mw 8.0), which is larger thanthe GCMT solution (6.7 × ­1020 Nm; Mw 7.8). The discrepancy between our estimates of the seismic moment
and those from the other inversion methods may be due
to a simplified seismic source model that may not be adequate to represent the complex source process. In the following, we will discuss how those rupture episodes relate
to the observed surface ruptures, to unravel the unprecedentedly complex rupture process of the 2016 Kaikoura
earthquake.

Page 7 of 12

As the initial strike-slip rupture propagated northeast during the first 10 s (Fig. 3a), the right-lateral nodal
planes of the potency-rate density tensors matched thestrike of the Humps fault (Langridge et al. 2016) (Fig. 4).
Hereafter, we interpret right-lateral or northwest-dipping
nodal planes as fault planes, because right-lateral strikeslip or northwest-directing dip slip is dominant in surface
faults reported in field surveys (Clark et al. 2017; Hamling et al. 2017; Stirling et al. 2017; Litchfield et al. 2018).
An oblique reverse rupture then propagated northeast
through the shallow part of the model plane. After 20 s,
the shallow rupture stagnated about 40 km northeast of
the epicenter, while the oblique reverse rupture continued to propagate deeper on the model plane (Fig. 3a).
The location where the shallow rupture stagnated corresponds to the gap in surface ruptures between the
Conway-Charwell and Manakau faults (Langridge et al.
2016) (Figs. 1 and 3a), and the deep oblique reverse slip
has also been identified by the finite-fault inversion ofgeodetic data (Hamling et al. 2017). Between 30 and 35 s,
oblique reverse rupture appeared near the ground surface about 70 km northeast of the epicenter, corresponding to the southwest end of the Manakau fault (Langridge
et al. 2016), and then propagated near the ground surface
until 45 s (Fig. 3a). Our results show that the initial rupture shifted deeper around the area of no surface rupture
between 20 and 30 s. However, because slips on multiple
fault planes are projected onto the single model plane in
our inversion, it is difficult to determine whether these
ruptures were connected at depth.
It is controversial how the plate interface contributed
to moment release in the 2016 Kaikoura earthquake (e.g.,
Lanza et al. 2019). Although the deep rupture between 20
and 30 s appeared at about 25 km depth, the resolved dip
angles (~ 40°) are steeper than those of the hypothesizedplate interface (e.g., Williams et al. 2013). Our solution isstill limited to be conclusive about a cause of the rupture
at depth to reach Kekerengu fault. This may be realized
by a listric geometry of faults connecting at depth (e.g.,
Xu et al. 2018), or the Kekerengu fault could be connected at depth via the Point Keen fault (e.g., Ulrich et al.
2019). We note here that our bilateral rupture scenario
would be useful to further constrain the rupture dynamics of the relevant deep faults by the independent studies
such as physics-based rupture simulation.
After 50 s, the main rupture appeared in the northeast part of the Kekerengu fault (Langridge et al. 2016)
and then propagated bilaterally until about 70 s, suchthat one end of the rupture appeared to propagate backward toward the epicenter (Fig. 3). Because we cannot
trace the rupture migration between 45 and 50 s, it isdifficult to determine how the initial rupture migrated
to the main rupture. The potency-rate density tensors

Ohara et al. Progress in Earth and Planetary Science(2023) 10:35

obtained at 50–55 s indicate both strike-slip and reverse
faulting, and the strikes of their right-lateral nodal planes
are consistent with that of the northeastern Kekerengu
fault (Langridge et al. 2016) (Fig. 4). For the northeastward strike-slip rupture, the strikes of the right-lateral
nodal planes match the orientation of the Needles fault(Langridge et al. 2016) (Fig. 4), and the dominance of
strike-slip faulting in the shallow part of the model plane(Fig. 3b) is consistent with other studies (Bai et al. 2017;
Cesca et al. 2017; Hollingsworth et al. 2017; Wang et al.
2018a, 2018b; Zheng et al. 2018; Xu et al. 2018; Mouslopoulou et al. 2019). For the backward rupture, thepotency-rate density tensors near the ground surface
show a transition from oblique strike-slip to obliquereverse faulting 80–110 km northeast of the epicenter(Fig. 3b), and the strikes of the right-lateral or northwestdipping nodal planes match those of the central Kekerengu fault and the Jordan Thrust (Langridge et al. 2016)(Fig. 4).
The potency-rate density tensors around the Jordan
Thrust and Papatea fault contain large non-double-couple components, reaching an 80% maximum from 60 to
66 s, that then rapidly decrease to less than 20% after 66 s
(Fig. 5). Our reproducibility tests also captured the time
variation of this component (see Additional file 1: Figs. S6and S7). The size of the non-double-couple component
from 60 to 66 s suggests that large non-double-couple
component is just apparent, and real rupture might havebeen 100% double-couple but slips occurred on multiple
faults with different orientations (e.g., Liu and Zahradník
2020; Turhan et al. 2023); this is consistent with reverse
faulting with sinistral strike-slip reported on the Papatea
fault (Hamling et al. 2017; Stirling et al. 2017; Litchfield
et al. 2018; Wang et al. 2018b; Xu et al. 2018), which isnearly perpendicular to the other surface ruptures (seeAdditional file 1: Fig. S11). Our result suggests that the

Page 8 of 12

backward rupture on the Kekerengu fault not only propagated into the Jordan Thrust, but also branched out andpropagated into the Papatea fault. Given that the rightlateral strike-slip rupture propagates in a southwest
direction along the Kekerengu fault, the Papatea fault issituated in the compressional quadrant; this suggests thatthe southwestward rupture along the Kekerengu fault
can better explain a trigger of the reverse faulting rupture along the Papatea fault than the northeastward rupture along the Jordan Thrust, which should require thePapatea fault to be located in the extensional quadrant.
Although we find it reasonable to explain the Papatea
rupture by our series of bilateral ruptures, more detailed
analyses and simulations incorporating the detailed
geometries of those faults will be required to testify
which of the scenarios is more favorable for the Papatea
rupture. Near the southwest end of the backward rupture of the Jordan Thrust and Papatea fault, the strikesof the right-lateral or northwest-dipping nodal planes
were about 10° clockwise from the model plane, whichis consistent with the strikes of the Upper Kowhai andManakau faults (Langridge et al. 2016) (Fig. 4).
The Papatea Fault has been a key element and its role inrupture propagation and deformation have been actively
discussed (e.g., Kilinger et al. 2018; Ulrich et al. 2019). In
particular, the seismological findings seem to have apparently diverging view of the rupture timing of the Papatea
faulting (e.g., Xu et al. 2018; Wang et al. 2018a, b). Our
study here is proposing that, based on our estimates ofrupture timing and the fault geometry, the rupture proceeds from the Kekerengu to Papatea faults as a part ofthe bilateral rupture propagation.
So far, an earthquake source modeling has often been
relying on a restricted degree of freedom, which has
been considered as a requirement for a plausible solution. However, the modeling employing fewer degrees

Fig. 5 Map views showing selected snapshots of potency-rate density tensors (lower-hemisphere projections) between 60 and 70 s in the top rowof the model plane 80–100 km northeast of the epicenter. The color of the beach ball symbols represents the potency-rate density. Above eachsymbol is shown the ratio of the apparent non-double-couple component that arises from mixing faulting types

Ohara et al. Progress in Earth and Planetary Science(2023) 10:35

of freedom might be easy to drop information that arerecorded in the observed data and critical to interpretthe source process (e.g., Shimizu et al. 2020), albeit thesolution derived from those modeling apparently looks
not bad. One of the advantages of employing a modelwith a high degree of freedom (e.g., this study) is thata solution is less susceptible to the modelers’ preconceptions. By estimating the potency tensor density distribution including the non-double-couple component,
we found that the backward rupture branched out andpropagated on the Papatea fault, which, to our best
knowledge, has not been reported in previous attemptsof the teleseismic body waves analyses.
Our analysis suggests the following scenario for themain rupture: It propagated bilaterally from the northeast part of the Kekerengu fault, the northeastward
rupture propagating along the Needles fault and thesouthwestward rupture propagating along the Kekerengu fault, Jordan Thrust, Papatea, Upper Kowhai,
and Manakau faults. We interpret the simultaneous
rupture events in the area around the Needles fault and
the Jordan Thrust noted in previous studies (Bai et al.
2017; Cesca et al. 2017; Hollingsworth et al. 2017) asbilateral rupture propagation. In addition, the backprojection image (Xu et al. 2018) shows that the seismic
wave radiation point moves toward the epicenter fromaround the south edge of the Papatea fault between 50
and 70 s, a finding consistent with backward rupture
propagation.
In the region of the backward rupture, multiple faults
may have ruptured during the initial rupture phase,because the aftershock region extends perpendicular to
the model plane and the focal mechanisms varied during
the initial rupture (Fig. 6). Because our model fault plane
may include projections of multiple independent ruptures, we cannot determine which faults participated in
the initial rupture. Therefore, we cannot say whether the
backward rupture was a re-rupture (Holden et al. 2017)or a rupture on a different fault, as in the 2010 El Mayor–
Cucapah earthquake (Yamashita et al. 2022a).
Back-propagating ruptures in seismic events are notso rare; they have been reported in the 2010 El Mayor–
Cucapah earthquake (Yamashita et al. 2022a), the 2011
Tohoku-Oki earthquake (Ide et al. 2011), the 2014
Iquique earthquake (Yagi et al. 2014), the 2016 Romanche transform-fault earthquake (Hicks et al. 2020), the2018 Peru earthquake (Hu et al. 2021), and the 2020
Elazig, Turkey earthquake (Gallovič et al. 2020). With
the exception of the 2011 Tohoku-Oki earthquake, where
the backward rupture followed an overshooting rupturenear the free surface (Ide et al. 2011), these earthquakes
have in common an initial weak rupture which triggers
a main rupture, at a point distant from the hypocenter,

Page 9 of 12

Fig. 6 Map views showing selected snapshots of potency-ratedensity tensors (lower-hemisphere projections) between 32 and 42 s
in the top row of the entire model plane. Black dots represent
aftershocks during the week after the mainshock (Lanza et al. 2019).
Gray, orange, blue, and green lines indicate surface faults (Langridge
et al. 2016)that involves a back-propagating rupture. It appears thatthe 2016 Kaikoura earthquake is another example of thiskind of event.
Our modeling approach requires few assumptions ofmodeling; that is, we solve for multiplicity of fault configuration and diverse rupture geometries on the flat single
model fault. This is still prone to non-uniqueness in the
Kaikoura rupture, primarily due to the limited spatial resolution of teleseismic records, but the rupture directions
and timing, involving back-rupture propagation resolvedin our model, in turn, can be useful for further inverse
and/or forward modeling using near-field datasets, which
contribute to converge to a realistic source model of the
Kaikoura earthquake.

6 Conclusions
We revealed the source process of the 2016 Kaikoura
earthquake by a potency density tensor inversion fromteleseismic P-waveform data, a method for which we did
not need to strictly define the fault geometry and rupture
directions. We found complex episodes including an initial unilateral and a delayed main bilateral rupture, andthe variations of the focal mechanisms are consistent
with the reported surface ruptures. The initial rupture
propagated northeastward at deep depths, when it passed
through a gap in reported surface ruptures. The main
rupture involved the southwestward back-rupture propagation, and it branched out and propagated into the Jordan Thrust and Papatea fault from the Kekerengu fault.
Our result suggests that teleseismic waveform data canresolve such a complex rupture process, and the potency
density tensor inversion approach of projecting slips on

Ohara et al. Progress in Earth and Planetary Science(2023) 10:35

multiple faults onto a single model plane, as opposed toan approach of prescribing fault planes, is useful for analyzing earthquakes with complex fault geometries.
Abbreviations
GCMT	Global Centroid Moment Tensor
PDTI	Potency density tensor inversion

Supplementary Information
The online version contains supplementary material available at https://​doi.​
org/​10.​1186/​s40645-​023-​00565-z.
Additional file 1: Tables S1 and S2. Velocity structure models used for
calculating the Green’s function. Figure S1. Comparison of waveforms
sampled at different intervals. Figure S2. Waveform fitting between
observed and synthetic waveforms. Figure S3. Full snapshots of potencyrate density evolution for the obtained source model of the 2016 Kaikoura
earthquake. Figure S4. Potency-rate density evolution along the strikedirection of the model plane. Figure S5. Waveform fitting between input
and output waveforms for the reproducibility test using the result of theanalysis for real waveforms as input. Figure S6. Snapshots of potency-ratedensity evolution before 55 s for reproducibility test using waveforms generated from the obtained source model as input. Figure S7. Snapshotsof potency-rate density evolution after 56 s for reproducibility test using
waveforms generated from the obtained source model as input. Figure
S8. Snapshots of potency-rate density evolution for sensitivity test using
1-D velocity structure near the epicenter determined based on CRUST1.0.
Figure S9. Comparison of snapshots of potency-rate density evolution
before 55 s for the preferred source model using the vertical model planesensitivity test under the assumption of the horizontal model plane.
Figure S10. Comparison of snapshots of potency-rate density evolution
after 56 s for the preferred source model using the vertical model planesensitivity test under the assumption of the horizontal model plane. Figure S11. Comparison of the obtained potency-rate density tensors during
60–66 s and visual summary of moment tensor transition between two
endmembers: PF and HP-JTKF-NF.
Acknowledgements
We would like to thank Editor František Gallovič and Reviewers Thomas Ulrich
and Jiří Zahradník for their thorough evaluation and constructive suggestions,which have led to significant improvements of our work. We thank Federica
Lanza for providing their relocated aftershock data (Lanza et al. 2019). The facilities of IRIS Data Services, and specifically the IRIS Data Management Center,
provided access to waveforms, related metadata, and derived products usedin this study. IRIS Data Services are funded through the Seismological Facilities
for the Advancement of Geoscience (SAGE) Award of the National Science
Foundation under Cooperative Support Agreement EAR-1851048. We used
FPSPACK (Gasperini and Vannucii 2003) software to handle the focal mechanism obtained in the inversion. The code for the PDTI was developed by Yuji
Yagi, Kosuke Shimizu, and Shinji Yamashita. All of the figures were generated
with Matplotlib (v3.4.0:https://​doi.​org/​10.​1109/​MCSE.​2007.​55) (Hunter 2007)
and Pyrocko (v2021.9.14:https://​doi.​org/​10.​5880/​GFZ.2.​1.​2017.​001) (Heimann
et al. 2017).
Author contributions
KO and YY designed this study, compiled the data, and performed theanalyses. All authors interpreted the research results. KO prepared figures andwrote the manuscript, which was revised by YY, SY, RO, SH, and YF. All authors
approved the manuscript. All authors agreed both to be personally accountable for their own contributions and to ensure that questions related to theaccuracy or integrity of any part of the work were appropriately investigatedand resolved and their resolution documented in the literature.

Page 10 of 12

Funding
This work was supported by Grant-in-Aid for Scientific Research (C) 22K03751
from the Ministry of Education, Culture, Sports, Science and Technology(MEXT).
Availability of data and materials
All seismic data were downloaded through the IRIS Wilber 3 system (https://​ds.​iris.​edu/​wilbe​r3/) or IRIS Web Services (https://​servi​ce.​iris.​edu/), including
the following seismic networks: (1) BDSN (https://​doi.​org/​10.​7932/​BDSN),(2) SCSN (https://​doi.​org/​10.​7914/​SN/​CI), (3) GEOSCOPE (https://​doi.​org/​10.​
18715/​GEOSC​OPE.G), (4) GEOFON (https://​doi.​org/​10.​14470/​TR560​404), (5)the Global Telemetered Seismograph Network (https://​doi.​org/​10.​7914/​SN/​GT), (6) the Hong Kong Seismograph Network, (7) the New China Digital
Seismograph Network (https://​doi.​org/​10.​7914/​SN/​IC), (8) the IRIS/IDA Seismic
Network (https://​doi.​org/​10.​7914/​SN/​II), and (9) the Global Seismograph Network (https://​doi.​org/​10.​7914/​SN/​IU). The CRUST1.0 and CRUST2.0 structural
velocity models are available from https://​igppw​eb.​ucsd.​edu/​~gabi/​crust1.​
html and https://​igppw​eb.​ucsd.​edu/​~gabi/​crust2.​html, respectively.

Declarations
Competing interests
The authors declare no competing interests. Correspondence and requests formaterials should be addressed to K.O. or Y.Y.
Received: 15 February 2023 Accepted: 23 June 2023

References
Akaike H (1980) Likelihood and the Bayes procedure. Trab Estad Investig Oper
31:143–166. https://​doi.​org/​10.​1007/​BF028​88350
Bai Y, Lay T, Cheung KF, Ye L (2017) Two regions of seafloor deformation generated the tsunami for the 13 November 2016, Kaikoura, New Zealandearthquake. Geophys Res Lett 44:6597–6606. https://​doi.​org/​10.​1002/​
2017G​L0737​17
Bassin C, Laske G, Masters G (2000) The current limits of resolution for surface
wave tomography in North America. EOS Trans AGU 81:F897
Bird P (2003) An updated digital model of plate boundaries. Geochem Geophys Geosyst. https://​doi.​org/​10.​1029/​2001G​C0002​52
Cesca S, Zhang Y, Mouslopoulou V, Wang R, Saul J, Savage M, Heimann S,
Kufner SK, Oncken O, Dahm T (2017) Complex rupture process of theMw 7.8, 2016, Kaikoura earthquake, New Zealand, and its aftershock
sequence. Earth Planet Sci Lett 478:110–120. https://​doi.​org/​10.​1016/j.​
epsl.​2017.​08.​024
Clark KJ, Nissen EK, Howarth JD, Hamling IJ, Mountjoy JJ, Ries WF, Jones K,
Goldstien S, Cochran UA, Villamor P, Hreinsdóttir S, Litchfield NJ, Mueller
C, Berryman KR, Strong DT (2017) Highly variable coastal deformation in
the 2016 MW 7.8 Kaikōura earthquake reflects rupture complexity along
a transpressional plate boundary. Earth Planet Sci Lett 474:334–344.
https://​doi.​org/​10.​1016/j.​epsl.​2017.​06.​048
DeMets C, Gordon RG, Argus DF, Stein S (1994) Effect of recent revisions to thegeomagnetic reversal time scale on estimates of current plate motions.
Geophys Res Lett 21:2191–2194. https://​doi.​org/​10.​1029/​94GL0​2118
Duputel Z, Rivera L (2017) Long-period analysis of the 2016 Kaikoura earthquake. Phys Earth Planet Inter 265:62–66. https://​doi.​org/​10.​1016/j.​pepi.​
2017.​02.​004
Dziewonski AM, Chou TA, Woodhouse JH (1981) Determination of earthquake
source parameters from waveform data for studies of global and regional
seismicity. J Geophys Res 86:2825–2852. https://​doi.​org/​10.​1029/​JB086​
iB04p​02825
Ekström G, Nettles M, Dziewoński AM (2012) The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes. Phys Earth Planet
Inter 200–201:1–9. https://​doi.​org/​10.​1016/j.​pepi.​2012.​04.​002

Ohara et al. Progress in Earth and Planetary Science(2023) 10:35

Fan W, Shearer PM (2015) Detailed rupture imaging of the 25 April 2015 Nepal
earthquake using teleseismic P waves. Geophys Res Lett 42:5744–5752.
https://​doi.​org/​10.​1002/​2015G​L0645​87
Frohlich C (2001) Display and quantitative assessment of distributions of earthquake focal mechanisms. Geophys J Int 144:300–308. https://​doi.​org/​10.​1046/j.​1365-​246x.​2001.​00341.x
Gallovič F, Zahradník J, Plicka V, Sokos E, Evangelidis C, Fountoulakis I, Turhan F(2020) Complex rupture dynamics on an immature fault during the 2020
Mw 6.8 Elazığ earthquake, Turkey. Commun Earth Environ. https://​doi.​
org/​10.​1038/​s43247-​020-​00038-x
Gasperini P, Vannucii G (2003) FPSPACK: a package of FORTRAN subroutines tomanage earthquake focal mechanism data. Comput Geosci 29:893–901.
https://​doi.​org/​10.​1016/​S0098-​3004(03)​00096-7
Hamling IJ, Hreinsdóttir S, Clark K, Elliott J, Liang C, Fielding E, Litchfield N,
Villamor P, Wallace L, Wright TJ, D’Anastasio E, Bannister S, Burbidge D,
Denys P, Gentle P, Howarth J, Mueller C, Palmer N, Pearson C, Power W,
Barnes P, Barrell DJA, Van Dissen R, Langridge R, Little T, Nicol A, Pettinga J,
Rowland J, Stirling M (2017) Complex multifault rupture during the 2016
Mw 7.8 Kaikōura earthquake, New Zealand. Science. https://​doi.​org/​10.​1126/​scien​ce.​aam71​94
Heimann S, Kriegerowski M, Isken M, Cesca S, Daout S, Grigoli F, Juretzek C,
Megies T, Nooshiri N, Steinberg A, Sudhaus H, Vasyura-Bathke H, Willey T,
Dahm T (2017) Pyrocko—an open-source seismology toolbox and library.
GFZ Data Serv. https://​doi.​org/​10.​5880/​GFZ.2.​1.​2017.​001
Hicks SP, Okuwaki R, Steinberg A, Rychert CA, Harmon N, Abercrombie RE,
Bogiatzis P, Schlaphorst D, Zahradnik J, Kendall JM, Yagi Y, Shimizu K, Sudhaus H (2020) Back-propagating supershear rupture in the 2016 M w 7.1
Romanche transform fault earthquake. Nat Geosci 13:647–653. https://​doi.​org/​10.​1038/​s41561-​020-​0619-9
Holden C, Kaneko Y, D’Anastasio E, Benites R, Fry B, Hamling IJ (2017) The 2016
Kaikōura earthquake revealed by kinematic source inversion and seismic
wavefield simulations: slow rupture propagation on a geometrically complex crustal fault network. Geophys Res Lett 44:11320–11328. https://​doi.​
org/​10.​1002/​2017G​L0753​01
Hollingsworth J, Ye L, Avouac JP (2017) Dynamically triggered slip on a splay
fault in the Mw 7.8, 2016 Kaikoura (New Zealand) earthquake. Geophys
Res Lett 44:3517–3525. https://​doi.​org/​10.​1002/​2016G​L0722​28
Hu Y, Yagi Y, Okuwaki R, Shimizu K (2021) Back-propagating rupture evolution within a curved slab during the 2019 Mw 8.0 Peru intraslab
earthquake. Geophys J Int 227:1602–1611. https://​doi.​org/​10.​1093/​gji/​ggab3​03
Hunter JD (2007) Matplotlib: a 2D graphics environment. Comput Sci Eng
9:90–95. https://​doi.​org/​10.​1109/​MCSE.​2007.​55
Ide S, Baltay A, Beroza GC (2011) Shallow dynamic overshoot and energetic
deep rupture in the 2011 Mw 9.0 Tohoku-Oki earthquake. Science
332:1426–1429. https://​doi.​org/​10.​1126/​SCIEN​CE.​12070​20/​SUPPL_​FILE/​12070​20-​IDE-​SOM.​PDF
Kennett BLN, Engdahl ER, Buland R (1995) Constraints on seismic velocities
in the Earth from traveltimes. Geophys J Int 122:108–124. https://​doi.​
org/​10.​1111/j.​1365-​246X.​1995.​tb035​40.x
Kikuchi M, Kanamori H (1991) Inversion of complex body waves—III. Bull
Seismol Soc Am 81:2335–2350. https://​doi.​org/​10.​1785/​BSSA0​81006​
2335
Klinger Y, Okubo K, Vallage A, Champenois J, Delorme A, Rougier E, Lei Z,
Knight EE, Munjiza A, Satriano C, Baize S, Langridge R, Bhat HS (2018)
Earthquake damage patterns resolve complex rupture processes. Geophys Res Lett 45:10279–10287. https://​doi.​org/​10.​1029/​2018G​L0788​42
Langridge R, Ries W, Litchfield N, Villamor P, Van Dissen R, Barrell D, Rattenbury M, Heron D, Haubrock S, Townsend D, Lee J, Berryman K, Nicol
A, Cox S, Stirling M (2016) The New Zealand active faults database.
NZ J Geol Geophys 59:86–96. https://​doi.​org/​10.​1080/​00288​306.​2015.​
11128​18
Lanza F, Chamberlain CJ, Jacobs K, Warren-Smith E, Godfrey HJ, Kortink M,
Thurber CH, Savage MK, Townend J, Roecker S, Eberhart-Phillips D(2019) Crustal fault connectivity of the ­Mw 7.8 2016 Kaikōura earthquake
constrained by aftershock relocations. Geophys Res Lett 46:6487–6496.
https://​doi.​org/​10.​1029/​2019G​L0827​80
Laske G, Masters G, Ma Z, Pasyanos M (2013) Update on CRUST10—A 1-degree
global model of Earth’s crust. Geophys Res Abstr 15:2658
Litchfield NJ, Villamor P, van Dissen RJ, Nicol A, Barnes PM, Barrell DJA, Pettinga
JR, Langridge RM, Little TA, Mountjoy JJ, Ries WF, Rowland J, Fenton C,

Page 11 of 12

Stirling MW, Kearse J, Berryman KR, Cochran UA, Clark KJ, Hemphill-Haley
M, Khajavi N, Jones KE, Archibald G, Upton P, Asher C, Benson A, Cox
SC, Gasston C, Hale D, Hall B, Hatem AE, Heron DW, Howarth J, Kane TJ,
Lamarche G, Lawson S, Lukovic B, McColl ST, Madugo C, Manousakis J,
Noble D, Pedley K, Sauer K, Stahl T, Strong DT, Townsend DB, Toy V, Williams J, Woelz S, Zinke R (2018) Surface rupture of multiple crustal faults
in the 2016 Mw 7.8 Kaikōura, New Zealand, earthquake. Bull Seismol Soc
Am 108:1496–1520. https://​doi.​org/​10.​1785/​01201​70300
Liu J, Zahradník J (2020) The 2019 MW 5.7 Changning Earthquake, Sichuan
Basin, China: a shallow doublet with different faulting styles. Geophys Res
Lett. https://​doi.​org/​10.​1029/​2019G​L0854​08
Mitchell JS, Mackay KA, Neil HL, Mackay EJ, Pallentin A, Notman P (2012)
Undersea New Zealand, 1:5,000,000. NIWA Chart Misc Ser 92
Mouslopoulou V, Saltogianni V, Nicol A, Oncken O, Begg J, Babeyko A, Cesca
S, Moreno M (2019) Breaking a subduction-termination from top to bottom: the large 2016 Kaikōura Earthquake, New Zealand. Earth Planet Sci
Lett 506:221–230. https://​doi.​org/​10.​1016/j.​epsl.​2018.​10.​020
Okuwaki R, Fan W (2022) Oblique convergence causes both thrust and strikeslip ruptures during the 2021 M 7.2 Haiti earthquake. Geophys Res Lett.
https://​doi.​org/​10.​1029/​2021G​L0963​73
Okuwaki R, Hirano S, Yagi Y, Shimizu K (2020) Inchworm-like source evolution
through a geometrically complex fault fueled persistent supershear
rupture during the 2018 Palu Indonesia earthquake. Earth Planet Sci Lett.
https://​doi.​org/​10.​1016/j.​epsl.​2020.​116449
Okuwaki R, Hicks SP, Craig TJ, Fan W, Goes S, Wright TJ, Yagi Y (2021) Illuminating a contorted slab with a complex Intraslab rupture evolution during
the 2021 Mw 7.3 East Cape, New Zealand. Earthquake. https://​doi.​org/​10.​1029/​2021G​L0951​17
Sato DSK, Fukahata Y, Nozue Y (2022) Appropriate reduction of the posterior
distribution in fully Bayesian inversions. Geophys J Int. https://​doi.​org/​10.​1093/​gji/​ggac2​31
Shimizu K, Yagi Y, Okuwaki R, Fukahata Y (2020) Development of an inversion
method to extract information on fault geometry from teleseismic data.
Geophys J Int 220:1055–1065. https://​doi.​org/​10.​1093/​gji/​ggz496
Stirling MW, Litchfield NJ, Villamor P, Van Dissen RJ, Nicol A, Pettinga J, Barnes
P, Langridge RM, Little T, Barrell DJA, Mountjoy J, Ries WF, Rowland J,
Fenton C, Hamling I, Asher C, Barrier A, Benson A, Bischoff A, Borella J,
Carne R, Cochran UA, Cockroft M, Cox SC, Duke G, Fenton F, Gasston C,
Grimshaw C, Hale D, Hall B, Hao KX, Hatem A, Hemphill-Haley M, Heron
DW, Howarth J, Juniper Z, Kane T, Kearse J, Khajavi N, Lamarche G, Lawson
S, Lukovic B, Madugo C, Manousakis I, Mccoll S, Noble D, Pedley K, Sauer
K, Stah T, Strong DT, Townsend DB, Toy V, Villeneuve M, Wandres A, Williams J, Woelz S, Zinke R (2017) THE M W 7.8 2016 Kaikōura earthquake:surface fault rupture and seismic hazard context. Bull N Z Soc Earthq Eng
50:73–84
Tadapansawut T, Okuwaki R, Yagi Y, Yamashita S (2021) Rupture process of the2020 Caribbean Earthquake along the oriente transform fault, involving
Supershear rupture and geometric complexity of fault. Geophys Res Lett.
https://​doi.​org/​10.​1029/​2020G​L0908​99
Tadapansawut T, Yagi Y, Okuwaki R, Yamashita S, Shimizu K (2022) Complex
rupture process on the conjugate fault system of the 2014 Mw 6.2
Thailand earthquake. Prog Earth Planet Sci. https://​doi.​org/​10.​1186/​
s40645-​022-​00484-5
Turhan F, Acarel D, Plicka V, Bohnhoff M, Polat R, Zahradník J (2023) Coseismic
faulting complexity of the 2019 Mw 5.7 Silivri Earthquake in the Central
Marmara seismic gap offshore Istanbul. Seismol Res Lett 94:75–86.
https://​doi.​org/​10.​1785/​02202​20111
Ulrich T, Gabriel AA, Ampuero JP, Xu W (2019) Dynamic viability of the 2016
Mw 7.8 Kaikōura earthquake cascade on weak crustal faults. Nat Commun. https://​doi.​org/​10.​1038/​s41467-​019-​09125-w
Wang D, Chen Y, Wang Q, Mori J (2018a) Complex rupture of the 13 November
2016 Mw 7.8 Kaikoura, New Zealand earthquake: comparison of highfrequency and low-frequency observations. Tectonophysics 733:100–107.
https://​doi.​org/​10.​1016/j.​tecto.​2018.​02.​004
Wang T, Wei S, Shi X, Qiu Q, Li L, Peng D, Weldon RJ, Barbot S (2018b) The 2016
Kaikōura earthquake: simultaneous rupture of the subduction interface
and overlying faults. Earth Planet Sci Lett 482:44–51. https://​doi.​org/​10.​1016/j.​epsl.​2017.​10.​056
Williams CA, Eberhart-Phillips D, Bannister S, Barker DHN, Henrys S, Reyners
M, Sutherland R (2013) Revised interface geometry for the Hikurangi

Ohara et al. Progress in Earth and Planetary Science(2023) 10:35

subduction zone, New Zealand. Seismol Res Lett 84:1066–1073. https://​doi.​org/​10.​1785/​02201​30035
Xu W, Feng G, Meng L, Zhang A, Ampuero JP, Bürgmann R, Fang L (2018)
Transpressional rupture cascade of the 2016 Mw 7.8 Kaikoura Earthquake,
New Zealand. J Geophys Res Solid Earth 123:2396–2409. https://​doi.​org/​
10.​1002/​2017J​B0151​68
Yabuki T, Matsu’ura M (1992) Geodetic data inversion using a Bayesian
information criterion for spatial distribution of fault slip. Geophys J Int
109(2):363–375. https://​doi.​org/​10.​1111/j.​1365-​246X.​1992.​tb001​02.x
Yagi Y, Fukahata Y (2011) Introduction of uncertainty of Green’s function
into waveform inversion for seismic source processes. Geophys J Int
186:711–720. https://​doi.​org/​10.​1111/j.​1365-​246X.​2011.​05043.x
Yagi Y, Okuwaki R, Enescu B, Hirano S, Yamagami Y, Endo S, Komoro T (2014)
Rupture process of the 2014 Iquique Chile Earthquake in relation withthe foreshock activity. Geophys Res Lett 41:4201–4206. https://​doi.​org/​10.​
1002/​2014G​L0602​74
Yagi Y, Okuwaki R, Enescu B, Lu J (2023) Irregular rupture process of the 2022
Taitung, Taiwan, earthquake sequence. Sci Rep. https://​doi.​org/​10.​1038/​
s41598-​023-​27384-y
Yamashita S, Yagi Y, Okuwaki R, Shimizu K, Agata R, Fukahata Y (2021) Consecutive ruptures on a complex conjugate fault system during the 2018
Gulf of Alaska earthquake. Sci Rep 11:5979. https://​doi.​org/​10.​1038/​
s41598-​021-​85522-w
Yamashita S, Yagi Y, Okuwaki R (2022a) Irregular rupture propagation andgeometric fault complexities during the 2010 Mw 7.2 El Mayor-Cucapah
earthquake. Sci Rep 12:4575. https://​doi.​org/​10.​1038/​s41598-​022-​08671-6
Yamashita S, Yagi Y, Okuwaki R, Shimizu K, Agata R, Fukahata Y (2022b) Potency
density tensor inversion of complex body waveforms with time-adaptive
smoothing constraint. Geophys J Int 231:91–107. https://​doi.​org/​10.​1093/​
gji/​ggac1​81
Zhang H, Koper KD, Pankow K, Ge Z (2017) Imaging the 2016 Mw 7.8 Kaikoura,
New Zealand, earthquake with teleseismic P waves: a cascading rupture
across multiple faults. Geophys Res Lett 44:4790–4798. https://​doi.​org/​10.​
1002/​2017G​L0734​61
Zheng A, Wang M, Yu X, Zhang W (2018) Source rupture process of the2016 Kaikoura, New Zealand earthquake estimated from the kinematic
waveform inversion of strong-motion data. Geophys J Int 212:1736–1746.
https://​doi.​org/​10.​1093/​gji/​ggx505

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

Page 12 of 12