Three-dimensional coseismic deformation of the 2016 MW7.8 Kaikuora, New Zealand

Jio Liu , Guohong Zhng ,,*, Jiqing Wng , Gungtong Sun , Yingfeng Zhng ,Ynzho Wng , Chunyn Qu , Xinjin Shn

a State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China

b Institute of Disaster Prevention, Yanjiao High Tech Zone, Sanhe City, Hebei Province 065201, China

c The First Monitoring and Application Center, China Earthquake Administration, Tianjin 300180, China

Keywords:The 2016 MW7.8 Kaikuora earthquake Three-dimensional InSAR Offset tracking

ABSTRACT

1. Introduction

Interferometry of SAR (InSAR) satellite images is commonly used to obtain coseismic deformation occurred during earthquakes.Compared with other geodetic means such as Global Positioning System (GPS) and Leveling, InSAR has advantages in high spatial sampling, large coverage of the epicentral area, and high vertical motion sensitivity. However, there are challenges with InSAR technology, such as easily being contaminated by water vapor or ionosphere disturbance. It is also difficult or even impossible for InSAR to monitor the surface displacements near faults during big earthquakes [1-5], because of severe loss of coherence caused by steep deformation gradients. Offset tracking of satellite images covered in the same scene with different sensing dates can be used to achieve the coseismic deformation based on SAR images and optical images. Offset tracking measurements usually have a precision of tens of centimeters[6],which is lower than that of InSAR measurements with several centimeters [3,7]. The very strong property of offset tracking method is the ability to acquire large deformation that occurred along near fault regions regardless of coherence,while high coherent situation is always the prerequisite for the application of InSAR.Surface measurements in regions with large deformation,especially for the near-source fault,are crucial to constraint seismic moment and fault slip distribution at depth.Thus, both InSAR and offset tracking techniques are important to obtain a full-coverage coseismic deformation, including near- and far-field of the epicentral area.In addition, it has long been recognized that InSAR coseismic measurements can detect the 1-dimensional (1D) deformation along Line-Of-Sight (LOS) direction effectively [8]. The offsets along Range and Azimuth directions in Range-Doppler plane or East-West and North-South motion components can be achieved by offset tracking method based on SAR images or optical images [6]. However, both InSAR and offset tracking methods could not be used to directly obtain the full 3-dimensional (3D) coseismic deformation. With multiple observing modes, including descending and ascending, together with applications of InSAR and offset tracking, it is promising to obtain the total 3D deformation of an earthquake. As noticed, 3D coseismic deformation is key to many issues related to the earthquake source, such as fault orientation identification and earthquake mechanism discrimination [9,10]. For instance, big earthquakes sometimes have spatially variable faulting mechanisms,as shown in the 2008 MW7.8 Wenchuan earthquake[4],and 3D coseismic deformation is essential to characterize such complexities of the seismogenic faults.

On 14th Nov 2016, the MW7.8 Kaikoura earthquake struck the northern part of south Island, New Zealand, causing two fatalities and dozens of injuries [11]. The Kaikoura earthquake was documented as one of the most complex events ever recorded on land[12-18]. For instance, the focal mechanism of the Kaikoura earthquake shows strong non-double-couple components[19,20],which signify dramatic variations of the earthquake rupture [12,13].Studies based on seismic measurements also claim that event's fault rupture was influenced by the heterogeneous crustal material properties of this region[21],and the nucleation during the rupture[22]. In addition, about 13 faults have been reported as coseismic ruptured and large displacements along surface ruptures found during the earthquake [23]. As for the tectonic background, the MW7.8 Kaikoura, New Zealand event is located within the active and complex Australia-Pacific plate boundary system (Fig. 1). To the northernmost part of the North Island of New Zealand, the Pacific plate subducts underneath the Australia plate almost perpendicular to the Hikurangi trench with a rate of 40 mm/a,and the subduction gradually becomes oblique to the northern part of the South Island, where the active Marlborough Fault System is generated to absorb the convergence between the two plates(Fig.1). To the southernmost of the South Island, on the contrary,the Australia plate subducts beneath the Pacific plate with an oblique rate of 34 mm/a along the Puysegur Trench. To the north,the subduction is mainly absorbed by the big dextral strike-slip fault, known as the Alpine Fault System, along the western coast of the southern part of the South Island, New Zealand. Due to the sharp contrasting subduction processes that happened between the North and the South, faults in the transitional region are very active, and the seismicity rate is very high, including the recent 2010 MW7.1 Darfield (Canterbury) earthquake and 2016 MW7.8 Kaikoura earthquake (Fig. 1). To sum it up, characterizing such a complicated event as the Kaikoura earthquake requires more comprehensive methods, including both InSAR and offset tracking techniques, to retrieve full-coverage and complete 3D surficial coseismic deformation. The previous studies focused on the calculation of the 3D coseismic deformation of 2016 MW7.8 Kaikoura earthquakes and the analysis of its mechanism [13,24,25].While the highlight of this study is accurate reconstruction of 3D deformation by add more dataset with high spatial resolution.Our study will capture a higher resolution than that of previous achievements (2 km-5 km), although no new conclusions about this event are drawn.

