# Irrigation-triggered landslides in a Peruvian desert caused by modern intensive farming

Pascal Lacroix, Amaury Dehecq, Edu Taipe

► To cite this version:

Pascal Lacroix, Amaury Dehecq, Edu Taipe. Irrigation-triggered landslides in a Peruvian desert caused by modern intensive farming. Nature Geoscience, 2019, 13 (1), pp.56 - 60. 10.1038/s41561-019-0500-x. hal-04945263

HAL Id: hal-04945263

https://hal.science/hal-04945263v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Irrigation-triggered landslides in a Peruvian desert caused by modern intensive farming

Pascal Lacroix ID <fnref n="1" />*, Amaury Dehecq ID ² and Edu Taipe ID ³

Intensification of agriculture leads to stress on the environment and subsequently can have strong societal and ecological impacts. In deserts, areas of very high sensitivity to land-use changes, these local-scale impacts are not well documented. On the arid southwestern coast of Peru, several vast irrigation programmes were developed in the 1950s on the flat detritic plateau surrounding narrow valleys to supply new farming areas. We document the long-term effects of irrigation on the erosion of arid deserts in the Vitor and Siguas valleys, south Peru, using 40 yr of satellite data. We demonstrate that irrigation initiated very large slow-moving landslides, affecting one-third of the valleys. Their kinematics present periods of quiescence and short periods of rapid activity, corresponding to landslide destabilization by their headscarp retrogression. This analysis suggests that the landslide motion continues long after their initiation by irrigations. Those landslides affect the fertile valley floors, leading to the destruction of villages and agricultural areas. We conclude that modern intensive farming can strongly impact traditional agriculture in desert areas where water management is particularly critical.

Agricultural intensification increases the pressure on landscapes<fnref n="1" />,². Land-use change negatively affects soil erosion³, leading to deterioration of onsite soil quality⁴ or excessive sedimentation in reservoirs. These environmental problems are well monitored in temperate areas³, but much less in deserts, despite the high sensitivity of desertic environments to climatic and anthropogenic changes⁵.

The landscape in hyper-arid deserts is in a geomorphic steady state over millennial timescales. Indeed, slowly eroding surfaces are found in the Atacama (0.1 m Myr⁻<fnref n="1" />)⁶ and Namib (0.5–1 m Myr⁻<fnref n="1" />)⁷ deserts. Over shorter timescales, despite highly weathered rocks, mass movements (an important contributor to erosion in slopes) are not often reported in deserts, due both to the low density of population⁸ and low precipitation that decrease the rate of landsliding. Low landslide rates are reported in desertic areas, mostly caused by strong earthquake shaking⁹, extreme climatic events at centennial timescales¹⁰⁻¹² and irrigation¹³,¹⁴. Irrigation-triggered landslides have been reported in the loess sediments of the semi-desertic area of central China¹⁵⁻¹⁷ and in the alluvial deposits of the Atacama hyper-arid desert in Peru¹³,¹⁴. The monitoring of these landslides realized over periods of a maximum of 6 yr suggests a large variability of the landslide kinematics from a few millimetres per year to a few metres per second¹⁴⁻¹⁷, depending on different landslide properties (material, mechanism and slope geometry) and forcing factors. The small number of case studies reported suggests either a low impact on populations of irrigation-triggered landslides in deserts or that the kinematics and failure mechanisms of such processes occur on timescales longer than years. Long-term analysis of the kinematics of irrigation-triggered landslides is therefore required to characterize the impact of irrigation on soil erosion.

On the western coast of Peru, a zone of extreme hyper-aridity¹⁸, agriculture was confined for more than 3,500 yr to deeply incised valleys where river flows coming from the Andes allow for the development of food-producing agriculture¹⁹. Between the 1950s and the 1990s, however, vast irrigation programmes started in the flat detritic plateaus surrounding those narrow valleys to supply

new farming areas. Those deeply incised valleys present steep walls with more than 200 m denivelation, making them prone to landsliding, as shown by the giant palaeo-landslides visible in the geomorphology¹³, certainly dating back from the valley incision before the mid Pliocene²⁰. We focus on this 2,000 km² area to study the impact of land-use change on the erosion of a desertic area over 40 yr (1978–2018).

## 1978–2016 changes in morphology reveal large landslides

