Investigating earthquake legacy effect on hillslope deformation using InSAR‐derived time series

Mountainous landscapes affected by strong earthquakes typically exhibit higher landslide susceptibility in post‐seismic periods compared to pre‐seismic conditions. This concept is referred to as the earthquake legacy effect, which needs to be better understood to develop an accurate post‐seismic landslide hazard assessment. The earthquake legacy effect is mainly assessed by monitoring either rapid landslide occurrences or slow‐moving landslides over time. However, landslide inventories provide discrete information both in spatial and temporal domains, whereas monitoring only selected slow‐moving landslides discards the response of hillslopes exposed to a given earthquake. Therefore, to provide a more comprehensive understanding of the concept, this research focuses on post‐seismic hillslope evolution by examining the deformation time series generated from the Interferometric Synthetic Aperture Radar technique over the area affected by the 2017 Mw 6.4 Nyingchi earthquake, China. We also analyse factors controlling these InSAR‐derived hillslope deformations. Our results show a high coherence between hydrologic conditions (i.e., precipitation and terrestrial water storage) and surface deformation in pre‐ and post‐seismic periods. The earthquake disturbs this strong correlation for a while (~2 years) right after the seismic tremor, and then, a seasonal deformation pattern depending on hydrologic conditions appears again. Our findings show that the average post‐seismic hillslope deformation is still higher than its pre‐seismic counterpart approximately four and a half years after the earthquake. These findings trigger further research questions regarding whether hillslopes could fully recover after a major earthquake or gain a new level of hillslope susceptibility caused by intense ground shaking.


| INTRODUCTION
Landslides pose serious threats to communities, especially in mountainous regions such as the Himalayan range where thousands of lives and billions of Euros are lost every year because of landslides (Upreti et al., 2001).For this mountainous landscape, earthquakes are a common triggering factor for landslides, causing disastrous cascading effects (e.g., Roback et al., 2018;Sato et al., 2007).Seismic shaking is not only responsible for co-seismic landslidesthe main secondary earthquake hazard (Daniell et al., 2017)-but can also be the reason for some long-term hillslope instability problems because of the intrinsic damage given to hillslope materials (e.g., Chen et al., 2020;Parker et al., 2015).In post-seismic periods, the earthquake legacy effect couples with hydrologic/anthropogenic disturbances and exacerbates hillslope instabilities (Kincey et al., 2021).Some argue that the landscape returns to pre-seismic landslide hazard level only after months or even years, depending on site-specific morphologic and hydrologic conditions typical of the area (e.g., Tian et al., 2020).
The process that any given landscape naturally restores to its preseismic slope stability conditions is commonly referred to as hillslope recovery.This recovery is mostly assessed through multi-temporal landslide inventories including various types of landslides including both shallow and deep-seated ones (Tanyaş, Kirschbaum, Görüm, et al., 2021).Specifically, pre-and post-landslide susceptibility levels are commonly identified based on the frequency of landslide occurrence, and the results are interpreted to infer the hillslope evolution processes in earthquake-affected areas.However, focusing only on landslide occurrences limits our observations with a small subset of the landscape under consideration.Shear strength perturbation could occur on any hillslopes affected by seismic shaking, regardless of landsliding (Tanyaş, Kirschbaum, Lombardo, et al., 2021).This also means that relatively higher hillslope deformation (HD) rates could be considered as a sign of relatively higher landslide susceptibility (i.e., the relative probability of landslide occurrence) compared to hillslopes with lower deformation rates.
InSAR-derived deformation time series could be used to investigate post-seismic hillslope reactions instead of multi-temporal landslide inventories, which provide discrete information in both spatial and temporal domains (e.g., Barth et al., 2020;Tang et al., 2016).InSAR techniques are widely used to provide a comprehensive picture of post-seismic hillslope response to a given earthquake (Bontemps et al., 2020;Lacroix et al., 2022;Song et al., 2022).For instance, Martino et al. (2020) examine hillslope reactions to the 2018 Mw 5.1 Montecilfone, Italy, earthquake using InSAR techniques and show activations and reactivations of landslides following the earthquake.Also, Cheaib et al. (2022) analyse slow-moving ancient landslides accelerated by the 2017 Mw 7.3 Sarpol Zahab, Iran, earthquake and demonstrate their kinematics using InSAR.In those contributions, hillslope observations are not bounded by rapidly occurring landslides but by continuous slope movements, albeit restricted to selected slow-moving landslides.
The earthquake legacy effect using surface deformations over the entirety of the area affected by a seismic tremor is yet to be explored.
And yet, our assumption here is that hillslope recovery analyses could be significantly enriched if observations are extended beyond failed hillslopes, encompassing even slight deformations (i.e., millimetre-scale) that take place over the entirety of seismically perturbed landscapes.InSAR is the only technique available to monitor such small-scale deformation rates over large areas for long periods.Therefore, following this hypothesis, here, we analyse HD patterns over the area affected by the 18 November 2017 Nyingchi earthquake over two periods: the 1.5 years before and the 4.5 years following the earthquake.
The manuscript will demonstrate the sudden increase in HD associated with the 2017 Nyingchi earthquake and a post-seismic overall landscape response where the HD still appears higher than its preseismic expression even after four and a half years from the mainshock.We will present our study site, methodology and results in Sections 2, 3 and 4, respectively.In Section 5, we will discuss our interpretation and share our concluding remarks.