In this study,we firstly obtain the 1D coseismic deformation along LOS using InSAR technique. Then offsets along azimuth and slang range in Range-Doppler coordinates, and components along East-West and North-South directions in geographical coordinates are all acquired by conducting offset tracking on SAR images and Landsat-8 optical images,respectively.Afterward,based on the coseismic offset tracking measurements, complete 3D coseismic deformation of the 2016 Kaikoura event is obtained using a weighted least square fitting.In order to verify the reliability of the results,we then compared them with GPS measurements and the horizontal offsets by Landsat-8 optical images.Finally,we analyze the rupture process based on the complexity of the coseismic deformation revealed by the full 3D coseismic deformation and the tectonic background.

2. Data description

In order to fully capture the coseismic deformation of the 2016 Kaikoura earthquake,we collect multiple kinds of satellite images,including two kinds of SAR images and one kind of optical images.The details of the data are listed in Table 1. Two different satellite sensors obtain the SAR satellite images, one is the Advanced Land Observing Satellite-2 (ALOS-2) PALSAR in L band (with a wavelength of 23.6 cm),and the other is the Sentinel-1A images in the C band(with a 5.6 cm wavelength).Both two kinds of SAR images are observed with descending and ascending modes. All the pairs's perpendicular baselines are within 150 m, while the maximum time interval is less than 120 days. The relatively short spatial baseline and time interval ensure the relative high coherence for this event to get a readable differential interferogram.The Landsat-8 optical images are also used to map the 2D horizontal surface displacements associated with this event. Note that, unlike Landsat-8 satellite's orthoscopic viewing, SAR satellite acquires data in slant direction. Both InSAR and offset tracking of SAR data get observations either along LOS direction or two orthogonal directions in radar slant plane (Fig. 2). However, offset tracking of Landsat-8 data obtains two orthogonal directions in the ground surface plane:East-west and North-south.Another thing worthy of notice, the acquisition dates of all the satellite images span the occurrence of the Kaikoura event,but also include short periods of inter-seismic and post-seismic phases. Possible inter-seismic and post-seismic deformation over such short periods of the time spans is likely negligible compared to the coseismic deformation,and we thus assume the measurements acquired by InSAR or offset tracking are records of only the coseismic deformation.

3. InSAR processing and 1D coseismic deformation along LOS

1D coseismic deformation along LOS of the 2016 Kaikoura earthquake is obtained using InSAR.We process both the Sentinel-1A and ALOS-2 data with GAMMA software.The traditional 2-pass DInSAR method is used to process the PALSAR-2 data,whose result presents excellent coherence in this earthquake. The processing procedure of Sentinel-1A TOPS mode data is quite different from that of the traditional strip map mode data, so it needs a more precise registration process. The new strategies for SLC coregistration have been incorporated into the new version of GAMMA software. The detailed data processing procedure for S1A TOPS mode data can be found in[26].After the SLC co-registration,the processing of the TOPS data can be followed by the traditional InSAR procedures. The 90 m SRTM is used to remove the phase related to ground topography. The unwrapped InSAR phase is detrended by a linear ramp estimated from a far-field area without seismic deformation to account for possible long-wavelength orbital errors and atmospheric errors. However, both ALOS-2 and Sentinel-1A state vectors are precise.The differential interferogram was unwrapped using the minimum cost flow (MCF) method. In addition, we masked out low coherence areas (coherence under 0.35), which were not be unwrapped. The resulting unwrapped phase was checked manually for unwrapping errors,especially for the regions covering the earthquake deforming zone. Afterward,the unwrapped phases were converted to the ground displacement,which was finally geocoded from radar slant coordinate to geographic coordinate to form (Fig. 3, Fig. 4).

