Source complexity of the 2016 MW7.8 Kaikoura (New Zealand) earthquake revealed from teleseismic and InSAR data

On November 13, 2016, an MW7.8 earthquake struck Kaikoura in South Island of New Zealand. By means of back‐projection of array recordings, ASTFs‐analysis of global seismic recordings, and joint inversion of global seismic data and co‐seismic InSAR data, we investigated complexity of the earthquake source. The result shows that the 2016 MW7.8 Kaikoura earthquake ruptured about 100 s unilaterally from south to northeast (~N28°–33°E), producing a rupture area about 160 km long and about 50 km wide and releasing scalar moment 1.01×1021 Nm. In particular, the rupture area consisted of two slip asperities, with one close to the initial rupture point having a maximal slip value ~6.9 m while the other far away in the northeast having a maximal slip value ~9.3 m. The first asperity slipped for about 65 s and the second one started 40 s after the first one had initiated. The two slipped simultaneously for about 25 s. Furthermore, the first had a nearly thrust slip while the second had both thrust and strike slip. It is interesting that the rupture velocity was not constant, and the whole process may be divided into 5 stages in which the velocities were estimated to be 1.4 km/s, 0 km/s, 2.1 km/s, 0 km/s and 1.1 km/s, respectively. The high‐frequency sources distributed nearly along the lower edge of the rupture area, the high‐frequency radiating mainly occurred at launching of the asperities, and it seemed that no high‐frequency energy was radiated when the rupturing was going to stop.


Introduction
It was reported by the U.S. Geological Survey (USGS)/National Earthquake Information Center (NEIC) that an M W 7.8 earthquake struck Kaikoura, New Zealand, on November 13, 2016 (the 2016 M W 7.8 Kaikoura earthquake), with epicenter location of 42.757°S, 173.077°E, and focal depth of 23 km (Figure 1). The aftershock seismicity and local tectonics indicate that the earthquake had ruptured a SW-NE fault. Considering the USGS W-phase moment tensor solution, the earthquake occurred on the fault plane of strike 219°, dip 38°, and rake 128°, and had centroid depth of 15.5 km. Considering the gCMT solution, the fault plane was of strike 226°, dip 33°, and rake 141°, and centroid depth was 18.8 km. Nevertheless, the solutions identically indicated that it was a right-lateral strike and thrust event and had a prominent nondouble-couple component (Duputel and Rivera, 2017). Hollinsworth et al., 2017;Kaiser et al., 2017;Zhang H et al., 2017). Different studies based on different datasets using different techniques reveal the different aspects of complexity. Kaiser et al. (2017) pointed out that the earthquake involved at least 12 separate faults extending over 150 km from the epicenter in northern Canterbury to near the Cook Strait. Hamling et al. (2017) concluded from field observations, InSAR, GPS, and seismic recordings that the rupture propagated northwards for more than 170 km along both mapped and unmapped faults before continuing offshore at its northeastern extent, and that the event involved at least 12 major faults, including possible slip along the southern Hikurangi subduction interface, extensive uplift along much of the coastline, and widespread anelastic deformation. By means of optical satellite imagery and seismic data, Hollingsworth et al. (2017) identified surface ruptures of two faults, Kekerengu Fault and Jordan Thrust, as well as another previously unrecognized fault, the Papatea Fault. Zhang H et al. (2017) obtained the tempospatial distribution of slip on two separate finite faults using teleseismic waveforms and the migration imagery of high-frequency sources jointly using two regional arrays. Duputel and Rivera (2017) simulated the long-period seismic trains using 4 doublecouple subevents, including one thrust-slip event. Obviously, the 2016 M W 7.8 Kaikoura earthquake has a very complex source; different techniques and/or datasets have revealed only some aspects of the complex source.
An additional consideration is that rupture characteristics of large earthquakes have been found to be strongly dependent on frequencies (e.g., Lay et al., 2012;Yao H J et al., 2013). For the 2016 M W 7.8 Kaikoura earthquake, several rupture features have been identified (Hamling et al., 2017;Hollingsworth et al., 2017;Kaiser et al., 2017;Zhang H et al., 2017). In this study, in order to present  (Langridge et al., 2016). The red beach ball denotes the W-phase moment tensor solution of the mainshock from USGS; the black one denotes the gCMT solution. The gray beach balls indicate the gCMT solutions of the aftershocks (M W ≥5.5). Inset: regional tectonic and relative movement of plates. The Pacific plate has been going under the Australian plate at a rate of 40-50 mm/a (DeMets et al., 2010). The purple and orange beach balls denote the gCMT solutions of earthquakes (M W ≥6.0) since 1979, where the purple ones are for the M W 6-7 earthquakes and the orange ones are for the M W 7-8 earthquakes.
another aspect of the complex source we configured several seismic arrays with stations across the world and selected the optimal one to image the sources radiating high-frequency signals (high-frequency sources), and selected 60 stations of broadband seismic recordings from the IRIS data center and three tracks of In-SAR data to invert for sources radiating low-frequency signals (low-frequency sources). These results will be not only helpful to the comprehensive understanding of the 2016 M W 7.8 Kaikoura earthquake itself, but also very significant to understanding rupture properties of other earthquakes.