| Study area
We examined an area located on the western edge of the Himalayan range (Figure 1).This area was shaken on 17 November 2017, by the M w 6.4 Nyingchi, China, earthquake, which occurred along a thrust fault.The same area was also affected by three large aftershocks (of magnitude greater or equal to 5.0) in the following months (USGS, 2023).Its direct effect on slope stability has been described by Zhao et al. (2019), with more than 1800 landslides mapped over a F I G U R E 1 The map of the study area overlaid by epicentres of the 2017, M w 6.4 Nyingchi earthquake and aftershocks (USGS, 2023).Glaciers (GLIMS Consortium, 2005;Raup et al., 2007) are indicated by green polygons.[Color figure can be viewed at wileyonlinelibrary.com] 530 km 2 area.After this seismic sequence, the same area was also affected by three subsequent earthquakes of magnitude 5.0 in 2019, 2020 and 2021 (Figure 1).This study thus focuses on a particularly complex system where the aforementioned earthquake sequence shook a rough mountainous landscape (80 km Â 95 km) with very large elevation gradients (up to $7 km over a relatively short distance).Despite the high mountainous system, the land surface temperature is above zero from March to December, and glacier bodies exclusively persist only across the high peaks (Figure 1).Anthropogenic impact is relatively limited in the area, with typically natural land cover (natural forest, grassland and barelands) and only a few patches being dedicated to agriculture (ESA Climate Change Initiative, 2022).

| Data
InSAR-derived deformations were generated for the same area to test the applicability of an interpretable model targeting HDs (He et al., 2023).However, in this research, we particularly looked at factors potentially influencing HD and analysed their variation over time.
To do so, we examined InSAR-derived surface displacements in relation to hydrologic variables, including precipitation (The Global Precipitation Measurement, GPM-Integrated Multi-satellite Retrievals, IMERG final product; Huffman et al., 2019) and terrestrial water storage (TWS, Li et al., 2019Li et al., , 2020)).IMERG final is derived based on an algorithm combining information from multiple sources, including satellite microwave, infrared precipitation estimates and ground gauge information.It provides precipitation with 0.1 spatial resolution and at a 30-min temporal resolution.TWS is a dataset (0.25 spatial resolution and 1-day temporal resolution) mainly generated from the Gravity and Recovery and Climate Experiment satellites, and it indicates the water storage as the sum of various factors including soil moisture, groundwater, surface waters, snow and ice, canopy interception and wet biomass (Li et al., 2019).Although precipitation is not one of the factors directly considered in TWS, it is a variable indirectly contributing to TWS because precipitation is essential for simulating groundwater storage changes (Li et al., 2015).
For InSAR time series analyses, we used 147 C-Band Sentinel-1 satellite SAR images acquired between May 2016 and July 2022, in the Vertical Transmit-Vertical Receive (VV) polarization channel.For the study area, Sentinel-1 ascending data are not available, and thus, we collected data in descending orbital direction, with Path 4 and Frame 491 (Table 1), which cover the area affected by the 2017 Nyingchi earthquake region (outlined by the black square in Figure 1).

