Long-period analysis of the 2016 Kaikoura earthquake

Z. Duputel¹ and L. Rivera¹

¹Institut de Physique du Globe de Strasbourg, UMR7516, Université de Strasbourg/EOST, CNRS, Strasbourg, France.

Abstract

The $M_w=7.8$ Kaikoura (New Zealand) earthquake involved a remarkably complex rupture propagating in an intricate network of faults at the transition between the Alpine fault in the South Island and the Kermadec-Tonga subduction zone. We investigate the main features of this complicated rupture process using long-period seismological observations. Apparent Rayleigh-wave moment-rate functions reveal a clear northeastward directivity with an unusually weak rupture initiation during 60 s followed by a major 20 s burst of moment rate. To further explore the rupture process, we perform a Bayesian exploration of multiple point-source parameters in a 3-D Earth model. The results show that the rupture initiated as a small strike-slip rupture and propagated to the northeast, triggering large slip on both strike-slip and thrust faults. The Kaikoura earthquake is thus a rare instance in which slip on intraplate faults trigger extensive interplate thrust faulting. This clearly outlines the importance of accounting for secondary faults when assessing seismic and tsunami hazard in subduction zones.

Keywords

New Zealand; Earthquake; Megathrust; Strike-slip; Surface waves

<footnote>Corresponding author: Z. Duputel, zacharie.duputel@unistra.fr</footnote>
**Highlights**

* We investigate the main features of the Kaikoura earthquake using long-period seismic waves

* The Kaikoura earthquake initiated as a small strike-slip rupture that propagated to the northeast.

* Slip on strike-slip faults triggered extensive thrust faulting under the eastern coast of the South Island

**1 Introduction**

