Precursory Motion and Time‐Of‐Failure Prediction of the Achoma Landslide, Peru, From High Frequency PlanetScope Satellites

Landslide time‐of‐failure prediction is crucial in natural hazards, often requiring precise measurements from in situ instruments. This instrumentation is not always possible, and remote‐sensing techniques have been questioned for detecting precursors and predicting landslides. Here, based on high frequency acquisitions of the PlanetScope satellite constellation, we study the kinematics of a large landslide located in Peru that failed in June 2020. We show that the landslide underwent a progressive acceleration in the 3 months before its failure, reaching at most 8 m of total displacement. The high frequency of satellite revisit allows us to apply the popular Fukuzono method for landslide time‐of‐failure prediction, with sufficient confidence for faster moving areas of the landslide. These results open new opportunities for landslide precursors detection from space, but also show the probable seldom applicability of the optical satellites for landslide time‐of‐failure prediction.

growth processes (Bruckl & Parotidis, 2001).This mechanism is characterized by a last kinematic phase with accelerating displacement until the rupture (Bhandari, 1988;Zvelebill & Moser, 2001).This accelerating pattern is a real opportunity for landslide prediction, and different methods have emerged to predict the time-of-failure of landslides based on this mathematical law.The most popular one is the Fukuzono method (Fukuzono, 1985) which relates the time-of-failure of the landslide to its inverse velocity thanks to a power-law coefficient α.For landslides, α is found to vary from 1.3 to 3.3 (Crosta & Agliardi, 2003;Federico et al., 2012;Intrieri et al., 2019).Assuming α = 2 leads to a linear relation between inverse velocity and time, which makes the method easy to apply to real data.
The inverse velocity method is widely used, even on an operational basis, mostly in mines (Rose & Hungr, 2007), based on time-series of landslide displacement from in situ measurements, coming from terrestrial laser scanning, extensometers, ground-based SAR, or global navigation satellite systems (GNSS).These measurements have the advantage of a high frequency rate and low uncertainties, required to derive velocity time-series.They also have the drawbacks of being punctuated, and require instrument maintenance.Most of all, these measurements require the landslide to be detected and instrumented before its failure, which is rarely the case.Remote sensing of slopes is thus of real interest for landslide prediction (Casagli et al., 2023), either from terrestrial cameras that scan a certain limited area (typically at a slope scale) or from airbornes or satellites.However, the noise of satellite remote sensing processing to derive ground surface displacement and their frequency of acquisition raise questions about the ability of satellites to detect these accelerations and predict their time-of-failure (Carlà et al., 2019;Intrieri et al., 2018;Lacroix et al., 2018;Moretto et al., 2017).
The development of time-series techniques and the increase of satellite revisit time has enabled the generation of high frequency time-series of landslide displacement (Intrieri et al., 2018;Lacroix et al., 2018).This spatially dense and increasingly temporally-resolved data set has enabled scientists to show the high variability of landslide kinematics at different time-scales, mostly focusing on non catastrophic landslides (Handwerger et al., 2013(Handwerger et al., , 2015;;Lacroix et al., 2019).However, at least 5 studies show the possibility for satellites to detect landslide pre-failure accelerations before their transition from slow to rapid motions (Carlà et al., 2019;Handwerger et al., 2019;Intrieri et al., 2018;Lacroix et al., 2018;Li et al., 2020) either using InSAR, or correlation of SAR or optical images.
Different limitations arise from the analysis of these different case-studies.First of all, the acceleration phase can produce very small displacements, hardly visible in the satellite image processing.For instance, on the Xianmo landslide the displacement observed by InSAR during the acceleration phase is only 4 cm (Intrieri et al., 2018).Then, the landslides can have quick accelerations relative to the sampling frequency.On the Xianmo landslide the acceleration is seen only over 6 Sentinel-1 acquisitions (Intrieri et al., 2018).On the Harmalière landslide, the acceleration, accomodating 4 m of displacement in 3 days, is seen only over 2 Sentinel-2 acquisitions (Lacroix et al., 2018).This low number of data points in the acceleration period prevents the prediction of the landslide time-of-failure.Very few case studies show substantial and long-lasting landslide accelerations seen with satellites.To our knowledge, this was observed only on the Baige landslide (Li et al., 2020), where the acceleration phase occurred over 1.5 years, with cumulative motions reaching 45 m.Usually, acceleration may occur over much shorter periods, with smaller magnitudes (Moretto et al., 2017).In this context, high frequency satellites like the PlanetScope cubesat constellation (Roy et al., 2021) offer a unique opportunity to increase the temporal sampling of these accelerations.Despite their sensitivity to cloud cover, one main advantage of optical image correlation compared to InSAR is its ability to resolve large motions over short period of times, and the possibility to obtain easily-interpretable 2D displacements.In this study we evaluate the potential of PlanetScope satellites to detect pre-failure accelerations and to predict landslide time-of-failure.For that purpose we back analyze a large landslide in Peru that failed catastrophically in June 2020.

Study Site
The Colca area, situated in South Peru (Figure 1), is a deep river incised valley inside lacustrine deposits from a paleo-lake of 30 km long.The terraces have allowed for villages and agriculture settlements for several millenia.The Colca river erodes the lacustrine deposits and trigger landslides that can then be remobilized by the seasonal rainfall infiltration (Zerathe et al., 2016) and M5+ earthquakes (Bontemps et al., 2020).Many rapid and slow-moving (m/yr) landslides were detected in the area, based on the comparison of DEMs (Zerathe et al., 2016) and the correlation of diachronic optical images (Bontemps et al., 2018;Lacroix et al., 2015;Pham et al., 2018).
The large Achoma landslide (picture in Figure S1 in Supporting Information S1), which failed the 18th of June 2020, is connected to the river.The first signs of activity, mostly vertical cracks, were reported by local inhabitants in May 2020.The landslide failure occurred during the night, but several seismic sensors situated between 8 and 70 km away detected the signal generated by the Achoma landslide failure (Vela Valdez & Chacaltana, 2020), dating the landslide occurrence to the 18th of June 2020 at 6h42 UTC.The rainy season stopped mid April 2020, so this failure occurred during the dry season.The triggering factor of this landslide is still an open question.
Digital Elevation Models were obtained before (2017/5/13) and right after (2020/06/19) the landslide, thanks to the processing of Pléiades (Lacroix et al., 2015) and drone acquisitions (see Figure 1b and Supporting Information S1), allowing the quantification of its characteristic size.The landslide is 800 × 500 m in size and the landslide headscarp, nearly vertical, is more than 100 m high.The landslide affects a zone of agriculture, and dammed the river, creating a lake upstream (Figure 1d).The risk of flooding posed by this dam was quickly solved by engineers digging inside the landslide mass and creating a path for the river flow.

PlanetScope CubeSat Images
In order to study the pre-failure motion of the landslide, we realized a time-series of ground displacement using the high-frequency PlanetScope satellites.The complete PlanetScope constellation, composed of approximately 180 satellites is able to image the entire land surface of the Earth every day (Roy et al., 2021).Each PlanetScope satellite is a CubeSat 3U form factor. Three generations of PlanetScope satellites have been launched since July 2014, with increasing footprints and multispectral capabilities.The third generation provides images after March 2020, so the images processed here are coming mostly from the two first generations.These two generations have similar characteristics (https://assets.planet.com/docs/Planet_Combined_Imagery_Product_Specs_letter_screen.pdf), namely a four-band frame imager, with resolution between 3 and 4.1 m resampled to 3 m.The data are provided at different levels of processing on the webserver (https://www.planet.com/explorer/).After several trials, we used the "PlanetScope Scene" product to avoid the mosaicking of different stripes of images acquired with different viewing angles on the same day.In this manner, we control the mosaicking of the images over the selected area and ease its correction.
We selected 79 images over an area of 3.3 × 3 km 2 (Figure 1), with a frequency of 1 per month during the low velocity period of the landslide from November 2017 to December 2019 to estimate the seasonal errors due to shadowing effects (Lacroix et al., 2019) and then all available cloud-free images in the period of interest from January to June 17 2020 (Table S2 in Supporting Information S1), corresponding to the last acquisition day prior to the landslide rupture.

Displacement Time-Series
A time-series of horizontal ground displacement is calculated following the approach developped in Bontemps et al. (2018) and adapted for the PlanetScope images (Figure S4 in Supporting Information S1): (a) correlation of all the pairs of images using Mic-Mac (Rupnik et al., 2017), which is well suited to images of low radiometric contrasts or small objects (Lacroix et al., 2020), (b) masking the low correlation coefficient values (CC < 0.7), (c) mosaicking correction, similar to stripe corrections (Bontemps et al., 2018), that we obtained by subtracting the median value of the stacked profile in the along-stripe direction, taking into account only stable areas, (d) least square inversion of the redundant system per pixel, weighted by the time separation between pairs (Bontemps et al., 2018), (e) correction of illumination effects (Lacroix et al., 2019), based on the 2 years of data between November 2017 and December 2019.

Time-Of-Failure Prediction
The most popular method for landslide time-of-failure prediction, the inverse velocity method, is based on a Fukuzono model (Fukuzono, 1985).This latter provides an analytical solution for the landslide velocity v to the common Voight law (Voight, 1989), that links the acceleration a and the velocity of any material in its final stage of acceleration under conditions of approximately constant stress, through two constants A and α:

𝑎𝑎(𝑡𝑡).𝑣𝑣(𝑡𝑡) −𝛼𝛼 = 𝐴𝐴
(1) With t f the time-of-failure.In the case of α = 2, 1/v becomes linear with t, and t f is therefore obtained by extrapolating the linear relation between the inverse-velocity of the landslide and time up to the zero inverse-velocity.This method is therefore applied on landslide velocity time-series, not displacements.Due to large uncertainties in the ground displacement values from optical satellites (see Section 3.1), the simple time derivation of the ground displacement time-series leads to unreliable velocities.Previous authors dealing with such problems first filtered the displacement time-series before their derivation, using different types of median filters or moving averages (Bozzano et al., 2018;Carlà et al., 2017;Intrieri et al., 2018;Sharifi et al., 2022;Zhang et al., 2021).These filters rely on a sufficiently high rate of observations, which can be common with in-situ measurements but not with satellite remote sensing.
Here we adapted the strategy previously used by Azimi et al. (1988), which consists of evaluating the velocity of the landslide over segments of equal displacement, to reduce the effect of noise in the displacement data.This idea takes advantage of the monotonic direction of the landslide displacement with time to provide more robust time-series.We adapt this strategy to consider the displacements not directly on the raw values but rather on linear fits of the raw values with time (see Figure S6 in Supporting Information S1).

Landslide Kinematics
The maps of cumulative displacements over the 2.5 years of measurements display a ground motion before the landslide rupture of between 3 and 8.1 m, located exactly on the mass that collapsed on the 18th of June 2020 (Figure 2).Stronger motions are encountered at the toe than at the headscarp location, suggesting either a fragmentation of the mass, or a rotational component of the landslide with a motion being mostly vertical at the headscarp and mostly horizontal at the toe.This last explanation is more plausible given the highly vertical headscarp of more than 100 m height observed after the landslide failure.
Around the landslide area, no clear signals of deformation are visible.The uncertainties (Figure S5 in Supporting Information S1), calculated on these stable areas, are between 0.2 and 0.9 m in EW and 0.2 and 1.2 m along the NS component, with a seasonal variation due to the seasonal variation of the shadows (Lacroix et al., 2019) in addition to a slight linear trend with time.Some dates clearly exceed this seasonal variation, so we remove these outliers (8 points out of 79 dates, see Figure S5 in Supporting Information S1 The behavior of all individual pixels on the landslide display a similar pattern with time, with a slow motion from November 2018 to March 2020, followed by a progressive acceleration over 3 months until the rupture on the 18th of June 2020 (Figure 2).

Landslide Acceleration Modeling
To better understand the physics behind the general pattern of acceleration, we first approximated the mathematical law of this acceleration using the classical Fukuzono model (see Equation 2).We fixed the time of failure t f , known from seismic sensors (Vela Valdez & Chacaltana, 2020), and inverted the other parameters of the model by minimizing the norm of the difference between the displacement time-series observed and the modeled one using a grid search over the parameters A and α.The Fukuzono model gives an analytical solution for the velocity, so that we integrate numerically the velocity to have a value for the modeled displacement.This inversion process leads to α = 2.14 ± 0.12 for the pixels with sufficiently large acceleration to be well resolved (Figures S7 and S8 in Supporting Information S1).We also find that α and A are not fully independent (see Figure S9 in Supporting Information S1).

Landslide Prediction
We now test the ability of this data set to predict the time-of-failure, by fixing α = 2 (despite we know that it is not exactly true), in order to use the classical inverse velocity model, and therefore invert t f .The time-series of displacement are highly noisy, with uncertainties of about 0.8 m for maximum pre-failure displacements of 8.1 m.Therefore, the question can be raised about their utility for landslide prediction.
We extract the velocity time-series from the displacement time-series by applying the algorithm presented in Section 2.3 to the pixels situated on the landslide.The time-series of velocity display an acceleration starting from the end of 2019 with values exceeding 1 cm/day around 1 March 2020 (Figure S6 in Supporting Information S1).
The inverse-velocity method applied to this data set is not strictly following the α = 2 regime, so that the linear fit over the last dates is contaminated by the previous measurements (Figure 3; Figure S10 in Supporting Information S1).To partly overcome this limitation, the linear regression of the inverse velocity with time is weighted by the inverse of the uncertainties of each velocity estimation, which are lower when the velocities are higher, that is close to the rupture.Using this method, a time-of-failure is obtained for every pixel of the landslide.Given the constraint that the time-of-failure occurs after the last satellite acquisition, the empirical probability density function (pdf) of the time-of-failure is positively skewed (Figure 3 and Figure S11 in Supporting Information S1).
We choose to characterize this pdf with two values: the median of the pdf (M pdf ), and the probability (P 3 ) to have a failure within 3 days after the last satellite acquisition (t sat ).The cumulative density function C(t f < t) is first calculated: P 3 is then obtained by calculating C(t f < t) for t = t sat + 3 days.These 3 days are chosen as it may be an acceptable time required to process the last acquired image, analyze the results, and eventually deliver an alert to the populations.
At the time of the last satellite acquisition (t sat = 17th June), the M pdf is found to be t sat (Figure 3b), due certainly to a value of α greater than 2. The probability P 3 to have a landslide failure within 3 days is found to be 0.72.
We can wonder how long in advance could we have predicted the time-of-failure of this event.For that purpose, we removed iteratively the last satellite acquisition and recalculated the time-series of velocities and the resulting t f prediction using the inverse velocity method.This process gives us a landslide time-of-failure estimate as a function of the last satellite acquisition t sat of the time-series (Figure S11 in Supporting Information S1).For each one of these time-series we also calculate M pdf and P 3 for groups of pixels that have different motion magnitudes (4-5 m, 5-6 m, 6-7 m, greater than 7 m) (Figure 4).The results show that for pixels of large motions (displacement greater than 7 m, 28% of the total pixels of the landslide), M pdf is fairly stable over the last 3 weeks before the failure, and estimates a time-of-failure within less than 8 days of the real rupture.For pixels with motions between 5 and 7 m, M pdf follows the same pattern but with less precision.For pixels with total motions less than 5 m, M pdf is always increasing with time and never converges.
The probability P 3 is always increasing with time, except for pixels with motion less than 5 m.This means that there is a minimum movement magnitude required to be able to predict a time-of-failure.This is also in agreement with the always increasing values of M pdf for pixels of total motion less than 5 m.For other groups of pixels, P 3 is increasing with time to reach the highest values close to the rupture time.The probability P 3 for motions larger than 7 m surpass 0.7 in the 3 days before the real rupture time.P 3 is also much lower for larger motions before the 10th of June.This means that the number of false alarms is considerably reduced by selecting only pixels with larger motions.We can evaluate the quality of the prediction by comparing P 3 in the first weeks of the acceleration and in the days just before the failure.For pixels of motion magnitudes greater than 7 m, P 3 ranges from 0.15 to 0.8, showing the ability of this predictor to evaluate the proximity of the landslide to the rupture.The satellite observations processed here, provide clear evidence that the landslide final stage of acceleration was not initiated during the dry season, but started at least in March during the rainy period.It must be noted that the 2020 rainy season was the rainiest season for at least 20 years, with clear impacts on the landslides of the area.
For instance, on the nearby slow-moving Maca landslide situated in the same geology (lacustrine sediments), the displacement reaches 10 m in 2020, more than the 7 m of the very active year 2012 (Zerathe et al., 2016).For that reason, and also due to its similarity with the Maca landslide, we can suspect the rain infiltration and/or the river erosion to be serious controlling factors of the Achoma landslide.Ground water storage might possibly play a role in the delayed acceleration of the landslide, as it is observed on many different slow-moving landslides (e.g., Bièvre et al., 2018) and in particular on the nearby Maca landslide (Zerathe et al., 2016).However in Maca, despite the motion continues for weeks after the rainy season, the velocity decreases with time.Here we observe the opposite behavior, with an acceleration toward the failure, discarding the ground water storage as the final landslide trigger.An analysis of the seismicity, obtained by a local seismic network operated by the IGP (Instituto Geofisico del Peru) with completude magnitude of Mℓ 3.5 doesn't show earthquakes greater than 4 in the vicinity of the landslide in the first 7 months of 2020.During that period, the largest earthquake recorded in the 50 km from the landslide was Mℓ 3.9, and the closest earthquake was 10 km away from the landslide (Figure S12 in Supporting Information S1).Therefore, it is highly unprobable that an earthquake produced this acceleration (Bontemps et al., 2020).
The pattern of progressive acceleration observed here follows a Voights law (Voight, 1989), typical of the last stage of tertiary creep, showing that the landslide failure is certainly linked to a process of progressive damage and fault maturation.This observation reinforces the hypothesis of a long process of landslide triggering, where the landslide acceleration is only the final stage of damage localization along a shear band (Zvelebill & Moser, 2001).The α coefficient of the Voight's law is spatially homogeneous (Figure S7 in Supporting Information S1), showing that the deformation occurs at depth, without much fracturation inside the moving mass.The variability of the motion magnitude is mainly caused by the rotational aspect of the failure surface that will change the ratio between the vertical and the horizontal motions.The exponent of the Voight's law is 2.14 ± 0.12, in the upper range of what is usually observed on other landslides (Intrieri et al., 2019).This high exponent may explain the long time-period of landslide acceleration.Indeed 3 months of landslide progressive pre-failure acceleration is in the upper limit of what is usually observed for this last stage of acceleration on other landslides (Federico et al., 2012;Scoppettuolo et al., 2020).Usually this phase can last between a few minutes to about 100 days for different types of rocks/soils and mechanisms (rockfalls, slides, toppling, …).
We can also wonder if ground movements were preceeding this final acceleration, as landslide kinematics can fluctuate at different time-scales (seasonally, interannualy) depending on the controlling factors of the landslide motion (e.g., seasonal rainfalls).On the Achoma landslide, the pattern of motion before March 2020 is not clear,

Time-Of-Failure Prediction
The process developed here allows us to predict the time-of-failure 3 weeks before the rupture with a precision of 8 days, for the points of largest motion magnitudes.The larger the motion magnitudes, the better the prediction.Using the probability of failure within 3 days, we show that the number of false alarms is highly reduced when selecting only the pixels of largest motion magnitudes, due to the large uncertainties of the displacement time-series from PlanetScope images.Despite the specific algorithm used to retrieve the velocities from the displacement time-series, the prediction is not possible for motions of magnitudes lower than 5 m, that is about six times the displacement field uncertainties obtained with PlanetScope.These encouraging results show the real advantage of the high frequency of PlanetScope image acquisition, which allows us to obtain a dense time-series of displacement, and therefore to improve the sampling of the pre-failure acceleration.This is important especially for the velocity estimate, which can be altered by the uncertainty of the displacement measurements from satellites.
A large number of predictions are providing a time-of-failure before the time of the last image acquisition, which is simply not possible.These false alarms are caused by the measurement uncertainties and by an α coefficient larger than 2 (mean of 2.14).For this reason the inverse velocity model should always bias the prediction toward lower time-of-failure.We partially solved this issue by adapting the linear fit of the inverse velocity model with time using a linear fit weighted by the inverse uncertainties of the velocities, thus providing more importance to the last (and fastest) measurements.To fully overcome this problem of α different than 2, other prediction methods should be tested.The method proposed by Hao et al. (2016), linear with α, seems appealing.However, this method requires the computation of not only the first but also the second derivative of the displacement time-series, which will be tricky to implement on the optical satellite measurements, due to their large uncer tain ties.Instead, we believe that further developments should try to find an analytical expression of the Fukuzono model in terms of displacement and not velocity.
Despite the relative success of this Achoma case-study for landslide prediction, we believe it is not representative of the majority of the landslides, and the time-of-failure prediction from space will be rather complex for the majority of other landslides.Indeed, first of all, the Achoma landslide was large enough to be easily detected and spatially covered by satellites, which enabled us to compute statistics.Then its location and time-of-failure, 2 months after the end of the rainy season, occurred in an area with little cloud cover, so that satellite acquisitions were cloud-free very often over the 2 last months before the failure.Finally its period of acceleration was fairly long (3 months), and rather slow acceleration as shown by the α coefficient greater than 2.These two factors enable a good sampling of the landslide acceleration period.For all these reasons, we believe Achoma was an "ideal" case-study for this time-of-failure prediction from optical satellite images.It is not obvious that existing SAR satellites can do better than optical for this purpose.Indeed, InSAR presents other challenges like phase unwrapping of small objects, with phase-ambiguities even more complex to fix with the recent loss in 2022 of the Sentinel-1B satellite.The future launch of the combined L-and S-bands NiSAR radar satellite, with 2-5 days revisit time and longer wavelength, will greatly improve the capabilities of SAR satellites for time-of-failure prediction.
Finally, the prediction from satellites, either SAR or optical, relies on time-series of ground deformation.This prediction has an interest mostly in cases where a time-series of ground displacement can be automatically generated and a slow-moving landslide automatically detected.In the case of PlanetScope images, the automatic generation of ground displacement suffers from the image selection process which is a manual step, involving selecting the images acquired with similar orientations and without cloud-cover.For landslide detection (and therefore prediction) from Planetscope images, it would be worth working on the automatic selection of the images.

Summary and Conclusions
Based on a 2.5 years time-series of ground displacement obtained from PlanetScope optical satellites, we document new observations of a landslide accelerating toward the rupture with an unprecedent view of the time and

Figure 1 .
Figure 1.Digital elevation Model of the area (a).Topography along the profile A-B before/after the event (b) measured using Pléiades and drone acquisitions (see Figure S3 in Supporting Information S1).PlanetScope image before/after the landslide occurrence (c).

Figure 2 .
Figure 2. Results of the time-series processing on the Achoma landslide: (a) map of cumulative displacement over 2.5 years before the failure, (b) time-series of cumulative displacement on one specific point (shown with a yellow dot on the map of subplot (a), (c) zoom over the last 3 months display in the progressive acceleration.The red dashed line represents the best fit of the Fukuzono model, fixing t f to its real time.

Figure 3 .
Figure 3. (a) Inverse velocity obtained for one specific point and linear fit (dashed black) leading to a prediction for t f of 16-June-2020 at 2:24.(b) PDF (white/black bars) and CDF (blue curve) of time-of-failure for all the points of the landslide using the whole time-series up to the last acquisition of 17th of June.The white bars in the PDF represent the t f prediction before the last satellite acquisition.

Figure 4 .
Figure 4. Time-of-failure median prediction (a) and probability to have a failure within 3 days (b) as a function of the last satellite image acquisition and for different motion magnitudes.