| Generating surface deformations through InSAR
We applied Persistent Scatterer Interferometry (PSI; Ferretti et al., 2001) to monitor InSAR-derived HD patterns for pre-and postseismic periods.While analysing the post-seismic period, we excluded periods in which large surface deformations might have happened.
Following the main shock, our study area was affected by subsequent three 5.0 Mw earthquakes in 2019, 2020 and 2021.Therefore, we split the post-seismic period into four time-stacks (TSs) with an average length of a year (Table 1).We used earthquake occurrence dates as a reference to determine the TSs.As a result, we exclude the coseismic period that could be associated with rapid-moving landslides, which may not be captured by InSAR techniques.This not only helps us to analyse pre-and post-seismic periods separately but also to examine shorter sequential TSs where coherence between acquisitions could be preserved to generate more Persistent Scatterer (PS) points separately during these five-time intervals.
The SNAP, SNAP2StaMPS and StaMPS software packages were used for the implementation of interferometric and time series analysis (Delgado Blasco et al., 2019;Foumelis et al., 2018;Hooper et al., 2012Hooper et al., , 2018)), with atmospheric correction conducted using the Generic Atmospheric Correction Online Service (GACOS) product (Yu et al., 2018).The master images are selected per TS in a way that we have the shortest possible temporal baseline between all the available images, which also helps in maintaining the coherence between master and slave images (Crosetto et al., 2016).Then, we removed the topographic and the flattened earth phases using the SRTM DEM (30 m) and applied 3D phase unwrapping to generate PS deformation time series estimation in line of sight (LOS) direction.We assume that the effect of the atmospheric phase on the deformation measurements is minimized by applying the GACOS product.Since the study area is affected by a sequence of earthquakes and owing to the unavailability of a known stable point, the value for the reference point is taken from the average value of all PSI points present in the entire region (Hooper et al., 2012).Details of the InSAR analyses are provided in the supporting information.
T A B L E 1 Sentinel-1 data used in the analyses

| Analysing and interpreting surface deformations
The PSI technique provides deformations in LOS directions as well as average annual velocities (V LOS ).Based on the relative position of hillslopes for the heading angle of the satellite (α, 190 for our case), V LOS could exhibit positive or negative signs (Colesanti & Wasowski, 2006;Notti et al., 2014).Deformations away from a given hillslope take positive values if the hillslope faces towards the SAR sensor.On the other hand, the same situation on hillslopes facing away from the sensor is identified with negative values.Because this research aims at identifying HD, we categorized PS points with positive or negative V LOS with different aspect values and only focused on the ones implying HDs away from a given hillslope.
We followed a three-step methodology (Figure 2).First, we masked flat areas (slope < 10 , e.g., Kritikos et al., 2015) to filter out deformation that may not be associated with hillslopes.To identify the threshold value of slope steepness, we performed visual checks to ensure that the majority of floodplains with gentle slopes were removed from the analyses.Second, we categorized the PS points into classes.Class (i) contains the PS points with an aspect facing towards to sensor (i.e., (α-180, α]), and Class (ii) includes the ones with an aspect away from the sensor (i.e., (α, α-180]).Ultimately, the PS points with positive V LOS from class (i) and negative ones from class (ii) were selected and combined.In the rest of the analyses, we took their absolute values and only used these PS points.
To visualize these HD systematically through the study area, we aggregated V LOS for each hillslope following Sadhasivam (2022).
There, we used the r.slopeunits software (Alvioli et al., 2016) to delineate landscape partitions, called slope units (SUs).These are characterized by similar aspect values and are mainly bounded by ridges between adjacent hillslopes.
To examine the link between hydrologic variables (precipitation and TWS) and surface deformation, we used Cross Wavelet Transform (XWT; Grinsted et al., 2004).XWT examines time-frequency domains and identifies the corresponding sections of time series carrying large common power with a consistent phase relationship to determine the coherence between analysed datasets.To perform this analysis, we also applied spline interpolation to our deformation time series to generate equally spaced time intervals (i.e., 12 days) where we can consistently compare the deformation time series with hydrologic variables (e.g., Schlögl et al., 2021;Tomás et al., 2016).
Between the five stacks, TS2 represents the one right after the 2017 M w 6.4 Nyingchi earthquake.TS2 is also the period with the lowest PS point density.This could be because of coherence loss associated with the earthquake.The first few months following a large seismic shock generally refer to a period where landsliding rates are still higher compared to pre-seismic conditions (Tanyaş, Kirschbaum, Görüm, et al., 2021).Therefore, the coherence loss associated with large postseismic hillslope failures, which cannot be captured using the PSI technique, could be the reason for low-density PS point distribution.

