Triggered afterslip on the southern Hikurangi subduction interface following the 2016 Kaikōura earthquake from InSAR time series with atmospheric corrections

The 2016 Mw 7.8 Kaikōura earthquake represents an extremely complex event involving over ten major crustal faults, altering conventional understanding of multi-fault ruptures. Although evidence for coseismic slip on the Hikurangi subduction interface is controversial, we present afterslip on the subduction zone beneath Marlborough using 13 months of Interferometric Synthetic Aperture Radar (InSAR) and Global Positioning System (GPS) observations. The spatially and temporally correlated atmospheric errors in SAR interferograms are problematic, and hence a new InSAR time series approach, combining the Generic Atmospheric Correction Online Service (GACOS) with a spatial-temporal Atmospheric Phase Screen (APS) filter to facilitate the InSAR time series analysis, is developed. For interferograms with over 250 km spatial extent, we achieve a 0.77 cm displacement RMS difference against GPS, improving 61% from the conventional InSAR time series method (TS). Comparisons between the overlapping region of two independent tracks show an RMS difference of 1.1 cm for the TS-GACOS-APS combined method, improving 54% from the TS method and 27% from using TS with an APS filter only. The APS filter reduces the short wavelength residuals substantially, but fails to remove the long wavelength error even after the ramp removal, revealing that the GACOS correction has played a key role in mitigating long wavelength atmospheric effects. The resultant InSAR displacements, together with the GPS displacements, are used to recover the time-dependent afterslip distribution on the Hikurangi subduction interface, which provides insights for reviewing the co-seismic slip sources, the present status of the subduction plate boundary and future seismic hazards.