Migration of the High-Frequency Sources
The back projection technique of seismic array recordings has been very popular in recent years in tracking high-frequency sources of large earthquakes (e.g. Du H L et al., 2009;Ishii et al., 2005;Kennett et al., 2014;Krüger & Ohrnberger, 2005;Yagi et al., 2012;Yao H J et al., 2012). This technique can be employed to estimate the velocity, duration, and dimension of rupture (e.g., Lay et al., 2009;Uchide et al., 2013;. In this study, we employed it firstly to image the tempo-spatial variation of the high-frequency sources, and then to estimate the rupture direction as well as the rupture velocity. Both configuring array (s) and choosing its (their) data are the keys to imaging high-frequency sources. A good array should have good array-response, and good data should have high signal-tonoise ratios. In order to have satisfying arrays with high-quality data, we configured arrays by choosing stations with epicenter distances of 30°-90° across the world at which the propagation paths of waves are simple. As a result, we configured four arrays as candidates, named Southeast Asia (SEA) array, Japan (JP) array, Australia (AU) array, and South America (SA) array, as shown in Figure 2 and Figure S1 in the supplementary material. After evaluating their array responses, we chose one of them, the SEA-array, because of its better array-response and recording quality. The array consists of 20 broadband stations with epicentral distances of 60°-90° and azimuth range of 280°-315°. Note: this array is nearly perpendicular to the rupture direction of the 2016 M W 7.8 Kaikoura earthquake, which allows us to objectively recognize the rupture details (Lay and Wallace, 1995). Having tried various frequencybands, we adopted 0.1-1.0 Hz as the final one.
It is well known that real arrival times are usually not the same as the theoretical ones due to differences between the real earth and the model earth; thus it is often necessary to calibrate the real arrival times to match the theoretical ones before back-projection.
Here we adopted the adaptive stacking technique (Rawlinson and Kennett, 2004) in calibrating the time differences. All the recordings (vertical components) shown in Figure 3a have been not only filtered with frequency-band of 0.1-1.0 Hz but also calibrated in arrival times.
Teleseismic recordings are known to be poor in depth resolution, so we had to give up information on depth variation of sources. According to the location information of the earthquake, mentioned above, we fixed all the source candidates at depth of 23 km (several kilometers of variation in depth will not have substantial impact on results. See Figure S2 in supplementary materials). Meanwhile, we chose a rectangular area of 170°E-176°E and 40°S-45°S to cover the possible rupture area, and gridded it into sub-areas of 0.02°×0.02°. We set a time-window of 15 s and kept a shift-step of 1 s after numerous tries based on our experiences (Du H L, 2007;Du H L et al., 2009).
Looking at the seismic recordings from the array (Figure 3a), the duration of high-frequency signals is at least 100 s. In order to guarantee a complete source process, we analyzed 150 s-long recordings following the first motions. As a result, we obtained 136 locations with the 15 s-wide time window and 1 s-long shift step.
In order finally to determine effective source locations and duration, we adopted the γ-criterion (i.e., γ=αβκ) (Du H L, 2007;Du H L et al., 2009), where α is the value of the correlation coefficient of the beamed wavelet with the reference one (from the reference station), β is the value of the amplitude ratio of the beamed wavelet to the reference one, and κ is the value of the amplitude-nor- malized beam-energy. The criterion demands that an effective or believable source should emit same-shaped signals to all the stations, the signal recorded at each station should have the same amplitude, and the signals recorded at all the stations should be strong enough. As shown in Figure 4, the γ value is getting close to zero around 90 s, and has reached zero at 100 s, suggesting the effective duration of the whole process was less than 100s. The γfunction clearly shows that there are 3 peaks during the process; the peak around 55 s is the highest while the one around 80 s is the lowest. Figure 5 directly presents the back-projected images for each of the time-windows, on which the migration of the highfrequency sources in the first 80 s can be clearly observed.
Peaking up the maximal energy-value for each of the time windows and determining their locations and relative strengths automatically resulted in an image as presented in Figure 6a, which shows the tempo-spatial evolution of the high-frequency sources. The direction of the evolution was estimated to be ~N33°E.
In order to determine the migration velocity of the high-frequency sources, we calculated distances of all the sources from the first one and put all of them on a time-distance diagram. As Figure 6b shows, the distance is generally increasing with running time, but the speed looks different in different time-windows. We divided the whole process into 5 stages labeled AB, BC, CD, DE and EF, respectively, and estimated the velocities for each.
As presented in Figure 6b, they were ~1.4 km/s, almost 0 km/s, 2.1 km/s, almost 0 km/s, and ~1.1 km/s, respectively. For the whole process, the average velocity was ~1.4 km/s, and the sources seemed to cover a distance of over 100 km.
Compared to the above result, by using local and regional recordings equal to and greater than 0.25 Hz Kaiser et al. (2017) obtained an amazing result showing that three belts of the sources seemed to be riding on three existed faults which were almost parallel to each other, and stated that the process lasted for at least 120 s, which were cut into 4 phases, 0-20 s, 20-40 s, 40-70 s and 70 s-end, respectively. The different phases had different apparent rupture velocities, for example, < 2 km/s in the first phase, 2.5 km/s in the second one, and ~2.5 km/s in the last one. Different from our work, Kaiser et al. (2017) employed local and regional recordings with higher frequencies, and obtained an image showing finer details. Zhang H et al. (2017), using dual high-frequency (0.4-2.5 Hz and 0.2-1 Hz) recordings from two regional arrays, Southeast Asia/ Oceania network (AU) and South America network (SA), obtained another result showing that the event involved at least two southwest-northeast striking faults, with an average rupture speed of 1.4-1.6 km/s and total duration of ~100 s. Actually, we had checked the arrays for their responses. As shown in Figure S1, the responses of the AU, JP, and SA are definitely worse than that of the SEA. But the apparent rupture velocity and duration reported by Zhang H et al. (2017) are similar to ours in spite of the migration images of the high-frequency sources.
Using the back projection method, Hollingsworth et al. (2017) concluded that the event had two main bursts of radiation, one around 20 s and near the epicenter, and the other around 70 s and close to the observed surface rupture. The rupture velocity estimated for the first burst was 2.0-2.5 km/s, and if the two bursts were considered together, the average velocity was ~1.5 km/s, which is similar to ours. Different from our work, they used the recordings of 0.5-2 Hz from the Australian seismic network, which was abandoned in our work due to its unsatisfying response, as shown in Figure S1. In all the work mentioned above, the different datasets were processed with the same method. All of the results revealed at least two common points, though they were not identical. The first is that the 2016 M W 7.8 Kaikoura earthquake was a unilateral rupture from southwest to northeast, and the second is that the event had an unexpectedly complex source.

Overall Rupture Azimuth
Analysis of apparent source time functions (ASTFs) is one of the important approaches to understanding general feature of earthquake rupture (e.g., Ammon et al., 2006aAmmon et al., , 2006bLay et al., 2009).
Like the back projection technique mentioned above, it can be used at least to obtain information on rupture directivity.
Several techniques have been used in retrieval of the ASTFs (e.g., Bertero et al., 1997;Chen Y T et al., 1991;Dreger, 1994;Helmberger and Wiggins, 1971;Kraeva, 2004;Piana & Bertero, 1997;Xu L S et al., 2002). Here we chose the PLD (Projected Landweber Deconvolution) technique for its advantages: causal relationships could be taken into account; non-negative constraints could be imposed; and it runs very efficiently Lanza et al., 1999;Piana and Bertero, 1997;Vallée, 2004;Zhang Y et al., 2009). This technique has been successfully applied in our previous work The distances change with time going of the sources from the initial point; the whole process is divided into 5 periods, AB, BC, CD, DE and EF, respectively. Different periods have different velocities, with a maximal value of ~2.1 km/s, but the average velocity over the whole process is ~1.4±0.1 km/s. (Xu L S et al., 2014;Zhang X and Xu L S, 2015;Zhang X et al., 2016).
Sixty stations with epicentral distances of 30°-90° and azimuth coverage of 360° were chosen from the IRIS data center ( Figure 7); they are well-distributed in space and their recordings have good signal-to-noise ratios. In computing Green's functions, we used the ak135 earth model (Kennett et al., 1995), the source parameters issued by USGS (epicenter, 42.757°S and 173.077°E; depth, 23.0 km; strike, 219°; dip, 38°; rake, 128°), and the computer codes developed by Wang R J (1999). The frequency band of both the observation data and the Green's functions was constrained to 0.005-0.05 Hz.
The ASTFs retrieved from the P waveform data are shown in Figure 8a. They suggest that the effective duration should be less than 100 s, with most of the energy being released from 50 s to 80 s. Looking at the time points of the maximal values on all of the AST-Fs, we are able to conclude that the rupture should have strong directivity. We noticed that the ASTFs retrieved from the Rayleigh waveforms (Duputel and Rivera, 2017) had very similar features.
In order to exhibit the azimuthal characteristics of the ASTFs more clearly, we projected the ASTFs onto a polar-coordinate system in which radial and tangential directions represent time and azimuth, respectively. As Figure 8b shows, a highlighted ring is chan-ging in shape with azimuth. The ring should have been like the dashed red circle if no directivity had existed. Obviously, the rupture azimuth should be NNE, as the white arrow shows.
In order to quantify the azimuth for the optimal rupture direction, we extracted the time data of the maximal values on the ASTFs and simulated them using the relationship , is azimuth of station, is azimuth of rupture, and is takeoff angle of the ray (e.g., Lay and Wallace, 1995;Park and Ishii, 2015). As Figure 8c shows, the optimal value is 28° when the reaches the minimum. In other words, the optimal rupture direction is at azimuth of 28°, which is very similar to the 33° estimated from the migration of the high-frequency sources. Figure 9 compares the observed waveforms with the synthetic ones obtained by convolving the retrieved ASTFs, as shown in Figure 8a, with the corresponding synthetic Green's functions. The average correlation coefficient reaches 0.91, which indicates that the ASTFs are reliable. However, at a few of the stations, such as FORT, COEN, CCD, PPTF, HOO and WRKA, the coefficients are below 0.80. This relatively poor fitting may be caused by the focal mechanism's being fixed.