Fig.1. Tectonic setting of the 2016 Kaikoura earthquake.

Table 1 The SAR and optical data information in this study.

Fig. 2. Geometric relationship between LOS and radar slant plane. ɑ is the radar look angle at the satellite, and θ is the incidence angle.

Figs. 3 and 4 give the rewrapped fringes and unwrapped displacement of the Kaikoura coseismic deformation obtained by InSAR, and one fringe denotes 11.8 cm deformation along LOS.Through the denseness of fringes,large gradient of deformation can be observed clearly,especially in regions from the eastern coast to the north part of the South Island, New Zealand. The fringes near the epicenter are enveloped and extend toward NE direction, and an enveloped center can also be found near the west side of the eastern coast of the South Island. There is an apparent incoherent zone bending inward to the island along the eastern coast,consistent with filed surveys [23]. According to the epicenter locations published by Geological and Nuclear Sciences Limited(New Zealand) (GNS) and GCMT, and the pattern of the fringes distribution, we can judge roughly that the earthquake source ruptures have involved multiple segments with an approximate striking of northeast, and may rupture at one-sided in multi-points. As revealed by the unwrapped measurements (Fig. 4), InSAR measurements show a peak relative motion of about 4 m both away from the satellite on the south-west side and towards the satellite on the north-east side (Fig. 4), which, if resolved into a purely vertical motion assuming an average incidence angle of 33°, is equivalent to an uplift of about 7.3 m. This is in good consistency with the maximum surface scarp of 9.0 m on the Kekerengu fault found by Morris et al. [27]. The difference between these two results may come from averaging the InSAR measurements within each pixel and the lack of near-fault observations of InSAR measurements.We find that ascending or descending tracks of different satellites show similar and comparable coseismic deformation(Fig.4A and C;Fig.4B and D)since they have only a slight difference between their observing incidence angles (Fig. 2; Table 1). Additionally, the far-field of the ascending ALOS-2 data seems contaminated by some ionospheric disturbance,masked to extrude further inversion investigations (Fig. 4A). An important feature worthy of notice is that,as shown in Figs.3 and 4,we obtain almost no InSAR measurements at near-surface rupture regions,which are commonly found in other InSAR earthquake studies [4,5,28].However, near-fault observations are crucial for further investigations, for instance, identifying surface ruptures. Hence, a technology that is not affected by decorrelation is essential for such a destructive earthquake.

4. Offset tracking and 2D coseismic deformation within radar slant plane

To complement InSAR measurements that are decorrelated at near-fault ruptures regions, offset tracking of SAR images and optical images has been implemented, which would help obtain 2D coseismic deformation along range and azimuth directions in radar slant plane, and along East-West and North-South in geographical horizontal plane respectively (Fig. 2). As also found by other studies, offset tracking is especially useful in providing near-fault measurements caused by destructive earthquakes [29-31]. In order to get redundant observations and mutual validations,both the S1 TOPS mode SAR images and ALOS2 SAR images are used to calculate the image offsets. During the process of interferometric pairs,the two SAR images acquired before and after the mainshock should be co-registered accurately at first, and then divided into many sub-images (patches). We follow a standard pixel tracking procedure to calculate the offsets between the corresponding patches using almost square search patches. To maintain similar pixel spacing for the two datasets, the offsets are estimated for every 20 range and 4 azimuth pixels in the S1 data,and for every 7 range and 35 azimuth pixels in the ALOS2 data.To reduce the noise further,an 8×8 window median filter is used,and finally,the range and azimuth offset measurements in radar slant plane are generated (Fig. 5).

Fig. 3. The coseismic interferometric fringes of the 2016 Kaikoura earthquake. The green and black stars indicate the epicenters published by GNS and GCMT, respectively.