Introduction
The 2016 Mw 7.8 Kaikōura earthquake, which struck the northern South Island of New Zealand on 13 November 2016, ruptured over ten major faults with up to 10 m surface displacements, generating a regional tsunami which peaked at ~7 m (Bai et al., 2017), and triggered numerous landslides (Massey et al., 2018). Geodetic and seismologic datasets have been widely used to constrain the complex multi-fault geometry and co-seismic slip distribution. Hamling et al. (2017) determined the fault geometry and inverted for the slip distribution, combining a range of datasets including field data, GPS and InSAR observations. In their models, at least 12 crustal faults slipped, which challenged the traditional assumptions about the degree to which earthquake ruptures are controlled by fault segmentation. Combining InSAR and seismologic data, Kaneko et al. (2017) reported that the rupture speed was overall slow (≤2 km/s), and the Conway-Charwell fault link aided the propagation of ruptures across the stepover from the Humps fault zone to the Hope fault (Xu et al., 2018). Holden et al. (2017) proposed kinematic models based on local strong-motion and high-rate GPS data, suggesting that the rupture propagated from south to north with half of the moment release occurring at the far north, 60 s after the origin time. However, none of these models has ruled-out possible slip along the southern Hikurangi subduction interface, which is still debated, and difficult to determine based on the existing coseismic observations (e.g. Holden et al., 2017;Hollingsworth et al., 2017;Xu et al., 2018).
After the 2016 Kaikōura event, there were considerable slips detected on the Hikurangi interface, which was previously found experiencing a slip rate deficit of approximately 10 mm/year (Wallace et al., 2012). Firstly, several Slow Slip Events (SSEs) were detected immediately after the mainshock beneath the North Island using a dense GPS network , including on the Kapiti coast a deep SSE accumulating up to 30 cm slip and a shallow (< 15 km), moderate (> 10 cm) east coastal SSE on the Hikurangi subduction interface (Wallace et al., 2018. Secondly, substantial afterslips were observed on the Hikurangi interface beneath the Marlborough Fault System (MFS) of the South Island (Wallace et al., 2018). Wallace et al. (2018) used 134 continuous GPS and semi-continuous GPS stations together with four descending Sentinel-1 interferograms from  to March 2017 to invert for the afterslip on the subduction interface. The ascending track with a longer time span (more than one year) and a large spatial extent was discarded due to substantial atmospheric disturbances. Therefore, they employed many more GPS data points (~55,000) than InSAR (~2500) in the time-dependent modelling. Jiang et al. (2018) used five months of data from GPS stations located in the east of the MFS and recovered a different but less widespread slip source localised in the southern portion of the slip distribution obtained by Wallace et al. (2018). Consequently, to further investigate the Hikurangi subduction interface and invert for its afterslip distribution, it is necessary to employ atmospheric-corrected, spatially and temporally continuous surface displacement observations around the MFS from the Sentinel-1 ascending track. This determination of the afterslip origin on the subduction interface and its timevarying distribution may subsequently be used to help resolve the debatable co-seismic component of any southern Hikurangi subduction interface slip, the present activity of the subduction plate boundary, and future seismic hazards.
To spatially and temporally constrain the afterslip model, InSAR interferograms across the entire MFS should be used and combined with any available GPS-estimated station displacements. The ascending track of Sentinel-1 (see Fig. 1) provides the largest and densest spatial-temporal coverage from any SAR satellites over the MFS, so is preferred and used here. However, large interferograms are subject to atmospheric disturbances (e.g. Parker et al., 2015), which have long been recognised as one of the major InSAR measurement error sources and must be sufficiently mitigated, otherwise actual tectonic displacements can be masked (e.g. Walters et al., 2013). Retrieving velocity maps with millimetric accuracy requires a stack of interferograms spanning as a time series to allow for the reduction of atmospheric errors whose temporal behaviour has previously been assumed random and separable from the ground movement (e.g. Cao et al., 2018;Fattahi and Amelung, 2015;Ferretti et al., 2001;González and Fernández, 2011;Hooper et al., 2007;Lauknes et al., 2011). For example, Ferretti et al. (2001) approximated atmospheric effects after removing a linear deformation component over small areas where the linear assumption and the constant velocity model held. Hooper et al. (2007) high-pass filtered the phase in time to isolate atmospheric contributions from deformation. Lauknes et al. (2011) modelled atmospheric errors as an additive Gaussian random process with a zero mean and 2-10 mm standard deviation. In summary, three fundamental assumptions are made by most current InSAR time series methods: (i) deformation is correlated in time and space, (ii) atmospheric effects are correlated in space but not in time, and (iii) noise is uncorrelated in space and time (e.g. Hooper et al., 2012). Nevertheless, despite the success of the abovementioned methods under specific situations, the displacement estimation may be biased as the atmospheric delay can be temporally correlated over a large spatial extent and therefore invalidates the three assumptions. First, water vapour can be highly correlated with its nearly invariant underlying topography (e.g. Li et al., 2006;Jolivet et al., 2011;Bekaert et al., 2015;Yu et al., 2017), implying a consistent phase-elevation correlation between independent interferograms collected on different dates. Moreover, a long wavelength error, caused by both atmospheric delay and orbital errors, can have a complicated spatial-temporal pattern which will eventually manifest as a bias on the velocity field. Removing a planar ramp from the observed phase, which has been done in many time series analysis methods (e.g. Biggs et al., 2007;Shirzaei and Walter, 2011;Wright et al., 2003;Zhao et al., 2016), is insufficient to reduce this effect, as will be discussed in Section 6.1, as well as having the potential undesirable effect of removing actual deformation signals, e.g. long wavelength inter-seismic and post-seismic motion.
To overcome these challenges, in this paper we propose a time series (TS) analysis method that takes advantage of the Generic Atmospheric Correction Online Service (GACOS; Yu et al. (2018)) to reduce the long wavelength and elevation dependent atmospheric errors from each interferogram, resulting in a more randomised error distribution of the phase measurements and which therefore helps satisfy the fundamental assumptions of InSAR time series analysis. Then a spatial-temporal Atmospheric Phase Screen (APS) filter is applied to further reduce the short wavelength atmospheric residual errors. The final velocity map and displacement time series are then estimated using least squares (e.g. Li et al., 2009;Mora et al., 2002). Our method (hereafter called TS-GACOS-APS) is validated across the MFS firstly against GPS displacement time series for discrete points, and then by comparing the results of an overlapping region covered by two different satellite tracks, on the basis that they are independent observations and have nearly identical look vectors. The more accurate displacement time series are then used to model the afterslip, in particular to determine the location of the afterslip origin on the Hikurangi interface.

Regional tectonic setting and geodetic observations
The oceanic Pacific plate obliquely converges under the continental Australian plate at various rates of 39-49 mm/year and causes active tectonics throughout New Zealand. Great earthquakes can occur, such as the 1855 Mw 8.2 Wairarapa event (Darby and Beanland, 1992), the 1976 Mw 8.2 Kermadec Island event (Habermann and Wyss, 1984), the 2009 Mw 7.8 Dusky Sound event (Beavan et al., 2010) and the 2016 Mw 7.8 Kaikōura event . In the northern South Island, the transition from the Hikurangi subduction zone to the strikeslip dominated Alpine Fault takes place via the MFS, which is a set of four large dextral strike-slip faults. Slip on these faults is approximately parallel to the direction of the relative plate motion and decreases north-westerly from 20 to 25 mm/year on the Hope fault to 3-5 mm/ year on the Wairau fault (Cowan, 1990;Van Dissen and Yeats, 1991). On the eastern side of the MFS, the fault swings anticlockwise by about 30 degrees, becoming the Jordan Thrust, which has a more northerly strike and a larger reverse slip component, and the Kekerengu fault with dominating dextral strike-slip (Van Dissen and Yeats, 1991), both of which ruptured during the 2016 Mw 7.8 Kaikōura earthquake (Fig. 1). Despite the underlying subduction interface at depths of 25-30 km (Williams et al., 2013), the shallow crustal faults accommodate the majority (> 75%) of the relative plate motion within the northern South Island according to Quaternary and geodetic observations (e.g. Holt and Haines, 1995;Norris and Cooper, 2001;Wallace et al., 2007). South of the MFS, at the latitude of northern Canterbury, the oblique plate convergence rate reaches 40 mm/year (DeMets et al., 1990) and is largely accommodated by a number of slowly deforming faults and folds, including the Humps and Hundalee faults (Pettinga et al., 2001).
The 2016 Mw 7.8 Kaikōura earthquake ruptured the Humps fault zone and the Hundalee fault, then propagated to the Hope fault through stepovers and splayed to further north on the Jordan Thrust, the Kekerengu fault and the Needles segment. At least 12 major faults were involved in this multi-fault rupture process with variable orientations and slip mechanisms, extending southwest-northeast for about 150 km, as shown in Fig. 1, following Hamling et al. (2017). It is much more complex than previously studied multi-fault rupture scenarios (e.g., the 2012 Sumatra event ruptured only three orthogonal strike-slip fault branches (Meng et al., 2012)). The stepovers between the Humps and Hope faults almost double the previously assumed limiting distance for halting a fault rupture , and transferred from more reverse faulting in the south to predominantly strike-slip in the north. The aftershocks first centred at the offshore end of the Hope fault with a broad northeast-southwest trend, then stepped approximately north following the Jordan Thrust and Kekerengu faults and clustered C. Yu, et al. Remote Sensing of Environment 251 (2020) 112097 at the Needles fault segment (near the Cape Campbell and Lake Grassmere area). Most of the aftershocks occurred at depths shallower than 30 km, with a mixture of reverse and strike slips according to the GeoNet Project, New Zealand, as shown in Fig. 1. Two ascending tracks of Sentinel-1 data (T52 and T154 -see Fig. 1) collected over the MFS during the period from November 2016 to December 2017 were used to investigate the post-seismic motion following the large earthquake. There were 33 acquisitions for T52 and 35 acquisitions for T154, with all the perpendicular baselines being no more than ~250 m (Fig. 2). Note that, as shown in Fig. 1, these two tracks have one overlapping swath spanning the same period but were observed on different dates, so their results can be compared for validation purposes. The descending track T73 has only 3 months of data (eight acquisitions), inconsistent with the ascending tracks, and therefore was not used for modelling. Nevertheless, for validation we will still process its time series using our proposed method and simulate its cumulative displacements based on our best-fit afterslip model. All the Sentinel-1 interferograms were generated using the GAMMA software (http://www.gamma-rs.ch) with the topographic phase contributions removed using a 3-arcsec (~90 m) Shuttle Radar Topography Mission digital elevation model (Farr et al., 2007). We also utilized the daily GPS time series from GeoNet (https://www.geonet.org.nz/), processed using the GAMIT/GLOBK software (the full processing strategy can be found at https://www.geonet.org.nz/data/supplementary/gnss_time_ series_notes). These GPS time series were detrended to remove secular inter-seismic deformation based on the inter-seismic velocity estimated from Beavan et al. (2016). The GPS displacement time series were used to validate the InSAR displacement time series (Section 4) and combined with InSAR to model the afterslip (Section 5).

InSAR time series and atmospheric correction: Methodology
Original interferograms may experience a mixture of topography correlated and turbulent atmospheric errors, manifesting as either short or long wavelength signals and degrading spatial-temporal filters when extracting deformation signals in InSAR time series analysis. Conventional InSAR time series analysis methods assume atmospheric errors are temporally random (e.g. Ferretti et al., 2001;Hooper et al., 2007;Lauknes et al., 2011), with only the spatial correlations considered. However, the atmospheric error temporal correlation should not be neglected, given that the tropospheric moisture content varies seasonally (e.g. fog is more prevalent at certain times of the year in coastal areas (Hooper et al., 2007)) and is correlated with topography, leading to an atmospheric error that can completely mask geophysical signals and introduce biases and unpredictable errors to the velocity estimates. Following the success of the reduction of long wavelength and topography related atmospheric errors on individual interferograms by GACOS (Yu et al., 2018), we now extend its usage to InSAR time series analysis. Our hypothesis is that the successful individual interferogram atmospheric correction will lead to a reduction of the atmospheric spatial-temporal correlation within a stack of interferograms. The residual short wavelength atmospheric errors are then filtered out by an APS filter and the displacement time series estimated by the TS algorithm.

GACOS atmospheric correction
GACOS incorporates model level surface pressure, temperature, and specific humidity from the High Resolution European Centre for Medium Resolution Weather Forecasts (HRES-ECMWF) numerical weather model (https://www.ecmwf.int/en/forecasts/datasets/set-i) to calculate atmospheric delays at each 0.125-degree grid point spatially (i.e. a spacing of approximately 9-12 km) and 6 h temporally (Yu et al., 2018). A simple linear temporal interpolation is applied as the acquisition time of the SAR image differs from the HRES-ECMWF data. The atmospheric delays are spatially densified to 90 m in accordance with the interferogram spacing using an iterative tropospheric decomposition spatial interpolator (Yu et al., 2017), which accounts for the elevation dependent and turbulent features of the troposphere. The interpolated zenith tropospheric delays are projected to the line-of-sight (LOS) direction by the Global Mapping Function (Boehm et al., 2006) and then subtracted from the unwrapped phases. Note that in fact the atmospheric error is a combined effect that comes from both the troposphere and the ionosphere. The ionospheric delay is generally moderate for short wavelength SAR satellites (e.g. the effect on C-band Sentinel-1 is ~5.5% of the effect on L-band satellites given a similar orbit geometry (Fattahi et al., 2017)) as it is inversely proportional to the frequency, at least over temperate zones, due to the shape of the Earth's magnetic field. We ignored the ionospheric effect in our study since 1) our study area (latitude 40-43S) is located in a temperate zone where the activity of the ionosphere is relatively quiet compared to polar and tropical regions; 2) the largest track we used only extends ~250 km (10 bursts per swath) compared to those in Liang et al. (2019) in which the ionospheric long wavelength effect is observed over much longer tracks (690 km long and 36 bursts per swath); and 3) the magnitude of the ionospheric delay on our interferograms (especially on track T154 which only covers 70 by 70 km) is much smaller compared to the tropospheric delay and has negligible effect on our subsequent geophysical modelling.
Substantial topography-related errors were observed on independent interferograms collected over different times of the 13month data set. Fig. 3 shows a sample of four from track T52 with notable such effects in which the measured displacement shows a positive correlation with elevation over the northeast part of the interferogram (denoted by the red box). If frequently observed throughout the time series, these spatial correlations could exert a temporal pattern on the phase measurement that could be wrongly interpreted as a ground deformation pattern, for example implying that high elevation areas are moving faster than low elevation areas. After GACOS atmospheric correction, it can be seen from Fig. 3 that the displacement-elevation dependence is reduced, with the fitted linear displacement-elevation functions approaching vertical lines and the average correlation coefficient reducing from 0.70 to 0.27, which implies their independence from elevation. As a result, the corrected interferograms ( Fig. 3b) tend to have much more random noise that better satisfies the basic time series analysis assumption of randomly distributed errors. They are therefore expected to be better handled by spatial-temporal filters. A sample of nine further interferograms from track T52 which exhibited substantial atmospheric effects is shown in Fig. 4. It can be seen that most of the original interferograms exhibit long wavelength signals along the northeast-southwest and northwestsoutheast directions with a maximum magnitude of over 10 cm, but which have been substantially reduced after applying GACOS atmospheric corrections.
In addition to the displacement-elevation correlation reductions showing the success of applying the GACOS atmospheric corrections, as an internal quality assessment we undertook cross-validation of the original 0.125-degree HRES-ECMWF atmospheric (zenith tropospheric) delays generated from GACOS. For each acquisition date of the two Sentinel-1 ascending tracks, we excluded one point from the whole HRES-ECMWF grid, determined its value from the remaining grid, and repeated for all grid points to obtain a cross-RMS difference between the interpolated and original values. The cross-RMS values, which Yu et al. (2018) show relate strongly with the atmospheric correction performance, are shown in Fig. 5 for each of the 33 track T52 and 35 track T154 image acquisition dates over the 13-month processing time span, and indicate a seasonal variation with a maximum in summer and minima in winter. If the cross-RMS for any image exceeded twice the standard deviation (computed using all 68 ascending track T52 and T154 cross-RMS values) then the GACOS corrections were deemed unsuitable for atmospheric correction. Based on this rejection threshold, 94.2% of the interferograms had GACOS corrections.

TS-GACOS-APS model
After applying the GACOS atmospheric corrections to those interferograms passing the cross-validation test, the conventional least squares based TS algorithm combined with a spatial-temporal APS filter was used to estimate time-dependent deformation (e.g. Li et al., 2009;Mora et al., 2002). For N interferograms from P identical dates, each map pixel complies with the following equation: where L is the original or atmospherically corrected phase observation on an interferogram with a master date of t m and a slave date of t s ; φ t is the cumulative displacement from the earthquake rupture time t 0 to t; B perp is the perpendicular baseline; D is the digital elevation model (DEM) error; r is the satellite-target distance (693 km for Sentinel-1); θ is the satellite incidence angle; ε accounts for the temporal decorrelation, orbital error, thermal noise effect and atmospheric error if not corrected. T is the coefficient matrix of the cumulative displacement and C is the coefficient matrix of the DEM error. If all the acquisitions are well connected, as was the case for our Sentinel-1 data, Eq. (1) can be well determined in a least squares sense.  For interferograms that have not had GACOS corrections applied because of failing the cross-validation test, their atmospheric errors are estimated by extending Eq. (1) to: where APS is the estimated atmospheric error for the uncorrected interferogram. Eq.
(2) is singular due to the correlation between the deformation and APS parameters. We therefore introduce a temporal deformation model as a constraint on the deformation parameter. For this postseismic study, a logarithmic deformation model may be used (Freed, 2007;Khoshmanesh et al., 2015): where φ is the phase change between time t k and the rupture time t 0 ; C and tau are the amplitude and the characteristic decay time, respectively, which are the parameters to be estimated. Substituting Eq. (3) into (2), we obtain: Eq. (4) can be determined with a well-connected acquisition network on a pixel-by-pixel basis. Notably, the temporal deformation model requires the APS parameter to be random, which can be largely satisfied in our situation as most of the interferograms have been corrected by GACOS and only a small portion of acquisitions require APS parameters (< 6%). Those limited number of failed GACOS corrections and those with a small improvement after correction will have a much lower chance of exerting systematic error patterns on the whole time series. To prevent unphysical oscillatory variations in the APS estimation, a spatial filter is performed on the APS parameters. Assuming that the atmospheric effect on each pixel within a given window W is identical, the final equation rearranged from Eq. (4) is obtained: where the DEM error for each pixel is introduced as an independent unknown parameter. Eq. (5) is an overdetermined system and can be easily solved in a least squares sense. Once the APS parameters are estimated, a whole network of interferograms corrected for atmospheric delays may be obtained.
A further refinement is done by solving Eq. (5) again using the atmospherically corrected interferograms (either by GACOS or the estimated APS). Instead of the total atmospheric delay, this time we assign a residual atmospheric delay as the APS parameter on each acquisition, provided that they are much smaller than the total values and more separable from the deformation signal. A step-by-step implementation of the proposed method is shown on the flowchart in Fig. 6.

InSAR time series and atmospheric correction: Validation
To evaluate the performance of the InSAR time series atmospheric corrections, four different methods were compared: (i) the conventional TS method without correcting or estimating the atmospheric delay (hereafter called TS); (ii) the conventional TS method after applying GACOS corrections for each interferogram (hereafter called TS-GACOS); (iii) the conventional TS method integrated with the APS model only (hereafter called TS-APS); and (iv) our proposed method, i.e. the conventional TS method integrated with the APS model and the GACOS correction for each interferogram (hereafter called TS-GACOS-APS). The major difference between the TS-APS and TS-GACOS methods is the APS parameter, for which TS-APS estimates the whole spatially and temporally correlated atmospheric delay for all acquisitions. Whereas for TS-GACOS-APS, the APS parameter only estimates the atmospheric delay residual, which can (to first order) be treated as random noise because a large portion of the atmospheric errors of its interferograms have been corrected. Fig. 7 shows the InSAR time series results for the two tracks T52 and T154 from the four abovementioned methods. The 3D GPS displacements were converted to LOS displacements according to Fialko et al. (2001). If the black circled area in Fig. 7a1 and b1 is considered, which forms part of the overlapping region of both tracks T52 and T154, an apparent long wavelength signal with a gradient from the northwest to southeast is seen on the uncorrected TS results (i.e. those without GACOS corrections) of track T52. However, the same region on track T154 shows a different pattern in the TS-only displacements, implying that this is not an actual deformation signal. After applying GACOS corrections, the atmospheric effect is largely reduced (see Fig. 7a2 and b2), with a further weakening of the short wavelength residual by applying the APS model apparent from Fig. 7a4 and b4. Fig. 7c shows detailed comparisons of displacement time series between InSAR and GPS at four sample stations where substantial displacements occurred. It can be seen that considerable reductions in the scatter arise for the TS-GACOS-APS time series over the TS-only time series. Without applying the GACOS correction, the displacement time series is either overestimated (e.g. station MUL1) or underestimated (e.g. stations LOK1 and TEN2). On station MAHA, where an obvious residual phase ramp is observed on the final displacement map (see Figs. 7a1, a3 and 8a3), the TS-APS method even provides a completely incorrect displacement time series due to the existence of the phase ramp. To quantify the reduction, GPS against InSAR displacement differences were computed using all epochs from all 34 stations which were located within the coverage of either the track T152 or T154 interferograms ( Fig. 7a1 and b1), and the two sets of values are shown as scatter plots in Fig. 7d for the four methods. A 0.77 cm RMS InSAR-GPS difference was obtained with the proposed TS-GACOS-APS method, improved from TS (1.95 cm) and TS-APS (0.88 cm), respectively, with an overall correlation of 0.80 between the TS-GACOS-APS derived deformation and GPS deformation time series, compared with 0.45 with TS alone. Note that these statistics only reflect the improvements at the 34 GPS stations, and do not indicate the method's performance for all parts of the interferograms. The RMS difference for TS-GACOS is smaller than the TS method, but both are greater than the TS-APS and TS-GACOS-APS methods because of the filtering of the short wavelength signals by the APS filter.

Validation with the overlapping swath
As the interferograms from tracks T52 and T154 were acquired on C. Yu, et al. Remote Sensing of Environment 251 (2020) 112097 different dates under different tropospheric conditions, their observations can be considered nearly independent and hence used for validation purposes. Before comparing their displacements, we computed the azimuth and elevation angles for the two tracks pixel-by-pixel and these are shown in Fig. 8c, producing a mean azimuth angle difference of 1.45 degrees and a mean incidence angle difference of −9.98 degrees. The magnitudes and variations of these differences are generally small (the standard deviations are 0.01 degrees for the azimuth and 0.29 degrees for the incidence), meaning that their displacements should be almost identical. A similar cumulative displacement pattern from the overlapping parts of the two tracks can be seen in Fig. 7a4 and b4 when using the TS-GACOS-APS method, with the red box denoting an area of particular similarity. To assess the similarity in more detail, differences in the cumulative displacement from the two tracks and the four correction methods were computed for all the overlapping region (42.4S-40.8S, 172.9E-174.5E) and plotted in Fig. 8. The smallest RMS of the differences across the region of 1.1 cm is obtained with the TS-GACOS-APS method, compared with the 2.4 cm RMS difference using the TS method. Also, the correlation coefficient between the displacements for all pixels increases from 0.29 for TS to 0.45 for TS-GACOS-APS. While the RMS difference using the TS-APS method is 1.5 cm and only slightly larger than the 1.1 cm obtained with TS-GACOS-APS, it can be seen from Fig. 8a3 that a long wavelength residual remains in the TS-APS solution which is not apparent in the TS-GACOS-APS solution and hence is an unmitigated atmospheric, not deformation signal. This long wavelength signal is also not apparent in the TS-GACOS differences shown in Fig. 8a2, but there are fairly large local, short wavelength variations still, with the RMS of the differences being 2.1 cm, but which reduces to 1.1 cm after applying the APS filter. Hence it has been shown that applying GACOS corrections leads to a reduction in the spatially and temporally correlated atmospheric errors, and that the APS filter is then able to reduce the short wavelength atmospheric residuals in the InSAR time series.

Afterslip modelling
We used the TS-GACOS-APS InSAR time series results from tracks T52 and T154, combined with the detrended GPS time series, for afterslip modelling. The cumulative displacements of InSAR and GPS were used to implement a static inversion and then the displacement time series of InSAR and GPS were used for a time-dependent inversion to recover the full history of the afterslip.
We first evaluated the spatial pattern of the observed deformation using Eq. (3) based on track T52's time series result with the TS-GACOS-APS model. It can be observed from Fig. 9d that the distribution of the amplitude (C in Eq. (3)) is similar to the distribution of the cumulative displacement (Fig. 7a4). The overall misfit of the logarithmic fitting is 0.36 cm, indicating the post-seismic surface deformation can be reasonably described by afterslip as suggested in Khoshmanesh et al. (2015).

Static afterslip modelling
We used the cumulative InSAR displacement maps of tracks T52 and T154 of the TS-GACOS-APS method. We calculated the GPS cumulative displacements from their time series by assuming a temporal logarithmic function (Eq. (3)) where the coefficients were constant parameters estimated by least squares (an example is shown in Fig. S1). The standard deviations of the GPS displacement misfits after fitting this equation were 0.42 cm (East), 0.51 cm (North), and 0.90 cm (Vertical), respectively. Considering the vast number of the spatially dense InSAR data points (~1500 compared to ~100), the relative weight between GPS and InSAR was determined as 5:1 to reduce the dominance of the InSAR data in the inversion based on the number of data points and their error variances (following Árnadóttir et al., 2004). A series of smoothing factors were employed in the slip inversions; the corresponding misfit RMS and slip distribution roughness are plotted in Fig. 9e. The optimal value was determined to be 5e5 when increasing the slip roughness will not be rewarded by considerable RMS reduction.
We utilized the fault geometry from Hamling et al. (2017), in which a subduction interface along with 19 crustal fault segments were included. To minimise the number of free parameters, we included only the five major crustal faults where there were the greatest co-seismic displacements (the Humps, Hope, Jordan Thrust, Kekerengu and Needles faults), together with the subduction interface, which was extended further north to include the entire northern South Island. We discretised the crustal fault planes into 2 by 2 km patches and the subduction interface into 4 by 4 km patches. For each patch, we estimated its strike and dip slip components. The estimation was a linear procedure in a least squares sense for rectangular dislocations in a uniform elastic half space (Okada, 1992).
The observations, modelled surface displacements and residuals for the two ascending InSAR tracks (T52 and T154) are shown in Fig. 9. The major deformation pattern is well explained, including the two major lobes and the northeast minor lobe on both tracks, with the RMS difference between the observed and predicted displacements of both tracks being less than 1 cm. The observed and modelled surface displacements at the GPS stations are shown in Fig. 10. The GPS displacements have both considerable horizontal and vertical components such as at stations CMBL, SEED and WITH, west of the Needles fault. The best-fit afterslip model reveals that the uplift at the GPS stations CMBL, SEDD and WITH is largely related to the oblique slip on the subduction interface and the horizontal displacement on GPS station CMBL may contribute to the right lateral striking slip on the Needles fault. There are small residuals along the Jordan Thrust's surface trace, suggesting the existence of small and shallow reverse slips on the shallow thrust. Northwest of the Jordan Thrust, the large residuals are probably due to the shallow (6.9 km) Mw 5.3 aftershock on 18 November 2016.
The fault slip distribution is shown in Fig. 11. The afterslip on the Needles fault, right-lateral dominating, is deeper (10-20 km) than the co-seismic slip (0-10 km), with considerable oblique slip found also on the subduction interface beneath the Needles fault and the north edge of the South Island. Another major oblique slip source on the subduction interface is located directly below the co-seismic slip, but with a smaller magnitude of maximum ~0.6 m (a moment release of Mw 6.9) compared with 5 m of co-seismic slip proposed by Hamling et al. (2017). Slight pure reverse slip is also observed on the shallow part of the interface to the east of station LRR1 (~30 km northeast of Kaikōura). We compared the afterslip distribution on the interface using all four of our time series analysis methods, as shown in Fig. 11b. The bestfit is given by TS-GACOS-APS, with an RMS difference of 0.7 cm and a correlation coefficient of 0.80 between the observed and modelled displacements, providing improvements over those from TS-APS of 59% for the RMS difference and 23% for the correlation coefficient. Our best-fit afterslip model has a smaller slipping area but follows a similar pattern compared with that of Wallace et al. (2018). As a result, we obtained a larger slip of about ~0.6 m on the interface but a smaller moment release of Mw 6.9, compared to the ~0.5 m interface slip of Wallace et al. (2018), corresponding to Mw 7.4. Nevertheless, our bestfit model fits reasonably well to the two InSAR tracks with an overall RMS difference less than 1 cm, including the descending track T73 (see Section 5.2). We have also conducted a checkerboard test on the subduction interface region by simulating a gridded slip pattern with a Gaussian random noise (variance 0.1 m: 10% of the simulated slip) added to the observations and plotted in Fig. 11c. Although the slip of the deeper parts of the fault (> 40 km) has a lower recovered resolution compared to the shallow parts as expected, it still explains the distribution of the gridded slip sources, with an RMS difference between the simulated and recovered surface displacements of 0.43 cm.

Time-dependent afterslip modelling
In order to investigate the full afterslip history, we performed a time-dependent afterslip modelling using the Principal Component Analysis-based Inversion Method (PCAIM) software (Kositsky and Avouac, 2010). The software relies on principal component analysis of the surface displacement time series without any adjustable parameters other than the number of components and the spatial smoothing factor. We used the TS-GACOS-APS derived InSAR time series and the detrended GPS time series in the modelling. The time step of the timedependent modelling was set to 6 days to reduce computation time and, more importantly, to comply with the InSAR time series of track T52 which is acquired every 6 days during the first two months after the main shock. In order to ensure all InSAR data points were located on the 6-day time step, we interpolated track T154 using the logarithmic Eq.
(3) onto track T52's acquisitions as was used in Kositsky and Avouac (2010). We expect that the error introduced by the interpolation is negligible as there is only a difference of one day between their acquisition times. The InSAR displacement maps were then down-sampled using a quad-tree quantisation algorithm to reduce the high spatial correlation. This resulted in ~7000 data points from GPS and ~ 110,000 data points from InSAR. We used the same relative weight between GPS and InSAR as in the static inversion determined by the number of data points and their error variances, which reduced the dominance of the InSAR data in the inversion (Árnadóttir et al., 2004), and employed the same fault segments, resulting in ~75,300 free parameters in total.
We determined the optimal spatial smoothing factor simultaneously with the number of principal components based on the reduced chisquare of all the observations (Fig. 12b). The reduced chi-squares with four (1.15) and five components (1.16) are similar but slightly smaller than those with three components (1.26), and therefore four components was selected as the optimal value. To avoid under-smoothing of the slip distribution, the optimal smoothing factor was determined as 10e5, which is the largest possible value as shown in Fig. 12b: increasing further will result in a significant jump of the reduced chisquare. Note that the smoothing factors for the static and the time-dependent inversions are different as the static smoothing factor is determined only by the one year cumulative displacement observations (Fig. 9e) whilst the time-dependent smoothing factor is determined by the whole displacement time series (Fig. 12b). Nearly 50% of the afterslip occurred in the first two months (November to December 2016) of the post-seismic period with a magnitude of up to ~0.35 m, which is also evidenced by the GPS displacement. The average slip magnitude reduced to less than 0.1 m between September and December 2017. The slip reaches a maximum at a depth of ~32 km which is a co-seismic slip deficit area (Fig. 12a), while its slip magnitude reduced to ~0.15 m between December 2016 and January 2017, to ~0.1 cm between March and June 2017, and to ~0.05 m between September and December 2017. Although the descending track T73 was not used in the afterslip inversion owing to its short and inconsistent time span with respect to the ascending tracks, we used it to validate the time-dependent afterslip model by simulating its corresponding cumulative displacement based on the cumulative afterslip between November 2016 and March 2017. The observed cumulative displacement of track T73 was computed in the same manner as the two ascending tracks, and the modelled and residual displacements are shown in Fig. 9c. Due to the different satellite geometries compared to the ascending tracks, the deformation signal on track T73 shows a different pattern but is well explained by our inverted time-dependent afterslip distribution, with an overall RMS difference between the modelled and observed displacement of 0.91 cm. The modelled and observed GPS displacements fit well with mean RMS differences of 0.31 cm, 0.28 cm and 0.65 cm for the east, north and up components respectively (full comparisons for all GPS stations are shown in Fig. S2). This provides further evidence of a reasonable modelling by the two atmospherically Fig. 9. Observed, modelled and residual interferograms for Sentinel-1 tracks T52 and T154 based on the best-fit static afterslip model, and track T73 based on the time-dependent afterslip model. The RMS differences indicated are between the modelled surface displacement and the observations. The black contours are the afterslip on the subduction interface with an interval of 0.1 m. (d) shows the spatial distributions of the post-seismic surface displacement amplitude, based on track T52's time series results. (e) shows the trade-off test to determine an optimal spatial smoothing factor for the static inversion. In summary, similar to the conclusions of previous studies (Jiang et al., 2018;Wallace et al., 2018), postseismic deformation following the Kaikōura earthquake provides clear evidence of afterslip on the southern Hikurangi subduction interface. It should be noted that the afterslip may be underestimated because of the SAR satellite geometry and absence of offshore observations, and we have neglected potential contributions from poroelastic rebound and viscoelastic relaxation, which we consider beyond the scope of this study. Poroelasticity is mainly visible as uplift and subsidence in the near-field with shallow processes and thus small spatial extent surface movement (Peltzer et al., 1996). This may relate to the small residuals seen near the Jordan Thrust trace in Fig. 9. After such a giant event, the crust would consist of an initial elastic rebound followed by a transient element of deformation controlled by the viscosity (Nur and Mavko, 1974), which means the viscous deformation may dominate over a decadal timescale, but is obscured by afterslip in the early stage of relaxation. To distinguish afterslip and viscoelastic deformation, further years (e.g. decades) of observations are required.

Discussion
The proposed TS-GACOS-APS InSAR method was validated with GPS displacements and by comparing the overlapping region of two independent tracks, with the small residuals of the afterslip model providing, to some extent, another validation of the retrieved cumulative displacement. Here we highlight two important aspects regarding our atmospheric correction model in InSAR time series analysis and the implication of the afterslip distribution for co-seismic slip models.

Long wavelength signal and orbital ramp deficiencies
The spatially and temporally correlated atmospheric errors are problematic in InSAR time series because APS filters, not only the one used in this study but many other filters used in conventional time series analysis (e.g. Ferretti et al., 2001;Hooper et al., 2007;Lauknes et al., 2011), are only valid with random noise. Shirzaei (2013) used time-frequency analysis rather than triangular and boxcar filters to reduce the effects of systematic artefacts, such as spatially correlated and temporally uncorrelated components of the atmospheric delay. The temporal deformation model, however, either linear or non-linear, requires the deformation signal to mix only with temporally random atmospheric noise. With our proposed method, the long wavelength and topography-related atmospheric errors are removed before filtering, reducing the spatial-temporal correlation of the atmospheric error and resulting in smaller and more random residuals. These largely satisfy the basic time series analysis assumptions and hence are more efficiently handled by filtering. The performance of our TS-GACOS-APS method mainly relies on the accuracy of the GACOS correction map and its quality control, to ensure only those corrections not deemed outliers are applied to interferograms.
In many InSAR time series analysis methods to date, a best-fit ramp has been removed from the original interferogram to compensate the long wavelength effect, caused by both orbital and atmospheric errors. However, the time series results of track T52 show that after the ramp removal, which is done in the TS and TS-APS methods, the long wavelength error still exists, suggesting that the atmospheric error can have more complicated spatial patterns than those of a ramp. Furthermore, even though a large portion of the atmospheric effect manifests as long wavelength signals, it has a completely different physical origin and spatial-temporal characteristics than those of an orbital ramp, implying that the best-fit ramp based on the original interferogram may be biased by both the atmospheric and the orbital long wavelength components (e.g. when they have opposite directions). We therefore suggest, as implemented here, that GACOS-based corrections are first applied to mitigate the long wavelength signal caused purely by atmospheric effects, leaving the orbital error induced signal, if any, dominant, and then to estimate the best-fit ramp from the corrected interferogram, as is done in the developed TS-GACOS-APS method.

Linkage between co-seismic slip and afterslip
The existence of co-seismic slip on the subduction interface, if any, remains controversial, with Hamling et al. (2017) suggesting two bestfit slip models, one without an interface slip and one including the interface and contributing ~10-30% of the total moment release. Hollingsworth et al. (2017) preferred a model that placed ~60% of the seismic moment release from the interface, however, with a centroid depth of ~22 km, shallower than those in Hamling et al. (2017). Holden et al. (2017) even found an opposite conclusion that the coseismic slip model which excludes the Hikurangi subduction interface explains the nearfield waveforms better. To summarise, the co-seismic slip on the subduction thrust is difficult to disentangle from crustal fault deformation measurements during the earthquake. However, the postseismic period may provide valuable information about the interface as it has been recognised that afterslip following large earthquakes is usually complementarily located with the co-seismic slip as compensations in magnitude and distribution, such as the 1999 Izmit earthquake (Reilinger et al., 2000) and the 2004 Sumatra-Andaman earthquake (Chlieh et al., 2007). By investigating our afterslip distribution, we find it agrees well with the co-seismic interface slip proposed by Hamling et al. (2017), with the major afterslip occurring at a co-seismic slip deficit area, distributing immediately below the main co-seismic slip area. This co-and after-slip distribution relationship provides indirect evidence that the subduction interface has been triggered during the mainshock, with a relatively small moment release compared with the crustal faults at a depth of ~30-35 km and is unceasingly moving afterwards.
The observations for the 2016 Kaikōura earthquake demonstrate that multiple faults can be triggered during a single event and the fault slip can propagate through fault stepovers and splays over a long distance (> 100 km, Fig. 1, e.g. Hamling et al., 2017;Xu et al., 2018). When considering the co-seismic slip on the interface, the far field subsidence and non-double-couple components of global moment tensors are better explained . This implies a much more complex event that undergoes slips along numerous faults with diverse orientations and directions, propagating from the hypocentre both horizontally to the adjacent fault segments and vertically to the underlying deep thrust.

Conclusions
Sentinel-1 line-of-sight range changes and GPS displacements were used to recover the time-dependent afterslip on the Hikurangi megathrust beneath the Marlborough Fault System following the Mw 7.8 Kaikōura earthquake. We detected displacements to an accuracy of about 1 cm through the development of a method (TS-GACOS-APS) to mitigate the spatially and temporally correlated atmospheric errors in InSAR time series. The method mitigates the long wavelength and topographic related atmospheric errors by applying HRES-ECMWF based atmospheric corrections using the GACOS service to those interferograms which pass a cross-validation test. This results in the atmospheric residuals in the corrected interferograms having a nearly random distribution, and then the deformation time series is extracted after applying an APS filter. The resultant cumulative InSAR displacements agreed with GPS to 0.77 cm RMS, considerably improved by 61% over the conventional TS method (1.95 cm RMS). The success of the InSAR displacement measurement was also confirmed from an RMS displacement difference of 1.1 cm for the pixels in the overlapping region of two ascending Sentinel-1 tracks using TS-GACOS-APS, compared with 2.4 cm using TS alone (54% improvement) and 1.5 cm with TS-APS (27% improvement). The APS filter reduced the short wavelength residuals substantially but, unlike when also including GACOS, failed to remove the long wavelength error. Improvements with the TS-GACOS-APS method were obtained on both Sentinel-1 tracks in a coastal environment, in particular on track T52 with a spatial coverage of over 250 km, suggesting that the method can be widely applied to studies with large spatial extents, where the long wavelength and topographic atmospheric errors are more variant and dominant.
In terms of geophysical modelling of afterslip using 13 months of ascending Sentinel-1 data and GPS station displacements, the inversion RMS between the observed and modelled surface displacement improved from 2.8 cm with no GACOS InSAR atmospheric corrections applied, to 1.7 cm with only an APS filter applied, and to 0.7 cm with the developed method using the combination of GACOS atmospheric corrections and an APS filter (TS-GACOS-APS). Compared with the previous studies that used short time span (4 months) descending Sentinel-1 interferograms and/or the GPS data, the wide coverage of the Sentinel-1 ascending interferograms, after careful atmospheric mitigation, provides more comprehensive constraints on the afterslip inversion and locates the spatial-temporal distribution of the afterslip at depths between 25 and 35 km. This location enables the investigation of the spatial relationship between the afterslip and co-seismic slip distributions on the Hikurangi megathrust, which implies that the interface has moved during the mainshock at a depth between 25 and 30 km, with the following afterslip occurring between 25 and 35 km. We also found considerable right lateral striking slip on the Needles   Fig. 12. (a) Cumulative afterslip distribution on the Hikurangi subduction interface from the time dependent inversion and the co-seismic slip distribution obtained from Hamling et al. (2017). The blue triangles are GPS stations and the red dotted lines indicate the modelled fault surface traces from Hamling et al. (2017). (b) is the trade-off test to determine the optimal number of components and the spatial smoothing factor. (c) is the afterslip of different periods within November 2016 and December 2017. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) fault, oblique slip on the subduction interface beneath the Needles fault and the northern edge of the South Island, and reverse slip on the shallow part of the interface located ~30 km northeast of Kaikōura.
The major benefits of the developed InSAR time series atmospheric correction method are: (i) the individual atmospheric corrections from GACOS are globally usable with a short delay (< 2 days); (ii) the long wavelength and topography related atmospheric errors are removed independently without contamination of other error sources, such as an orbital ramp, and without assuming the atmospheric signal can be fitted linearly to the height or simply by a ramp in a sufficiently small window; (iii) it is suitable for both small and large areas, which is especially helpful to construct deformation maps at national or continental scales.

Declaration of Competing Interest
None.