Tempo-Spatial Distribution of the Slip
The joint inversion technique has been proven effective in determining tempo-spatial distribution of slip on fault plane (e.g., Avouac et al., 2015;Delouis et al., 2010;Simons et al., 2011;Wald et al., 1996;Yue H et al., 2013;Zhang Y et al., 2012). This technique allows multiple types of observation data to be jointly inverted to obtain the optimal result, and has been successfully applied to many studies of dynamic rupture process. In this study, we jointly inverted the direct P-waves of teleseismic recordings and the near-source Line-of-Sight (LOS) coseismic displacements, using the joint inversion technique developed by ourselves (Zhang X, 2016).
The teleseismic data were the same as those used in the ASTFs analysis. But here we used 160 s-long P-wave recordings, including 10 s before the P first arrival, with a frequency band of 0.005-0.1 Hz. Similarly, we used the ak135 earth model (Kennett et al., 1995), the same source and fault parameters issued by USGS as above, and the computer codes developed by Wang R J (1999) in computing the Green's functions.
The LOS coseismic displacement data came from three tracks of the coseismic interferograms from Sentinel-1 and ALOS2 (Table 1). Thanks to the excellent imaging ability of Sentinel-1, the first postseismic TOPS (Terrain Observation with Progressive Scans) SAR acquisition in track 52 covering the epicentral area was made just two days after the mainshock. With a minimum time interval of 6 days, the T52 coseismic interferometric pair of 20161103-20161115 (Table 1) maintains very good coherence. The decorrelation in the T52 pair resulted from large deformation gradients caused by the coseismic rupture along the surface rupture where the deformation gradients were captured in another two StripMap L-band interferograms, T102 and T103 (Table 1). The Lband SAR data have better spatial resolution (2 m by 3 m) than Sentinel-1 TOPS SAR data. The later has pixel spacing ~5 m by 20 m in range and azimuth, respectively. The near-field measurements in the L-band interferograms also suffered from the large deformation gradients in a traditional InSAR processing strategy. To fully recover the near-field coseismic displacements in the two L-band interferograms, a sub-pixel offset compensation strategy was applied during the InSAR processing, in which the residual offsets The illustration of how to obtain the optimal rupture azimuth using the peak times. The vertical axis is arrival times of the peaks and the horizontal one is azimuths of the stations. The gray shaded area shows the error range of the measurements and the red discontinuous curve is the optimal fit to the measured values (blue dots).
between the master and resampled slave SAR data (that had been resampled into the master geometry) were used to redo the co-registration for the slave acquisitions. Due to limited spatial resolution of the TOPS SAR data, the correction method did not work well for the T52 pair. The whole InSAR processing was implemented with an automated InSAR processing package, gInSAR, that is being developed at CCEMO for the Canadian government (Feng W P et al., 2017). In computing the Green's functions for the InSAR coseismic displacements, we used computer codes EDGRN/ED-CMP developed by Wang R J et al. (2003) and the ak135 earth model (Kennett et al., 1995).
The initial finite fault model was set up according to strike of 219°a nd dip of 38° based on the USGS focal mechanism solution, 300 km long in the strike direction and 70 km wide in the down-dip direction. The whole fault plane was divided into sub-planes (faults) of 10 km×10 km. All the sub-faults had the same strike and dip, but their rakes were allowed to change by at most ±45°a round 128°, the given rake.
The slip functions of sub-faults were allowed to have different shapes and to be automatically determined in the inversion process. In addition, in order to obtain physically reasonable results, smoothing constraints in both time and space (Yagi et al., 2004;Zhang Y, 2008) were imposed, as well as the constraint of minim-izing the scalar moment (Antolik and Dreger, 2003;Hartzell and Heaton, 1983;Zhang Y, 2008).
Although the slip-function shapes of sub-faults were variable, the slip duration of each sub-fault had to be limited. As we know, there is trade-off between slip duration and rupture velocity, and inversion results will be strongly dependent on the trade-off (e.g., Gusman et al., 2015;Lay et al., 2010). Therefore, we first adopted the maximal rupture velocity 2.1 km/s, which had been obtained by means of the back projection of high-frequency signals, as the maximum in the whole rupture process, and then we determined the optimal slip duration by trial and error. As Figure 10a shows, the misfit value obviously changes at 40 s, suggesting that 40 s should be the optimal choice.
Relative weights of seismic and geodetic data also strongly affect final results. In this study, we determined the relative weight of the InSAR data with respect to the seismic data by evaluating the fitness of both the seismic data and the InSAR data. Different weights resulted in different fitness, as Figure 10b shows, but the weight value of 2 made the two types of data fit equally well. Therefore, we chose 2 to be the optimal weight.
The inverted STF (Figure 11a) shows that the rupture lasted a little bit less than 100 s and released scalar moment of 1.01×10 21 Nm; the dominant moment was released in the first 80 s. The static slip  Figure 9. Comparison of the observed waveforms with the synthetic ones calculated using the retrieved ASTFs. On the left side of each subplot are station name, epicentral distance (°) and azimuth (°), respectively, and on the right side are correlation coefficient and phase name, respectively. The average correlation coefficient reaches 0.91 as shown at the bottom right.