Fig.4. The coseismic deformation fields associated with the 2016 Kaikoura earthquake.The mapped major faults are denoted by black thin lines. WF,Wairau fault;AWF, Awatere fault; CF, Clarence fault; KF, Kekerengu fault; JT, Jordan Thrust; HFZ, Humps fault zone; HF, Hope fault; HDF, Hundalee fault.

Fig.5. The range and azimuth offset measurements in radar slant plane by ALOS-2 and Sentinel-1.The green thick lines show the location of surface ruptures mapped referenced from [13,24,25] and inferred boundaries. PF, Papatea fault; CCF, Conway-Charwell fault zone.

As expected, the offsets based on SAR images clearly show continuous observations,especially at near surface rupture regions(Fig. 5). There is a good consistency between the result of ALOS-2 SAR images and that of S1 TOPS SAR images, while the former with high SNR (Signal Noise Ratio) is more reliable than the latter(Fig. 5A-D). Note that the range offset has the same vector direction with the LOS (Line of sight) measurement acquired by InSAR for a specific SAR satellite,and different names are just used to distinguish measurements obtained by different techniques.We find remarkable consistencies between range offset results and LOS measurements (e.g., Figs. 4B and 5B; Figs. 4D and 5D; Figs. 4C and 5F). In addition, we observe some boundary lines from all the 6 offset components, which may delineate the earthquake fault surface ruptures. Most of the surface ruptures are detected by the offset tracking results and are comparable with the results found by field investigations and other previous studies [13,24,25].

In addition to the SAR dataset,we use the COSI-Corr package to calculate the 2D horizontal surface displacement with the Landsat 8 images. We set the initial and final search window size as 32 pixels×32 pixels,step as 4 pixels×4 pixels,and the corresponding ground resolution as 60 m×60 m.To mitigate the effects of highfrequency noises on the accuracy of correlation results, we set the frequency mask threshold as 0.9 and robust iteration as 2. The 2D surface displacements along East-west and North-south are shown in Fig.6A and B.The footprint of Landsat-8 images covers only the northern part of the epicentral area,and the SNR is even worse than the S1A data,thus,we abandon the optical offset results in the later calculation of 3D deformation filed. Nonetheless, the 2D surface displacements by Landsat-8 offset tracking provide an independent data source to evaluate the resolved 3D coseismic deformation.

As mentioned above, coseismic measurements obtained by the offset tracking method always have a lower precision of 15 cm or worse than that of InSAR measurements (2-3 cm);such a low precision would always bias the inversion of fault slip distribution too much.For instance,a previous study shows that 5 cm error of surface deformation may lead to about 1 m uncertainty of fault slip at depth[32].Sowe donot use theoffset results asadirectinput tothefaultslip inversion in this study,instead,we use them as inputs of resolving 3D deformation and constraints to simplify the fault ruptures[33].

5. Full 3D coseismic deformation in the ground surface by inversion

Neither InSAR nor offset tracking measurements cannot characterize the real true 3D surface deformation,and cause difficulties in recognizing fault movement characteristics, and cannot be compared with GPS or field investigation measurements directly.In other words, for the same earthquake event, the LOS-oriented deformation fields may have opposite characteristics due to different observation modes,which may lead to misunderstanding for readers. For example, the same areas in Fig. 4C and D always exhibit opposite deforming characteristics shown as different colors.As for the 2016 Kaikoura earthquake,81 GPS stations detect the 3D surface coseismic measurements [13], but sparsely distributed within the South Island of New Zealand, has greatly limited the capability of characterizing the coseismic deformation. To achieve the full 3D coseismic deformation field with the high resolution,we thus use 6 components of offset measurements consist of azimuth offsets and range offsets from ALOS-2 descending track and Sentinel-1 descending and ascending tracks as the inputs measurements(Fig.5). The relation between offset measurements and the real 3D surface deformation can be expressed as follow:

