Detection of slow-moving landslides through automated monitoring of surface deformation using Sentinel-2 satellite imagery

Landslides are one of the most damaging natural hazards and have killed tens of thousands of people around the world over the past decade. Slow-moving landslides, with surface velocities on the order of 10 (cid:1) 2 – 10 2 m a (cid:1) 1 , can damage buildings and infrastructure and be precursors to catastrophic collapses. However, due to their slow rates of deformation and at times subtle geomorphic signatures, they are often overlooked in local and large-scale hazard inventories. Here, we present a remote-sensing workflow to automatically map slow-moving landslides using feature tracking of freely and globally available optical satellite imagery. We evaluate this proof-of-concept workflow through three case studies from different environments: the extensively instrumented Slumgullion landslide in the United States, an unstable lateral moraine in Chilean Patagonia and a high-relief landscape in central Nepal. This workflow is able to delineate known landslides and identify previously unknown areas of hillslope deformation, which we consider as candidate slow-moving land-slides. Improved mapping of the spatial distribution, character and surface displacement rates of slow-moving landslides will improve our understanding of their role in the multi-hazard chain and their sensitivity to climatic changes and can direct future detailed localised investigations into their dynamics.

Landslides have killed over 50 000 people in the 21st century alone and continue to cause several thousand deaths each year (Froude & Petley, 2018;Petley, 2012).Monitoring and mitigating the impacts of landslide hazard is challenging due to the extent of potentially susceptible areas-virtually any mountainous region might be predisposed to landsliding-and the wide range of potential triggering mechanisms.
Landslides can be triggered by climatic drivers such as rainfall (Caine, 1980;Carey et al., 2019), earthquakes (e.g., Fan et al., 2018;Harp & Jibson, 1996;Roback et al., 2018) and anthropogenic activity such as road construction (Amatya et al., 2019;Petley et al., 2007;Tanyaş et al., 2022).Due to their wide range of predisposing and triggering factors, developing landslide inventories is key to improving our understanding of the underlying processes and for defining landslide susceptibility, hazard and ultimately risk.
Each remote-sensing and ground-based method is associated with different spatiotemporal resolutions, uncertainties and intrinsic limitations.To date, large-scale landslide mapping has primarily been implemented through manual mapping or automated change detection in optical and radar imagery (Ambrosi et al., 2018;Burrows et al., 2019Burrows et al., , 2020;;Kincey et al., 2021;Scheip & Wegmann, 2021).
Optical change detection-based landslide mapping relies on identifying a distinctive alteration of Earth-surface properties, most commonly a reduction in NDVI as a proxy for the loss of vegetation cover (e.g., Scheip & Wegmann, 2021).This is well suited for identifying the footprints of rapid, shallow landslides, but is less successful in identifying slow-moving landslides that might not remove vegetation or significantly change the spectral character of the ground surface (Figure 1).These landslides, with typical rates of surface displacement of 10 À2 to 10 2 m/year, can cause chronic damage or destruction of infrastructure and have major impacts on land use or livelihoods.
While some slopes might only experience slow but continuous deformation, others can experience increases in velocity and rapid failure (Berg et al., 2018;Lacroix et al., 2020;Mansour et al., 2011;Palmer, 2017).In some specific settings, monitoring of landslide velocities can be used to forecast the time of failure (e.g., Manconi, 2021;Moretto et al., 2017;Petley et al., 2017), but to date, this has primarily been applied retrospectively on known landslides.Slow-moving landslides are an important but often overlooked population in landslide inventories and associated hazard assessments, although they may be detected using alternative techniques.
Because slow-moving landslides often cannot be detected through changes in surface properties on a scale observable by satellites, alternative methods are necessary to identify them (Figure 1).Three main techniques have been used: topographic change detection, interferometric synthetic aperture radar (InSAR) and feature tracking (Dai et al., 2020;Dille et al., 2021;Lacroix et al., 2022;Shugar et al., 2021;Stumpf et al., 2014;Stumpf et al., 2017;Van Wyk de Vries, Bhushan, et al., 2022).Topographic change detection can be used to identify vertical change associated with slope deformation F I G U R E 1 Range of common spatial scales and deformation velocities for slow-moving landslides.Four different example landslides are shown, along with possible displacement monitoring techniques: satellite-based interferometric synthetic aperture radar (InSAR), feature tracking and ground-based methods such as terrestrial laser scanning.Assessments of elevation change are not shown here but can also be used to monitor landslides.Velocity classification is from Cruden and Varnes (1996) and Varnes (1978).
through ground-based, airborne or spaceborne LiDAR or photogrammetry (e.g., subsidence or bulging; Bernard et al., 2021;Dai et al., 2020).InSAR, using either ground-based or spaceborne instruments, can be used to map line-of-sight displacements, providing a combination of horizontal and vertical movement depending on the viewing geometry (Ferrigno et al., 2017;Wasowski & Bovenga, 2014).
Feature tracking measures the horizontal displacement of coherent patterns between a pair of images, which can be derived from optical or radar sensors (Bickel et al., 2018).These methods have been widely used for historical and ongoing monitoring of known slow-moving landslides (Casagli et al., 2023;Delacourt et al., 2007;Ferrigno et al., 2017;Wasowski & Bovenga, 2014), but low signal-to-noise ratios complicate their use for automated detection of previously unknown landslides, particularly in steep mountainous terrain (Bekaert et al., 2020;Lacroix, 2016;Lacroix et al., 2022;Rosi et al., 2018;Van Wyk de Vries, Bhushan, et al., 2022).
The majority of studies of slow-moving landslides have focussed on known unstable slopes, using remote-sensing techniques either alone or in concert with ground monitoring to better understand the spatial and temporal changes in slope velocity (Dai et al., 2020;Dille et al., 2021;Hu et al., 2020;Lacroix et al., 2018;Van Wyk de Vries et al., 2022).Very large landslides, known as deep-seated gravitational slope deformation (DSGSD) have been inventoried on a regional scale in the European Alps (Agliardi et al., 2013;Clague & Stead, 2012;Crosta et al., 2013), although these inventories cannot always distinguish inactive DSGSDs from ones with present-day motion.Other studies have leveraged the displacement-tracking capabilities of either InSAR (Dini, Daout, et al., 2019;Dini, Manconi, & Loew, 2019;Solari et al., 2019) or optical feature tracking (Zhang et al., 2022) to inventory slow-moving landslides on a local or regional scale.Dini et al. (2019) for instance use InSAR in the Bhutan Himalaya to identify almost 700 previously unknown unstable slopes in an area with very limited field data.
Here, we develop and implement a proof-of-concept automated workflow for the detection of deformation typical of slow-moving landslides using optical satellite imagery.Through a series of case studies, we show that our landslide detection workflow is able to automatically delineate known and identify previously unknown slowmoving landslide candidates across a wide range of environments.We then discuss how slow-moving landslide identification, without the requirement for prior knowledge of landslide locations or manual processing, can assist with the construction of regional or global inventories.Finally, we discuss the value of including slow-moving landslides in multi-hazard inventories and possible interactions with other Earth-surface processes.