Earth and Planetary Physics
doi: 10.26464/epp2018029 317 distribution, as shown in Figure 11b, suggests that the slip mainly concentrated on two asperities, with one close to the initial rupture point, relatively weaker, while the other far away and much stronger; the whole rupture was unilateral. Details of the rupture propagation can be seen clearly in Figure 11c, which displays the rupture times and the slip functions of all the sub-faults.
Considering the spatial features of the slip distribution, we divided the whole process into two sub-events A and B, as shown in Figure 11b. Sub-event A lasted 65 s, covering ~70 km along strike, producing a maximal slip of ~6.9 m, and releasing ~30% of the total moment. It is equivalent to an event of M W 7.6. Sub-event B started 40 s later and also lasted about 40 s, covering ~90 km along strike, producing a maximal slip of ~9.3 m, and releasing 70% of the moment. It is equivalent to an event of M W 7.8. It is worth stressing that event A was nearly thrust while event B had both strike and thrust components.
We also noticed that sub-event A had not ended before sub-event B started (Figure 11a). In order to illustrate the tempo-spatial complexity of the rupture process, we made snapshots of slip-rates as shown in Figure 12. In the first 10 s, the slip was very weak and had no observable directivity, which can be considered to be bilateral. Afterwards, the slipping area started moving obviously toward the left, and its slipping amplitude was gradually increasing till about 80 s later, which was absolutely unilateral.
We computed the synthetic waveforms using the inverted dynamic rupture model in order to compare with the observed ones. As shown in Figure 13, the correlation coefficient at each station is above 0.8, and the average reaches 0.93. Similarly, we computed the LOS displacements using the same model as above, and compared them with the observed ones, which are exhibited in Figure 14. The variance reductions (Kim and Dreger, 2008) for the LOS displacements of Track 52, Track 102 and Track 103 are 93%, 71% and 81%, respectively, which means that the main part of the observed information has been explained. Looking at the pictures shown in Figure 14, fitness along the coastal line is quite poor, even unacceptable. We think there are two possible causes. One is that the InSAR observation itself has been in question because of its special geographical location, and the other is that our simple plane model is incapable of explaining the InSAR observa-   13. Comparison of the observed waveforms with the synthetic ones calculated using the inverted model. On the left side of each subplot are station name, epicentral distance (°) and azimuth (°), respectively, and on the right side are correlation coefficient, phase name and channel, respectively. The average correlation coefficient reaches 0.93 as shown at the right bottom.