where O is the offset observation, a 6 × 1 vector, representing 6 components of range offset or azimuth offset obtained by ALOS-2 or Sentinel-1A SAR images (Fig. 5), δ is the data uncertainty, a 6 × 1 vector,S is derived from the relation between offset measurements and surface deformation,which are introduced in detail in previous studies[9],U is a 3×1 vector,U =(dE,dN,dU)T,dU、dEand dNare the vertical, east-west and north-south displacement component,respectively. Since we have six components of observations but only three unknowns for each pixel, the equation is overdetermined, and can be solved using the least square fitting method to derive the reliable 3D surface deformation components[9,10,34]. The relative weights between different datasets were determined by based on their uncertainties,the calculation method is the same as that in Ref. [25,35].

Fig. 6. The 2D surface displacements along East-west and North-south by Landsat-8 images.

Fig.7. The 3D surface displacements of the 2016 Kaikoura earthquake.Subplot A/B/C denotes the east-west,north-south,and vertical deformation components,subplot D is the 3D deformation components from GPS stations [13], and inversion results, blue and red arrows denote vertical displacement components, purple and yellow arrows represent the horizontal velocity. The mapped major faults are denoted by thin black lines. The thick green lines show the location of surface ruptures mapped referenced from [13,24,25] and inferred boundaries.

The 3D coseismic surface deformation of the 2016 Kaikoura earthquake is shown in Fig.7;a very complex feature of coseismic surface deformation field is presented as a whole. Despite the complexity of coseismic deformation, it is still clear that the apparent deformation zones are closely related to the faults distributed on the ground surface spatially,which indicates that at least more than ten faults were ruptured during the earthquake,including some faults that were previously thought to have unclear characteristics of tectonic activity [12,36,37]. These faults across two seismotectonic zones with significant differences in activity modes and intensity and generally trended northeast and distributed in a range of about 160 km in length and 30 km in width.The intense horizontal faulting occurs in the northern section of the eastern coast of the South Island, mainly concentrated near the Kekerengu fault,the Papatea fault,the Jordan Thrust,and the Hope fault, while the faulting in the southern part of the rupture is relatively weak. The maximum horizontal deformation of about 10 m occurs along the Kekerengu fault with the vertical deformation up to 2 m and is consistent with the field investigation measurements [23] and the geodetic results obtained by Wang et al.[38]. Unlike the horizontal field, the vertical deformation exhibits concentrated and discontinuous features,while three uplift centers appear in the area defined by the Hope fault,the Humps fault zone and the Conway-Charwell fault zone, the intersection of Papatea fault and Jordan Thrust, and the region surrounded by the Hope fault,the Hundalee fault and the Stone Jug fault,and the maximum uplifts of about 3 m are observed along the Conway-Charwell fault zone and the Papatea fault. These can tell us that rupture in these sections experienced a more significant shortening. In addition,there is a large area with minor magnitude subsidence near the Awatere fault,and the emergence of such a phenomenon will be an issue to be discussed.

6. Discussions

In order to verify the solution result, we choose the offset measurements by Landsat-8 images and the 3D deformation components by GPS as independent observations quoted from Hamling et al., 2017 [13]. First, the east-west and north-south components obtained separately by inversion and Landsat-8 images offset tracking have enormous consistency (Figs. 7A and 6A;Figs. 7B and 6B). Both show similar displacement characteristics along the Hundalee fault,Papatea fault,and Jordan Thrust.We have also extracted two profiles across the ruptured faults to demonstrate the confidence of the inversion results(Fig.8).All the profiles along East-West and North-South components show good consistency,although the standard deviations of P1-P1′,P2-P2’in East-West and North-South components are 1.0 m and 0.82 m, 0.91 m,and 0.71 m.Moreover,the GPS measurements are also comparable to the resolved deformation field by InSAR, which however still contains some noise signals,especially in the horizontal directions(shown by vectors in Fig. 7D). With redundant azimuth offset measurements,even the north-south movement(which is the least insensitive to the SAR satellite) has been successfully resolved. It can also be seen that the East-West component calculated by GPS is larger than that calculated by inversion and Landsat-8 images. We also calculate the standard deviations of differences between our results and GPS data. Because of the location difference between the GPS station and its corresponding pixel in inversion results,we apply the nearest neighbour interpolation to the inversion results,and then get the 3D components at the GPS sites. Given the incoherent areas in inversion results,not all the sites have valid values.The standard deviation of differences between our results and GPS in EW,NS and UP direction is 0.76 m,0.56 m,0.83 m,respectively.In view of the large-scale deformation caused by the event and the interpolation errors effects, the deviations are still reasonable.However, the 3D coseismic surface deformation of the 2016 Kaikoura earthquake we obtained has much better continuous coverage of the epicentral region and much higher spatial sampling than GPS. This allows us to take a new insight into the ruptured faults.

