A new method to detect changes in displacement rates of slow-moving landslides using InSAR time series

Slow-moving landslides move downslope at velocities that range from mm year−1 to m year−1. Such deformations can be measured using satellite-based synthetic aperture radar interferometry (InSAR). We developed a new method to systematically detect and quantify accelerations and decelerations of slowly deforming areas using InSAR displacement time series. The displacement time series are filtered using an outlier detector and subsequently piecewise linear functions are fitted to identify changes in the displacement rate (i.e., accelerations or decelerations). Grouped accelerations and decelerations are inventoried as indicators of potential unstable areas. We tested and refined our new method using a high-quality dataset from the Mud Creek landslide, CA, USA. Our method detects accelerations and decelerations that coincide with those previously detected by manual examination. Second, we tested our method in the region around the Mazar dam and reservoir in Southeast Ecuador, where the time series data were of considerably lower quality. We detected accelerations and decelerations occurring during the entire study period near and upslope of the reservoir. Application of our method results in a wealth of information on the dynamics of the surface displacement of hillslopes and provides an objective way to identify changes in displacement rates. The displacement rates, their spatial variation, and the timing of accelerations and decelerations can be used to study the physical behavior of a slow-moving slope or for regional hazard assessment by linking the timing of changes in displacement rates to landslide causal and triggering factors.


Introduction
Landslides are a major natural hazard that occur worldwide and cause high economic losses and a large number of fatalities annually (Guzzetti et al. 2003;Petley 2012;Papathoma-Köhle et al. 2015;Froude and Petley 2018). To reduce and understand the effects of these hazards, many investigations have been focused on causes, triggers, and predictions of rapid or catastrophic landslides (Xu et al. 2017;Bogaard and Greco 2018;Monsieurs et al. 2019;Wang et al. 2020). These landslides have velocities that range from m day −1 to m s −1 (Cruden and Varnes 1996;Hungr et al. 2014;Lacroix et al. 2020b) and are characterized by the rapid downward motion of slope material. On the other hand, slow-moving landslides move downslope at velocities that range from mm year −1 to m year −1 for months to hundreds of years (Cruden and Varnes 1996;Lacroix et al. 2020b). Investigations of slow-moving landslides include landslide inventories (Lu et al. 2012;Borrelli et al. 2018;Zhang et al. 2018), long-term monitoring (Parise et al. 2003;Macfarlane 2009;Kavoura et al. 2020), and case studies that define, characterize, and determine the conditions of slope movement (Tomás et al. 2016;Kang et al. 2017;Bounab et al. 2021;Dille et al. 2021;Li et al. 2021;Jacquemart and Tiampo 2021).
Slow-moving landslides are usually present in areas with mechanically weak soils, with intermittent layers of clay between weathered sedimentary or metamorphic rocks (Mainsant et al. 2012;Schulz et al. 2018). Their constant movement over the years make them a key factor that influences the evolution of landscapes in mountainous areas (Roering et al. 2009;Mackey and Roering 2011;Simoni et al. 2013). Although they are not characterized by high velocities or catastrophic failure, their occurrence affects infrastructure (Soto et al. 2017;Dille et al. 2019;Nappo et al. 2019) and agriculture (Lacroix et al. 2020a;Garcia-Chevesich et al. 2021), and alters the normal development of communities. Slow-moving landslides often display a seasonal behavior, which means that their rate of movement is influenced by an external factor such as heavy rainfall (Handwerger et al. 2015(Handwerger et al. , 2019aBayer et al. 2018;Finnegan et al. 2021).
Both remote sensing and in situ approaches can be used to monitor slow-moving landslides. Remote sensing approaches include light detection and ranging (lidar) (Mackey and Roering 2011;Pirasteh et al. 2018;Jaboyedoff and Derron 2020), interferometric synthetic aperture radar (InSAR) (Strozzi et al. 2005;Handwerger et al. 2013;Bayer et al. 2018;Dai et al. 2020), and optical remote sensing (Bennett et al. 2016;Lacroix et al. 2020a). In situ approaches include the global navigation satellite system (GNSS) Notti et al. 2020), terrestrial laser scanners (Rosser et al. 2007;Aryal et al. 2012;Booth et al. 2018;Huang et al. 2019), geophysical methods (Whiteley et al. 2019b, a), accelerometers (Bagwari et al. 2021), slope deformation sensors (Askarinejad and Springman 2017), inclinometers (Lollino et al. 2018), extensometers (Klimeš 2018), and electronic distance measurement (EDM) (Petley et al. 2005;Pecoraro et al. 2019). While the highest data quality comes from in situ measurements, these are limited to single locations within a landslide, can be challenging to install and maintain (especially in remote regions), and as a result fail to capture large-scale spatial and temporal changes in landslide behavior. Therefore, satellite-based data are a better approach to identify and monitor large regions of active slow-moving landslides (Lu et al. 2012;Bianchini et al. 2018;Del Soldato et al. 2019;van Natijne et al. 2020).
Satellite-based InSAR data have been used to monitor slowmoving landslides around the world for several decades. Long-term monitoring provides the opportunity to use cumulative displacement time series to detect changes in the motion of a landslide (Cigna et al. 2012;Berti et al. 2013;Raspini et al. 2018) or to detect landslides over broad areas (Bordoni et al. 2018). To better understand landslide processes, these prior studies have focused on displacement time series evaluation, in comparison to external triggering factors such as rainfall. Our new approach similarly focuses on long-term kinematic changes and, in addition, expands on prior work by incorporating the spatial variation and timing of such changes at a regional scale. This is a step forward for the regional evaluation of landslides, in particular, more broadly constraining landslide dynamics, physical behaviors, and trigger responses.
Here, we present a new method to detect, quantify, and inventory changes in the surface displacement rate of slowly deforming areas, such as landslides, across regional scales. Our method uses InSAR displacement time series to identify slowly deforming areas and detect the moment that a deforming area begins to accelerate or decelerate. All identified accelerations and decelerations are analyzed and inventoried to determine the timing and location of changes in the displacement rate of potential unstable areas. We first test and refine our new approach as a proof of concept at the well-studied and analyzed Mud Creek landslide, located on the Big Sur coast, CA, USA. Then, we apply our method to a regular, unscreened data set along a reservoir upstream of the Mazar Dam, Ecuador.

InSAR processing for Big Sur, CA
For the California case study, we examined published InSAR time series from Handwerger et al. (2019b). This time series was made using data acquired by the Copernicus Sentinel-1 A/B (S1) satellites. These data are freely available and are provided by the European Space Agency (Desnos et al. 2014). The S1 satellites operate with a C-band (5.6 cm) radar wavelength and acquire data with a minimum acquisition interval of 6 days at a given location. Data are collected in both ascending (flying north and looking east) and descending (flying south and looking west) flight geometries. Handwerger et al. (2019b) processed data from descending track 42 between March 2015 and May 2017 and applied corrections to their InSAR data by using a scalable deformation model to reduce and correct unwrapping errors. They also manually removed poorquality interferograms prior to performing the time series inversion. These two steps were important for creating a high-quality InSAR dataset that was used to reveal complex landslide motions. Yet, this type of data correction is time-consuming and challenging and is therefore infeasible for regional investigations that may consist of tens to hundreds of landslides and hundreds or thousands of interferograms. For the full InSAR processing details and analyses, please see Handwerger et al. (2019b).

InSAR processing for Southeast Ecuador
To identify and monitor active landslides near the Mazar Dam, Ecuador, we constructed differential interferograms from InSAR data collected by S1 satellites. We processed the S1 data acquired in the Interferometric Wide (IW) swath mode, which has a 250km wide swath and a pixel spacing size of ~ 2.3 m in the looking (i.e., range) direction and ~ 15.6 m in the along-flight (i.e., azimuth) direction. We processed 966 interferograms using the Jet Propulsion Laboratory (JPL) InSAR Scientific Computing Environment (ISCE) software package (Rosen et al. 2012). Our processing strategy was to construct interferogram pairs with two nearest neighbors. We processed 495 interferograms on ascending track 18 (T18A) and 471 interferograms on descending track 142 (T142D). All of the interferogram pairs used in this study are listed in Online Resources (ESM) 1 and 2. To geocode the data and remove topographic phase contributions, we incorporated a ~ 30-m DEM from the Shuttle Radar Topography Mission (SRTM) into our processing (Farr et al. 2007). To reduce noise, we multi-looked the interferograms by taking 9 looks in the range direction and 2 looks in the along-flight direction and applied a standard power spectral filter with a value of 0.5 (Goldstein and Werner 1998).
Finally, we quantified the time-dependent behavior of active landslides by constructing time series with the open-access Miami InSAR time-series software in Python (MintPy) (Yunjun et al. 2019). More specifically, we used the Small Baseline Subset (SBAS) technique (Berardino et al. 2002) weighted by the inverse of phase covariance (Tough et al. 1995;Guarnieri and Tebaldini 2008;Yunjun et al. 2019). We applied a coherence threshold mask and dropped noisy pixels with coherence less than 0.4. We also corrected for tropospheric delay using ERA5 data from the European Center for Medium-Range Weather Forecasts (ECMWF) (Jolivet et al. 2011(Jolivet et al. , 2014. To reduce long-wavelength noise, we selected a local stable reference point near the active landslides. The additional InSAR processing steps (i.e., unwrapping error corrections) performed by Handwerger et al. (2019b) for the case of the Mud Creek landslide were not implemented for the Ecuador case study because our goal was to develop and test a InSAR processing strategy that does not require individual corrections such that it can be applied to regional landslide detection and monitoring. The final result is a time series of cumulative displacements measured along the satellite line-of-sight (LOS) for each pixel.

Methodology for the detection of accelerations and decelerations
We developed a method to identify and quantify changes in the displacement rate over time of slowly deforming areas by evaluating the InSAR cumulative deformation time series of each pixel in our study areas. We assume that the deforming areas identified with InSAR are slope movement, but they could also be related to deforming structures in the area. Ultimately, we do not expect InSAR signals based on anything else than surface displacement. Our method consists of four steps (Fig. 1). First, we select pixels in the area of interest that show a significant movement (the "Pixel selection" section). Second, we perform outlier detection on the time series of each selected pixel (the "Outlier detection" section). Third, we fit a piecewise-linear function model to each selected cumulative displacement time series to identify accelerations and decelerations (the "Model fitting, evaluation, and selection" section). Fourth, we perform a spatial analysis on the detected accelerations and decelerations by identifying neighboring pixels with similar accelerations and decelerations (the "Spatial analysis: detection of accelerations and decelerations" section). Our final result is a monthly inventory of the change points in the displacement time series corresponding to accelerations and decelerations. This information can be used to identify and monitor active slowmoving landslides and other localized ground deformations.

Pixel selection
The pixel selection is performed by analyzing the InSAR data in the area of interest. To identify the areas that most likely represent slope movement, we define our selection criterion based on the magnitude of the movement of pixels. We select pixels that exhibit the largest magnitude of displacement, above a specified percentile, for further analysis. It is recommended that the selected threshold percentile includes the largest displacement magnitude pixels without including noisy pixels (i.e., isolated pixels that are not likely representing slope movement). A lower threshold may include such isolated pixels, while a larger threshold would avoid noisy pixels, but it would also exclude pixels that belong to potential deforming areas. A preliminary inspection of the data results in a balanced threshold selection. The selected pixels for the case studies in California and Ecuador are presented in the "Pixel selection, case 1" and "Pixel selection, case 2" sections, respectively.

Outlier detection
The cumulative displacement time series from the InSAR data may contain outliers, and we use the Hampel method (Pearson 2005(Pearson , 2011 to identify them. The Hampel method uses a sliding window that scans the data and identifies an outlier when a data point differs from the median in the sliding window by a specified number of standard deviations. The sliding window is the window size (W) on each side of the evaluated point, so the total sliding window size is 2 W + 1. The value of the window size is based on the temporal sampling of the InSAR data. Datasets with a lower temporal sampling require a smaller window size to avoid including more than one season in the sliding window. A higher temporal sampling allows a larger window size. A lower number of standard deviations result in a stricter filter and in the identification of more outliers, while a higher number of standard deviations result in a coarser filter and fewer outliers. All identified outliers are removed from the time series. Outlier detection is applied to the case studies in California and Ecuador in the "Outlier detection, case 1" and "Outlier detection, case 2" sections, respectively.

Model fitting, evaluation, and selection
After the outliers are identified and removed, a piecewise linear function is fitted to each cumulative displacement time series. A piecewise linear function consists of a number of straight segments where the slope of each segment represents a period of movement at a constant velocity. The specific time at which the slope (i.e., velocity) of the linear segment changes is called a breakpoint, and these breakpoints represent the timing of a change in velocity resulting from an acceleration or deceleration. We apply the PWLF Python package (Jekel and Venter 2019), which was developed to fit continuous piecewise linear functions, provided that the number of breakpoints is specified.
We fit multiple piecewise linear function models to each time series. Each model has m breakpoints, where m ranges from 1 to the maximum number of breakpoints. Model m has 2 m + 2 parameters: m breakpoints, m + 1 slopes, and an intercept. The maximum number of breakpoints may be set based on the timespan of the dataset and the expected maximum number of accelerations and decelerations in a given time frame. For example, in regions where landslides have documented seasonal velocity changes related to wet and dry seasons, we expect two breakpoints (1 acceleration and 1 deceleration) per year (e.g., Handwerger et al. 2019b;Bayer et al. 2018).
The InSAR displacement time series indicate deformation in the LOS direction, and can be positive or negative (depending on the direction of motion relative to the satellite look direction). In order to simplify the analysis, we converted the negative LOS values to positive LOS values since our objective is to detect changes in the time series. Yet, the piecewise linear fit may return sections with negative slopes that correspond to motion with a LOS direction that is opposite a landslides' downslope motion. For landslides, we assume that they are always moving in the same downslope direction within the period of study and there is no obvious physical explanation as to how the sign of the LOS can switch from positive to negative (or vice versa) during a short sliding period. Therefore, we assume that negative slopes are a result of inversion or unwrapping errors and remove linear fits with negative slopes unless it is the first or last segment. In this latter case, slopes with LOS opposite of the downslope direction are likely an artifact of a limited number of data points when breakpoints occur near the beginning or end of the time series.
Next, each estimated breakpoint is evaluated using two criteria: the uncertainty of the breakpoint, referred to as the breakpoint criterion, and the estimated confidence intervals of the slopes of the segments on the two sides of each breakpoint, referred to as the slope criterion. The breakpoint criterion is assessed by evaluating the standard error of the estimated timing of the breakpoint. The standard error of the estimated breakpoint must be lower than a predefined threshold, which is set based on the temporal sampling of the InSAR data. The slope criterion considers the confidence intervals of the slopes of two consecutive segments in a model. We estimated the 95% confidence interval of a particular slope as ± 1.96 times the standard error of the estimated slope. A change in slope is considered significant when the confidence intervals of two consecutive slopes do not overlap.
All models that meet both the breakpoint criterion and the slope criterion are further evaluated using the Akaike information criterion (AIC) to determine the optimal number of breakpoints. The AIC criterion is used to assess the overall fit of a model and penalizes for the number of estimated parameters, which prevents overfitting (Burnham and Anderson 2002). The AIC is computed as follows: where SSR is the sum of squared residuals, n is the number of data points, and k is the number of parameters (Burnham and Anderson 2002). A smaller AIC indicates a better model fit. The model with the optimal number of breakpoints has the lowest AIC. Model fitting, evaluation, and selection are applied to the case studies in California and Ecuador in the "Model fitting, evaluation, and selection, case 1" and "Model fitting, evaluation, and selection, case 2" sections, respectively.

Spatial analysis: detection of accelerations and decelerations
After the breakpoints (i.e., points of acceleration/deceleration) are identified, a monthly spatial analysis is performed. Because of the uncertainty in the estimated timing of a breakpoint, a breakpoint is partly counted in the month of the estimated timing and partly in the months before and after. The distribution across the 3 months is based on the estimated standard error of the breakpoint as follows. A breakpoint is assigned to the middle of a month and the probability that it occurs in that month is estimated using a Normal distribution with the estimated standard error as the standard deviation. The remaining probability (i.e., the probability that the breakpoint does not occur on the estimated month) is distributed equally over the months before and after. As a result, the number of accelerations or decelerations in a month is not an integer and the total number of accelerations and decelerations sums to the total number of detected breakpoints. Finally, a spatial analysis is performed to identify neighboring pixels exhibiting similar time series behaviors. We used the Python implementation of the density-based spatial clustering of applications with noise (DBSCAN) algorithm (Ester et al. 1996;Schubert et al. 2017) to identify pixels that belong to clusters with similar behavior (accelerations or decelerations). The algorithm uses two parameters: the maximum distance between pixels in a cluster and the minimum number of pixels in a cluster (Boeing 2018). These two parameters are set based on the InSAR data density. More details in the "Spatial analysis: detection of accelerations and decelerations, case 1" and "Spatial analysis: detection of accelerations and decelerations, case 2" sections.
Our final product is an inventory of the timing of the changes in the displacement rate of pixels within a cluster of pixels that show similar behavior. The inventory is accompanied by multitemporal maps of grouped pixels corresponding to areas with similar behavioral patterns, likely representing slope movement. Note that in this paper, we focus on identifying and quantifying the accelerations and decelerations that we find within pixels that are part of a cluster per month. The grouped pixels are indicators of deforming areas.

Case study 1: Mud Creek landslide
We tested and refined our new method to detect the timing of accelerations and decelerations from InSAR time series on data of the Mud Creek landslide. Mud Creek was a landslide that moved slowly for at least 8 years (likely much longer) before accelerating rapidly and failing catastrophically on 20 May 2017 (Warrick et al. 2019;Handwerger et al. 2019b;Jacquemart and Tiampo 2021). This landslide had a pre-catastrophic failure area of approximately 0.23 km 2 (Handwerger et al. 2019b) and a mean slope angle of 38 degrees (Warrick et al. 2019). The landslide's bedrock geology is composed of the Franciscan mélange rock, which is characterized by a clayey granular matrix with highly sheared sandstone, siltstone, metasandstone, shale, serpentinite, and blueschist (McLaughlin et al. 1982(McLaughlin et al. , 2000. The average precipitation is around 1000 mm/year and occurs primarily between October and May. The Mud Creek landslide experienced periods of extreme drought and extreme rainfall during our study period. A historic drought lasted from 2012 to 2016, while 2017 was one of the wettest years on record. Previous work by Handwerger et al. (2019b) and Jacquemart and Tiampo (2021) used InSAR time series to detect changes in landslide motion and found that seasonal rainfall drives these changes. This landslide was selected as a proof of concept for our new method due to its high-quality time series (as described in the "InSAR processing for Big Sur, CA" section) and its documented seasonal behavior. We examined the period of slow landslide motion captured by the S1 InSAR time series between 2015 and 2017.

Pixel selection, case 1
For the Mud Creek landslide, we selected the 2% of pixels with the largest (absolute) displacement (147.4 mm) resulting in 1124 of 93,590 pixels. This 2% threshold selected pixels that are part of the deforming area (Fig. 2). In the case of Mud Creek, a higher percentile (> 2%) would include pixels that may represent noise, and not real displacement. A lower percentile (< 2%) would leave out important pixels that may be part of an unstable area. A preliminary evaluation of the pixels is advised to select the threshold that captures most pixels within moving areas without noisy pixels. All selected pixels are entirely within the boundaries of the precatastrophic polygon mapped by Handwerger et al. (2019b).

Outlier detection, case 1
We used the Hampel method to identify outliers as described in the "Outlier detection" section with a sliding window size of 7 data points (representing ~ 84 days) and 2 as the number of standard deviations. Very few outliers were detected and removed from this high-quality data set. In total, 335 outliers were found and removed from 244 time series (0.47% of all data points), with a maximum of 4 outliers in one time series. In Fig. 3, we show some examples of the application of the outlier filter to the displacement time series of the Mud Creek landslide.

Model fitting, evaluation, and selection, case 1
We identified the number of breakpoints in each filtered time series using the method described in the "Model fitting, evaluation, and selection" section. We used a threshold standard error value of 30 days for the breakpoint criterion. The maximum number of breakpoints in a time series is set to four because the displacement time series are available for 1 year and 9 months, where we can observe two rainy seasons (October to May each year). Based on the typical landslide behavior in coastal California (e.g., Handwerger et al. 2019a, b), we expect at most two accelerations and two decelerations. A total of 1120 out of the 1124 time series were fitted successfully. We identified 2967 breakpoints in total: 121 time series (i.e., pixels) with 1 breakpoint, 315 time series with 2 breakpoints, 520 time series with 3 breakpoints, and 164 time series with 4 breakpoints (see Fig. 4).

Spatial analysis: detection of accelerations and decelerations, case 1
We detected the timing of accelerations and decelerations of the fitted time series, as described in the "Spatial analysis: detection of accelerations and decelerations" section. We determined the location of the pixels of the 1120 fitted time series and we used the DBSCAN algorithm to identify and select pixels within a cluster, using 12 m as the maximum distance between pixels and 4 pixels as the minimum number of pixels that form a cluster. The maximum distance between pixels is the space between the edges of 2 pixels, and is set to 12 m, which is the spacing of the digital elevation model that was used to geocode the interferograms (Handwerger et al. 2019b). We then compiled an inventory of the total number of accelerations and decelerations detected from pixels within a cluster per month. For Mud Creek, the number of breakpoints after clustering remains the same. All pixels are part of a cluster during the studied period. In We compared our breakpoint detection inventory to the local precipitation patterns known to have controlled the behavior of Mud Creek. We found very little activity between July and October 2015 (i.e., there are no accelerations and only a few decelerations). This period is at the end of the dry season of a historic drought period. The landslide behavior changed when the 2015-2016 rainy season began and we detected many acceleration breakpoints between October 2015 and February 2016. Figure 4b shows that a large portion of the landslide was accelerating in February 2016 (rainy season). This period of acceleration was followed by a period of deceleration during the 2016 dry season. Figure 4c shows the spatial distribution of deceleration points within the landslide occurring in June 2016. Comparison with the November 2016 map shows that there are spatial differences in the timing of accelerations and decelerations within the landslide.
The landslide then started to accelerate again shortly after the onset of seasonal rainfall in the wet season of 2016-2017. We found that the largest number of accelerations recorded in a month occurred earlier (2 months after the onset of the rainy season) and was higher than the previous 2015-2016 wet season (Fig. 4a). This change in behavior was presumably driven by the large changes in rainfall that occurred during the study period and our findings agree with the detailed analysis presented by Handwerger et al. (2019b). We found that almost all pixels within the landslide are accelerating in November 2016 (Fig. 4d). Accelerations occurred during the entire wet season until the catastrophic failure in May 2017.
We also compared our inventory of accelerations and decelerations to the velocity time series of Handwerger et al. (2019b) for 36 pixels within an area of 60 × 60 m in the landslide (Fig. 2). The timing of accelerations and decelerations found with our method matches those found by Handwerger et al. (2019b). Accelerations detected with our method correspond to periods of increasing velocities, while decelerations were detected during the period of decreasing velocities (see Fig. 5). We observe that 32 out of 36 pixels accelerated from December 2016 to March 2017, 31 out of 36 pixels decelerated from May 2016 to August 2016, and 35 out of 36 pixels accelerated again from October 2016 to January 2017. Only 2 pixels decelerated in this last period. The green box indicates a representative area of 60 × 60 m used by Handwerger et al. (2019b) to derive landslide velocities

Case study 2: Mazar Region
Our second study site is the region surrounding the Mazar hydroelectric power plant and its reservoir in southeast Ecuador (Fig. 6). Here, a major hydroelectric complex extends from the Andes to the Amazonian region. The area is known to have several deep-seated landslides (e.g., Nicole 2015; Urgilez Vinueza et al. 2020). Around the reservoir, seventeen slow-moving landslides have been identified by the electrical corporation of Ecuador (CELEC). CELEC identified the unstable areas during the construction of the Mazar dam and has been monitoring them because they are a threat to sustainable hydropower generation. The lithology of the landslide area is composed of two geological units: Alao Paute and El Pan, characterized by metamorphic rocks, overlain by colluvium deposits ranging from 2 to 28 m (Nicole 2015). The precipitation in the area is around 3000 mm/ year and occurs primarily between April and August.

Pixel selection, case 2
For the Mazar region, we selected an area of 211 km 2 around the reservoir to examine the InSAR displacement time series from October 2016 to August 2020. Our pixel selection resulted in 3230 pixels with an absolute displacement value above the 98th percentile (99.9 mm). Twenty-eight percent of the selected pixels fall within the boundaries of the ground-based landslide inventory carried out by CELEC, while 72% fall outside the boundaries of the identified unstable areas.

Outlier detection, case 2
We used the Hampel method as described in the "Outlier detection" section with a sliding window size of 7 data points and 2 as the number of standard deviations. Factors such as the vegetated soil cover and atmospheric disturbances resulted in noisy displacement time series. Additionally, and contrary to the case of the Mud Creek landslide, a prior quality control of the data was not conducted for the Mazar landslides. One of the main objectives of our new detection method is to process large quantities of data at a regional spatial scale and in a relatively fast and semi-automated manner. In total, 25,860 outliers were found in 2268 time series, each time series with 1 to 24 outliers. A total of 4.7% of the data points were identified as outliers and removed from the time series. In Fig. 7, we show two examples of the application of the outlier filter to the displacement time series of the pixels in the Mazar region.

Model fitting, evaluation, and selection, case 2
We identified the breakpoints in each time series using a threshold standard error value of 30 days for the breakpoint criterion. In the case of the Mazar region, we set the maximum number of breakpoints to eight due to the expected number of accelerations and decelerations that can occur over the course of 3 years and 6 months under the influence of the rainy season (April-August). Our method was able to fit 2273 out of the 3230 time series successfully. We identified 3397 breakpoints in total: 1155 time series with 1

Spatial analysis: detection of accelerations and decelerations, case 2
We performed the spatial analysis using the pixels of the 2273 fitted time series per month, using 30 m (i.e., DEM pixel spacing) as the maximum distance between pixels and 4 pixels as the minimum number of pixels that create a cluster. In the case of the Mazar region, the number of breakpoints after clustering diminished. There are 2793 of 3397 breakpoints in total: 561 time series with 1 breakpoint, 566 time series with 2 breakpoints, 186 time series with 3 breakpoints, 90 time series with 4 breakpoints, 34 time series with 5 breakpoints, and 2 time series with 6 breakpoints. In Fig. 8, we present a spatio-temporal inventory of the number of detected acceleration and deceleration points within a cluster, and the location of the corresponding pixels for four periods. Figure 8a and b show a seasonal and yearly distribution of the monthly number of accelerations and decelerations, respectively. Figure 8a shows that the number of decelerations is higher than the number of accelerations in the first months of the year. There is a modest increase in the number of accelerations and a modest decrease in the number of decelerations once the wet season starts. However, we observe that by the end of the wet season and after, both accelerations and decelerations occur. Figure 8b indicates that the number of accelerations increases throughout the study period, while the number of decelerations decreases. In the year 2020, the number of accelerations is higher than the number of decelerations. Figure 8c-f shows the location of pixels with a mild to high probability of occurrence of accelerations and decelerations.
By examining the spatial variability of accelerations and decelerations over the area around the Mazar reservoir, we find that most of the activity occurs on the south side of the reservoir, where two reservoir tributaries meet. Some activity is observed near the dam, on the north side, and on the central-east side of the reservoir. Our inventory reveals that accelerations and decelerations occur throughout the year and are sparse around the reservoir. These are concentrated in specific locations at the end of the study period. This variability of accelerations and decelerations occurs within groups of pixels as well as between groups of pixels.

Discussion
In this study, we developed a systematic method to detect, quantify, and inventory changes in the surface deformation rate of slowly deforming areas at a local and regional scale to investigate their temporal and spatial dynamics. Slow-moving landslides have been studied prior to this work using satellite data to identify ground motion areas and shifts in the displacement time series (Cigna et al. 2012;Berti et al. 2013;Raspini et al. 2019). Our method differs from previous work in that our InSAR detection analysis provides an objective way to construct multitemporal maps of unstable areas and an inventory of the timing of changes in deformation rate of unstable areas.
Due to the nature of the InSAR data, the time series we used for the analysis often contained outliers. These outliers are usually related to specific data acquisitions in the time series and impede the fitting of the piecewise linear functions. For the Mud Creek landslide, there were very few outliers because analyses of a single landslide allow for more in-depth quality control measures. For Mazar, there were many outliers because we did not carefully inspect individual interferograms or perform unwrapping error corrections. This was intentional as one of our main goals is to develop a method to analyze large quantities of regional slope deformation data where it is infeasible to inspect and/or make corrections to thousands of interferograms. Therefore, we opted for the Hampel method as a filter routine to identify and remove the outliers, while at the same time, avoiding the exclusion of important data.
The optimization of the Hampel parameters was carried out considering the temporal sampling of the InSAR data. Smaller window sizes will not detect the short-term outliers, while larger window sizes fail to identify outliers within a small portion of the window due to a higher median value. In our case, we used a window size of seven data points, representing a period of ~ 90 days and the number of standard deviations was set to two. Moreover, we found that using a standard deviation smaller than 2 tends to over-smooth the time series, while using a standard deviation over 3 did not identify all outliers.
We evaluated the uncertainty of the timing of the breakpoints using the breakpoint criterion and decided on a threshold standard error value of 30 days. A smaller threshold value leads to a stricter algorithm, so that very few breakpoints are accepted. On the contrary, a higher threshold value allows more breakpoints to be identified, but then, the time frame when accelerations and decelerations occur becomes too large, and becomes meaningless with respect to slow-moving landslide dynamics. Any threshold that yields a time frame larger than the duration of a (wet) season will not give useful information about the temporal dynamics of the slow-moving slopes. Therefore, we used a threshold standard error of 30 days and a monthly temporal resolution.
We selected our pixels considering the 98th percentile of the absolute cumulative displacement of the InSAR data. Some of the selected pixels that comply with this condition were isolated and did not have an immediate neighboring pixel that also showed significant displacement. We assumed that isolated pixels do not Fig. 7 Examples of displacement time series and the identified outliers for pixels P18 (11 outliers) and P1410 (12 outliers) in a and b and the fitted breakpoints of P18 (3 breakpoints) and P1410 (4 break-points) in c and d in the Mazar region. Pink points are the initial and final data points of the time series and are not considered as breakpoints correspond to landslides. Several studies have shown that clusters of pixels with relatively high LOS displacement can be used to identify active landslides (e.g., Bekaert et al. 2020). Therefore, we selected pixels that belong to a cluster that showed activity (acceleration or deceleration) in the same month. We followed this approach in order to achieve spatial consistency and temporal persistency (pixels with a significant change in displacement rate) .
Our method identified breakpoints that show clear changes in deformation velocity that can be related to seasonal rainfall. In the Mud Creek landslide, accelerations take place during the rainy season, while decelerations occur during the dry season, as has been shown by hundreds of landslides in coastal California (Handwerger et al. 2019a). Previous work on the Mud Creek landslide showed that the slope dynamics are directly related to large changes in seasonal rainfall (Warrick et al. 2019;Handwerger et al. 2019b;Jacquemart and Tiampo 2021). Our change detection approach captured the seasonal kinematics of Mud Creek and allowed us to explore spatial trends and accelerations and decelerations by fully utilizing the rich information provided by InSAR. During the time period between Feb and May 2017, the landslide was likely moving faster than InSAR is able to detect, as the landslide approached catastrophic collapse. This causes phase-bias, an additional unwrapping challenge, that obscures the true deformation rate and was not possible to manually correct. However, it is encouraging that our method detected the overall change in the signal of the data, dominated by accelerations, and did not capture the small apparent deceleration in the months prior to failure (Fig. 5). For the Mazar region, we observed more complex behavior resulting from analyses of numerous large and spatially variable deforming areas. We found that the accelerations and decelerations occurred during the entire period of study and are distributed over the area around the reservoir. In Fig. 8a and b, we showed that both accelerations and decelerations occurred during the entire period, and that in 2020, the number of accelerations was higher than the number of decelerations. The high number of unstable areas that were identified using our method may have caused this somewhat less predictable behavior, which can be related to lagged responses of deep-seated landslides in the area, as well as to creep behavior of more surficial landslides. Local factors such as slope, distance to the reservoir, specific land use, irrigation, and local geomorphology may influence these different behaviors and the occurrence of accelerations and decelerations at different times. This behavior is also captured at Mud Creek, where even a single landslide can show spatial variation on the timing of accelerations and decelerations. The overall behavior of the Mazar region is complex and needs further in-depth analysis, such as a hydro-meteorological and geotechnical analysis of the larger Mazar region, which is out of the scope of this paper and is part of our next follow up study.

Conclusions
We developed an objective and systematic method for the detection of accelerations and decelerations of slowly deforming areas from InSAR data. Our method consists of InSAR time series analyses corresponding to the selected pixels (with the highest cumulative displacement). These time series are filtered, and breakpoints are detected using piecewise linear functions fitted to the time series. These breakpoints represent the times when the displacement rate changes significantly. We analyzed the spatial distribution of the successfully modeled pixels and inventoried the accelerations and decelerations that showed similar spatial behavior.
We tested our method on the high-quality InSAR dataset of the Mud Creek landslide, CA. Our method successfully detected the timing of accelerations and decelerations at Mud Creek that were driven by changes in seasonal rainfall. Next, we investigated a landslide prone area impacting a hydropower area in southeast Ecuador. Although the time series data were of significantly lower quality (compared to Mud Creek), we identified deforming areas with complex patterns of accelerations and decelerations, both within and between groups of pixels that did not always coincide with wet and dry seasons.
We conclude that our method is able to identify changes in the ground surface displacement rate of deforming areas that can be used to examine this behavior, and inventory these changes in an objective and straightforward manner. The ability to determine the temporal and spatial variation of velocity changes is a step forward in the large-scale interpretation of the physical behavior of slowmoving deforming areas. Ultimately, our inventory of accelerations and decelerations can be used as a tool to shed light on the dynamics of slow-moving landslides at both sub-landslide and regional scales with high spatial and temporal resolution.