Earth and Planetary Physics
doi: 10.26464/epp2018029 319 tion along the coastal line, which previous studies have revealed to be largely dominated by local tectonics (Hamling et al., 2017;Cesca et al., 2017;Lo et al., 2018) .   aspects obtained by varieties of techniques from varieties of datasets carrying varieties of frequency components.

Discussion
The 2016 M W 7.8 Kaikoura earthquake has an extremely complex source process. Different studies have given different source models due to differences in techniques and/or datasets (Hamling et al., 2017;Hollingsworth et al., 2017;Kaiser et al., 2017;Zhang H et al., 2017). In our study, the back-projection method, the ASTF-analysis method, and the joint inversion method are adopted, and the high-frequency (0.1-1 Hz) P-wave recordings from the SEA array, the low-frequency (0.005-0.1 Hz) P-wave recordings from globally distributed stations, and the InSAR co-seismic displacement data were used, in order to obtain information on source complexity.
It is interesting that, as Figure 15 shows, spatially the high-frequency sources situate at deeper locations, forming the lower edge of the rupture area instead of overlapping on it. A similar phenomenon was reported in previous studies (e.g., Kiser and Ishii, 2011;Koper et al., 2011;Lay et al., 2012;Uchide et al., 2013;Yao H J et al., 2013). Especially for earthquakes between oceanic block and continental block, for example, most of the low-frequency signals are radiated by rupturings at depths of 15-35 km while most of the high-frequency ones are produced at depths of 35-55 km Yao H J et al., 2013).
It is more interesting that, as shown in Figure 16, temporally the high-frequency radiation occurs before the low-frequency radiation. Building source time functions for the asperities A and B, respectively ( Figure 11a and Figure 16c), and comparing them with the normalized-power time function shown in Figure 16a, we get a clear impression that the first peak of the high-frequency power appears around 13 s while the peak of the asperity A's slip function occurs around 40 s, and the second peak of the high-frequency power appears around 50 s while the peak of asperity B's slip function takes place around 70 s. On the normalized power function does the third peak appear, but no more peak follows the asperity B's peak on the source time function, it is likely related with healing of the asperity B as shown in Figure 16f.
The rupture of asperity B starts when the slip of asperity A is still continuing. As Figure 16e shows, the much stronger high-frequency energy is being radiated while asperity B is starting its rupture, and no high-frequency energy is being radiated while asperity A is powering off.
In addition, we notice that the high-frequency sources always appear in the "backyard" of the low-frequency ones, whose slip directions denote their "foreyard". Calculating average rakes for asperities A and B, and projecting onto the ground surface, as the red arrows show in Figure 15, we find that the two asperities have different slip directions, but the high-frequency sources always stay behind the low-frequency sources, implying that the high- Having identified the temporal and spatial nature of the high-frequency sources, we can easily understand why the high-fre-quency sources have different spatial distribution from the lowfrequency ones, and why the times at which the high-frequency signals are the strongest are different from those at which the low-frequency signals are the strongest, and we can also easily understand why no high-frequency sources appear on the Northern and Northwestern edge of asperity B.
The STF obtained by inverting the teleseismic data and InSAR data indicates that the 2016 M W 7.8 Kaikoura earthquake lasted about 80 s while the time function of the high-frequency power shows about 100 s. We prefer the 100 s to be the duration time of the event, because something was very likely happening around the high-frequency sources while asperity B was "powering off". In other words, the high-frequency radiation still belongs to the rupture process though the low-frequency radiating is getting weaker and weaker. Duputel and Rivera (2017)  The slip direction of asperity A is found to be different from that of asperity B, which has a greater strike-slip component. We think this can only be considered as a general feature of the whole event, for we adopted a simple planar model with fixed dip and strike. In fact, at least 12 mapped and unmapped faults with varieties of focal mechanisms were involved in this event (Hamling et al., 2017;Kaiser et al., 2017), and most of them were of strike-slip types. Duputel and Rivera (2017) employed three strike-slip events and one thrust-slip event in order to explain the main fea-ture of ground motion, and thought the first event around the epicenter should be of the strike-slip type and the thrust-slip event should be at the northern end; but we feel strongly that the purely thrust event should be under all the strike-slip events, that more strike-slip events should be on northern area (around asperity B) than on the southern area (around asperity A), and that a possible model should be like the one shown in Figure 17, which is similar to the model proposed by Lo et al. (2018). Of course, more comprehensive investigation and study will be required in order to clarify this.

