# Source rupture process of the 2016 Kaikoura, New Zealand earthquake estimated from the kinematic waveform inversion of strong-motion data

Source rupture process of the 2016 Kaikoura, New Zealand earthquake estimated from the kinematic... Summary On 2016 November 13, an Mw 7.8 earthquake occurred in the northeast of the South Island of New Zealand near Kaikoura. The earthquake caused severe damages and great impacts on local nature and society. Referring to the tectonic environment and defined active faults, the field investigation and geodetic evidence reveal that at least 12 fault sections ruptured in the earthquake, and the focal mechanism is one of the most complicated in historical earthquakes. On account of the complexity of the source rupture, we propose a multisegment fault model based on the distribution of surface ruptures and active tectonics. We derive the source rupture process of the earthquake using the kinematic waveform inversion method with the multisegment fault model from strong-motion data of 21 stations (0.05–0.35 Hz). The inversion result suggests the rupture initiates in the epicentral area near the Humps fault, and then propagates northeastward along several faults, until the offshore Needles fault. The Mw 7.8 event is a mixture of right-lateral strike and reverse slip, and the maximum slip is approximately 19 m. The synthetic waveforms reproduce the characteristics of the observed ones well. In addition, we synthesize the coseismic offsets distribution of the ruptured region from the slips of upper subfaults in the fault model, which is roughly consistent with the surface breaks observed in the field survey. New Zealand, Waveform inversion, Body waves, Theoretical seismology INTRODUCTION A great earthquake with Mw 7.8 occurred in Kaikoura, New Zealand, at 11:02:56 AM (UTC), 2016 November 13, which was the most powerful earthquake experienced in the past more than 150 yr in the area. Shaking was felt almost throughout New Zealand, and serious damages were found along the northeastern part of the South Island, which significantly impacted the natural and social environments there. Apart from numerous buildings and structures destroyed by the great shaking, the earthquake also generated severe secondary disasters. Helicopter survey indicated this earthquake triggered tens of thousands of landslides (Dellow et al.2017), and caused the destruction of many roads and railways. Moreover, the coseismic uplifts were observed in areas of the East Coast in the upper South Island, which even led to tsunami waves. The Institute of Geological and Nuclear Sciences in New Zealand reported the earthquake epicentre was at 42.69°S, 173.02°E, with a depth of about 15 km. The earthquake occurred within the tectonically active and complex Australian-Pacific plate boundary region, which controls the tectonics of New Zealand. As a result of the oblique convergence between two plates, the Pacific plate is subducting beneath the Australian plate along the Hikurangi subduction zone in the east of the North Island. While in the south, the tectonics is dominated by the Alpine fault system with purely dextral motion. The Marlborough fault system in the central region of the northern South Island is the transition from subduction zone to right-lateral faults, including four major faults, the Wairau, Awatere, Clarence and Hope faults. Numerous smaller active faults between major faults with diverse orientations form a complicated fault network across this area (Shi et al. 2017). Based on the analysis of focal mechanisms of aftershocks, the earthquake rupture region showed a mixture of reverse and dextral faulting. The moment tensor solution of the main shock released by the Global CMT Project (GCMT, Ekström et al.2012) was also consistent with the feature, which had a rake angle of 143°. It had been revealed that the Kaikoura earthquake was one of the most complex earthquakes ever recorded. Ground deformation in conjunction with field observation demonstrated that at least 12 fault segments ruptured, including several that might be not mapped previously (Ref. GeoNet, http://info.geonet.org.nz/pages/viewpage.action?pageId=20971550, last accessed at: 2017 March 25). The earthquake seemed to initiate near the Humps fault zone, in the north of the epicentre, and then propagated northeastward to the Leader fault, the Conway-Charwell fault and the Hundalee fault. After that, the rupture jumped northward through the Upper-Kowhai fault, the Jordan thrust, the Kekerengu fault and the Needles fault. Surface displacement was observed on the Fidget and the Papatea fault as well. Among these ruptured faults, the Hundalee fault and the Needles fault extended across the coast and offshore. The number of ruptured faults and the rupture propagating pattern showed an extreme complexity, implying that we cannot investigate the source process following an ordinary perspective. In source process analysis routines, we are accustomed to using a single-fault plane as the source model which is commonly adapted from a moment tensor solution or the distribution of active tectonics. However, it is apparent that a source model derived from this method is not enough to describe the complicated rupture process of the Kaikoura earthquake. The observation of surface rupture and the change in moment tensor solutions of aftershocks suggested that the source model should have transformations of focal mechanism in the rupture procedure of the event. In this study, we proposed a multisegment fault model based on the spatial distribution of surface rupture (Litchfield et al.2016), active faults (Litchfield et al.2014) and aftershocks. Using the multisegment fault model, we estimated the source process of the Mw 7.8 with strong-motion data. In addition, we compared the inversion result with that derived from a single-segment fault model to demonstrate the significance of using the multisegment fault model. METHOD AND DATA In this paper, we used the multitime-window linear waveform inversion method (Olson & Apsel 1982; Hartzell & Heaton 1983; Sekiguchi et al.2000, 2002; Asano et al.2005) to estimate the kinematic source rupture process of the 2016 Kaikoura earthquake. In this method, the observation (displacement, velocity and acceleration) is discretized in space and time based on the representation theorem (Aki & Richards 2002). For discretization in space, a finite-fault plane is divided into small subfaults, and each subfault is regarded as a point source. For discretization in the time domain, a moment release history on a subfault is represented by several time-window functions. Thus, the observation equation is described as follows:   $$\dot{U}\left( t \right) = \sum\limits_{if}^{nf} {\sum\limits_{itw}^{ntw} {\sum\limits_{is}^{ns} {{m_{if,itw,is}}T\left( {t - \Delta {t_{{\rm{trig}}}}} \right)*{{\dot{G}}_{if,is}}\left( t \right)} } }$$ (1)with   \begin{equation*}1 \le if \le nf,\quad 1 \le itw \le ntw,\quad 1 \le is \le ns\left( {ns = 2} \right),\end{equation*} where mif, itw, is is the seismic moment, T(t) is a smoothed ramp function used as the time-window function and Δttrig = R/Vr + Δtw(itm − 1), the first time window is triggered by a constant velocity Vr, R is the distance between a subfault and the rupture starting point, Δtw is the time-window interval, Gif, is(t) is the Green's function, $$\dot{U}( t )$$ is the synthetic velocity waveform and the subscripts if, itw, is are indexes for subfault, time window and slip direction, respectively. Arbitrary directed moment release on each subfault is realized by two orthonormal unit vectors, thus the number of slip directions ns is 2. Therefore, we obtain the observation equation in a vector form as follows:   $$\left[ {\begin{array}{@{}*{1}{c}@{}} {{\boldsymbol A}}\\ {\lambda {{\boldsymbol S}}} \end{array}} \right]{{\boldsymbol m}} = \left[ {\begin{array}{@{}*{1}{c}@{}} {{\boldsymbol d}}\\ {{\boldsymbol 0}} \end{array}} \right]$$ (2)where matrix A is composed of the derivative Green's function with time-window functions; S and λ are, respectively, a matrix of a smoothing constraint posed on the model parameters to stabilize the inversion problem and a scalar coefficient to control the smoothing strength; m is a model vector composed of the seismic moment and d is a data vector composed of the observed waveforms. Data and Green's functions for each station were normalized by the maximum amplitude of the observed waveforms for three components. To determine the proper value of λ, we used Akaike's Bayesian Information Criterion (ABIC, Akaike 1980). Yoshida (1989) and Ide et al. (1996) introduced ABIC into source modeling. We applied non-negative constraint to prohibit back slip, which is limiting the rake angle variation within a certain range of ±45° centred at the prescribed slip angle of assumed source mechanisms. The non-negative least-squares algorithm developed by Lawson & Hanson (1974) was used to solve the observation equation. We used three-component strong-motion acceleration waveform records at 21 stations of GeoNet strong-motion database (http://info.-geonetorg.nz/display/appdata/Strong-Motion+Data, last accessed at: 2016 November 23, Van Houtte et al.2017). The location distribution of stations used is shown in Fig. 1. The epicentral distance of these stations ranged from 8 to 259 km. The original data records should have been pre-processed before inversion calculation. At first, the acceleration waveforms were numerically integrated in time domain into velocity. Next, the velocity records were bandpass filtered between 0.05 and 0.35 Hz, and resampled at 2 Hz. Because of interest in the S-wave portion of the records, we used the length of 100 s data starting from 1 s before the S-wave arrival, which was identified by visual inspection and adjustment. Green's functions were calculated by the discrete wave number method (Bouchon 1981) and the reflection/transmission coefficient matrix method (Kennett & Kerry 1979), assuming a horizontally layered velocity structure. The 1-D velocity structure model extracted from the New Zealand wide 3-D velocity model (Eberhart-Phillips et al.2010) is available in the Supporting Information. Figure 1. View largeDownload slide Location of the strong-motion stations (triangle) used in this study and the main shock (star) with its aftershocks (circle) happened in 24 hr. (a) 13 solid lines donates the projection trace of the multisegment fault model (the strike side close to the surface), and focal mechanisms of inversion result and GCMT are shown in the bottom right (b). The rectangle donates the projection of the single-segment fault model, and the solid line is the strike side close to the surface. Figure 1. View largeDownload slide Location of the strong-motion stations (triangle) used in this study and the main shock (star) with its aftershocks (circle) happened in 24 hr. (a) 13 solid lines donates the projection trace of the multisegment fault model (the strike side close to the surface), and focal mechanisms of inversion result and GCMT are shown in the bottom right (b). The rectangle donates the projection of the single-segment fault model, and the solid line is the strike side close to the surface. Fault model As a general rule, the source rupture process is described by a model containing a single-fault plane based on the moment tensor solution or the epicentre distribution of aftershocks. It is conventional to examine the fault geometry which should cover the range of surface ruptures, aftershocks, or active faults. So, we constructed a single-fault-plane source model to conduct a preliminary inversion. The length and width of the planar fault plane were set as 204 and 60 km, respectively. The strike (226°), dip angle (33°) and the central rake angle (143°) were followed the moment tensor solution determined by GCMT. The projection of the fault model is shown in Fig. 1(b). The fault plane was divided into 34 subfaults along the strike direction and 10 subfaults along the dip direction, amount to 340 subfaults, each with a size of 6 km × 6 km. The slip time history of each subfault was constrained into seven time windows, each with a width of 2.6 s, and a lag of 1.3 s. Thus, the allowed slip duration for each subfault was 10.4 s. The first time-window starting time was defined as the time prescribed by a circular rupture propagation with the constant velocity Vr. We performed inversions using several combinations of Vr and smoothing strength coefficient λ. Fig. 2 shows the distribution of slips inverted from the single-segment fault model, in which large slips concentrated on the top of the plane, and the maximum was close to 17 m. The waveform fit coefficient, Wxc, was included as a criterion for evaluating the effect of waveform fit,   \begin{eqnarray}{W_{xc}} \!=\! \frac{1}{{stn}}\sum\limits_{ist \!=\! 1}^{stn} {\frac{{\sum\limits_{icmp \!=\! 1}^{ncmp} {\sum\limits_{it \!=\! 1}^{nt} {{D_{ist,it,icmp}} \cdot {S_{ist,it,icmp}}} } }}{{{{\left[ {\sum\limits_{icmp \!=\! 1}^{ncmp} {\sum\limits_{it \!=\! 1}^{nt} {{D_{ist,it,icmp}}^2} } \cdot \sum\limits_{icmp \!=\! 1}^{ncmp} {\sum\limits_{it \!=\! 1}^{nt} {{S_{ist,it,icmp}}^2} } } \right]}^{1/2}}}}}.\nonumber\\\end{eqnarray} (3) Figure 2. View largeDownload slide Slip distribution derived from the single-segment fault model. The star indicates the rupture staring point, and arrows indicate the slip vector. The contour interval is 2 m. Figure 2. View largeDownload slide Slip distribution derived from the single-segment fault model. The star indicates the rupture staring point, and arrows indicate the slip vector. The contour interval is 2 m. Here, stn is the number of used stations, nt is the length of the waveform record, ncmp is the number of components of the record, Dist, it, icmp represents the observed waveform record and Sist, it, icmp represents the synthesized one. 0 ≤ Wxc ≤ 1, and the larger Wxc is, the better the waveform fit between observed waveforms and synthetic waveforms is. Fig. 3 shows the comparison of observed waveforms and synthetic waveforms and Wxc is 0.53. The synthetic waveforms did not fit the observed data well, and the peak ground velocity of every component at all stations was smaller than the observed one. In order to improve the waveform fit, we increased the amplitude of time window. However, even if the corresponding magnitude of the released energy had exceeded the measured magnitude, the waveform fit was still not improved. This indicated that the source model, the single-segment fault, was not enough to explain the observed data. Although the spatial projection range of the fault plane covered most epicentres of aftershocks, it only meant that the geometric dimension of the model was close to that of ruptured faults. Hence, a single focal mechanism was far from the realistic situation, and such a model was not enough to extract information about the source rupture process. Figure 3. View largeDownload slide Comparison between the synthetic waveforms obtained from the inversion with the single-segment fault model and the observed ones (velocity). The maximum values of the observed and synthetic waveforms are shown at the upper right of each waveform, respectively. Figure 3. View largeDownload slide Comparison between the synthetic waveforms obtained from the inversion with the single-segment fault model and the observed ones (velocity). The maximum values of the observed and synthetic waveforms are shown at the upper right of each waveform, respectively. Referring to the epicentre distribution of 24-hr aftershocks, the seismicity of the Kaikoura earthquake was mainly focused in the northeastern region of South Island. The hypocentre of the main shock was located in the southwest of this region, thus the rupture started in the southwest of the whole earthquake rupture area, and then propagated northeastward until to Cape Campbell. It is obvious that the rupture propagating direction was changed near Kaikoura with a contrast of about 30° inferred from the distribution of aftershocks in Fig. 1(a). As stated in the Introduction, the rupture area is in the transition zone between Hikurangi subduction and Alpine fault, and the sense of slip on faults in this area are diverse from each other, such a single-fault plane is certainly not adequate to represent the change in the focal mechanism. Observed surface breaks suggested that faults ruptured along a snake-shape pattern in the south with branches in the northern part. Based on the distribution of active tectonics and the trace of surface breaks, we constructed a fault model consisted of 13 segments with different orientations and dip angles. The size of subfaults was changed into 5 km × 5 km. Each segment was not simply ruptured in order from the starting point to its spatial position. We calculated the triggering time of the first time window on each subfault to approach the complicated rupture pattern of the Kaikoura earthquake. On account of the existence of surface breaks, each fault plane of this model has a width of 30 km and a common top depth of approximately 5 m. The multisegment fault model represented the probable seismogenic tectonics, which consisted of at least 12 ruptured faults. Detailed information of each segment is listed in Table 1, including strike direction, dip angle, length and width and rake angle variation. We estimated the source rupture process using the multisegment fault model, and inversion parameters were the same with those in the preliminary inversion by the single-segment model. As for Vr and λ, the solution was considered as the most appropriate when the ABIC value was minimum. Five values of Vr−2.10, 2.29, 2.48, 2.67 and 2.87 km s−1 (corresponding to 55, 60, 65, 70 and 75 per cent of the S-wave velocity at the depth of the hypocentre) and ten values of λ (0.001–0.01, with an interval of 0.001) were tested. These two parameters were selected as 2.29 km s−1 and 0.004, respectively (Fig. 4). Figure 4. View largeDownload slide Trade-off curves of ABIC with Vr (first time-window triggering velocity) and λ (smoothing strength coefficient). Figure 4. View largeDownload slide Trade-off curves of ABIC with Vr (first time-window triggering velocity) and λ (smoothing strength coefficient). Table 1. Parameters of the multisegment fault model. Number of segment  1  2  3  4  5  6  7  8  9  10  11  12  13  Strike  226°  242°  215°  157°  255°  245°  181°  220°  157°  250°  180°  265°  216°  Dip  80°  60°  37°  60°  60°  60°  60°  55°  60°  70°  60°  60°  60°  Length (km)  35  20  25  20  25  20  15  20  20  10  20  35  25  Width (km)  30  30  30  30  30  30  30  30  30  30  30  30  30  Rake angle variation  143° ± 45°  143° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  143° ± 45°  Number of segment  1  2  3  4  5  6  7  8  9  10  11  12  13  Strike  226°  242°  215°  157°  255°  245°  181°  220°  157°  250°  180°  265°  216°  Dip  80°  60°  37°  60°  60°  60°  60°  55°  60°  70°  60°  60°  60°  Length (km)  35  20  25  20  25  20  15  20  20  10  20  35  25  Width (km)  30  30  30  30  30  30  30  30  30  30  30  30  30  Rake angle variation  143° ± 45°  143° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  143° ± 45°  View Large RESULTS AND DISCUSSION Fig. 5 shows the resulting slip distribution derived from the multisegment fault model. The obtained source model had a released seismic moment of 6.34 × 1020 N · m (Mw 7.8). Large slips mainly concentrated at shallow depth on the Kekerengu fault (segments 2 and 3) and the Jordan thrust (segment 4), on which the observed surface ruptures were also large. The maximum slip was 19.0 m. The slips on segments 4, 7, 9 and 11 were close to reverse slip, consistent with the rupture interpretation based on Sentinel-1 SAR and ALOS-2 data (Shi et al.2017). Other segments performed a right-lateral strike or oblique reverse slip motion. Fig. 6 compares the observed velocity waveforms and the synthetic waveforms. The synthetic waveforms reproduced the characteristics of the observed waveforms well, and the fit of almost every station was getting greatly improved compared with the result obtained by a single-segment fault model as shown in Fig. 3. The waveform fit coefficient derived from the multisegment model is 0.61, which is increased by about 15 per cent compared to that using the single-segment model. Figure 5. View largeDownload slide Slip distribution derived from the multisegment fault model. The star indicates the rupture staring point, and arrows indicate the slip vector. The contour interval is 2 m. Figure 5. View largeDownload slide Slip distribution derived from the multisegment fault model. The star indicates the rupture staring point, and arrows indicate the slip vector. The contour interval is 2 m. Figure 6. View largeDownload slide Comparison between the synthetic waveforms obtained from the inversion with the multisegment fault model and the observed ones (velocity). The maximum values of the observed and synthetic waveforms are shown at the upper right of each waveform, respectively. Figure 6. View largeDownload slide Comparison between the synthetic waveforms obtained from the inversion with the multisegment fault model and the observed ones (velocity). The maximum values of the observed and synthetic waveforms are shown at the upper right of each waveform, respectively. Fig. 7 shows the rupture progression of the source model. The total duration of the rupture was approximately 100 s. The rupture started at the hypocentre on segment 13, and then propagated toward segments 12 and 11. After the rupture front arrived on segment 3, the rupture also propagated toward the other two branches, segments 4 and 5. As shown by the moment rate function for each subfault in Fig. 8, most part of the seismic moment was released by the Kekerengu fault and the Jordan thrust. These two faults ruptured about 40 s after the onset of the whole progression, and lasted for approximately 40 s. That was also obvious in the moment rate functions of individual fault segments and the entire fault plane presented in Fig. 9. We synthesized the focal mechanism from the cumulative moment tensor of the inversion result, and showed it in Fig. 1(a) (‘beach ball’). We used MoPaD (a moment tensor analysis tool, Krieger & Heimann 2012) to analyse the cumulative moment tensor and the GCMT, and they both included a non-double-couple (non-DC) component. The cumulative moment tensor had a 95 per cent double-couple (DC) component, which indicated a mixture of reverse and strike slip faulting, and a 5 per cent CLVD (compensated linear vector dipoles) component. The GCMT had a 68 per cent DC component and a strong CLVD component about 32 per cent. In case of decomposing the moment tensor, the CLVD could be considered as a sum of several DC sources, so the 5 per cent CLVD component of the cumulative moment tensor might be due to near-simultaneous ruptures on nearby fault segments with different geometries, but it was negligible compared to the CLVD component of the GCMT. Hence, our model indicated that the strong-motion records could be explained well without a non-DC component. Considering the rupture area was in the Hikurangi subduction zone, which could imply the existence of the possible slip on the subducting slab interface (Hamling et al.2017). While the fault model used in the study was mainly located in the crust, it was hard to verify the possible slip. In other words, similar to the geodetic measurements (Hamling et al.2017), the bandlimited (0.05–0.35 Hz) near source strong-motion records might also be insensitive to the slip on the interface, particularly if the rupture had long rise time. Figure 7. View largeDownload slide The rupture process derived from the multisegment fault model. Figure 7. View largeDownload slide The rupture process derived from the multisegment fault model. Figure 8. View largeDownload slide Moment rate functions of subfaults obtained from the multisegment fault model. The background is the final slip distribution on the fault plane. Figure 8. View largeDownload slide Moment rate functions of subfaults obtained from the multisegment fault model. The background is the final slip distribution on the fault plane. Figure 9. View largeDownload slide Moment rate functions of individual segments and the entire fault obtained from the multisegment fault model. Panels are arranged in the rupture propagating order of segments. Figure 9. View largeDownload slide Moment rate functions of individual segments and the entire fault obtained from the multisegment fault model. Panels are arranged in the rupture propagating order of segments. Field investigation and ground deformation data indicated the existence of surface breaks, so the fault rupture of this earthquake should have reached the ground surface. Permanent offsets associated with surface ruptures were mainly caused by the coseismic displacement. The depth of the top of the source model used in this study, which was constructed based on the distribution of surface ruptures, was almost close to the ground surface. The inversion result of the earthquake showed that subfaults with large slips concentrated on the top of the fault model as well. Thus, slips on the uppermost subfaults could be thought to generate ground ruptures with coseismic offsets. The horizontal and vertical components of surface offsets were synthesized from the projection of slips on the uppermost subfaults, whose distribution is shown in Fig. 10. The horizontal offsets of the Jordan thrust and the Kekerengu fault were very large, while relatively small in the southern fault array. Vertical offsets were produced on the most of segments, suggesting the uplift of the Australian plate in the side of the hangingwall. This is consistent with the feature of the convergent plate boundary in the area. In the numerical comparison, the result shown in Fig. 10 is just roughly similar with the observed surface offsets. As the uppermost subfaults were close to the surface, but each subfault was treated as a point source at its centre, such inverted slip or released moment was actually the sum of faulting in its geometrical dimension. The length and width of subfaults in the fault model were both 5 km. Therefore, the resolution scale resulted from the subfault size was not quite fit for the realistic variation of surface ruptures, but of certain significance from the perspective of the spread of earthquake damages. Figure 10. View largeDownload slide The distribution of coseismic offsets synthetized from slips of shallow subfaults. (a) Horizontal offsets. (b) Vertical offsets. Figure 10. View largeDownload slide The distribution of coseismic offsets synthetized from slips of shallow subfaults. (a) Horizontal offsets. (b) Vertical offsets. The source rupture of the 2016 Kaikoura earthquake contained extreme complexity, in which the most important was the focal mechanism, including numerous fault tectonics, with different orientations and variation of the slip direction. When we are building up a source model, parameters of the model are set according to available geological or geophysical data, or adjusted in an inversion procedure. The synthetic waveforms could fit the observed record pretty well, but this result is only the ‘best solution’ based on the prior condition and the selected data, which might be not close to the nature. The multitime-window linear waveform inversion method used in the study, in which the propagation triggering velocity of the first time window is prescribed as a constant. Affected by the propagation velocity, the rupture front arrival time of each subfault is fixed in a duration range determined by subjectivity. In the assumed circular rupture pattern, the rupture time of a subfault in a source model depends on the distance between the subfault and the rupture starting point or the hypocentre. However, there seems to be stepovers in the Kaikoura earthquake. Faults in the northeast of the South Island, including the Needles fault, the Kekerengu fault, the Jordan thrust and the Upper Kowhai fault, join together and form a complete fault array. The fault array is not directly linked with the faults in the south around the epicentre, instead of a gap zone with a length more than 10 km. Thus, in the propagation of the fault rupture from south to north, how the rupture passes over the gap is worth considering carefully. In addition, the distribution of faults around the epicentre is unusual, and the hypocentre is unsurely located in a defined fault. Moreover, how the rupture propagates among these faults is not certainly determined. As if there are stepovers in the rupture process, it would be difficult to constrain the rupture time of subfaults by the inversion method used in the paper, unless through other technique or data. In the conventional kinematic inversion method, we used to set the source model as a single-fault plane, in which the rupture propagates successively in a circular or bilateral pattern. However, discontinuous rupture usually occurs in the complicated multifault earthquake, for example, the stepover mentioned above, so the conventional method may no longer meet the actual situation. In summary, the fact implies that we should reconsider the conventional research method of source rupture and earthquake hazard assessment in the region with complex fault tectonics from a new perspective. CONCLUSIONS The source process of the 2016 Kaikoura earthquake was estimated by the multitime-window linear waveform inversion method from strong-motion data. In order to analyse the rupture process of the Mw 7.8 event, we proposed a fault model consisted of 13 segments. The inversion result showed a south-to-north rupture, starting at the epicentre and terminated near the Cook Strait. The duration of the whole rupture process was about 100 s, and most large slips concentrated at the shallow depth. An apparent asperity zone occurred along the Kekerengu fault and the Jordan thrust, and the rupture duration of the asperity zone was approximately 40 s. The distribution of slips on the fault model suggested the mechanism of the event was a mixture of right-lateral strike and reverse slip. The total released seismic moment was 6.34 × 1020 N · m (Mw 7.8), and the maximum slip was 19.0 m. A comparison with the result obtained by using a single-fault-plane model demonstrated that the use of the multisegment fault model led to improved waveform fit at most stations. The coseismic offset synthesized from slips of shallow subfaults were roughly consistent with the distribution of surface breaks observed in the field investigation. ACKNOWLEDGEMENTS Strong-motion data and hypocentre catalogue used here were released by GeoNET. We also used SRTM DEM data (Jarvis et al.2008) and Generic Mapping Tools (Wessel & Smith 1995) to draw figures in the paper. This work is supported by the National Natural Science Foundation of China (grants 41674056 and 41374105) and the CAS/CAFEA international partnership Program for creative research teams (grant KZZD-EW-TZ-19). The comments from the reviewer Nori Nakata (University of Oklahoma), another anonymous reviewer and the editor, Prof Jean Virieux are helpful in improving the manuscript. Here, we would like to express our gratitude together. REFERENCES Akaike H., 1980. Likelihood and the Bayes procedure, in Bayesian Statics , pp. 143–166, eds Bernardo J.M., DeGroot M.H., Lindlely D.V., Smith A.F.M., University Press, Valencia, Spain. Aki K., Richards P.G., 2002. Quantitative Seismology , 2nd edn, University Science Books, Sausalito, CA. Asano K., Iwata T., Irikura K., 2005. Estimation of source rupture process and strong ground motion simulation of the 2002 Denali, Alaska, earthquake, Bull. seism. Soc. Am. , 95( 5), 1701– 1715. https://doi.org/10.1785/0120040154 Google Scholar CrossRef Search ADS   Bouchon M., 1981. A simple method to calculate Green's functions for elastic layered media, Bull. seism. Soc. Am. , 71( 4), 959– 971. Dellow S. et al.  , 2017. Landslides caused by the Mw 7.8 Kaikoura earthquake and the immediate response, Bull. N. Z. Soc. Earthq. Eng. , 50( 2), 106– 116. Eberhart-Phillips D., Reyners M.E., Bannister S., Chadwick M.P., Ellis S.M., 2010. Establishing a versatile 3-D seismic velocity model for New Zealand, Seismol. Res. Lett. , 81( 6), 992– 1000. https://doi.org/10.1785/gssrl.81.6.992 Google Scholar CrossRef Search ADS   Ekström G., Nettles M., Dziewoński A.M., 2012. The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes, Phys. Earth planet. Inter. , 200, 1– 9. https://doi.org/10.1016/j.pepi.2012.04.002 Google Scholar CrossRef Search ADS   Hamling I.J. et al.  , 2017. Complex multifault rupture during the 2016 Mw 7.8 Kaikoura earthquake, New Zealand, Science , 356( 6334), doi:10.1126/science.aam7194. https://doi.org/10.1126/science.aam7194 Hartzell S.H., Heaton T.H., 1983. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake, Bull. seism. Soc. Am. , 73, 1553– 1583. Ide S., Takeo M., Yoshida Y., 1996. Source process of the 1995 Kobe earthquake: determination of spatio-temporal slip distribution by Bayesian modeling, Bull. seism. Soc. Am. , 86, 547– 566. Jarvis A., Reuter H.I., Nelson A., Guevara E., 2008. Hole-filled SRTM for the globe Version 4, Available from the CGIAR-CSI SRTM 90 m Database: http://srtm.csi.cgiar.org, last accessed 30 November 2016. Kennett B.L.N., Kerry N.J., 1979. Seismic waves in a stratified half space, Geophys. J. Int. , 57( 3), 557– 583. https://doi.org/10.1111/j.1365-246X.1979.tb06779.x Google Scholar CrossRef Search ADS   Krieger L., Heimann S., 2012. MoPaD–moment tensor plotting and decomposition: a tool for graphical and numerical analysis of seismic moment tensors, Seismol. Res. Lett. , 83( 3), 589– 595. https://doi.org/10.1785/gssrl.83.3.589 Google Scholar CrossRef Search ADS   Lawson C.L., Hanson R.J., 1974. Solving Least Squares Problems , Vol. 161, Prentice-Hall. Litchfield N.J. et al.  , 2014. A model of active faulting in New Zealand, N. Z. J. Geol. Geophys. , 57( 1), 32– 56. https://doi.org/10.1080/00288306.2013.854256 Google Scholar CrossRef Search ADS   Litchfield N.J. et al.  , 2016. 14th November 2016 M 7.8 Kaikoura earthquake. Preliminary surface fault displacement measurements, Version 2, GNS Science , doi:10.21420/G2J01F. Olson A.H., Apsel R.J., 1982. Finite faults and inverse theory with applications to the 1979 Imperial Valley earthquake, Bull. seism. Soc. Am. , 72, 1969– 2001. Sekiguchi H., Irikura K., Iwata T., 2000. Fault geometry at the rupture termination of the 1995 Hyogo-ken Nanbu earthquake, Bull. seism. Soc. Am. , 90( 1), 117– 133. https://doi.org/10.1785/0119990027 Google Scholar CrossRef Search ADS   Sekiguchi H., Irikura K., Iwata T., 2002. Source inversion for estimating the continuous slip distribution on a fault-introduction of Green's functions convolved with a correction function to give moving dislocation effects in subfaults, Geophys. J. Int. , 150( 2), 377– 391. https://doi.org/10.1046/j.1365-246X.2002.01669.x Google Scholar CrossRef Search ADS   Shi X., Wang Y., Liu-Zeng J., Weldon R., Wei S., Wang T., Sieh K., 2017. How complex is the 2016 M w 7.8 Kaikoura earthquake, South Island, New Zealand?, Sci. Bull. , 62( 5), 309– 311. https://doi.org/10.1016/j.scib.2017.01.033 Google Scholar CrossRef Search ADS   Van Houtte C., Bannister S., Holden C., Bourguignon S., McVerry G., 2017. The New Zealand strong motion database, Bull. N. Z. Soc. Earthq. Eng. , 50( 1), 1– 20. Wessel P., Smith W.H.F., 1995. New version of the generic mapping tools, EOS, Trans. Am. geophys. Un. , 76( 33), 329– 329. https://doi.org/10.1029/95EO00198 Google Scholar CrossRef Search ADS   Yoshida S., 1989. Waveform inversion using ABIC for the rupture process of the 1983 Hindu Kush earthquake, Phys. Earth planet. Inter. , 56( 3–4), 389– 405. https://doi.org/10.1016/0031-9201(89)90172-6 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. supplementary data.zip Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Source rupture process of the 2016 Kaikoura, New Zealand earthquake estimated from the kinematic waveform inversion of strong-motion data

, Volume 212 (3) – Mar 1, 2018
11 pages

/lp/ou_press/source-rupture-process-of-the-2016-kaikoura-new-zealand-earthquake-y0XLCjrHDn
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx505
Publisher site
See Article on Publisher Site

### Abstract

Summary On 2016 November 13, an Mw 7.8 earthquake occurred in the northeast of the South Island of New Zealand near Kaikoura. The earthquake caused severe damages and great impacts on local nature and society. Referring to the tectonic environment and defined active faults, the field investigation and geodetic evidence reveal that at least 12 fault sections ruptured in the earthquake, and the focal mechanism is one of the most complicated in historical earthquakes. On account of the complexity of the source rupture, we propose a multisegment fault model based on the distribution of surface ruptures and active tectonics. We derive the source rupture process of the earthquake using the kinematic waveform inversion method with the multisegment fault model from strong-motion data of 21 stations (0.05–0.35 Hz). The inversion result suggests the rupture initiates in the epicentral area near the Humps fault, and then propagates northeastward along several faults, until the offshore Needles fault. The Mw 7.8 event is a mixture of right-lateral strike and reverse slip, and the maximum slip is approximately 19 m. The synthetic waveforms reproduce the characteristics of the observed ones well. In addition, we synthesize the coseismic offsets distribution of the ruptured region from the slips of upper subfaults in the fault model, which is roughly consistent with the surface breaks observed in the field survey. New Zealand, Waveform inversion, Body waves, Theoretical seismology INTRODUCTION A great earthquake with Mw 7.8 occurred in Kaikoura, New Zealand, at 11:02:56 AM (UTC), 2016 November 13, which was the most powerful earthquake experienced in the past more than 150 yr in the area. Shaking was felt almost throughout New Zealand, and serious damages were found along the northeastern part of the South Island, which significantly impacted the natural and social environments there. Apart from numerous buildings and structures destroyed by the great shaking, the earthquake also generated severe secondary disasters. Helicopter survey indicated this earthquake triggered tens of thousands of landslides (Dellow et al.2017), and caused the destruction of many roads and railways. Moreover, the coseismic uplifts were observed in areas of the East Coast in the upper South Island, which even led to tsunami waves. The Institute of Geological and Nuclear Sciences in New Zealand reported the earthquake epicentre was at 42.69°S, 173.02°E, with a depth of about 15 km. The earthquake occurred within the tectonically active and complex Australian-Pacific plate boundary region, which controls the tectonics of New Zealand. As a result of the oblique convergence between two plates, the Pacific plate is subducting beneath the Australian plate along the Hikurangi subduction zone in the east of the North Island. While in the south, the tectonics is dominated by the Alpine fault system with purely dextral motion. The Marlborough fault system in the central region of the northern South Island is the transition from subduction zone to right-lateral faults, including four major faults, the Wairau, Awatere, Clarence and Hope faults. Numerous smaller active faults between major faults with diverse orientations form a complicated fault network across this area (Shi et al. 2017). Based on the analysis of focal mechanisms of aftershocks, the earthquake rupture region showed a mixture of reverse and dextral faulting. The moment tensor solution of the main shock released by the Global CMT Project (GCMT, Ekström et al.2012) was also consistent with the feature, which had a rake angle of 143°. It had been revealed that the Kaikoura earthquake was one of the most complex earthquakes ever recorded. Ground deformation in conjunction with field observation demonstrated that at least 12 fault segments ruptured, including several that might be not mapped previously (Ref. GeoNet, http://info.geonet.org.nz/pages/viewpage.action?pageId=20971550, last accessed at: 2017 March 25). The earthquake seemed to initiate near the Humps fault zone, in the north of the epicentre, and then propagated northeastward to the Leader fault, the Conway-Charwell fault and the Hundalee fault. After that, the rupture jumped northward through the Upper-Kowhai fault, the Jordan thrust, the Kekerengu fault and the Needles fault. Surface displacement was observed on the Fidget and the Papatea fault as well. Among these ruptured faults, the Hundalee fault and the Needles fault extended across the coast and offshore. The number of ruptured faults and the rupture propagating pattern showed an extreme complexity, implying that we cannot investigate the source process following an ordinary perspective. In source process analysis routines, we are accustomed to using a single-fault plane as the source model which is commonly adapted from a moment tensor solution or the distribution of active tectonics. However, it is apparent that a source model derived from this method is not enough to describe the complicated rupture process of the Kaikoura earthquake. The observation of surface rupture and the change in moment tensor solutions of aftershocks suggested that the source model should have transformations of focal mechanism in the rupture procedure of the event. In this study, we proposed a multisegment fault model based on the spatial distribution of surface rupture (Litchfield et al.2016), active faults (Litchfield et al.2014) and aftershocks. Using the multisegment fault model, we estimated the source process of the Mw 7.8 with strong-motion data. In addition, we compared the inversion result with that derived from a single-segment fault model to demonstrate the significance of using the multisegment fault model. METHOD AND DATA In this paper, we used the multitime-window linear waveform inversion method (Olson & Apsel 1982; Hartzell & Heaton 1983; Sekiguchi et al.2000, 2002; Asano et al.2005) to estimate the kinematic source rupture process of the 2016 Kaikoura earthquake. In this method, the observation (displacement, velocity and acceleration) is discretized in space and time based on the representation theorem (Aki & Richards 2002). For discretization in space, a finite-fault plane is divided into small subfaults, and each subfault is regarded as a point source. For discretization in the time domain, a moment release history on a subfault is represented by several time-window functions. Thus, the observation equation is described as follows:   $$\dot{U}\left( t \right) = \sum\limits_{if}^{nf} {\sum\limits_{itw}^{ntw} {\sum\limits_{is}^{ns} {{m_{if,itw,is}}T\left( {t - \Delta {t_{{\rm{trig}}}}} \right)*{{\dot{G}}_{if,is}}\left( t \right)} } }$$ (1)with   \begin{equation*}1 \le if \le nf,\quad 1 \le itw \le ntw,\quad 1 \le is \le ns\left( {ns = 2} \right),\end{equation*} where mif, itw, is is the seismic moment, T(t) is a smoothed ramp function used as the time-window function and Δttrig = R/Vr + Δtw(itm − 1), the first time window is triggered by a constant velocity Vr, R is the distance between a subfault and the rupture starting point, Δtw is the time-window interval, Gif, is(t) is the Green's function, $$\dot{U}( t )$$ is the synthetic velocity waveform and the subscripts if, itw, is are indexes for subfault, time window and slip direction, respectively. Arbitrary directed moment release on each subfault is realized by two orthonormal unit vectors, thus the number of slip directions ns is 2. Therefore, we obtain the observation equation in a vector form as follows:   $$\left[ {\begin{array}{@{}*{1}{c}@{}} {{\boldsymbol A}}\\ {\lambda {{\boldsymbol S}}} \end{array}} \right]{{\boldsymbol m}} = \left[ {\begin{array}{@{}*{1}{c}@{}} {{\boldsymbol d}}\\ {{\boldsymbol 0}} \end{array}} \right]$$ (2)where matrix A is composed of the derivative Green's function with time-window functions; S and λ are, respectively, a matrix of a smoothing constraint posed on the model parameters to stabilize the inversion problem and a scalar coefficient to control the smoothing strength; m is a model vector composed of the seismic moment and d is a data vector composed of the observed waveforms. Data and Green's functions for each station were normalized by the maximum amplitude of the observed waveforms for three components. To determine the proper value of λ, we used Akaike's Bayesian Information Criterion (ABIC, Akaike 1980). Yoshida (1989) and Ide et al. (1996) introduced ABIC into source modeling. We applied non-negative constraint to prohibit back slip, which is limiting the rake angle variation within a certain range of ±45° centred at the prescribed slip angle of assumed source mechanisms. The non-negative least-squares algorithm developed by Lawson & Hanson (1974) was used to solve the observation equation. We used three-component strong-motion acceleration waveform records at 21 stations of GeoNet strong-motion database (http://info.-geonetorg.nz/display/appdata/Strong-Motion+Data, last accessed at: 2016 November 23, Van Houtte et al.2017). The location distribution of stations used is shown in Fig. 1. The epicentral distance of these stations ranged from 8 to 259 km. The original data records should have been pre-processed before inversion calculation. At first, the acceleration waveforms were numerically integrated in time domain into velocity. Next, the velocity records were bandpass filtered between 0.05 and 0.35 Hz, and resampled at 2 Hz. Because of interest in the S-wave portion of the records, we used the length of 100 s data starting from 1 s before the S-wave arrival, which was identified by visual inspection and adjustment. Green's functions were calculated by the discrete wave number method (Bouchon 1981) and the reflection/transmission coefficient matrix method (Kennett & Kerry 1979), assuming a horizontally layered velocity structure. The 1-D velocity structure model extracted from the New Zealand wide 3-D velocity model (Eberhart-Phillips et al.2010) is available in the Supporting Information. Figure 1. View largeDownload slide Location of the strong-motion stations (triangle) used in this study and the main shock (star) with its aftershocks (circle) happened in 24 hr. (a) 13 solid lines donates the projection trace of the multisegment fault model (the strike side close to the surface), and focal mechanisms of inversion result and GCMT are shown in the bottom right (b). The rectangle donates the projection of the single-segment fault model, and the solid line is the strike side close to the surface. Figure 1. View largeDownload slide Location of the strong-motion stations (triangle) used in this study and the main shock (star) with its aftershocks (circle) happened in 24 hr. (a) 13 solid lines donates the projection trace of the multisegment fault model (the strike side close to the surface), and focal mechanisms of inversion result and GCMT are shown in the bottom right (b). The rectangle donates the projection of the single-segment fault model, and the solid line is the strike side close to the surface. Fault model As a general rule, the source rupture process is described by a model containing a single-fault plane based on the moment tensor solution or the epicentre distribution of aftershocks. It is conventional to examine the fault geometry which should cover the range of surface ruptures, aftershocks, or active faults. So, we constructed a single-fault-plane source model to conduct a preliminary inversion. The length and width of the planar fault plane were set as 204 and 60 km, respectively. The strike (226°), dip angle (33°) and the central rake angle (143°) were followed the moment tensor solution determined by GCMT. The projection of the fault model is shown in Fig. 1(b). The fault plane was divided into 34 subfaults along the strike direction and 10 subfaults along the dip direction, amount to 340 subfaults, each with a size of 6 km × 6 km. The slip time history of each subfault was constrained into seven time windows, each with a width of 2.6 s, and a lag of 1.3 s. Thus, the allowed slip duration for each subfault was 10.4 s. The first time-window starting time was defined as the time prescribed by a circular rupture propagation with the constant velocity Vr. We performed inversions using several combinations of Vr and smoothing strength coefficient λ. Fig. 2 shows the distribution of slips inverted from the single-segment fault model, in which large slips concentrated on the top of the plane, and the maximum was close to 17 m. The waveform fit coefficient, Wxc, was included as a criterion for evaluating the effect of waveform fit,   \begin{eqnarray}{W_{xc}} \!=\! \frac{1}{{stn}}\sum\limits_{ist \!=\! 1}^{stn} {\frac{{\sum\limits_{icmp \!=\! 1}^{ncmp} {\sum\limits_{it \!=\! 1}^{nt} {{D_{ist,it,icmp}} \cdot {S_{ist,it,icmp}}} } }}{{{{\left[ {\sum\limits_{icmp \!=\! 1}^{ncmp} {\sum\limits_{it \!=\! 1}^{nt} {{D_{ist,it,icmp}}^2} } \cdot \sum\limits_{icmp \!=\! 1}^{ncmp} {\sum\limits_{it \!=\! 1}^{nt} {{S_{ist,it,icmp}}^2} } } \right]}^{1/2}}}}}.\nonumber\\\end{eqnarray} (3) Figure 2. View largeDownload slide Slip distribution derived from the single-segment fault model. The star indicates the rupture staring point, and arrows indicate the slip vector. The contour interval is 2 m. Figure 2. View largeDownload slide Slip distribution derived from the single-segment fault model. The star indicates the rupture staring point, and arrows indicate the slip vector. The contour interval is 2 m. Here, stn is the number of used stations, nt is the length of the waveform record, ncmp is the number of components of the record, Dist, it, icmp represents the observed waveform record and Sist, it, icmp represents the synthesized one. 0 ≤ Wxc ≤ 1, and the larger Wxc is, the better the waveform fit between observed waveforms and synthetic waveforms is. Fig. 3 shows the comparison of observed waveforms and synthetic waveforms and Wxc is 0.53. The synthetic waveforms did not fit the observed data well, and the peak ground velocity of every component at all stations was smaller than the observed one. In order to improve the waveform fit, we increased the amplitude of time window. However, even if the corresponding magnitude of the released energy had exceeded the measured magnitude, the waveform fit was still not improved. This indicated that the source model, the single-segment fault, was not enough to explain the observed data. Although the spatial projection range of the fault plane covered most epicentres of aftershocks, it only meant that the geometric dimension of the model was close to that of ruptured faults. Hence, a single focal mechanism was far from the realistic situation, and such a model was not enough to extract information about the source rupture process. Figure 3. View largeDownload slide Comparison between the synthetic waveforms obtained from the inversion with the single-segment fault model and the observed ones (velocity). The maximum values of the observed and synthetic waveforms are shown at the upper right of each waveform, respectively. Figure 3. View largeDownload slide Comparison between the synthetic waveforms obtained from the inversion with the single-segment fault model and the observed ones (velocity). The maximum values of the observed and synthetic waveforms are shown at the upper right of each waveform, respectively. Referring to the epicentre distribution of 24-hr aftershocks, the seismicity of the Kaikoura earthquake was mainly focused in the northeastern region of South Island. The hypocentre of the main shock was located in the southwest of this region, thus the rupture started in the southwest of the whole earthquake rupture area, and then propagated northeastward until to Cape Campbell. It is obvious that the rupture propagating direction was changed near Kaikoura with a contrast of about 30° inferred from the distribution of aftershocks in Fig. 1(a). As stated in the Introduction, the rupture area is in the transition zone between Hikurangi subduction and Alpine fault, and the sense of slip on faults in this area are diverse from each other, such a single-fault plane is certainly not adequate to represent the change in the focal mechanism. Observed surface breaks suggested that faults ruptured along a snake-shape pattern in the south with branches in the northern part. Based on the distribution of active tectonics and the trace of surface breaks, we constructed a fault model consisted of 13 segments with different orientations and dip angles. The size of subfaults was changed into 5 km × 5 km. Each segment was not simply ruptured in order from the starting point to its spatial position. We calculated the triggering time of the first time window on each subfault to approach the complicated rupture pattern of the Kaikoura earthquake. On account of the existence of surface breaks, each fault plane of this model has a width of 30 km and a common top depth of approximately 5 m. The multisegment fault model represented the probable seismogenic tectonics, which consisted of at least 12 ruptured faults. Detailed information of each segment is listed in Table 1, including strike direction, dip angle, length and width and rake angle variation. We estimated the source rupture process using the multisegment fault model, and inversion parameters were the same with those in the preliminary inversion by the single-segment model. As for Vr and λ, the solution was considered as the most appropriate when the ABIC value was minimum. Five values of Vr−2.10, 2.29, 2.48, 2.67 and 2.87 km s−1 (corresponding to 55, 60, 65, 70 and 75 per cent of the S-wave velocity at the depth of the hypocentre) and ten values of λ (0.001–0.01, with an interval of 0.001) were tested. These two parameters were selected as 2.29 km s−1 and 0.004, respectively (Fig. 4). Figure 4. View largeDownload slide Trade-off curves of ABIC with Vr (first time-window triggering velocity) and λ (smoothing strength coefficient). Figure 4. View largeDownload slide Trade-off curves of ABIC with Vr (first time-window triggering velocity) and λ (smoothing strength coefficient). Table 1. Parameters of the multisegment fault model. Number of segment  1  2  3  4  5  6  7  8  9  10  11  12  13  Strike  226°  242°  215°  157°  255°  245°  181°  220°  157°  250°  180°  265°  216°  Dip  80°  60°  37°  60°  60°  60°  60°  55°  60°  70°  60°  60°  60°  Length (km)  35  20  25  20  25  20  15  20  20  10  20  35  25  Width (km)  30  30  30  30  30  30  30  30  30  30  30  30  30  Rake angle variation  143° ± 45°  143° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  143° ± 45°  Number of segment  1  2  3  4  5  6  7  8  9  10  11  12  13  Strike  226°  242°  215°  157°  255°  245°  181°  220°  157°  250°  180°  265°  216°  Dip  80°  60°  37°  60°  60°  60°  60°  55°  60°  70°  60°  60°  60°  Length (km)  35  20  25  20  25  20  15  20  20  10  20  35  25  Width (km)  30  30  30  30  30  30  30  30  30  30  30  30  30  Rake angle variation  143° ± 45°  143° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  90° ± 45°  143° ± 45°  143° ± 45°  View Large RESULTS AND DISCUSSION Fig. 5 shows the resulting slip distribution derived from the multisegment fault model. The obtained source model had a released seismic moment of 6.34 × 1020 N · m (Mw 7.8). Large slips mainly concentrated at shallow depth on the Kekerengu fault (segments 2 and 3) and the Jordan thrust (segment 4), on which the observed surface ruptures were also large. The maximum slip was 19.0 m. The slips on segments 4, 7, 9 and 11 were close to reverse slip, consistent with the rupture interpretation based on Sentinel-1 SAR and ALOS-2 data (Shi et al.2017). Other segments performed a right-lateral strike or oblique reverse slip motion. Fig. 6 compares the observed velocity waveforms and the synthetic waveforms. The synthetic waveforms reproduced the characteristics of the observed waveforms well, and the fit of almost every station was getting greatly improved compared with the result obtained by a single-segment fault model as shown in Fig. 3. The waveform fit coefficient derived from the multisegment model is 0.61, which is increased by about 15 per cent compared to that using the single-segment model. Figure 5. View largeDownload slide Slip distribution derived from the multisegment fault model. The star indicates the rupture staring point, and arrows indicate the slip vector. The contour interval is 2 m. Figure 5. View largeDownload slide Slip distribution derived from the multisegment fault model. The star indicates the rupture staring point, and arrows indicate the slip vector. The contour interval is 2 m. Figure 6. View largeDownload slide Comparison between the synthetic waveforms obtained from the inversion with the multisegment fault model and the observed ones (velocity). The maximum values of the observed and synthetic waveforms are shown at the upper right of each waveform, respectively. Figure 6. View largeDownload slide Comparison between the synthetic waveforms obtained from the inversion with the multisegment fault model and the observed ones (velocity). The maximum values of the observed and synthetic waveforms are shown at the upper right of each waveform, respectively. Fig. 7 shows the rupture progression of the source model. The total duration of the rupture was approximately 100 s. The rupture started at the hypocentre on segment 13, and then propagated toward segments 12 and 11. After the rupture front arrived on segment 3, the rupture also propagated toward the other two branches, segments 4 and 5. As shown by the moment rate function for each subfault in Fig. 8, most part of the seismic moment was released by the Kekerengu fault and the Jordan thrust. These two faults ruptured about 40 s after the onset of the whole progression, and lasted for approximately 40 s. That was also obvious in the moment rate functions of individual fault segments and the entire fault plane presented in Fig. 9. We synthesized the focal mechanism from the cumulative moment tensor of the inversion result, and showed it in Fig. 1(a) (‘beach ball’). We used MoPaD (a moment tensor analysis tool, Krieger & Heimann 2012) to analyse the cumulative moment tensor and the GCMT, and they both included a non-double-couple (non-DC) component. The cumulative moment tensor had a 95 per cent double-couple (DC) component, which indicated a mixture of reverse and strike slip faulting, and a 5 per cent CLVD (compensated linear vector dipoles) component. The GCMT had a 68 per cent DC component and a strong CLVD component about 32 per cent. In case of decomposing the moment tensor, the CLVD could be considered as a sum of several DC sources, so the 5 per cent CLVD component of the cumulative moment tensor might be due to near-simultaneous ruptures on nearby fault segments with different geometries, but it was negligible compared to the CLVD component of the GCMT. Hence, our model indicated that the strong-motion records could be explained well without a non-DC component. Considering the rupture area was in the Hikurangi subduction zone, which could imply the existence of the possible slip on the subducting slab interface (Hamling et al.2017). While the fault model used in the study was mainly located in the crust, it was hard to verify the possible slip. In other words, similar to the geodetic measurements (Hamling et al.2017), the bandlimited (0.05–0.35 Hz) near source strong-motion records might also be insensitive to the slip on the interface, particularly if the rupture had long rise time. Figure 7. View largeDownload slide The rupture process derived from the multisegment fault model. Figure 7. View largeDownload slide The rupture process derived from the multisegment fault model. Figure 8. View largeDownload slide Moment rate functions of subfaults obtained from the multisegment fault model. The background is the final slip distribution on the fault plane. Figure 8. View largeDownload slide Moment rate functions of subfaults obtained from the multisegment fault model. The background is the final slip distribution on the fault plane. Figure 9. View largeDownload slide Moment rate functions of individual segments and the entire fault obtained from the multisegment fault model. Panels are arranged in the rupture propagating order of segments. Figure 9. View largeDownload slide Moment rate functions of individual segments and the entire fault obtained from the multisegment fault model. Panels are arranged in the rupture propagating order of segments. Field investigation and ground deformation data indicated the existence of surface breaks, so the fault rupture of this earthquake should have reached the ground surface. Permanent offsets associated with surface ruptures were mainly caused by the coseismic displacement. The depth of the top of the source model used in this study, which was constructed based on the distribution of surface ruptures, was almost close to the ground surface. The inversion result of the earthquake showed that subfaults with large slips concentrated on the top of the fault model as well. Thus, slips on the uppermost subfaults could be thought to generate ground ruptures with coseismic offsets. The horizontal and vertical components of surface offsets were synthesized from the projection of slips on the uppermost subfaults, whose distribution is shown in Fig. 10. The horizontal offsets of the Jordan thrust and the Kekerengu fault were very large, while relatively small in the southern fault array. Vertical offsets were produced on the most of segments, suggesting the uplift of the Australian plate in the side of the hangingwall. This is consistent with the feature of the convergent plate boundary in the area. In the numerical comparison, the result shown in Fig. 10 is just roughly similar with the observed surface offsets. As the uppermost subfaults were close to the surface, but each subfault was treated as a point source at its centre, such inverted slip or released moment was actually the sum of faulting in its geometrical dimension. The length and width of subfaults in the fault model were both 5 km. Therefore, the resolution scale resulted from the subfault size was not quite fit for the realistic variation of surface ruptures, but of certain significance from the perspective of the spread of earthquake damages. Figure 10. View largeDownload slide The distribution of coseismic offsets synthetized from slips of shallow subfaults. (a) Horizontal offsets. (b) Vertical offsets. Figure 10. View largeDownload slide The distribution of coseismic offsets synthetized from slips of shallow subfaults. (a) Horizontal offsets. (b) Vertical offsets. The source rupture of the 2016 Kaikoura earthquake contained extreme complexity, in which the most important was the focal mechanism, including numerous fault tectonics, with different orientations and variation of the slip direction. When we are building up a source model, parameters of the model are set according to available geological or geophysical data, or adjusted in an inversion procedure. The synthetic waveforms could fit the observed record pretty well, but this result is only the ‘best solution’ based on the prior condition and the selected data, which might be not close to the nature. The multitime-window linear waveform inversion method used in the study, in which the propagation triggering velocity of the first time window is prescribed as a constant. Affected by the propagation velocity, the rupture front arrival time of each subfault is fixed in a duration range determined by subjectivity. In the assumed circular rupture pattern, the rupture time of a subfault in a source model depends on the distance between the subfault and the rupture starting point or the hypocentre. However, there seems to be stepovers in the Kaikoura earthquake. Faults in the northeast of the South Island, including the Needles fault, the Kekerengu fault, the Jordan thrust and the Upper Kowhai fault, join together and form a complete fault array. The fault array is not directly linked with the faults in the south around the epicentre, instead of a gap zone with a length more than 10 km. Thus, in the propagation of the fault rupture from south to north, how the rupture passes over the gap is worth considering carefully. In addition, the distribution of faults around the epicentre is unusual, and the hypocentre is unsurely located in a defined fault. Moreover, how the rupture propagates among these faults is not certainly determined. As if there are stepovers in the rupture process, it would be difficult to constrain the rupture time of subfaults by the inversion method used in the paper, unless through other technique or data. In the conventional kinematic inversion method, we used to set the source model as a single-fault plane, in which the rupture propagates successively in a circular or bilateral pattern. However, discontinuous rupture usually occurs in the complicated multifault earthquake, for example, the stepover mentioned above, so the conventional method may no longer meet the actual situation. In summary, the fact implies that we should reconsider the conventional research method of source rupture and earthquake hazard assessment in the region with complex fault tectonics from a new perspective. CONCLUSIONS The source process of the 2016 Kaikoura earthquake was estimated by the multitime-window linear waveform inversion method from strong-motion data. In order to analyse the rupture process of the Mw 7.8 event, we proposed a fault model consisted of 13 segments. The inversion result showed a south-to-north rupture, starting at the epicentre and terminated near the Cook Strait. The duration of the whole rupture process was about 100 s, and most large slips concentrated at the shallow depth. An apparent asperity zone occurred along the Kekerengu fault and the Jordan thrust, and the rupture duration of the asperity zone was approximately 40 s. The distribution of slips on the fault model suggested the mechanism of the event was a mixture of right-lateral strike and reverse slip. The total released seismic moment was 6.34 × 1020 N · m (Mw 7.8), and the maximum slip was 19.0 m. A comparison with the result obtained by using a single-fault-plane model demonstrated that the use of the multisegment fault model led to improved waveform fit at most stations. The coseismic offset synthesized from slips of shallow subfaults were roughly consistent with the distribution of surface breaks observed in the field investigation. ACKNOWLEDGEMENTS Strong-motion data and hypocentre catalogue used here were released by GeoNET. We also used SRTM DEM data (Jarvis et al.2008) and Generic Mapping Tools (Wessel & Smith 1995) to draw figures in the paper. This work is supported by the National Natural Science Foundation of China (grants 41674056 and 41374105) and the CAS/CAFEA international partnership Program for creative research teams (grant KZZD-EW-TZ-19). The comments from the reviewer Nori Nakata (University of Oklahoma), another anonymous reviewer and the editor, Prof Jean Virieux are helpful in improving the manuscript. Here, we would like to express our gratitude together. REFERENCES Akaike H., 1980. Likelihood and the Bayes procedure, in Bayesian Statics , pp. 143–166, eds Bernardo J.M., DeGroot M.H., Lindlely D.V., Smith A.F.M., University Press, Valencia, Spain. Aki K., Richards P.G., 2002. Quantitative Seismology , 2nd edn, University Science Books, Sausalito, CA. Asano K., Iwata T., Irikura K., 2005. Estimation of source rupture process and strong ground motion simulation of the 2002 Denali, Alaska, earthquake, Bull. seism. Soc. Am. , 95( 5), 1701– 1715. https://doi.org/10.1785/0120040154 Google Scholar CrossRef Search ADS   Bouchon M., 1981. A simple method to calculate Green's functions for elastic layered media, Bull. seism. Soc. Am. , 71( 4), 959– 971. Dellow S. et al.  , 2017. Landslides caused by the Mw 7.8 Kaikoura earthquake and the immediate response, Bull. N. Z. Soc. Earthq. Eng. , 50( 2), 106– 116. Eberhart-Phillips D., Reyners M.E., Bannister S., Chadwick M.P., Ellis S.M., 2010. Establishing a versatile 3-D seismic velocity model for New Zealand, Seismol. Res. Lett. , 81( 6), 992– 1000. https://doi.org/10.1785/gssrl.81.6.992 Google Scholar CrossRef Search ADS   Ekström G., Nettles M., Dziewoński A.M., 2012. The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes, Phys. Earth planet. Inter. , 200, 1– 9. https://doi.org/10.1016/j.pepi.2012.04.002 Google Scholar CrossRef Search ADS   Hamling I.J. et al.  , 2017. Complex multifault rupture during the 2016 Mw 7.8 Kaikoura earthquake, New Zealand, Science , 356( 6334), doi:10.1126/science.aam7194. https://doi.org/10.1126/science.aam7194 Hartzell S.H., Heaton T.H., 1983. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake, Bull. seism. Soc. Am. , 73, 1553– 1583. Ide S., Takeo M., Yoshida Y., 1996. Source process of the 1995 Kobe earthquake: determination of spatio-temporal slip distribution by Bayesian modeling, Bull. seism. Soc. Am. , 86, 547– 566. Jarvis A., Reuter H.I., Nelson A., Guevara E., 2008. Hole-filled SRTM for the globe Version 4, Available from the CGIAR-CSI SRTM 90 m Database: http://srtm.csi.cgiar.org, last accessed 30 November 2016. Kennett B.L.N., Kerry N.J., 1979. Seismic waves in a stratified half space, Geophys. J. Int. , 57( 3), 557– 583. https://doi.org/10.1111/j.1365-246X.1979.tb06779.x Google Scholar CrossRef Search ADS   Krieger L., Heimann S., 2012. MoPaD–moment tensor plotting and decomposition: a tool for graphical and numerical analysis of seismic moment tensors, Seismol. Res. Lett. , 83( 3), 589– 595. https://doi.org/10.1785/gssrl.83.3.589 Google Scholar CrossRef Search ADS   Lawson C.L., Hanson R.J., 1974. Solving Least Squares Problems , Vol. 161, Prentice-Hall. Litchfield N.J. et al.  , 2014. A model of active faulting in New Zealand, N. Z. J. Geol. Geophys. , 57( 1), 32– 56. https://doi.org/10.1080/00288306.2013.854256 Google Scholar CrossRef Search ADS   Litchfield N.J. et al.  , 2016. 14th November 2016 M 7.8 Kaikoura earthquake. Preliminary surface fault displacement measurements, Version 2, GNS Science , doi:10.21420/G2J01F. Olson A.H., Apsel R.J., 1982. Finite faults and inverse theory with applications to the 1979 Imperial Valley earthquake, Bull. seism. Soc. Am. , 72, 1969– 2001. Sekiguchi H., Irikura K., Iwata T., 2000. Fault geometry at the rupture termination of the 1995 Hyogo-ken Nanbu earthquake, Bull. seism. Soc. Am. , 90( 1), 117– 133. https://doi.org/10.1785/0119990027 Google Scholar CrossRef Search ADS   Sekiguchi H., Irikura K., Iwata T., 2002. Source inversion for estimating the continuous slip distribution on a fault-introduction of Green's functions convolved with a correction function to give moving dislocation effects in subfaults, Geophys. J. Int. , 150( 2), 377– 391. https://doi.org/10.1046/j.1365-246X.2002.01669.x Google Scholar CrossRef Search ADS   Shi X., Wang Y., Liu-Zeng J., Weldon R., Wei S., Wang T., Sieh K., 2017. How complex is the 2016 M w 7.8 Kaikoura earthquake, South Island, New Zealand?, Sci. Bull. , 62( 5), 309– 311. https://doi.org/10.1016/j.scib.2017.01.033 Google Scholar CrossRef Search ADS   Van Houtte C., Bannister S., Holden C., Bourguignon S., McVerry G., 2017. The New Zealand strong motion database, Bull. N. Z. Soc. Earthq. Eng. , 50( 1), 1– 20. Wessel P., Smith W.H.F., 1995. New version of the generic mapping tools, EOS, Trans. Am. geophys. Un. , 76( 33), 329– 329. https://doi.org/10.1029/95EO00198 Google Scholar CrossRef Search ADS   Yoshida S., 1989. Waveform inversion using ABIC for the rupture process of the 1983 Hindu Kush earthquake, Phys. Earth planet. Inter. , 56( 3–4), 389– 405. https://doi.org/10.1016/0031-9201(89)90172-6 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. supplementary data.zip Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations