Slow‐Moving Landslides Triggered by the 2016 Mw 7.8 Kaikōura Earthquake, New Zealand: A New InSAR Phase‐Gradient Based Time‐Series Approach

Earthquake‐triggered slow‐moving landslides are not well studied mainly due to a lack of high‐resolution in‐situ geodetic observations both in time and space. Satellite‐based interferometric synthetic aperture radar (InSAR) has shown potential in landslides applications, however, it is challenging to detect earthquake‐triggered slow‐moving landslides over large areas due to the effects of post‐seismic tectonic deformations, atmospheric delays, and other spatially propagated errors (e.g., unwrapped errors caused by decorrelation noises). Here, we present a novel InSAR phase‐gradient‐based time‐series approach to detect slow‐moving landslides that triggered by the 2016 Mw 7.8 Kaikōura earthquake. Twenty‐one earthquake‐triggered large (>0.1 km2) slow‐moving landslides are detected and studied. Our results reveal decaying characteristics of the temporal evolutions of these landslides, that averagely 3.9 years after the earthquake, their post‐seismic velocity will decay by 90% to reach approximately the pre‐seismic level. Our study opens new perspectives for understanding mass balance of earthquakes and helps reduce associated hazards.

• A new interferometric synthetic aperture radar phase-gradient based time-series approach used to detect earthquake-triggered slow-moving landslides • Twenty-one slow-moving landslides triggered by the 2016 Mw 7.8 Kaikōura earthquake are detected and monitored • Our results reveal decaying characteristics of post-seismic velocity of the earthquake-triggered slow-moving landslides