| METHODOLOGY AND METHODS
We develop a workflow for identifying the location and extent of slow-moving landslides without the requirement for any prior landslide location information.This workflow requires a minimum of two optical satellite images and a digital elevation model (DEM) of the location of interest, the latter being downloaded automatically through an application programming interface (API).The optical satellite image time series is first pre-processed, then used to calculate pairwise displacement maps and finally post-processed into a single median velocity map (Figure 2).The final median velocity map and DEM are used to generate three products that describe potential landslides: a binary landslide map, a map of landslide velocities and a table of attributes for each coherent of deformation or landslide (e.g., elevation and area; Figure 2).The details of this workflow are described in more detail in Sections 2.1 and 2.2, and Sections 2.3-2.5 describe three case studies used to evaluate the workflow.

| Optical feature tracking
We adapt an optical feature-tracking workflow from the ice flow speed mapping toolbox Glacier Image Velocimetry (GIV; Van Wyk de Vries & Wickert, 2021).Two key differences between glacier flow and landslide deformation are (i) in general, the displacement related to glacier flow is on the order of 10 1 -10 3 of m a À1 , while slow-moving landslide deformation is very slow to moderate at 10 À2 -10 2 m a À1 (Cruden & Varnes, 1996;Varnes, 1978); (ii) glacier surfaces can change appearance rapidly due to high flow speeds and surface ablation, while areas of slow-moving landsliding more commonly retain a constant surface character over several years, notwithstanding seasonal ground surface changes.The theoretical detection limit of feature tracking depends on the resolution of the images used, the number of individual velocity maps stacked (each obtained from a unique image pair) and the type of surface patterns.It is usually in the range of 0.1-1 m for Sentinel-2 imagery (e.g., Millan et al., 2019;Stumpf et al., 2017), corresponding to a landslide speed detection limit of 0.13-1.3m a À1 for a 9-month temporal baseline and 0.05-0.5 m a À1 for a 2-year temporal baseline.The relatively low velocity and increased persistence of surface features over time lead to three changes to the feature-tracking workflow: 1. Longer temporal baselines (time intervals) between images can be used for landslides than is generally possible for glaciers.
2. The scale of displacement on a landslide are generally less than the spatial resolution of the imagery used in tracking, so the accuracy relies upon subpixel displacement algorithms.
3. Relatively low velocities associated with landslides mean that the stacking of multiple individual velocity maps might be required to distil actual deformation from background noise.
We apply feature tracking to Band 8 (833 nm wavelength; near infrared) of Sentinel-2 L1C imagery, which has a 10 m spatial resolution (Van Wyk de Vries & Wickert, 2021).Starting with a stack of all Sentinel-2 images (2015-2023) over each study region, we filter out all images with greater than 10% scene-wide cloud cover and manually filter out remaining images with any degree of local cloud cover.
We manually edit Randolph Glacier Inventory glacier polygons to cover the full current extent of any glaciers and mask these out (Pfeffer et al., 2014).We then pre-process the imagery by applying an orientation filter to enhance local spectral contrasts (Van Wyk de Vries & Wickert, 2021) and standardize the resulting orientation images.
We construct all possible image pairs (termed 'multiple-pairwise image correlation', Stumpf et al., 2017) with a temporal baseline greater than 9 months and run the feature-tracking algorithm on these.Setting a minimum temporal baseline of 9 months increases the chance that a detectable level of displacement has occurred for typical slow-moving landslide velocities of $10 0 m a À1 .We do not apply a maximum bound on the temporal separation between images as we expect most surfaces to remain coherent over the $6 year time series used, with the possible exception of seasonal changes (e.g., vegetation loss and snow cover) or the incidence of shallow, rapid landslides.We use a single pass chip-wise cross-correlation algorithm with a window size of 16 pixels and 50% window overlap, for an equivalent ground resolution of 80 m.This parameter combination maximises the benefits of larger window sizes for background noise reduction while retaining a final resolution of less than 100 m per pixel.We choose a single pass algorithm instead of the GIV default three-tier multipass algorithm as integer pixel displacements are expected to be $0, reducing the benefit of prior integer displacement information of initial coarse-chip searches (Van Wyk de Vries & Wickert, 2021).The effective resolution of output velocity maps is 80 m, which balances the benefits of high spatial resolution with the increased noise levels of small window sizes.We linearize the chipwise cross-correlation such that it is performed across the entire image in a single operation instead of requiring a computationally expensive double loop.We also implement an ensemble of subpixel estimation methods instead of a single algorithm, which allows us to quantify the subpixel-level uncertainty and propagate this into future calculations.Modifications to the GIV cross-correlation algorithm are described in more detail in the supporting online-only information.
We post-process each individual displacement map prior to merging the stack into a single median velocity map.We first multiply by the image resolution and divide by the temporal baseline (in years) to convert displacements to velocity vector.We then filter out all pixels with a cross-correlation peak ratio-the ratio of the cross-correlation maxima to the second largest correlation peak-of less than 1.5 (Van Wyk de Vries & Wickert, 2021).Filtering by peak ratio removes false matches in areas where decorrelation has occurred between the images.We apply a local gap-filling algorithm in all areas with less than 15 empty pixels within a 5 Â 5 pixel neighbourhood, with any larger gaps left unfilled.Finally, we remove the scene-wide median velocity in the E-W and N-S directions to correct for any systematic georeferencing errors between images.We compute median velocity from the E-W and N-S velocity components separately to avoid transforming the noise.