We quantify the changes in both land use and morphology over the period 1978–2016 based on satellite images from the Hexagon spy satellite and recent SPOT6/7 images all acquired in stereoscopic mode. The comparison of the images from the two periods shows that the irrigation programmes of Vitor and Siguas led to large expansions of agricultural crops close to valley borders from 1978 to 2016, covering in 2016 an area of 105 and 76 km² respectively, as estimated by a precise mapping of the fields on the satellite images (Fig. 1a).

The stereo acquisition of these images enabled us to quantify the changes in morphology based on the computation of two digital elevation models (DEMs) in 1978 and 2016 (see Methods). The difference of the two DEMs highlights specific areas situated on the valley walls with similar patterns: an elevation loss close to the crest, and an elevation gain in the valley floor (Fig. 1a). This pattern reveals a transportation of sediments from the crest line of the valley walls towards the valley floor. Major changes of the morphology are thus found on 9 km over the 25-km length of the irrigated parts of the Vitor Valley (Fig. 1c,d), and 8 km over the 20-km length of the Siguas Valley (Fig. 1b). Those mass movements are encountered only on the banks of the valleys where irrigation programmes occur. The prominent changes in morphology occurred over 12 distinct areas of size between 0.10 and 1.82 km² (Supplementary Table 2 and Supplementary Fig. 2), including two areas where the elevation differs by more than 50 m: on the Siguas Valley, where a major slow-moving landslide (no. 1 in Fig. 1b) of at least 40 Mm³ initiates in 2005¹⁴ (Supplementary Fig. 2a), and on the Vitor Valley in an area called Punillo Sur (no. 7 in Fig. 1c). This last area is affected by a

<footnote n="1">¹ISTerre - Université Grenoble Alpes, IRD, CNRS, IFSTTAR, Université Savoie Mont Blanc, Grenoble, France. ²Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA. ³INGEMMET, Lima, Peru. *e-mail: pascal.lacroix@univ-grenoble-alpes.fr</footnote>
a

<figure><img src="image_2.png" /></figure>

b

<figure><img src="image_3.png" /></figure>

c

<figure><img src="image_4.png" /></figure>

d

<figure><img src="image_5.png" /></figure>

**Fig. 1 |** Elevation change between 1978 and 2016 in the Siguas and Vitor valleys in southern Peru. **a**, A map showing the study area. The areas outlined in red are magnified in **b–d**. The green areas contour the irrigated fields in 2016 (the dates indicate the origin time of the irrigation programme). Major active landslide scarps are delineated in black. **b–d**, Elevation changes (dZ), indicated according to the colour scale, between 1978 and 2016 as measured by the difference of two DEMs generated from the KH9 spy satellite (1978) and the SPOT6/7 satellite (2016).

large slow-moving landslide of 1.82 km² (Supplementary Fig. 2b) and a volume estimated to 75 M³ (see Supplementary Table 2).

## Irrigation triggered landslides

We process satellite data to follow the kinematics of all these land-slides (see Methods), based on image correlation techniques²¹⁻²³,

combining KH9 spy images (1978), Landsat5 images for the period 1984–2011 and Landsat8 images for the period 2013–2018, all resa-mpled at 30-m resolution. This dataset shows that the 12 detected areas have significant horizontal displacements over the 1978–2018 period, highlighting the presence of active slow-moving landslides. Those landslides are activated later in the Siguas Valley (after 2003,
<figure><img src="image_6.png" /><figcaption>**Fig. 2** | Horizontal ground displacement of different landslides over 40 yr (1978–2018). **a–c**, The slow-moving landslides from the northern Vitor (a), southern Vitor (b) and Siguas (c) valleys. The irrigation onsets in these three different areas are indicated with a black arrow. The width of the bands indicates the uncertainty of the measurements (Supplementary Fig. 15).</figcaption></figure>

Fig. 2c) than in the southern (between 1988 and 1995, Fig. 2b) and northern (before 1984, Fig. 2a) parts of the Vitor Valley. These different initiation times coincide with different onsets of the irrigation programmes in these different areas (Fig. 1a), the irrigation programmes being developed earlier in the northern Vitor Valley (1953), than in the southern Vitor (1974) and Siguas (1983) valleys.