| Aggregating HDs for SUs
The row PS points and associated V LOS values presented in Figure 3 were filtered out to extract only the HD away from hillslopes (see Overall, deformations in TS3 are still higher than in TS1 but with a slight decrease compared to TS2 (Figure 4g).From TS3 onwards, the elevated HD noticed in TS2 begins to fade away.And, although deformations are overall still above the pre-seismic conditions, they are not as high as TS2.
Also, the spatial distribution of HD only matches with the epicentral locations of the 2017 M w 6.4 Nyingchi earthquake and its nearby aftershocks (Figure 4b).However, from TS3 to TS5, other earthquakes of M w 5.0 (see the north-western corner of the study area in Figure 4) do not show any ground motion pattern that could be visibly linked to that of the estimated HD.We should stress that this is not conclusive evidence given the limited observation at the border of the study area.On the contrary, relatively high deformations still concentrate around the epicentre of the Nyingchi earthquake (Figure 4c-e).This shows that from TS2 to TS5, the overall post-seismic HD is mainly derived from the legacy effect of the 2017 Nyingchi earthquake.
This information is summarized in Figure 5, where the mean average V LOS for TS1 to TS5 is represented by boxplots.It shows a statistically significant increase in HD right after the main earthquake followed by a decreasing trend until TS5.For instance, the respective HD median values are 11.7, 20.4,15.9, 14.2 and 14.0 for TS1 to TS5.
Note that this suggests that post-seismic HD is still higher than the pre-seismic level 4.5 years after the event.

| Time series analyses
We also examined post-seismic deformation patterns using deformation time series instead of V LOS .To carry out this analysis for the whole study area, we calculated mean deformation values and compared them with respect to precipitation and TWS.The visual comparison shows that TWS is high in wet seasons as one can expect (Figure 6a,b).Because the study area receives most of its precipitation in these seasons, TWS is high in these periods, and deformations appear to be high as well.This is to highlight that hillslopes demonstrate a seasonal deformation pattern in TS1, TS3, TS4 and TS5.Nevertheless, right after the 2017 Nyingchi earthquake in TS2, the correlation between both hydrologic proxies and HD seems to be negative (Figure 6a,b).
To numerically examine the correlation, we used XWTs of precipitation and HD (Figure 6c) as well as TWS and HD (Figure 6d).Results show a significant coherence between the time series (i.e., common features in the wavelet power of the two-time series) in both cases with approximately 1 year period (see black polygons in Figure 6c,d

| DISCUSSION
Our results show that average V LOS values are still higher than the pre-seismic period, approximately four and a half years after the earthquake.This implies that the earthquake legacy effect has still been playing a role in surface deformations (see Figure 5).However, it the apparent recovery of the HD to pre-seismic levels could be dependent on a lower water content rather than on processes of hillslope recovery.
Testing this hypothesis requires decoupling the water content signal from the HD time series because minimizing the hydrologic contribution could provide a better and less biassed insight into post-seismic hillslope evolution processes.To accomplish this task, we used TWS, which is a variable representing the delayed effect of precipitation on hillslopes.We initially standardized both TWS and HD time series subtracting their respective means and dividing each time point by their respective standard deviation.This procedure, commonly referred to as mean-zero, unit-variance ensures that both time series are rescaled to the same unitless range (Figure 7a).
In a second step, acknowledging that HD should be better analysed by minimizing the influence of the TWS signal, we normalized once more both time series, anchoring their respective intra-seismic distributions to the pre-seismic one (Figure 7b).In other words, we calculated the mean TWS values for each TS and shifted the TWS time series in TS2-TS5 to the level shown in TS1.We also applied the same shifts in the corresponding HD time series, for TS2-TS5 (Figure 7a).Results show the HD time series normalized as a function of TWS, for each TS (Figure 7b).With this normalization, we assume that we minimized the contribution of TWS as a factor causing surface deformations.By decreasing the contribution of the hillslope water storage, we observe that HDs are still higher than the pre-seismic level, something we observed in the raw HD data plotted in Figure 5. Therefore, our analyses confirmed that postseismic deformations are higher than their pre-seismic counterparts regardless of hydrologic conditions.This also implies that the earthquake legacy effect still influences HD approximately four and a half years after the earthquake.
However, we should stress that the hillslope recovery time we mention here is not comparable to the landslide recovery time often discussed in the literature where authors exploited discrete landslide information over time to assess it (e.g., Fan et al., 2018;Kincey et al., 2021;Tang et al., 2016).In this context, examples of landslide recovery could be plausible for periods longer than four and half years only for a few cases corresponding to quite strong earthquakes, such as 1999 M w 7.7 Chi-Chi, 2005 M w 7.6 Kashmir, 2008 M w 7.9 Wenchuan or 2015 M w 7.8 Gorkha (Tanyaş, Kirschbaum, Görüm, et al., 2021).However, this consideration cannot be easily generalized for the number of observations we can rely upon is very limited.For instance, even among the very few long-term studies, significant differences do exist.In some exceptional cases, authors hypothesize persisting earthquake legacies over multiple decades whereas in other situations this is confined within a few years.Parker et al. (2015), for instance, argue that the legacy effect of the 1929 M w 7.7 Buller earthquake still persists in the landslide distribution associated with the 1968 M w 7.1 Inangahua earthquake.

| CONCLUSIONS
This research focuses on the evolution of post-seismic HD, a concept mostly studied using exclusively the visible deformation recorded in multi-temporal landslide inventories or some targeted slow-moving landslides.The novel contribution in this work is the use of InSAR-derived HD to study earthquake legacies.To do so, we applied PSI to monitor InSAR-derived HDs in the area affected by the 2017 M w 6.4 Nyingchi, China earthquake.We also analysed hydrologic variables (i.e., precipitation and terrestrial water storage) influencing surface deformations.Our findings indicated a strong correlation between these hydrologic proxies and HDs in pre-and post-seismic periods.Although this correlation was disturbed for approximately 2 years following the earthquake, a seasonal deformation pattern, which was also observed in the pre-seismic phase, appeared again after this 2-year long gap.We also demonstrated that the average post-seismic HD is still higher than its pre-seismic counterpart approximately four and a half years after the earthquake.Aside from the specific example at hand, the most generic argument is that higher HD should be associated with higher landslide susceptibility/hazard.In this context, the landslide recovery time identified from landslide inventories might not reflect the actual postseismic hillslope conditions because not all the ground motion disturbance translates into an actual failure.Slopes that were on the brink of instability certainly may have a higher chance of failure, but it is also true that slopes that were previously stable may be brought close to failure without it manifesting.This is the reason why exclusively focusing on landslide inventories for the estimation of recovery times may largely underestimate or at least provide strongly biassed information related to the earthquake legacy effect on hillslope stability.
This discussion points out some further research questions that still need to be addressed.Specifically, the link between hillslope recovery and landslide susceptibility/hazard should require further analysis to numerically express how the variation in hillslope recovery influences landslide susceptibility level.Further, a more robust identification of post-seismic hillslope strength could also help us improve landslide susceptibility and hazard assessment after strong earthquakes.

AUTHOR CONTRIBUTIONS
Hakan Tanyas conceptualized the research idea.InSAR Section 3).After aggregating absolute values of V LOS for SUs, we visualized those HDs, which are colour-coded as a function of the mean average annual V LOS in Figure4a-e.Post-seismic evolution of F I G U R E 2 Schematic drawing showing the methodology applied to identify HD away from hillslopes.[Color figure can be viewed at wileyonlinelibrary.com]hillslope displacements shows a rapid increase in deformations right after the 2017 Nyingchi earthquake (TS2), with mean average annual V LOS reaching up to 50 mm/year.Figure4falso shows the same high deformation rates over the entirety of the study area.In the majority of SU, the variation in HD shows an increasing trend from TS1 to TS2.
), though TS2 shows lack of coherence as well as slightly different period and phase.Arrows indicate the phase difference (i.e., the lag time) between time series.Overall, arrows pointing to the right indicate a positive correlation and arrows pointing left represent a negative correlation.XWT of precipitation and HD shows arrows pointing out upper right: $45 from the horizontal axis in the pre-seismic period, which indicates approximately 45 days (1/8 of a cycle) lagtime between precipitation and HD (Figure6c).Also, arrows pointing out right in the post-seismic periods (i.e., TS3-TS5) show a strong inphase correlation indicating that precipitation and HD are coincidental in the post-seismic periods (Figure6c).A similar observation is also valid for the XWT of TWS and HD (Figure6d).From pre-seismic to post-seismic periods, the lag time between TWS and HD gradually decreases, and in TS5, the two time series become almost coincidental in time.Notably, precipitation gives a relatively shorter lag time compared to TWS because the former refers to a process feeding the latter with a time lag.Results indicate a strong correlation with varying time lagsbetween the two-time series except for TS2.There, HD might still be connected to hydrologic variables, but they do not appear as the ones F I G U R E 3 Spatial distribution of PS points from pre-seismic (a) to post-seismic (b-e) periods.[Color figure can be viewed at wileyonlinelibrary.com]dominating the overall deformation.Also, some local variations, for instance, the water table fluctuation, might not be well represented in this global dataset of TWS.Notably, strong earthquakes could cause changes in groundwater level and its recovery may take several hours to several months depending on tectonic and lithological conditions(Liu et al., 2018).