On November 13, 2016, a large earthquake struck the north-east coast of the South Island in New Zealand (GeoNet hypocenter: latitude=-42.69°, longitude=173.02°, depth=14 km, O.T.=11:02:56 UTC; Global CMT $M_w=7.8$). This earthquake occurred in the Marlborough Fault system, an intricate network of right lateral strike-slip faults connecting the Alpine fault in the South Island to the Hikurangi subduction zone (cf., Fig. 1). A tsunami swept onto the coastlines with wave-heights of 2.5 m at Kaikoura (https://www.geonet.org.nz/tsunami). This earthquake is the largest event in the region since a magnitude 7.5 earthquake that occured 100 km to the northeast in October 1848 [Grapes et al., 1998; Mason and Little, 2006]. The 1848 earthquake ruptured ~100 km of the Awatere Fault near Blenheim with horizontal displacements as large as 6 m causing significant damages in Wellington and the Awatere valley.

Several facts indicate that the 2016 Kaikoura earthquake involved a complex rupture. Long-period moment tensor solutions indicate an oblique thrust focal mechanism with a large non-double couple component. Using the definition of Hara et al. [1996], the non-double couple component of Global CMT (GCMT) and USGS W-phase solutions are $\epsilon = -0.12$ and $\epsilon = -0.21$ respectively. This suggest that the mainshock is not well represented by a single fault plane. Both GCMT and W-phase solutions have large centroid time-delay $\tau_c$ ~ 57 s that indicate an anomalously long rupture duration, more than 2.5 times longer than what is expected from standard scaling laws [Duputel et al., 2013]. In addition, preliminary fields reports indicate that multiple faults were involved with surface strike-slip offsets as large as 10 m across the Kekerengu fault and coastal uplift between 2 and 5 m northeast of Kaikoura [Litchfield et al., 2016].

In this study, we investigate the mainshock rupture using long-period records available at teleseismic distances. Using this dataset, we conduct a directivity analysis using apparent
moment rate functions and perform a multiple-point-source inversion accounting for 3-D
Earth structures.

## 2 Rayleigh-wave moment rate functions

To study the time-history of the rupture and investigate possible directivity effects visible at long-period for the $M_w=7.8$ Kaikoura earthquake, we compute apparent Rayleigh-wave moment rate functions (MRFs). The dispersive wave-propagation effects are removed by deconvolving the data by point-source synthetic seismograms. To reduce biases in Rayleigh-wave MRFs due to unaccounted lateral heterogeneities, we use broadband (10-600 s) SEM syntethics computed for a 3D Earth model (S362ANI and CRUST2.0) using the spectral element code SPECFEM3D_GLOBE [Komatitsch and Tromp, 2002; Kustowski et al., 2008]. Deconvolution is performed using the projected Landweber deconvolution method [Bertero et al., 1999; Lanza et al., 1999] imposing causality and positivity [Duputel et al., 2016].

The MRFs shown in Fig. 2 indicate an unusual initiation with very small moment rate in the first 60s. Following this slow initiation phase, the moment rate rises abruptly during ~20s. This main energy burst shows clear azimuth-dependent time-shifts that are consistent with unilateral rupture propagation to the northeast. Such directivity is consistent with both GCMT and W-phase centroid location that lie ~120 km northeast of the epicenter. To study this apparent directivity, we image the spatial distribution of long-period seismic wave radiation. To do so, the MRFs are averaged in 10° azimuthal windows and back-projected over a gridded region around the epicenter assuming an average phase-velocity of 4 km/s. The resulting peak stacked amplitude shown in Fig. 1 indicates that this main moment-rate burst emanated from a region including the Kekerengu fault and the east coast of the upper South Island.

## 3 Multiple point source analysis of the 2016 Kaikoura earthquake

### 3.1 Multiple CMT inversion approach

We employ a strategy similar to Duputel et al. [2012a] where multiple moment tensor sources are inverted simultaneously using W-phase waveforms. Here, we extend this approach to surface waves using a larger time-window in the period range of 100-450s. The W-phase being mainly sensitive to first-order source parameters, incorporating surface waves and extending our pass-band to shorter periods improve our ability to capture more details of
the rupture process. While most of the W-phase energy propagates into the mantle and are therefore not strongly affected by shallow structures, fundamental mode surface waves are sensitive to shallow lateral heterogeneities such as those associated by the oceans and continents. As for Rayleigh-wave MRFs (section 2), such 3D effects are accounted for using SEM Green's functions that are computed for the 3D Earth model S362ANI.

Using a Bayesian approach, we invert for multiple double-couple parameters (i.e., the strike, dip, rake, seismic moment and centroid time of each sub-event), while sub-event locations are fixed based on field and satellite observations [e.g., Litchfield et al., 2016]. From Bayes' theorem, the posterior probability density function (PDF) is given by:

$$p(\mathbf{m}|\mathbf{d}_{\text{obs}}) \propto p(\mathbf{m}) p(\mathbf{d}_{\text{obs}}|\mathbf{m}), \quad (1)$$

where $p(\mathbf{d}_{\text{obs}}|\mathbf{m})$ is the likelihood function representing the ability of a source model $\mathbf{m}$ to fit the observations $\mathbf{d}_{\text{obs}}$ and $p(\mathbf{m})$ is the a priori PDF, which describes our prior knowledge about the source. We define the likelihood function as

$$p(\mathbf{d}_{\text{obs}}|\mathbf{m}) \propto \exp\left(-\frac{1}{2}[\mathbf{d}_{\text{obs}} - \mathbf{g}(\mathbf{m})]^T \mathbf{C}_{\chi}^{-1} [\mathbf{d}_{\text{obs}} - \mathbf{g}(\mathbf{m})]\right) \quad (2)$$

where $\mathbf{d}_{\text{obs}}$ and $\mathbf{g}(\mathbf{m})$ are the observed and predicted waveforms for a multiple-point-source model $\mathbf{m}$ and $\mathbf{C}_{\chi}$ is the misfit covariance matrix. We assume a block-diagonal $\mathbf{C}_{\chi}$, which is given by the following expression for a station n:

$$(\mathbf{C}_{\chi})_{n}^{ij} = (\sigma_{n})^{2} \exp(-|\Delta t^{ij}|/t_{0}). \quad (3)$$

In the previous equation, $\sigma_n$ is the data uncertainty at station n, $t_0$ represents a characteristic correlation duration and $\Delta t^{ij}$ is the time difference between data samples i and j [Duputel et al., 2012b]. As discussed in Minson et al. [2013] and Duputel et al. [2014], model prediction error is expected to be roughly proportional with data amplitude. We use a conservative estimate of 20% uncertainty (i.e., $\sigma_n = 0.2 \times \max(\mathbf{d}^n)$, where $\mathbf{d}^n$ is the long-period waveform recorded at station n) that is consistent with the variability in RMS misfits obtained using different 3D Earth models [e.g., S362ANI and S40RTS in Duputel et al., 2016]. The correlation duration $t_0$ is estimated to match the main lobe of the auto-correlation of the residual vector $\mathbf{r} = \mathbf{d}_{\text{obs}} - \mathbf{g}(\tilde{\mathbf{m}})$, where $\tilde{\mathbf{m}}$ is the maximum a posteriori model obtained from a first inversion. After different tests using different number of sub-events, we found that $t_0 = 30$ sec is appropriate for the dataset used in this study (cf., Fig. S1b).

Our solution is thus not a single "optimum" model but the full posterior PDF $p(\mathbf{m}|\mathbf{d}_{\text{obs}})$ describing the ensemble of source models that are consistent with observations and prior as-
sumptions. With this purpose, we employ a Monte Carlo strategy that allows us to generate
an ensemble of models that are statistically distributed according to $p(m|d_{\text{obs}})$. This is done
using an adaptive Metropolis Method [Brooks et al., 2011], which is described in the supple-
mentary material.

3.2 Results

We use records from 60 broadband stations within an epicentral distance of 90°. Be-
fore inversion, waveforms are deconvolved to displacement and band-pass filtered in the pe-
riod range of 100-450 s. Using this dataset, we perform a three-point source inversion with
sub-events placed along the main ruptures planes reported by field and satellite observations:
(1) Sub-event E1 is located in the Humps fault zone close to the epicenter, (2) Sub-event E2
is put along the Hope fault and Jordan thrust, (3) Sub-event E3 is placed along the Kekerengu
fault [Litchfield et al., 2016, http://comet.nerc.ac.uk, http://www.gsi.go.jp/cais/topic161117-
index-e.html]. Despite a limited sensitivity with respect to source depths, we observe smaller
RMS misfits when using shallow focal depths for sub-events E1 and E2 (~10 km) and a
deeper source for sub-events E3 (~20 km; cf., Table S2). The resulting solution shown in
Fig. 3a present oblique strike-slip mechanisms with strike angles that are globally consistent
with fault orientations. There is a dominant moment contribution from sub-event E3 with a
centroid time-shift of 73 s, in good agreement with the time and location of the main burst of
moment-rate seen in Fig. 1 and Fig. 2. Comparison between data and predictions in Fig. 3b
indicate that fundamental mode rayleigh-wave amplitudes are slightly overestimated for sta-
tions located north of the epicenter (in the Pacific Ocean) and underestimated for stations
located to the south (in Antarctic and the Southern Ocean). This suggest that a stronger di-
rectivity toward the north is needed to better match the observed waveforms.

To investigate the possibility that the rupture triggered slip along the Hikurangi megath-
rust, we add a fourth sub-event to our multiple point-source model. This additional point
source is also motivated by field reports and InSAR observations indicating significant ver-
tical displacement southeast of the Kekerengu fault that cannot be explained by our three-
point source model [3-5 m coastal uplift measured by Litchfield et al., 2016, up to 5m uplift
according to http://www.gsi.go.jp/cais/topic161117-index-e.html]. We thus perform a four
double-couple inversion adding sub-event E4 where large uplift is observed. The solution
presented in Fig. 4 show oblique strike-slip mechanisms for sub-events E1, E2 and E3 (quite
similar to our three double-couple solution in Fig. 3) and a pure thrust mechanism for sub-

110

111

112

113

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141
event E4. Waveform fits are significantly improved, which suggest that this additional thrust source fulfill the northward directivity that was missing in our three point-source model. Moreover, the composite mechanism in Fig. 4a obtained by summing moment tensors of the four sub-events is remarkably consistent with GCMT, while the composite three sub-event solution show larger differences (cf., Fig. 3a).

# 4 Discussion and conclusion

Long-period analysis of the Kaikoura earthquake reveals a rather unique source process. Apparent MRFs indicate a northeastward directivity with an unusually weak radiation during the first 60 s followed by a major 20 s burst of moment rate around the Kekerengu fault. This slow rupture initiation results in an anomalously long rupture duration, with a centroid time-shift more than 2.5 times longer than what is expected from standard scaling laws [Duputel et al., 2013].

In addition, our multiple point source analysis reveals a remarkably complex rupture process. The first sub-event (E1) suggests that the Kaikoura earthquake initiated as a small strike-slip rupture consistent with field reports of ~1 m horizontal offsets in the humps fault zone close to the epicenter [Litchfield et al., 2016]. As the rupture propagated to the north-east with large strike-slip motion (sub-event E3), significant thrust-slip was also triggered on a deeper shallow-dipping thrust fault (sub-event E4). This model explains W-phase and Global CMT single point-source solutions involving a significant non-double component with a large thrust slip motion (Fig. 4). This also agrees with field and satellite observations showing large horizontal offsets along the Kekerengu and Papatea faults accompanied by significant coastal uplifts along with 2.5 m tsunami waves recorded at Kaikoura (https://www.geonet.org.nz/tsunami). Notice that sub-event E2 involves a non-negligible thrust component suggesting that dip-slip motion might have started during the early stages of the rupture process. We are aware that this model only describe the first order features of the overall rupture process, which in detail might involve additional complexities (e.g., slip on the Papatea fault that is nearly orthogonal to the Kekerengu fault). We leave to further investigations the determination of a refined source model providing a detailed mapping of slip from available near-field geodetic and seismic data.

The strike-slip and thrust focal mechanisms of the two largest sub-events (E3 and E4) can be interpreted as slip on the Kekerengu fault splaying from a deeper thrust fault underlying
ing the Marlborough coastline. This previously undocumented shallow-dipping fault might correspond to the southwest extension of the Hikurangi Megathrust or some other thrust fault within the forearc wedge [Barnes and Audru, 1999]. The strike of sub-event E4 thrust mechanism is consistent with the azimuth of the Kekerengu Bank fault, although our probabilistic results yield a relatively large uncertainty on the strike angle. Activation of secondary faults are common during large megathrust earthquakes [Li et al., 2014]. Triggering of splay faults have been identified on many subduction margins and represent significant seismic and tsunami threats [Park et al., 2002; Wendt et al., 2009; Melnick et al., 2012]. However, the Kaikoura earthquake is a rare instance clearly showing that megathrust ruptures can also be triggered by slip on secondary faults. Another example of such interactions with intraplate earthquakes was observed during the 2009 Samoa-Tonga sequence in which an outer-rise event triggered extensive interplate faulting along the subduction interface [Lay et al., 2010]. Those rare instances reverse the typical pattern between megathrust and secondary faults and expand seismic and tsunami hazard in subduction zones.

## Acknowledgments

This work was supported by the Initiative d'Excellence (IDEX) funding framework (Université de Strasbourg). The research described herein used seismological data from various global networks available through the IRIS DMC (including Geoscope and IRIS/USGS stations in Antarctica and the Southern Ocean). We are grateful to the operators of these networks for ensuring the high quality of the data and making them publicly available. We also acknowledge the New Zealand GeoNet project and its sponsors EQC, GNS Science and LINZ, for providing data/images used in this study. We also thank COMET and GSI for providing InSAR images. This study contributed from fruitful discussions with Lingling Ye and James Hollingsworth. We thank the Editor, Vernon F. Cormier, and an anonymous reviewer for his constructive comments, which helped improve this manuscript.

## References

Barnes, P. M., and J. C. Audru (1999), Recognition of active strike-slip faulting from high-resolution marine seismic reflection profiles: Eastern Marlborough fault system, New Zealand, *Geol Soc America Bull*, 111(4), 538-559.

Bertero, M., D. Bindi, P. Boccacci, M. Cattaneo, C. Eva, and V. Lanza (1999), Application of the projected Landweber method to the estimation of the source time function in seismol-
204 ogy, *Inverse Problems*, **13**(2), 465–486, doi:10.1088/0266-5611/13/2/017.

205 Brooks, S., A. Gelman, G. L. Jones, and X.-L. Meng (2011), *Handbook of Markov Chain Monte Carlo*, Chapman & Hall/CRC, New York.

207 Duputel, Z., H. Kanamori, V. C. Tsai, and L. Rivera (2012a), The 2012 Sumatra great earthquake sequence, *Earth Planet. Sci. Lett.*, **351-352**, 247–257.

209 Duputel, Z., L. Rivera, Y. Fukahata, and H. Kanamori (2012b), Uncertainty estimations for seismic source inversions, *Geophys. J. Int.*, **190**(2), 1243–1256.

211 Duputel, Z., V. C. Tsai, L. Rivera, and H. Kanamori (2013), Using centroid time-delays to characterize source durations and identify earthquakes with unique characteristics, *Earth Planet. Sci. Lett.*, **374**, 92–100, doi:10.1016/j.epsl.2013.05.024.

214 Duputel, Z., P. S. Agram, M. Simons, S. E. Minson, and J. L. Beck (2014), Accounting for prediction uncertainty when inferring subsurface fault slip, *Geophys. J. Int.*, **197**(1), 464–482, doi:10.1093/gji/ggt517.

217 Duputel, Z., J. Vergne, L. Rivera, G. Wittlinger, V. Farra, and G. Hetényi (2016), The 2015 Gorkha earthquake: A large event illuminating the Main Himalayan Thrust fault, *Geophys. Res. Lett.*, **43**(6), 2517–2525, doi:10.1002/2016GL068083.

220 Grapes, R., T. Little, and G. Downes (1998), Rupturing of the Awatere Fault during the 1848 October 16 Marlborough earthquake, New Zealand: historical and present day evidence, *New Zealand Journal of Geology and Geophysics*, **41**(4), 387–399.

223 Hara, T., K. Kuge, and H. Kawakatsu (1996), Determination of the isotropic component of deep focus earthquakes by inversion of normal-mode data, *Geophys. J. Int.*, **127**, 515–528.

225 Komatitsch, D., and J. Tromp (2002), Spectral-element simulations of global seismic wave propagation—II. Three-dimensional models, oceans, rotation and self-gravitation - Komatitsch - 2002 - Geophysical Journal International - Wiley Online Library, *Geophys. J. Int.*

229 Kustowski, B., G. Ekström, and A. M. Dziewonski (2008), Anisotropic shear-wave velocity structure of the Earth’s mantle: A global model, *J. Geophys. Res.*, **113**, B06,306.

231 Lanza, V., D. Spallarossa, M. Cattaneo, D. Bindi, and P. Augliera (1999), Source parameters of small events using constrained deconvolution with empirical Green’s functions, *Geo-phys. J. Int.*, **137**(3), 651–662, doi:10.1046/j.1365-246x.1999.00809.x.

234 Lay, T., C. J. Ammon, H. Kanamori, L. Rivera, K. Koper, and A. R. Hutko (2010), The 2009 Samoa-Tonga great earthquake triggered doublet, *Nature*, **466**(7309), 964–968.
236 Li, S., M. Moreno, M. Rosenau, D. Melnick, and O. Oncken (2014), Splay fault triggering by great subduction earthquakes inferred from finite element models, *Geophys. Res. Lett.*, 41(2), 385-391, doi:10.1002/2013GL058598.

239 Litchfield, N. J., D. J. A. Barrell, J. G. Begg, A. Benson, K. R. Berryman, K. J. Clark, U. A. Cochran, S. C. Cox, C. Gasston, P. J. Glassey, D. W. Heron, J. D. Howarth, R. M. Langridge, T. Little, B. Lukovic, A. Nicol, J. Pettinga, W. F. Ries, R. J. M. W. Stirling, D. B. Townsend, P. Upton, R. J. Van Dissen, and P. Villamor (2016), 4th November 2016 Kaikoura Earthquake. Preliminary earthquake geology observations, *GNS Science*, doi: 10.21420/G2MW2D.

245 Mason, D. P. M., and T. A. Little (2006), Refined slip distribution and moment magnitude of the 1848 Marlborough earthquake, Awatere Fault, New Zealand, *New Zealand Journal of Geology and Geophysics*, 49(3), 375-382, doi:10.1080/00288306.2006.9515174.

248 Melnick, D., M. Moreno, M. Motagh, M. Cisternas, and R. L. Wesson (2012), Splay fault slip during the Mw 8.8 2010 Maule Chile earthquake, *Geology*, 40(3), 251-254, doi: 10.1130/G32712.1.

251 Minson, S. E., M. Simons, and J. L. Beck (2013), Bayesian inversion for finite fault earthquake source models I—theory and algorithm, *Geophys. J. Int.*, 194(3), 1701-1726, doi: 10.1093/gji/ggt180.

254 Park, J.-O., T. Tsuru, S. Kodaira, P. R. Cummins, and Y. Kaneda (2002), Splay Fault Branching Along the Nankai Subduction Zone, *Science*, 297(5584), 1157-1160, doi: 10.1126/science.1074111.

257 Wendt, J., D. D. Oglesby, and E. L. Geist (2009), Tsunamis and splay fault dynamics, *Geo-phys. Res. Lett.*, 36(15), n/a-n/a, doi:10.1029/2009GL038295.
<figure><img src="image_1.png" /><figcaption>**Figure 1.** The 2016 Kaikoura earthquake sequence - Red mechanisms correspond to our preferred four-point-source model obtained using body waves and surface waves assuming a 3-D Earth model. Red colors indicate peak-stacked amplitude in the source region from moment rate functions (MRF) backprojected relative to the main shock epicentral location. Blue mechanisms are the Global CMT solutions obtained for $M_W \ge 5$ aftershocks (2016/11/14 to 2016/11/22). Red lines are faults traces from the New Zealand Active Faults Database (GNS Science, https://data.gns.cri.nz/af). Yellow dashed line shows the approximate trench location of the Hikurangi megathrust.</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Rayleigh-wave moment rate functions - Apparent moment rate functions (MRFs) are shown as a function of azimuth and colored by epicentral distance. These represent seismic moment as a function of time observed at different stations. The black curve corresponds to the maximum MRF peak amplitude, showing the predicted arrival time of energy radiated from this location (cf., Fig. 1).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Three double-couple inversion - (a) Red mechanisms correspond to our three-point-source model obtained from Bayesian inversion of long-period seismic waves assuming a 3-D Earth model. Grey lines indicate the posterior population of double-couple mechanisms and the black lines correspond to the posterior mean model. Top left inset shows the timing and moment-rate function of each sub-event. Bottom right inset presents a comparison of the Global CMT (GCMT) solution and a composite mechanism corresponding to the sum of moment tensors for sub-events E1, E2 and E3. (b) Comparison between data (black) and synthetic (red) waveforms for representative stations. The part of the signal used for multiple point source determination is delimited by red dots. Yellow circles in right insets show the global distribution of stations used for the inversion. The station azimuth ($\phi$) and epicentral distance ($\Delta$) is indicated on top of each seismic trace.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>(b) Model performance</figcaption></figure>

<figure><img src="image_5.png" /></figure>

**Figure 4.** **Four double-couple inversion - (a)** Red mechanisms correspond to our four-point-source model obtained from Bayesian inversion of long-period seismic waves assuming a 3-D Earth model. Grey lines indicate the posterior population of double-couple mechanisms and the black lines correspond to the posterior mean model (cf., Table S1 and S3). Top left inset shows the timing and moment-rate function of each sub-event. Bottom right inset presents a comparison of the Global CMT (GCMT) solution and a composite mechanism corresponding to the sum of moment tensors for sub-events E1, E2, E3 and E4. **(b)** Comparison between data (black) and synthetic (red) waveforms for representative stations. The part of the signal used for multiple point source determination is delimited by red dots. Yellow circles in right insets show the global distribution of stations used for the inversion. The station azimuth ($\phi$) and epicentral distance ($\Delta$) is indicated on top of each seismic trace. Additional waveform comparisons for this inversion are shown in the supplementary material.