The spatial distribution of the landslides, all situated on the side of the valleys where irrigation programmes were developed, together with the different initiation timing in the two valleys clearly indicates that the landslides are triggered by irrigation. Our long-term dataset provides clear evidence and timing constraints on this previously hypothesized triggering, based on local testimonies¹³.

We also note that the initiations of the landslides all occur about 20 yr after the irrigation onsets (Fig. 2). This observation can be explained by two main factors. First, the groundwater amount reaching the slopes can be small for many years, depending on the soil permeability and the evolution of the water supply with time. Second, the maturation of a shear band at the landslide basal interface is a time-consuming process, caused by progressive damage from subcritical crack growth²⁴,²⁵. The displaced mass is the consequence of a shear band localization, and strain rate acceleration as observed during tertiary creep. The timing between shear band initiation and landslide motion initiation has sometimes been measured on deglaciation-triggered landslides, showing a large range of time between zero and several thousands of years²⁶,²⁷. Modelling suggests that this delay is a strong function of the mechanical parameters of the slope²⁸. For the highly weathered sedimentary rocks found here, a time delay of 20 yr seems therefore reasonable.

## Headscarp retrogression accelerates the landslide motion

The motion rate of those landslides is highly variable, ranging from 50 m to about 400 m over the 1978–2018 period (Fig. 2). All of the landslides display nonlinear behaviour with time, with periods of

quiescence and short periods of faster motions. For instance, the Punillo Sur landslide underwent several periods of acceleration in 1989–1990, 2000–2001 and 2016 (Fig. 2b). All of the landslides present such patterns, but with different timing of the acceleration periods (Fig. 2). This observation highlights that no common regional factor (neither heavy precipitation nor earthquakes) can explain the landslide kinematics.

The random timing of the landslide accelerations raises questions about the controlling factors of these landslides. To investigate these factors, we focus on a major acceleration period of the Punillo Sur landslide occurring in 2016, and retrieve a high-frequency time series of ground displacement²²,²³ using the images from the Sentinel2 satellite launched in 2015 (10 days of revisit time) over the 2015–2017 period. The displacement time series shows that the landslide accelerated suddenly between June and September 2016 to reach 48 ± 1.5 m in 3 months in the southern part of the landslide, 17 ± 1.5 m in its northern part (Fig. 3b), one-sixth of its total displacement over the 1978–2018 period. This spatial pattern of movement is also confirmed by high-resolution SPOT6/7 satellite images (Fig. 3a). The observed transient motion lasts 3 months, with an initiation between 19 June and 9 July 2016 and a velocity peak of 1 m d⁻¹ calculated by the mean velocity between 9 July and 29 July, decreasing slowly to a steady-state end of September 2016.

The initiation and peak of velocity coincide with a retrogression of the headscarp on the rear part of the landslide (Fig. 3a). Analysis of its headscarp, based on both Landsat8 and Sentinel2 images, shows headscarp retrogression lasting between 12 June and 9 July 2016 (Fig. 3b). The total volume of this retrogression, estimated from a difference of SPOT6/7 DEMs in 2014 and 2016, is 1 Mm³ (see Supplementary Fig. 17). The sudden acceleration of the landslide motion occurs just after the initiation of the headscarp activity, clearly indicating that the sediment supply triggered the landslide motion. This interaction between mass movements has already
<figure><img src="image_7.png" /><figcaption>**Fig. 3** | Horizontal ground displacement of the Punillo Sur landslide over the 2014–2017 period. **a,** The displacement field between 2014 and 2016 estimated from SPOT6/7 stereo images. The white and yellow dots mark the location of the two time series of displacement shown in **b**. The area affected by the headscarp retrogression of 1Mm³ in June 2016 is represented by the blue patch. **b,** The white and yellow lines show the time series of landslide displacement over the December 2015–May 2017 period, estimated from Sentinel2 satellites. The error bars are estimated by the mean standard deviation of east–west and north–south components over the whole area. The blue bars show the headscarp retrogression area, which corresponds to a measurement of the loss of material from the cliff at the rear of the landslide.</figcaption></figure>

been observed<sup>14,29,30</sup> and explained by either a mechanism where overloading induces an undrained increase of pore water pressure<sup>31</sup> or/and a mechanism of titling and pushing of the landslide body by the retrogressive headscarp<sup>14</sup>. The singularity of the present case study is the small sediment supply compared to the total landslide volume (total volume of sediment supply is less than 1 Mm³ for a landslide volume of about 75 Mm³), which raises questions about the factors that favour this acceleration.

Unexpectedly, the landslide motion is greater far from the headscarp retrogression, and less just below the headscarp (Fig. 3a). We can also see that the area of stronger velocities coincides with the presence of many ponds, showing the higher humidity of the area. Despite the lack of in situ measurements of pore water pressure, the temporal and spatial pattern of movement suggests that the landslide kinematics could be due to variations of pore water pressure over space and time, causing variations in the effective normal stress at the landslide basal interface<sup>32</sup>. This transient variation of pore water pressure can be created by the pushing of the advancing landslide body induced by the headscarp retrogression<sup>14</sup>. The observed temporal behaviour could then be explained by undrained loading of the landslide body, where pore pressure immediately increases due to the loading by the sediment supply, and water cannot be directly squeezed out due to the low permeability of the porous media<sup>30,31,33</sup>. This poro-elastic effect can also explain the landslide motion in the months after the headscarp retrogression. Indeed, water is progressively evacuated from the loaded porous media, and the water pore pressure returns to its previous steady state, or drained conditions, after a few months. This analysis shows that two factors impact the landslide motion pattern and velocities: the headscarp activity (the frequency of the retrogressions), and the water content in the landslide mass. Differences in the water content over the different landslides can be caused by both differences in water supply and slight differences in lithology.

## Societal implications

Water infiltration from irrigation, due to overwatering of crops, is observed everywhere in the monitored valleys, and many water springs can be detected gushing from the cliffs from impermeable layers (see Supplementary Fig. 3). Those water springs are associated either with slow-moving landslides (as described here) or with rockfalls (see Supplementary Fig. 3), sometimes affecting villages (see Supplementary Fig. 4). These gravitational motions have large impacts on the traditional agricultural crops. The loss of territory due to landsliding between 1978 and 2016 is 1.86 km² that is, 7% of the total valley surface). Furthermore, new large irrigation programmes (Majes-Siguas II) are starting on the eastern side of the Siguas Valley, and new landslides are thus doubtless expected in the near future. We are therefore witnessing the progressive disappearance of traditional agricultural crops in the fertile valleys, due to intensive irrigation on the above plateau. A paradox arises with regard to the irrigation solutions used. Indeed, valley floors are easily irrigated thanks to the presence of a natural river, whereas irrigation on the plateau requires waterworks, drainage and intensive maintenance. Furthermore, some of the canals built are now affected by the retrogression of the cliffs, and even abandoned<sup>34</sup>. This case study therefore highlights the competition between intensive and traditional agriculture in a desert area, where hydraulic stress is predominant, and raises sensitive political questions for local development.

To face this erosion problem, vast cactus fields have been planted on the plateau to farm cochineal (Supplementary Fig. 1b), an insect providing a natural colourant used in the food industry. Those crops require less irrigation than previous lucerne crops, resulting in less water infiltration. This adaptation strategy is however not always economically sustainable, due to large fluctuations in the cochineal world market prices<sup>35</sup>. Despite this agronomic answer, the landslide motions are still going on, showing the long-term impact
of the irrigation on the landslide kinematics. As such this study raises important questions about the long-term impacts of irrigation in desert regions and more generally in regions submitted to poor water management¹⁶.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-019-0500-x.

Received: 9 May 2019; Accepted: 6 November 2019;
Published online: 9 December 2019

## References

1. Matson, P. A., Parton, W. J., Power, A. G. & Swift, M. J. Agricultural intensification and ecosystem properties. *Science* **277**, 504–509 (1997).

2. Reganold, J. P., Elliott, L. F. & Unger, Y. L. Long-term effects of organic and conventional farming on soil erosion. *Nature* **330**, 370–372 (1987).

3. Bakker, M. M. et al. The response of soil erosion and sediment export to land-use change in four areas of Europe: the importance of landscape pattern. *Geomorphology* **98**, 213–226 (2008).

4. Li, X. G. et al. Changes in soil organic carbon, nutrients and aggregation after conversion of native desert soil into irrigated arable land. *Soil Tillage Res.* **104**, 263–269 (2009).

5. Jickells, T. et al. Global iron connections between desert dust, ocean biogeochemistry, and climate. *Science* **308**, 67–71 (2005).

6. Dunai, T. J., López, G. A. G. & Juez-Larré, J. Oligocene–Miocene age of aridity in the Atacama desert revealed by exposure dating of erosion-sensitive landforms. *Geology* **33**, 321–324 (2005).

7. Der Wateren, F. M. V. & Dunai, T. J. Late Neogene passive margin denudation history—cosmogenic isotope measurements from the central Namib desert. *Glob. Planet. Change* **30**, 271–307 (2001).

8. Petley, D. Global patterns of loss of life from landslides. *Geology* **40**, 927–930 (2012).

9. Lacroix, P., Zavala, B., Berthier, E. & Audin, L. Supervised method of landslide inventory using panchromatic SPOT5 images and application to the earthquake-triggered landslides of Pisco (Peru, 2007, Mw8.0). *Remote Sens.* **5**, 2590–2616 (2013).

10. Carretier, S. et al. Slope and climate variability control of erosion in the Andes of central Chile. *Geology* **41**, 195–198 (2013).

11. Sepúlveda, S. A., Moreiras, S. M., Lara, M. & Alfaro, A. Debris flows in the Andean ranges of central Chile and Argentina triggered by 2013 summer storms: characteristics and consequences. *Landslides* **12**, 115–133 (2015).

12. Wilcox, A. C. et al. An integrated analysis of the March 2015 Atacama floods. *Geophys. Res. Lett.* **43**, 8035–8043 (2016).

13. Hermanns, R. L. et al. Landslides in the Andes and the need to communicate on an interandean level on landslide mapping and research. *Rev. Asoc. Geol. Argent* **69**, 321–327 (2012).

14. Lacroix, P., Araujo, G., Hollingsworth, J. & Taipe, E. Self entrainment motion of a slow-moving landslide inferred from Landsat-8 time-series. *J. Geophys. Res. Earth Surf.* **124**, 1201–1216 (2019).

15. Xu, L. et al. Landslides in a loess platform, north-west China. *Landslides* **11**, 993–1005 (2014).

16. Zhang, D., Wang, G., Luo, C., Chen, J. & Zhou, Y. A rapid loess flowslide triggered by irrigation in China. *Landslides* **6**, 55–60 (2009).

17. Zhao, C. et al. Small-scale loess landslide monitoring with small baseline subsets interferometric synthetic aperture radar technique—case study of Xingyuan landslide, Shaanxi, China. *J. Appl. Remote Sens.* **10**, 026030 (2016).

18. Houston, J. Variability of precipitation in the Atacama desert: its causes and hydrological impact. *Int. J. Climatol.* **26**, 2181–2198 (2006).

19. Hesse, R. & Baade, J. Irrigation agriculture and the sedimentary record in the Palpa valley, southern Peru. *Catena* **77**, 119–129 (2009).

20. Thouret, J.-C., Gunnell, Y., Jicha, B., Paquette, J.-L. & Braucher, R. Canyon incision chronology based on ignimbrite stratigraphy and cut-and-fill sediment sequences in SW Peru documents intermittent uplift of the western Central Andes. *Geomorphology* **298**, 1–19 (2017).

21. Leprince, S., Barbot, S., Ayoub, F. & Avouac, J.-P. Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements. *IEEE Trans. Geosci. Remote Sens.* **45**, 1529–1558 (2007).

22. Bontemps, N., Lacroix, P. & Doin, M. Inversion of deformation fields time-series from optical images, and application to the long term kinematics of slow-moving landslides in Peru. *Remote Sens. Environ.* **210**, 144–145 (2018).

23. Lacroix, P., Bièvre, G., Pathier, E., Kniess, U. & Jongmans, D. Use of Sentinel-2 images for the detection of precursory motions before landslide failures. *Remote Sens. Environ.* **215**, 507–516 (2018).

24. Atkinson, B. K. Subcritical crack growth in geological materials. *J. Geophys. Res.* **89**, 4077–4114 (1984).

25. Bruckl, E. & Parotidis, M. Estimation of large-scale mechanical properties of a large landslide on the basis of seismic results. *Int. J. Rock Mech. Min. Sci.* **38**, 877–883 (2001).

26. Ballantyne, C. K. Paraglacial geomorphology. *Quat. Sci. Rev.* **21**, 1935–2017 (2002).