Conclusions
The 2016 M W 7.8 Kaikoura earthquake has been found to exhibitthe most complex source process ever recorded, as revealed by the latest studies (Duputel and Rivera, 2017;Hamling et al., 2017;Hollingsworth et al., 2017;Kaiser et al., 2017;Zhang H et al., 2017). According to back-projection of the high-frequency (0.1-1 Hz) P-wave recordings from a deliberately chosen and configured array, quantitative analysis of the ASTFs from 60 globally distributed stations, and joint inversion of the low-frequency (0.005-0.1 Hz) P-waveform data from the 60 stations and co-seismic InSAR displacement data from 3 tracks of satellite interferograms, the event lasted ~100 s, produced a rupture area about 160 km long along strike direction and about 50 km wide in down-dip direction and released a scalar moment of 1.01×10 21 Nm. Spatially, the rupture area consisted of two asperities, one close to the initial rupture point having a maximal slip value ~6.9 m and the other far away in the north having a maximal slip valuẽ 9.3 m. Temporally, the first asperity slipped for about 65 s and the second initiated 40 s later and lasted for 40 s, with both of them slipping simultaneously for about 25 s. As to focal mechanism, the first asperity had a nearly thrust slip while the second had both thrust and strike slip components. Based on the analysis of the high-frequency signals, the ASTFs, and the joint inversion, the rupture propagated generally in a unilateral way at an average speed of 1.4 km/s, but in details, as revealed by the high-frequency analysis, the rupture velocities were changing, with 1.4 km/s, 0 km/s, 2.1 km/s, 0 km/s and 1.1 km/s corresponding to the 5 stages of the whole process, respectively. Most of the highfrequency energy was generated on the lower edge of the rupture area, and the stronger high-frequency energy was radiated by the launching of each asperity; thus the high-frequency energy released some time before the low-frequency energy, and no high-frequency energy seemed to be radiated when the asperities' rupturing was coming to its end.