F
I G U R E 4 Evolution of HD from pre-to post-seismic periods.Panels from (a) to (e) show mean average V LOS from TS1 to TS5, respectively.Panels (f)-(i) indicate the difference in V LOS between each successive time stuck, namely Δ [V LOS ]. [Color figure can be viewed at wileyonlinelibrary.com] also appears evident that V LOS values are following a decreasing trend.When we focus on the deformation time series, the pre-seismic deformation level seems to be reached in TS3, approximately 2 years after the Nyingchi earthquake (see Figure 6a).There, a discussion point should be raised to further investigate whether this recovery in hillslopes is fully associated with the earthquake legacy effect or if some other external factors may play a role in this recovery.For instance, if the hillslopes received less precipitation in the postseismic periods compared to their pre-seismic counterparts, this might have also caused relatively smaller HD in the post-seismic period.The decreasing HD trend shown from TS3 onwards largely matched the decreasing trend in TWS during the same period.For this reason, F I G U R E 5 Evolution of mean average V LOS from TS1 to TS5.Outliers were removed from the analyses.[Color figure can be viewed at wileyonlinelibrary.com] F I G U R E 6 Panels showing variation in hillslope deformation (HD) time series in relation to (a) precipitation and (b) total terrestrial water storage (TWS).Panels (c) and (d) show cross wavelet transforms (XWT) of precipitation and HD as well as TWS and HD, respectively.In panels (a) and (b), the grey-shaded areas indicate the range of pre-seismic HD.In panels (c) and (d), colour pallets from blue to yellow indicate increasing similarities between common patterns in the examined time series.The 5% significance level against red noise is indicated by the black contour lines, and the cone of influence where edge effects might distort the picture is shown as a lighter shade (Grinsted et al., 2004).[Color figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 7 Time series of total terrestrial water storage (TWS) and mean hillslope deformations (HD) after (a) mean-zero and unit-variance normalization and (b) re-normalization of HD with respect to TWS in pre-seismic phase.The grey-shaded area indicates the range of pre-seismic HD.Blue dashed lines indicate mean TWS values for each time stack.[Color figure can be viewed at wileyonlinelibrary.com] analyses were carried out by Kun He, Ling Chang, Nitheshnirmal Sadhasivam and Hakan Tanyas.Hakan Tanyas, Luigi Lombardo and Kun He wrote the manuscript, whereas Ling Chang, Nitheshnirmal Sadhasivam, Xiewen Hu, Zhice Fang, Islam Fadel, Ashok Dahal and Gang Luo provided feedback on it.