27. ElBedoui, S., Guglielmi, Y., Lebourg, T. & Pérez, J.-L. Deep-seated failure propagation in a fractured rock slope over 10,000 years: the La Clapiere slope, the south-eastern French Alps. *Geomorphology* **105**, 232–238 (2009).

28. Lacroix, P. & Amitrano, D. Long-term dynamics of rockslides and damage propagation inferred from mechanical modeling. *J. Geophys. Res. Earth Surf.* **118**, 2292–2307 (2013).

29. Budetta, P. Rockfall-induced impact force causing a debris flow on a volcanoclastic soil slope: a case study in southern Italy. *Nat. Hazards Earth Syst. Sci.* **10**, 1995–2006 (2010).

30. Booth, A. M. et al. Transient reactivation of a deep-seated landslide by undrained loading captured with repeat airborne and terrestrial lidar. *Geophys. Res. Lett.* **45**, 4841–4850 (2018).

31. Hutchinson, J. & Bhandari, R. Undrained loading, a fundamental mechanism of mudflows and other mass movements. *Geotechnique* **21**, 353–358 (1971).

32. Iverson, R. M. & Major, J. J. Rainfall, ground-water flow, and seasonal movement at Minor Creek landslide, northwestern California: physical interpretation of empirical relations. *Geol. Soc. Am. Bull.* **99**, 579–594 (1987).

33. Sassa, K., Fukuoka, H., Wang, G. & Ishikawa, N. Undrained dynamic-loading ring-shear apparatus and its application to landslide dynamics. *Landslides* **1**, 7–19 (2004).

34. Araujo, E. Evolucion, Dinamica, Implicancias y Monitoreo del Deslizamiento de Siguas, Arequipa. PhD thesis, Univ. Nacional San Antonio Abad del Cusco (2017).

35. Müller-Maatsch, J. & Gras, C. in *Handbook on Natural Pigments in Food and Beverages* (eds Carle, R. & Schweiggert, R.) 385–428 (Woodhead Publishing, 2016).

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2019
## Methods