| Landslide identification
Landslide identification directly from a median velocity map is challenging, as background noise levels are commonly of a similar magnitude to the landslide displacement signal.Therefore, direct thresholding or classification of the velocity commonly leads to many type 1 errors (false positive identifications), including physically implausible scenarios of upslope movement.To reduce some of this complexity, we include additional topographic information in our workflow.We automatically download the 30 m Copernicus DEM of the study areas using the OpenTopography API (www.opentopography.com)and use these to calculate the local slope Diagram showing the main steps in our automated slow-moving landslide detection algorithm.direction using TopoToolbox (Schwanghart & Scherler, 2014).We then calculate the angle difference between feature tracking-derived velocity direction and slope direction as projected onto a horizontal plane.
We create a binary landslide mask by thresholding the median velocity map using the 95th percentile.This is an empirically determined threshold, which we explore below through a sensitivity test in the first case study with comparison to available field datasets.The 95th percentile provides a balance between a very low fraction of false positives (necessary for application in areas where a minority of the study area is actually moving) and sufficiently high true positive rate (60% to 90% depending on the comparison used; see supplement).We then compute regions retained by this velocity mask that meet a set of further conditions: an angle between feature trackingderived velocity direction and slope direction of less than 45 , a slope greater than 5 , a peak ratio greater than 1.5 and a signal-to-noise ratio greater than 5 (Van Wyk de Vries & Wickert, 2021).The slope direction threshold excludes areas with physically implausible apparent surface motion not in the downhill direction.The slope threshold is used to exclude lakes and other low slope regions which might otherwise interfere with the results.The feature-tracking peak ratio threshold is used to exclude areas with poor coherence, such as persistent snow cover or river channels.
We then identify clusters of adjacent thresholded pixels within the remaining mask and remove clusters smaller than a minimum landslide area cut-off-set by default to 10 contiguous pixels in our workflow.The choice of pixel threshold is a trade-off between effective random noise removal and detection of small landslidesthe exact threshold choice is again explored in a sensitivity test and comparison to available ground data in the first case study.Connected component analysis allows the workflow to distinguish isolated significant pixels (likely related to random noise) from clusters of multiple significant pixels (likely related to a real signal).Following this, we bridge significant pixels that have two unconnected neighbours, buffer the resulting polygon by one pixel, filter out isolated gaps in polygons using a 3 Â 3 pixel neighbourhood majority vote and finally erode the polygon boundary by one pixel.This sequence of morphological operations smooths and closes gaps between immediately adjacent clusters of significant pixels.The removal of significant pixel clusters smaller than a threshold prior to this operation prevents merging of landslide regions with adjacent noise.We note that any area without motion is not considered a landslide according to this workflow, such that a 'landslide' with motion observed only in some sectors separated by inactive zones would be considered separate incidences.Any clusters retained following this procedure are considered candidate slow-moving landslides and their attributes (mean velocity, size and elevation range) are calculated accordingly.
In the following sections, we describe the three case studies from a range of different environments: (1) the well-studied Slumgullion landslide, USA, which provides an evaluation of the landslide delineation workflow; (2) a region around Glaciar Oriental, Chile, through which we investigate the stability of lateral moraines around a rapidly thinning glacier in an area with frequent snow cover; and (3) central Sindhupalchok district, Nepal, a densely populated area combining high relief, extreme seasonal rainfall, dense vegetation and recent seismic activity.We explore these three very different case studies to demonstrate the wide applicability of our workflow.

| Case Study 1: Slumgullion landslide
The Slumgullion landslide is a $4-km-long and $0.5-km-wide earthflow located in the San Juan Mountains, Colorado, USA (Coe et al., 2003;Parise & Guzzi, 1992).The active portion of the landslide, with an estimated volume of 2.0 Â 10 7 m 3 (Parise & Guzzi, 1992), is nested within an even larger inactive flow extending an additional 2 km downslope and damming Colorado's largest natural lake, Lake San Cristobal.The Slumgullion landslide is an ideal test case for this workflow as (i) the landslide has been extremely well studied, with several decades of remote-sensing and ground-based measurements; (ii) it is located within complex mountainous terrain with variable vegetation cover and a lake, with potential for feature-tracking errors; and (iii) the active portion of the landslide is nested within a larger inactive flow, and the boundary between the two cannot easily be distinguished without velocity measurements.We run the workflow over the region surrounding the Slumgullion landslide using 12 cloud-free Sentinel-2 L1C images ranging from October 2016 to September 2022, for a total of 49 image pairs.

| Case Study 2: Glaciar Oriental, Chile
Glaciar Oriental is a large lake-terminating glacier located at the northeastern corner of the Southern Patagonian Icefield (SPI) in Chile.Similar to much of the SPI, Glaciar Oriental has undergone rapid thinning and frontal retreat in the 21st century (Abdel Jaber et al., 2019;Foresta et al., 2018;Malz et al., 2018;Minowa et al., 2021;Willis et al., 2012).
This rapid retreat has exposed large expanses of loose, unconsolidated sediment at the glacier flanks, which have the potential to form large landslides.Ice-marginal landslides have been associated with tsunamis (Winocur et al., 2015) and changes in ice dynamics (Van Wyk de Vries, Wickert, et al., 2022) in Patagonia and other similar environments (Dai et al., 2020).We choose the region around Glaciar Oriental as a case study due to its complex environment (high relief, frequent snow cover, presence of fast-flowing glaciers and large lakes).Due to the immediately adjacent large glacier, we carefully verify that any ground displacements are not affected by measurements of glacier flow.We run the workflow over the region surrounding the Glaciar Oriental using 19 cloud-free Sentinel-2 L1C images from January 2016 to January 2023, for a total of 127 image pairs.

| Case Study 3: Central Sindhupalchok district, Nepal
Nepal, located on the southern margin of the Himalaya, has one of the steepest topographic gradients in the world.Combined with extreme seasonal rainfall and frequent seismic activity, this high relief makes much of the country highly susceptible to landslides.Largescale surveys of shallow landslides have been compiled (Adhikari et al., 2017;Ambrosi et al., 2018;Kincey et al., 2021;Rosser et al., 2021), but to date, studies of slow-moving landslides have focussed on relatively small spatial extents (Bekaert et al., 2020;Lacroix et al., 2022).Here, we run the landslide detection workflow over a 165 km 2 region of central Sindhupalchok district, including the Bhote Koshi valley and adjacent areas.This valley hosts one of the main Nepal-China cross-border roads and was affected by a large number of coseismic and post-seismic landslides associated with the 2015 Gorka earthquake (Kincey et al., 2021).A previous study has identified a number of slow-moving landslides in this region using feature tracking on very high-resolution Pleiades imagery and lower resolution Sentinel-1 SAR imagery (Lacroix et al., 2022).We run the slow-moving landslide detection workflow over this region using 11 cloud-free Sentinel-2 L1C images over the post-earthquake period between January 2016 and March 2022, for a total of 35 image pairs.As this area is densely populated, we compare the footprints of candi-  The sensitivity test shows that the active portion of the Slumgullion earthflow is detected in all cases, with an area ranging from 1.08 to 1.35 km 2 .Additional small candidate landslides are detected when the minimum pixel threshold is reduced close to zero, particularly when the cut-off threshold is also lowered (Figure S1), but our confidence in these is low.There is a wide range of parameters in which the key landslide is identified with very little change, and the recommended parameters in this workflow effectively suppress less reliable landslide candidates for automated usage over large areas.
Our assessment of false positive and false negative rates (Figures S2 and S3) further shows that these parameters greatly suppress false positive area detection (<0.5% of the Slumgullion landslide area regardless the minimum area cut-off chosen) while retaining a moderate to high true positive area detection (62% to 85% depending on which prior extent is used) for the Slumgullion earthflow.Altogether, these show that the workflow successfully detects most of the active landslide while falsely detecting little to none of the inactive landslide area or surrounding stable terrain.
Further details on these tests are discussed in the Supporting Information.

| Glaciar Oriental
The Glaciar Oriental study area has a median velocity magnitude of À0.0050 m a À1 (IQR = À0.203 to 0.134) in the E-W direction and 0.0015 m a À1 (IQR = À0.149 to 0.122) in the N-S direction.This excludes all signal related to the flow of Glaciar Oriental and marginal glaciers, the velocity of which reaches several hundred m a À1 on the main branch.Our automated landslide detection workflow identifies two candidate slow-moving landslides with areas of 0.65 and 6.16 km 2 .The candidate landslides have median E-W velocities of À4.0 and À5.8 m a À1 and median N-S speeds of À3.9 and À4.0 m a À1 , respectively (Figure 4).The 6.16 km 2 candidate landslide complex identified is composed of multiple lobes across the majority of Glaciar Oriental's eastern lateral moraine.This feature is inferred to be a composite landslide that initiates at an elevation of $1330 m and has an elevation range of approximately 1000 m from crest to toe.This landslide has a maximum E-W speed of 26.9 m a À1 and a N-S speed of 18.9 m a À1 .The $300 m.The two candidate landslides might be a part of the same landslide complex, with the intermediary area deforming at a velocity slower than the detection threshold of this method or temporarily stabilised.

| Central Sindhupalchok district, Nepal
The central Sindhupalchok district study area has a median velocity magnitude of À0.0076 m a À1 (IQR = À0.138 to 0.121) in the E-W direction and 0.0040 m a À1 (IQR = À0.127 to 0.121) in the N-S direction.Our automated landslide detection workflow identifies 20 candidate slow-moving landslides with areas ranging from 0.14 1.72 km 2 (mean = 0.50 km 2 ).The largest candidate landslide (Figure 5, Object 4) has a median speed of À0.11 and À1.46 m a À1 in the E-W and N-S directions, respectively.A number of small landslide scars are visible in optical satellite imagery at the toe of this landslide.
The attributes of all 20 landslides are summarized in Table 1.within this radius.Similarly, 11 landslides have at least one building within the zone of active deformation (55%), with two having more than 50 buildings (10%).These buildings are potentially at high risk of destabilization due to gradual slope motion.All landslides have buildings with a 2 km radius, with an average of 811 buildings within this zone.Assessment of the hazard resulting from these landslides would require a further in depth analysis, but their proximity to populated areas suggests that they should be closely monitored.

| DISCUSSION
We show, through a detailed analysis of three case studies, that it is possible to automatically detect and map the geographic extent of slow-moving landslides using openly available, frequently collected Sentinel-2 optical satellite imagery.In a test case on the welldocumented Slumgullion earthflow, USA, our workflow was able to effectively delineate the active portion of the earthflow and produced velocity magnitudes consistent with both ground-based instrumentation and previous remote-sensing estimates.One particular benefit of optical feature tracking over alternate techniques such as InSAR is its capacity to measure velocities in the range of 10 0 -10 2 m/year and, therefore, might identify whether any landslides transition to more hazardous rapid motion.In another test case on the margin of the SPI, Chile, we identified a very large (>6 km 2 ) composite landslide in the eastern lateral moraine of Glacier Occidental and, conversely, demonstrated that the western moraine of this glacier appears to be stable within the bounds of this technique.Finally, we identified 20 candidate slow-moving landslides in a 165 km 2 area around central Sindhupalchok district, Nepal, some of which have previously been mapped and underlie or are in close proximity to densely populated areas.These results, using three case studies from very different geographical environments, demonstrate the effectiveness and wide applicability of our proposed workflow.
Our workflow is an advance on previous optical satellite imagery landslide identification methods in two main ways.Firstly, we identify the location and extent of landslides using Sentinel-2 imagery, instead of the airphotos or very high-resolution commercial satellite imagery used in most previous studies (e.g., Lacroix et al., 2022;Stumpf et al., 2017).Sentinel-2 imagery is freely available at a high revisit interval (10 days at most and typically 2-3 days in mid-latitude locations) for any location on Earth (Lacroix et al., 2018;Provost et al., 2022).Zhang et al., 2022 describe an alternative method for blind slow-moving landslide detection at a site in SW China using Sentinel-2 data, with noise reduction performed using time series analysis methods.Second, our workflow produces an inventory of slow-moving landslides in an automated manner, without the need for prior location information, image georeferencing or manual data processing.Taken together, this means that our workflow can be applied to any region on Earth, regardless of the availability or reliability of prior information.In addition, eliminating the need for time-consuming manual intervention and computationally efficient feature-tracking implementation facilitates the upscaling of this method for national to continental-scale slow-moving landslide surveys.
While the candidate slow-moving landslides from this workflow can be used as a final product, further validation of their location, extent and displacement magnitude using ground-based surveys or alternative remote-sensing approaches is always recommended.There are situations in which our workflow will either fail to identify true T A B L E 1 Attributes of the Sindhupalchok district landslides.slow-moving landslides (false negative) or falsely identify nonlandslide areas (false positive).True landslides will be missed in two main situations: either their movement is too slow or too intermittent to be detected by a feature tracking-derived median velocity map or rapid displacement and/or other surface change caused rapid surface decorrelation.Landslides with velocities too low for detection using feature tracking might be identifiable using alternative methods such as InSAR (Ferrigno et al., 2017;Wasowski & Bovenga, 2014) or ground-based monitoring (Khan et al., 2021).This workflow may work less well in areas of the globe with conditions commonly leading to decorrelation, most commonly due to persistent snow cover, but this is likely to be the case with any optical feature-tracking workflow.
Another rare case in which real landslides would be missed is in the event that their motion is not approximately oriented downslope, as may be the case for lateral spreads.The use of median velocities may also result in failure to detect landslides with movement over too small a fraction of the time period covered, particularly where this is close the beginning or end of the observation period.Future work will determine whether time series analysis can identify some of these false negative cases without greatly increasing the false positive rate.
False positive landslide identifications can currently occur in situations with image distortions leading to large, topographically correlated error or where other real, non-landslide surface displacements are occurring.Glaciers are a common non-landslide source of surface displacements, and any errors in ice masks will introduce artefacts into landslide maps.Other geophysical surface displacements, such as coseismic slip, can also lead to false positive identifications.In summary, we can place bounds on the type of landslides which may be detected using this technique: Their area must be greater than $0.05 km 2 , their velocity must be greater than $0.1 m/year, their surface must retain coherence between images, their motion must be approximately downslope and they must not be underneath a glacier (which are masked out).We, therefore, recommend that any resulting candidate landslides or large-scale inventories be interpreted with these limitations in mind.
The three examples investigated in this study highlight the importance of considering slow-moving landslides in multi-hazard assessments.In the case of the Slumgullion earthflow, the toe of the landslide has dammed a >1 km 2 lake.Monitoring of the Slumgullion landslide can therefore inform about any changes in the dam and provide advance warning of associated flooding or dam instability.In this case, we see that the active portion of the landslide does not extend into the dam, a finding confirmed by an extensive in situ monitoring network (Coe et al., 2003;Parise & Guzzi, 1992).In the case of the Glaciar Occidental lateral moraine, a very large and possibly paraglacial landslide has formed adjacent to a rapidly retreating glacier.
The median total displacement of this landslide over the past 7 years (49 m) is of the same scale as the lateral retreat of the glacier over the same period (40-80 m), suggesting that ice loss and slope instability might be coupled.Close monitoring of glacier margins can provide insight into whether unstable areas are likely to collapse onto the ice during retreat, collapse into the fjord or restabilize once the glacier has retreated.Both collapses onto the ice and into the fjord can be hazardous (Dai et al., 2020;Gardner & Hewitt, 1990;Van Wyk de Vries, Wickert, et al., 2022), and ongoing monitoring of this area is recommended to understand whether this landslide poses a glacial lake outburst flood hazard to communities around Lago O'Higgins/ San Martin or Rio Pascua.For the central Sindhupalchok district study area, Lacroix et al., 2022 showed that some slow-moving landslides in the area accelerated in response to the 2015 Gorka earthquake.This region is also extremely susceptible to shallow landslides, which are found on at least 75% (15/20) of the 20 candidate slow-moving landslides and have killed more than 100 people in Sindhupalchok district alone over the past decade (source: http://drrportal.gov.np/).The high population densities of this region mean that any feedback with other hazards would be particularly impactful.
Slow-moving landslides are sensitive to changes in climate, both directly through changes in their hydrology and water availability (Ardizzone et al., 2023;Gariano & Guzzetti, 2016;Handwerger et al., 2019Handwerger et al., , 2022) ) and indirectly through changes in their boundary conditions (Dai et al., 2020;Mititelu-Ionus¸et al., 2011;Van Wyk de Vries, Bhushan, et al., 2022;Yang et al., 2021).The close match between the average displacement of the Oriental lateral moraine and lateral retreat rate of the glacier is an example of this second category.
Rapid climate warming is a driving factor for the current high rates of ice loss (Abdel Jaber et al., 2019;Willis et al., 2012), which in turn impact the stability of surrounding slopes.Other changes in local conditions such as vegetation cover, snow accumulation, rivers and landslide dam or glacial outburst floods can also affect the long-term motion of slow-moving landslides (Cook et al., 2018;Mititelu-Ionus ȩt al., 2011;Van Wyk de Vries, Bhushan, et al., 2022;Yang et al., 2021).In many instances, slow-moving landslide motion is driven by the availability of water which increases internal pore-water pressure and decreases the frictional strength of slopes (Bogaard & Greco, 2016;Lacroix et al., 2020).Climate change is projected to cause a complex array of hydrological changes around the globe, with increases in precipitation in some regions and drying in others (Eyring et al., 2016;John et al., 2022;Pörtner et al., 2022;Tebaldi et al., 2021).In Nepal specifically, forecasts of the Asian monsoon change predict an increase in both the mean and intensity of precipitation over the 21st century (Chevuturi et al., 2018;Wu et al., 2022), although there remains considerable uncertainty in the magnitude of these changes.Basic information about the abundance and geographic distribution of slow-moving landslides is a requirement for evaluating their sensitivity to future climate change but is currently lacking in most regions of the globe.
Systematic mapping of slow-moving landslides across varied environments can increase our sample size of multi-hazard interactions between slow-moving landslides and other hazards-including earthquakes, shallow landslides, debris flows and landslide-dammed lakesand raises the possibility of identifying new interactions that had previously been overlooked or unresolvable due to a lack of data.
Large-scale, automated slow-moving landslide mapping can be used to identify new locations for detailed field investigations, either to improve our understanding of slow-moving landslide processes, to quantify of the sensitivity of these landslides to future climate change, or to monitor particularly dangerous landslides.

| CONCLUSIONS
Slow-moving landslides cause gradual deformation of the Earth's surface and have traditionally been more difficult to detect due to this limited surface signature.Here, we describe a proof-of-concept date slow-moving landslides to a population density map (Meta High Resolution Population Density data, with 30 m cell size) and a building location map (OpenStreetMap), which provide two different estimates of local exposure to hazards.
The 17 km by 10 km Slumgullion study area has a median velocity magnitude (speed) of À0.0055 m a À1 (interquartile range[IQR] = À0.132 to 0.128) in the E-W direction, À0.0088 m a À1 (IQR = À0.238 to 0.218) in the N-S direction.Our automated landslide detection workflow identifies three candidate slow-moving landslides with areas of 1.16, 0.22 and 0.18 km 2 .The three candidate landslides have median E-W speeds of 0.17, À1.4 and À0.18 m a À1 and median N-S speeds of À0.75, 0.33 and À0.89 m a À1 , respectively (Figure3shows a zoom on the main landslide).The largest candidate slow-moving landslide also has the most rapid motion, with a maximum E-W speed of 3.8 m a À1 and N-S speed of 2.55 m a À1 .The spatial extent of this candidate slowmoving landslide closely matches the independently mapped active section of the Slumgullion earthflow, and the maximum velocity we calculate (4.40 m a À1 , or 12.1 mm per day) is within the range of previous ground and satellite-based estimates (3.65-5.84m a À1 , Coe et al., 2003; $5.10 m a À1 , Hu et al., 2020).Visual inspection of satellite imagery suggests that two other candidate slow-moving landslides might correspond to rock glaciers, but we cannot confirm this without additional ground-based observations.We conduct both a sensitivity test and an assessment of false positive and false negative matches relative to prior knowledge on the extent(Hu et al., 2020;Parise & Guzzi, 1992) of the landslide.F I G U R E 3 Results of the Slumgullion automated landslide detection.(a) The active portion of the Slumgullion earthflow (extents from Hu et al., 2020; Parise & Guzzi, 1992) is clearly delineated while the inactive portion is not.(b,c) Pixel frequency histograms of the E-W and N-S velocity components across the entire study area, respectively.The red line represents the outline of the active earthflow (with dashed areas uncertain) while the white line represents the extent of the inactive earthflow (Hu et al., 2020).
second, smaller, candidate landslide is located in the same eastern lateral moraine and separated from the main composite landslide by only F I G U R E 4 Results of the Glaciar Oriental automated landslide detection.(a) A large portion of the eastern lateral moraine of this glacier is unstable, while the western lateral moraine is stable.(b,c) Pixel frequency histograms of the E-W and N-S velocity components across the entire study area, respectively.
Only one of the candidate landslides we identify directly matches a landslide shown byLacroix et al. (2022): 20 or d at Pokhan.This is consistent withLacroix et al.'s (2022) finding that landslides b, e, g and h stabilized post-2015, while motion continued on d.Lacroix et al.'s (2022) landslide c partially overlaps with Landslide 13 from this work but with a different overall geometry, possibly representing the activation of a different part of a larger unstable slope.Lacroix et al.'s (2022) study however covers a different time period (2014-2017) to our own (2016-2022) and observes that landslide velocities were (i) substantially modified, in the short term, by the 2015 Gorka earthquake and (ii) vary through time even outside of this event.The methods used by Lacroix et al. (2022), feature tracking of very highresolution Pleiades imagery and Sentinel-1 range displacement, also differ from our own, as do the number of individual velocity maps combined prior to landslide identification.Therefore, despite considering overlapping areas, the two studies are measuring different expected surface deformation fields and differences in landslide identification and maximum velocity are expected.The close match between the Pokhan landslide (20/d) in this study and Lacroix et al. (2022) confirms that landslides with continued motion can be detected with both methods.A comparison of our candidate landslides with a high-resolution population density map shows that the region of active deformation is inhabited for 10 of these landslides (50%), with over 50 people living on five of these (25%).All landslides have people living within 2 km of the zone of active deformation, with an average of 2460 people living F I G U R E 5 Results of the Sindhupalchok district automated landslide detection.(a) A number of large deforming slopes are present across this high-relief terrain.The dashed lines show outlines of landslides b-g identified by Lacroix et al. (2022).(b,c) Pixel frequency histograms of the E-W and N-S velocity components across the entire study area, respectively.