Supporting Information:
Supporting Information may be found in the online version of this article.
2 of 11 tend to be debris and rock avalanches and flows (e.g., Keefer, 1984; C. I. Massey et al., 2020). The movement patterns of coherent (Keefer, 1984), co-seismic landslides during and after an earthquake are also difficult to monitor using in-situ geodetic techniques as the location of where the landslide is often unknown prior to the earthquake, which precludes the instrument being deployed in advance. Even if a landslide is pre-existing and monitoring equipment is installed on it, to capture its reactivation during an earthquake, requires an earthquake with ground shaking large enough to trigger movement to occur during the lifetime of the monitoring. Remote sensing data sets combined with digital image correlation and deep learning techniques can be used to derive landslide movement patterns from pre-and post-landslide movement, optical image and lidar surveys (e.g., Senogles et al., 2022), however, such techniques lack the temporal resolution to identify the changes in landslide movement patterns over time, therefore not allowing patterns of movement to be linked to their causes. As a result of these constraints: (a) the movement patterns at high temporal resolutions of coherent landslides, in response to, and after strong ground shaking remain largely unknown; and (b) the mechanisms governing such displacements are poorly understood. Most coherent slides and slumps are typically assumed to displace via "sliding block" style mechanisms (e.g., Newmark, 1965). Such mechanisms assume sliding occurs along a plane, with or without internal deformation of the sliding mass (e.g., Makdisi & Seed, 1978). In such analyses, permanent slope displacement occurs when the acceleration of the slope, caused by the earthquake inertia forces, is larger than the yield acceleration of the slope. Movement of the landslide is assumed to stop when the earthquake accelerations are below the yield acceleration. However, movement of the landslide could continue if for example, the coseismic displacement causes a change in the strength of the moving mass and/or a change in pore-water pressures within the mass. Earthquake-induced landslides, could therefore in theory, continue to move after an earthquake.
Satellite-based interferometric synthetic aperture radar (InSAR) techniques are increasingly being used to monitor slow-moving landslides (e.g., Bekaert et al., 2020;Schlögel et al., 2015;Sun et al., 2015;Y. Xu et al., 2022), over large areas (e.g., ∼250 km) with high spatial resolution (e.g., ∼10 m). However, it is still challenging to use current InSAR techniques to detect earthquake-triggered slow-moving landslides because of at least three reasons: (a) the locations of these landslides are often unknown and their distributions can cover large areas (e.g., up to 100 km away from the epicenter); (b) the weak (e.g., several centimeters per year) and local (e.g., ∼1 km) landslide deformations would be largely contaminated by atmospheric delays (e.g., Cao et al., 2021; and post-seismic deformation; and (c) possible unwrapping errors mainly caused by decorrelation noises particularly over dense vegetation regions (e.g., Yunjun et al., 2019), usually cause large uncertainties or even wrong values in InSAR-derived time-series results (e.g., velocity map or time-series displacements), especially when the pixel of interest is far away from the reference pixel due to the spatially propagated characteristics of the unwrap errors.
To address the above-mentioned problems for hunting and monitoring earthquake-triggered slow-moving landslides over large areas, here we propose a new InSAR phase-gradient based time-series approach (details in Section 2). Using the new method, we identify and investigate slow-moving landslides triggered by the 14 November 2016 Mw 7.8 Kaikōura earthquake, New Zealand, which is one of the most complex earthquakes ever recorded (e.g., Hamling et al., 2017;W. Xu et al., 2018). More than 29,000 mapped landslides were triggered by the Kaikōura earthquake over an area of about 10,000 km 2 . Most were mapped using aerial imagery, optical satellite images, and digital elevation model (DEM) differences, and most comprised rapid debris and rock avalanches (C. Massey et al., 2018;C. I. Massey et al., 2020). In contrast to these landslides, many large coherent landslides were already present on the slopes, which were strongly shaken by the 2016 Kaikōura earthquake. Most of these landslides can be classified as rock slides or slumps (Hungr et al., 2013), but their activity rates before, during and after the Kaikōura earthquake were largely unknown. Given the "coherent" morphology of these landslides it could be assumed that they move relatively slowly (Hungr et al., 2013), or rapidly over very short time periods, thus giving the impression that they are either dormant or creeping, and/or move episodically in response to strong earthquake shaking. In this paper, we focus on these large coherent rock slides/slumps, which are assumed to be slowly moving and could be either first-time failures (initiated by the 2016 Kaikoura earthquake) or reactivated by the earthquake. Hereafter, for simplicity, we define these earthquake-triggered slow-moving landslides as eSMLs.
CAO ET AL.

InSAR Phase Gradient Based Time-Series Approach
InSAR phase gradients calculated from wrapped interferograms can be written as: where ∆ ( ) and ∆ ( ) are the phase gradients (see Figure S1 in Supporting Information S1) at location ( ) along the "Range" direction (i.e., line-of-sight, LOS) and "Azimuth" direction (i.e., flight direction), respectively, { } denotes the wrapped phase, { } is the operator of calculating phase gradient from wrapped phase. Here, we consider that the topographic phase in the interferograms has been removed from the wrapped phases using an external DEM. Note that the "phase gradient" in this study means the "adjacent-pixel phase difference," and we do not consider the pixel-spacing of the SAR images.
Building on the basic idea of small baseline subset (SBAS) InSAR (Berardino et al., 2002) and using multi-temporal short-baseline interferograms, we can calculate and get multi-temporal of phase-gradients at each pixel along the two directions, respectively, which can be written as where M means the number of the used short-baseline interferograms. Relationships between multi-temporal phase gradients and time-series of displacement gradients can be described by where is the design matrix determined by the multi-temporal interferometric network (Z. Li et al., 2019). Note that because the "Range" phase gradients and the "Azimuth" phase gradients have the same interferometric network, the design matrixes are the same. Time-series solutions of the gradients at the two directions can be written as: where N denotes the number of SAR acquisitions. ∆ ( ) and ∆ ( ) represent the estimated displacement gradients of the ith SAR acquisition at position of ( ) along the two directions, respectively. We consider all of the interferograms are included in one network, and we solve time-series of the gradients using a weighted least-square approach that is the weighting matrix of the multi-temporal phase gradients which can be written as where 2 ∆ denotes the variance of the kth phase gradients. By assuming the decorrelation noises at the two adjacent pixels are not correlated, we calculate the phase-gradient variance as below based on the variance-covariance propagation theory where Var{⋅} means the operator of calculating variance, Equations 7a and 7b represent the variances of the phase-gradient along "Range" and "Azimuth," respectively. Here, we model the variance of interferometric phase according to its coherence and the used multi-looks (Tough et al., 1995). We evaluate the quality of the time-series solutions of the displacement-gradients using temporal coherence (Pepe & Lanari, 2006), which can be calculated as below: where is the temporal coherence that varies between (0, 1) with quality changes from "bad" to "good," ∆ denotes the measurement (i.e., observation) of the kth phase gradient ( = 1, 2, . . . , ), and ∆ denotes the phase-gradient that recovered according to the time-series gradient solutions. We should note that, for regular time-series InSAR (TS-InSAR), the temporal coherence would be affected by both unwrapping errors and decorrelation noises, and for that of the gradient-based approach, it would be mainly affected by decorrelation noises.
In comparison to regular TS-InSAR (e.g., SBAS), advantages of the new method include: (a) effects of those spatial correlated large-scale signals (e.g., atmospheric delays and post-seismic deformations) are minimized; and (b) reliability of the gradient results would be independent of the size of research area, because it would be not affected by spatially propagated unwrapping errors, while that of regular TS-InSAR would be often related with the distance between the pixel of interest and the reference pixel, due to possible unwrap errors particularly over vegetated regions. Also compared with stacking-based gradient approach (e.g., Fu et al., 2022;Hu et al., 2020), the advantages of our new approach include: (a) quality of the time-series products (e.g., gradient velocities, or time-series solutions of displacement-gradients) can be well evaluated based on the temporal gradient-phase closures; and (b) temporal evolutions of the landslides can be analyzed directly from time-series of the gradient. But here, we should also note that, our new method cannot resolve the results over vegetated areas or the areas where are strongly affected by decorrelation noises.

Application to Monitor Slow-Moving Landslides Triggered by the 2016 Mw 7.8 Kaikōura Earthquake
To detect the slow-moving landslides that triggered by the 2016 Mw 7.8 Kaikōura earthquake, we use 166 Sentinel-1 images acquired from ascending track during October 2014 to September 2021, and we totally generate 776 short-baseline interferograms that are all included in one network (see Figure S2 in Supporting Information S1). We simulate the topographic phase and the flatten earth phase using the 30 m Copernicus DEM (https://spacedata.copernicus.eu/), and we remove the simulated components from the raw interferogram to get the differential interferogram, and we apply 10 and 2 multi-look numbers to the range and the azimuth direction, separately, which derives a spatial resolution of around 30 m for each pixel. To further mitigate the effects of decorrelation noises, we apply an improved Goldstein filter (Z. W. Li et al., 2008) to the differential interferogram, and we calculate the phase gradients from the filtered differential interferogram.
We estimate time-series of the gradients using the weighted least squares method (see Equation 5), and we calculate the temporal coherence of the time-series solution to evaluate its quality and remove low-quality ( 0.6 ) pixels that are seriously affected by decorrelation noises in further analysis. We correct DEM errors of the raw gradient time-series based on the spatial baseline history (Fattahi & Amelung, 2013 (Figures 1a-1d). We calculate the differences in the gradient velocities before and after the earthquake to detect those earthquake-triggered local small-scale deformations, and considering gradient velocities have different sensitivities along the "Range" and the "Azimuth" directions, we fuse the two directions' velocity differences as follow: where ∆ is the fused gradient velocity of the two directions that can be used to describe changes of the local signals before and after the earthquake; and ∆ and ∆ are the velocity differences of the displacement gradients before and after the earthquake along the range and the azimuth directions, respectively.
From Figure 1, we can find that those earthquake-triggered local signals are mainly distributed along the ruptures (red line) or active faults (e.g., Clarence fault and London Hill fault) which are close to the main fault ruptures. But we also can see that not all of local signals are associated with eSMLs, with some being caused by short-wavelength afterslip or postseismic rebound (e.g., signals near London Hill fault). By using the fused gradient velocity (Figure 1e), the gradient velocities before and after the earthquake (Figures 1a-1d), and jointly considering the geomorphology (based on DEM and Google Earth) and time-series evolutions of the displacement gradients, we pick up earthquake triggered slow-moving landslides (gray closed curves in Figure 1f, Figure S4 in Supporting Information S1), and totally 21 eSMLs are detected within an area of about 80 km × 80 km. Note here we use the gradient results only for detecting the locations and approximate boundaries of the eSMLs. After figuring out the locations of these eSMLs, we dig out the subset interferograms for each landslide, and we unwrap the local subset interferograms using the minimum cost flow algorithm (Werner et al., 2002). Note here the subset interferograms, that associated with only one landslide, only cover a very small area (e.g., ∼5 km 2 ), as a result, we can consider that the spatially propagated unwrap errors can be largely constrained compared with that of processing the full interferogram. For each subset interferograms (i.e., each landslide), we select one reference point from stable and good-coherence area (i.e., an area with no or limited deformation with good interferometric coherence). Then we use the regular SBAS approach to estimate the spatio-temporal displacements for each landslide. Workflow of our methods that built on the new phase-gradient approach used for detecting eSMLs caused by the 2016 Mw 7.8 earthquake can be summarized into seven steps (see Figure S3 in Supporting Information S1).
We calculate the pre-and post-seismic velocity patterns of the detected eSMLs based on the displacement time-series solutions that estimated case by case using a local reference point, and we find the pre-seismic velocities (October 2014 to November 2016) of these eSMLs are all around zero or can be ignored, while the post-seismic velocities (November 2016 to January 2018) all increased significantly (e.g., ∼6 cm/year). Adopting the scheme of Hungr et al. (2013), we find most of the landslides (eSMLs) comprise compound planar or rotational slides in rock (see Table S1 in Supporting Information S1). Given their size (source area), it suggests that their basal slide surfaces could be relatively deep-seated (C. I. Massey et al., 2020). Most are relict and reactivated in response to the Kaikōura Earthquake and thus their basal slide surfaces are likely to be either fully formed or at least in part formed. All of the detected eSMLs are larger than 0.1 km 2 (varies from 0.19 to 3.31 km 2 , see Table S2 in Supporting Information S1), however, we should note that some small eSMLs could be not detected in our study, as the signals of those small landslides are easier to be polluted by other source of errors or signals.

Temporal Evolutions and Models of the Detected eSMLs
To further analyze temporal evolutions of these eSMLs, we calculate time-series of the gradients and displacements of each eSML. Considering spatial heterogeneities of the displacement of each eSML, we estimate the time-series by averaging a 5 × 5 window (location of window center presented in Table S2 in Supporting Information S1) of the raw results in space at each SAR acquisition, and we select the center point for each eSML from the central deforming area with temporal coherence larger than 0.7. To pick apart the likely mechanisms driving the displacements, here we have converted the averaged LOS displacements to downslope displacements by assuming the moving direction is approximately along the slope. We can see that, before the earthquake, both the gradients (Figures 2a and 2b) and the displacements (Figure 2c) of the eSMLs are quite stable and fluctuate around zero (or evolve with a tiny velocity), after that, we can find significant changes (e.g., increase or decrease significantly) in the gradients and the displacements, which indicates that these changes are caused by the 2016 Kaikōura earthquake. After about 5 years since the earthquake, the accumulated post-seismic range-and azimuth-gradients reach to ∼8 mm (e.g., the 21st eSML) and ∼7 mm (e.g., the 4th eSML), respectively, and that of the downslope displacements reach up to about 25 cm (e.g., the 21st eSML).
Figure 2 also clearly shows that, temporal evolutions of the gradients and the displacements all show significant "decaying" characteristics after the earthquake, for example, the absolute velocities (gradient or displacement) increase significantly and reach to a maximum value just after the earthquake, while that decrease gradually in the years after the earthquake. To quantitatively evaluate the temporal evolutions, here, we propose to use the following piecewise function to model time-series of the downslope displacements  where pre( ) is a linear function that denotes the modeled pre-seismic process, and , are to be estimated parameters; post ( ) is a decay-shaped exponential function that represents the modeled post-seismic process, and , , 0 are to be estimated parameters. and will determine the decay shape, and 0 will determine the shift at the start-point associated with the co-seismic displacement of the eSML and the InSAR observation errors (e.g., decorrelation noise, atmospheric delay). We use a least-squares method to estimate the unknown parameters based on the time-series estimations. The best fitted models are illustrated as solid gray lines in Figure 2c, and the averaging value of the correlations between measurements and the models reaches 0.95 (see Table S3 in Supporting Information S1), which indicates that the used models are suitable for these eSML cases. From the models, we find that the averaged downslope velocity of the 21 eSMLs increased significantly from 4 mm/year before the earthquake to 34 mm/year after the earthquake, and the averaged co-seismic displacement, that estimated by co = post ( EQ) − pre( EQ) , is about 1.4 cm. Note that here we can consider that the co-seismic tectonic displacement, which can be up to meters (Hamling et al., 2017), has been removed from the co-seismic landslide movements, as we use a local reference point for each landslide and the residual co-seismic tectonic displacement should be very limited. From the displacement time-series (Figure 2c), we also can find weak seasonal variations that could be caused by pore pressure changes on the sliding slope due to seasonal variations of the rainfalls. We estimate the seasonal components using a simple annual sine function based on the model residuals, but we do not find clear correlations between the mean monthly rainfall and the seasonal amplitude (about 1.4 mm averagely) of the displacement (see Table S2 in Supporting Information S1), this indicates that the seasonal variations may be dominated by multiple factors (e.g., geological unit, sliding type).
To quantitatively evaluate the decaying characteristics of the eSMLs, here we calculate the temporal velocity evolutions of each eSMLs based on the first-derivative of the estimated postseismic model. Instead of analyzing the full sliding velocity, we calculate and compare the normalized post-seismic velocities only, which can be obtained by where post ( ) is the normalized post-seismic velocities of the eSML, that decays from 1 to 0, and means the time since the earthquake, is the same parameter as that in the displacement model, which controls the decay speed of the eSML. Post-seismic models of the normalized velocities (Figure 3a) clearly show that the landslide velocities 9 of 11 decay rapidly spanning years since the earthquake, and different eSMLs show different decaying speeds, for example, landslide velocities of the 10 th eSML decay to around zero within 1 year since the earthquake, whereas that of the 3 rd eSML experiences more than 8 years. The decay-timings (Figure 3b) of the post-seismic velocities reveal that after 6 years since the earthquake, the landslide velocities of the eSMLs all decay at least 85%, and within less than 10 years after the earthquake, all of the eSMLs would decay 95%. Averagely, post-seismic velocity of the eSML decay 90% within around 3.9 years since the earthquake. Spatial distribution of the post-seismic rainfalls over the regions of the eSMLs (gray histograms in Figure 3b) indicate that, the relationship between the rainfalls and the decay-timings is limited.

Discussion and Conclusions
We developed a new InSAR-phase-gradient based time-series approach to detect locations of earthquake-triggered slow-moving landslides over large areas, that are difficult to hunt using regular InSAR techniques (e.g., D-InSAR, Stacking-, or SBAS-InSAR, see Figures S5 and S6 in Supporting Information S1) due to the effects of post-seismic tectonic deformations, atmospheric delays, and possible unwrapped errors. In total, 21 slow-moving landslides triggered by the 2016 Mw 7.8 Kaikōura earthquake are detected adding to the catalog of co-seismic landslide failures previously identified (e.g., C. Massey et al., 2018;C. I. Massey et al., 2020).
Our results reveal decaying characteristics of the temporal evolutions of these landslides, that averagely after 3.9 years since the earthquake, their post-seismic velocity will decay 90% and recover to a level of near pre-seismic, and this may indicate self-healing of the earthquake-triggered slow-moving landslides in some degree. We suggest that this decay was likely related to a combination of internal deformation of the landslides, as well as slip along their basal slide surfaces, in response to the earthquake-triggered displacement, as post-earthquake landslide displacements did not correlate to rainfall or any other monitored factor. Time-series evolutions of the displacements of each eSML (Figure 2c) is likely to be a function of the local site "susceptibility" conditions, mainly geological (structure and geomechanical properties) and morphological (slope angle and local slope relief), along with the amplitude and duration of shaking each site experienced during the earthquake (e.g., Brain et al., 2021;. Given the landslide types (reactivations), their sizes (depths), and eSML time-series displacements, we hypothesis that the co-seismic changes in the velocities of the gradients and the displacements are likely to relate mainly to basal sliding mechanisms along with some distributed shear within the overlying slide mass (e.g., Bray & Travasarou, 2007;Makdisi & Seed, 1978). We also suggest that the post-earthquake time series displacements, given their patterns of displacement, are more likely to relate mainly to distributed shear within the slide mass material, possibly related to its consolidation in response to coseismic deformation. This is because the exponential decay in displacement, post-earthquake, appears unrelated to seasonal changes in stress caused by rainfall-induced fluctuations in pore-water pressures or subsequent earthquake loading, or any other environmental factor that could initiation movement (C. I. Massey et al., 2013). Such post-seismic deformation might relate to the consolidation of the slide mass, which will be further studied in the future.
Only one ascending track of Sentinel-1 data sets were used in this study, as that of the descending track is not well covered with enough data sets, thus some of the eSMLs signals may not be detected, due to the geometry limitations of InSAR (e.g., shortening and shading) and the 1-D displacement imaging limitation (e.g., InSAR cannot detect those landslides whose sliding direction is vertical to the LOS direction). In addition, due to the dense vegetation coverage over the Kaikōura region which could cause strong decorrelation noises in InSAR observations, those small eSMLs (e.g., less than 100 m × 100 m ) may not be well detected in our study. To better separate eSML signals from other non-landslide signals in the gradient-results, we may also consider to use some artificial intelligence techniques in future applications.