This study relies on measurements of vertical and horizontal ground displacements from satellite images all available on the earthexplorer (https://earthexplorer.usgs.gov/) and Copernicus (https://scihub.copernicus.eu/dhus/) repositories.

**Topography changes.** The elevation change between 1978 and 2016 was calculated by construction and differencing of two DEMs from stereo optical images. The 1978 DEM was derived from a pair of 6-m-resolution Hexagon (KH-9) images acquired on 16 May 1978. The declassified surveillance satellite images were stored on photographic film and scanned by the United States Geological Survey in two segments at 7 µm (3,600 dpi) with about 1 cm of overlap, requiring heavy preprocessing to detect image fiducial markers, correct for distortions and merge the two image segments³⁶⁻³⁹. During acquisition, images were overlaid with a reseau grid providing the opportunity to reconstruct the image geometry at the time of exposure and correct for distortions created during development, duplication and storage of the film³⁶. The reseau grid consists of 47 × 23 black cross-like markers spaced by 10 mm. The markers were detected in each half-image with subpixel accuracy by convoluting the images with a cross-like pattern⁹. The distortion of each marker to a regular grid was computed, interpolated at each pixel using a thin plate spline interpolation and corrected, providing distortion-free images. Matching features between both image halves were then searched using normalized cross-correlation to stitch together the two scans. The image edges that were not exposed during acquisition were cropped on the basis of the distance to the reseau grid and fiducial markers were filled with white noise to prevent erroneous matches during the stereo processing stage. All of the preprocessing was performed with the OpenCV (https://openCV.org/) and scikit-image (http://scikit-image.org) Python libraries. To estimate the satellite position at the time of acquisition (information that remains classified), we manually identified ground control points in the images and the SPOT DEM. A 24-m-resolution DEM was generated from the stereo pair using the Ames Stereo Pipeline¹⁰, following the methodology of ref. ³⁹.

The 2016 DEM was generated from a stereo pair of SPOT7 images acquired on 9 December 2016, with the method described by ref. ⁴¹ using the Ames Stereo Pipeline software¹⁰. We first orthorectify the images using the ASTER-GDEM with 30-m resolution, to remove the long-wavelength stereo component. Correlation of the resulting ortho-images allows the refinement of the low-resolution DEM. This iterative process also helps to reduce artefacts on steep slopes. Similarly, a DEM was also generated using a tri-stereo pair from SPOT6 images acquired on 24 May 2014. All of the DEMs were aligned¹² on the 2016 DEM, before being subtracted, leading to maps of elevation variations at 6-m resolution between 2014 and 2016, and 24-m resolution between 1978 and 2016. A detriping process¹² is applied to the 1978–2016 difference map to correct for stripes due to errors in the 1978 photographic film scan. The resulting uncertainty for the 1978–2016 elevation change is 4.8 m as estimated by the standard deviation of the difference map over the whole stable area (Supplementary Fig. 9). This uncertainty is similar to those obtained in previous studies³⁶,³⁷,³⁹.

**Time series of horizontal ground displacement.** The KH9 and SPOT6/7 satellite images were orthorectified on their concomitant DEM with respective resolutions of 6 m and 1.5 m. The 2014–2016 displacement field has been produced by correlation of the two orthorectified SPOT6/7 images, using the COSI-corr software²¹.

Then, two different time series of horizontal ground displacement were produced based on two datasets: one long-term (1978–2018) dataset at 30-m resolution using KH9-Landsat5-Landsat8 images over the whole area of study; and one high-frequency, short-term time series at 10-m resolution using Sentinel2 images between December 2015 and May 2017 over the Vitor Valley alone. We first selected images of the Landsat5 satellite with approximately 1-yr separation between 1984 and 2011. A total of 24 images were selected on track 071 in June or July, and 26 images on track 072 in March or April, which together allow a full coverage of the study area. Similar months of acquisition were chosen to reduce the illumination effects that affect the estimation of the ground displacement from optical images¹⁴. Six images from between 2013 and 2018 were then selected on the same tracks from the Landsat8 satellite. The same spectral band 0.35–0.45 µm was chosen for both satellites (that is, band B3 for Landsat5 and B4 for Landsat8) to optimize the similarity between images. All of the Landsat and KH9 images were downscaled at 30-m resolution and stacked to cover the same area. All of the pairs of images from a stack were correlated using MicMac, which provides better results than other correlators (for example, COSI-corr) with those low-resolution images¹³, leading to a horizontal displacement field split into east-west and north-south components. Post-processing of the displacement fields includes a removal of a mean bias between pairs²¹, masking of the points of low quality of the correlation⁴¹, and along- and across-track stripe removal⁴⁵,⁴⁶. Finally, the ground displacement time series was derived following an inversion strategy of the redundant system of displacement pairs²⁵, using a weight between pairs depending on their time intervals. The time series of east-west and north-south displacement are then projected along the mean direction of motion for each pixel to obtain the represented time series from Fig. 2 over the different landslides (see also Supplementary Figs. 10–13). The use of a 30-m dataset is new for ground motion estimation, so we decided to validate the results with in situ observations. A comparison of the displacement time series is realized with 24 markers measured

by the Global Positioning System on the Siguas landslide twice, in 2015 and 2017¹⁴. Satellite time series are first interpolated at the same Global Positioning System dates using a linear interpolation. This last comparison shows that the uncertainty is much lower (and even negligible) for points selected in the centre of the landslide mass than at the borders (Supplementary Fig. 14), which is a classical artefact of correlation methods²². This conclusion drives the selection of the points shown in Fig. 2. The uncertainty of the ground displacement was then estimated with the standard deviation of each east-west and north-south component over the stable areas (Supplementary Fig. 15). We found uncertainties of between 6 and 10 m.

We followed the same strategy (except the across-track stripe removal, non-existent in the Sentinel2 images) to compute a high-frequency ground displacement between December 2015 to May 2017 using 20 cloud-free images from the Sentinel2 satellite at 10-m resolution (Supplementary Fig. 16). Band B2 was used, as already used in a previous study²³. The displacement field uncertainty (estimated as the root-mean-square error) of the medium-resolution satellites has been shown to be close to 1–1.5 m on the Siguas landslide¹⁴, but might differ from one site to the other and also depends on the satellite resolution. The uncertainty of the ground displacement is therefore estimated in a conservative way by the standard deviation of each east-west and north-south component over the stable area (excluding the landslide areas). We found uncertainties of between 0.4 and 2.8 m.

**Geomorphological mapping.** The activity of the Punillo Sur main headscarp has been estimated by a precise mapping of the crest line in each of the available Sentinel2 (10-m resolution multispectral) and Landsat 8 (15-m resolution panchromatic) images in 2016. The individual mapping of the crest line is shown in Supplementary Fig. 18.

## Data availability

The satellite images are available on the earthexplorer (https://earthexplorer.usgs.gov/) and Copernicus (https://scihub.copernicus.eu/dhus/) repositories. Landsat-5, Landsat-8 and Sentinel-2 images are available under https://doi.org/10.5066/F7N015TQ, https://doi.org/10.5066/F71835S6 and https://doi.org/10.5066/F76W992G, respectively. Any additional data can be requested by e-mailing the corresponding author.

## Code availability

The Ames Stereo Pipeline code for DEM processing is available at https://ti.arc.nasa.gov/tech/asr/groups/intelligent-robotics/ngt/stereo/. The code for processing time series of ground displacement from optical images is available via svn (https://sourcesup.cru.fr/scm/?group_id=3177).

## References

36. Surazakov, A. & Aizen, V. Positional accuracy evaluation of declassified Hexagon KH-9 mapping camera imagery. *Photogramm. Eng. Remote Sens.* **76**, 603–608 (2010).

37. Pieczonka, T., Bolch, T., Junfeng, W. & Shiyin, L. Heterogeneous mass loss of glaciers in the Aksu-Tarim catchment (central Tien Shan) revealed by 1976 KH-9 Hexagon and 2009 SPOT-5 stereo imagery. *Remote Sens. Environ.* **130**, 233–244 (2013).

38. Holzer, N. et al. Four decades of glacier variations at Muztagh Ata (eastern Pamir): a multi-sensor study including Hexagon KH-9 and Pléiades data. *Cryosphere* **9**, 2071–2088 (2015).

39. Maurer, J. & Rupper, S. Tapping into the Hexagon spy imagery database: a new automated pipeline for geomorphic change detection. *ISPRS J. Photogramm. Remote Sens.* **108**, 113–127 (2015).

40. Shean, D. E. et al. An automated, open-source pipeline for mass production of digital elevation models (dems) from very-high-resolution commercial stereo satellite imagery. *ISPRS J. Photogramm. and Remote Sens.* **116**, 101–117 (2016).

41. Lacroix, P. Landslides triggered by the Gorkha earthquake in the Langtang valley, volumes and initiation processes. *Earth Planets Space* **68**, 46 (2016).

42. Nuth, C. & Kääb, A. Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change. *Cryosphere* **5**, 271–290 (2011).

43. Rosu, A.-M., Pierrot-Deseilligny, M., Delorme, A., Binet, R. & Klinger, Y. Measurement of ground displacement from optical satellite image correlation using the free open-source software MicMac. *ISPRS J. Photogramm. Remote Sens.* **100**, 48–59 (2015).

44. Berthier, E. et al. Surface motion of mountain glaciers derived from satellite optical imagery. *Remote Sens. Environ.* **95**, 14–28 (2005).

45. Torres, J. & Infante, S. O. Wavelet analysis for the elimination of striping noise in satellite images. *Opt. Eng.* **40**, 1309–1315 (2001).

46. Leprince, S., Muse, P. & Avouac, J. P. In-flight CCD distortion calibration for pushbroom satellites based on subpixel correlation. *IEEE Trans. Geosci. Remote Sens.* **46**, 2675–2683 (2008).
**Acknowledgements**

We acknowledge a careful reading of the manuscript by J. Palmer and E. Berthier. This work has been supported by a grant from ESA through the Alcantara project 'Monitoring and Detection of Landslides from optical Images time-Series' (ESA 15/P26). We also acknowledge the CNES support through the ISIS programme that provided the SPOT6/7 images.

**Author contributions**

P.L. coordinated the study, processed and analysed the image correlation data and the SPOT6/7 DEMs, and wrote the drafts of the manuscript. A.D. processed the KH9 DEM, cross-examined the observations and results, and revised the manuscript. E.T. realized the field measurements and discussed the content of the paper.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41561-019-0500-x.

**Correspondence and requests for materials** should be addressed to P.L.

**Peer review information** Primary Handling Editor(s): Melissa Plail; Heike Langenberg.

**Reprints and permissions information** is available at www.nature.com/reprints.