Fig.8. The profiles of surface east-west,north-south movement components from inversion result and Landsat-8.The locations of 2 profiles(P1-P1’,P2-P2’)are indicated in Fig.7.The blue and black dots represent components extracted from the inversion result and Landsat-8 offsets.

The epicenter released by GNS shows that the initial rupture point is near the Humps fault [13]. However, relocation of the mainshock is supported by the evidence that the earthquake began to rupture from the Humps fault [23]. The 3D deformation clearly shows that the rupture mainly occurs in the north area of the epicenter. The fringes and LOS measurements mentioned above also show the same displacement, indicating the earthquake initiated on the Humps fault. The complexity of surface rupture zones in mega-earthquakes often reflects the complex correspondence between surface active faults and deep seismogenic structures [13]. If the dip angle of the fault plane does not change, the surface rupture zone should be approximately a straight strip.The fringes show before both have a curved incoherent strip with wide ends and narrow middle, which means that the dip angle of the fault varies considerably along the strike.

If we observe the coseismic faulting from the epicenter location along the NE direction,we can find that the Humps fault zone and the Hope fault near the epicenter show mainly dextral strike-slip,the Conway-Charwell fault zone and the Hundalee fault are mainly thrusting,while Jordan Thrust and Papatea fault have both dextral strike-slip and thrusting, and Kekerengu fault changes to dextral strike-slip again.This may indicate that the earthquake may be a rupture process of strike-slip, compression, transpressional rupture, and strike-slip in space along the NE direction, showing the mainshock is a multi-segment faulting. This may be caused by the shearing action transforming into compression when the strike-slip at depth encounters the blocking of NS-trending faults in its propagation along the NE direction. The Kekerengu fault in the northeastern segment was solely responsible for tectonic deformation,while multi-faults were involved in structural deformation in other parts of the rupture. As for the subsidence area near the Awatere fault, it may be caused by the pull-apart action of strikeslip faults on both sides, or a barrier may exist underground,which hinders the strike-slip movement along the North-East direction.The earthquake also tells us to focus not only on the known highly active areas,but also on the areas of the previously thought inactive faults.

7. Conclusions

Using the multi-perspective InSAR and offset-tracking measurements,we geodetically mapped and characterized the complex spatial pattern of the coseismic surface deformation produced by the 2016 Kaikoura earthquake. The rupture of at least ten active faults, including both the previously mapped and tectonically inactive faults,should be responsible for the complex deformation.With the precise 3-D coseismic deformation field resolved by the least-squares algorithm, we find that the largest coseismic displacement occurred along the Kekerengu fault, with a vertical displacement of up to 2 m and a horizontal displacement of ～10 m.Our 3-D deformation field is also consistent with the coseismic GPS measurements. Based on these observations and study, we think the mainshock was multi-segment faulting with a rupture process of strike-slip,compression,transpressional rupture,and strike-slip in space along the NE direction.

Author statement

Jiao Liu and Guohong Zhang designed the study, analyzed the results and revised the manuscript; Jiaqing Wang performed data processing, analyzed the results and revised the manuscript; Jiao Liu performed original draft and revised the manuscript; Guangtong Sun,Yingfeng Zhang,YanZhao Wang,Chunyan Qu,and Xinjian Shan participated in data processing and inversion work. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest

The authors declare that there is no conflicts of interest.

Acknowledgements

This work was co-supported by the National Key Research and Development Program of China (Grant No. 2019YFC1509204), the National Nonprofit Fundamental Research Grant of China,Institute of Geology, China Earthquake Administration (Grant No.IGCEA2005 and No. IGCEA2014), and the National Science Foundation of China (Grant No. 41631073). Figures are generated by Generic Mapping Tools [39].