Remote analysis of an open-pit slope failure: Las Cruces case study, Spain

Slope failures occur in open-pit mining areas worldwide, producing considerable damage in addition to economic loss. Identifying the triggering factors and detecting unstable slopes and precursory displacements —which can be achieved by exploiting remote sensing data— are critical for reducing their impact. Here we present a methodology that combines digital photogrammetry, satellite radar interferometry, and geo-mechanical modeling, to perform remote analyses of slope instabilities in open-pit mining areas. We illustrate this approach through the back analysis of a massive landslide that occurred in an active open-pit mine in southwest Spain in January 2019. Based on pre- and post-event high-resolution digital elevation models derived from digital photogrammetry, we estimate an entire sliding mass volume of around 14 million m3. Radar interferometry reveals that during the year preceding the landslide, the line of sight accumulated displacement in the slope reached − 5.7 and 4.6 cm in ascending and descending geometry, respectively, showing two acceleration events clearly correlated with rainfall in descending geometry. By means of 3D and 2D stability analyses we located the slope instability, and remote sensing monitoring led us to identify the likely triggers of failure. Las Cruces event can be attributed to delayed and progressive failure mechanisms triggered by two factors: (i) the loss of historical suction due to a pore-water pressure increase driven by rainfall and (ii) the strain-softening behavior of the sliding material. Finally, we discuss the potential of this methodological approach either to remotely perform post-event analyses of mining-related landslides and evaluate potential triggering factors or to remotely identify critical slopes in mining areas and provide pre-alert warning.


Introduction
Mass movements and surface deformations, caused by ground subsidence or slope instability processes, are very frequent phenomena in both abandoned and active mining areas (Bell and Donnelly 2006). In open-pit mining, ground movements potentially lead to slope failures entailing risks for personnel, equipment, and infrastructures, in addition to disrupting mine scheduling and increasing production costs (Paradella et al. 2015). Yet the absence of failure in an open-pit mine could be a sign of overconservative slope design (Hoek et al. 2000). Conducting comprehensive slope monitoring and modeling programs is essential for the open-pit mining industry to ensure safety. The literature describes large-scale landslides in open-pit mines whose mobilized mass volumes range from 200,000 m 3 to 65 million m 3 . These studies focus on slope monitoring (Brawner and Stacey 1979;Pankow et al. 2014;Carlà et al. 2018), slope modeling (Voight and Kennedy 1979;Tutluoglu et al. 2011;Ozbay and Cabalar 2014), and slope stabilization (Seegmiller 1979).
In open-pit mining, slope monitoring usually involves in situ techniques such as topographic surveying (total stations/prisms, GPS, etc.) and subsurface techniques (inclinometers, extensometers, etc.) (Vaziri et al. 2010). The application of remote sensing is recent and its use in the mining sector remains scarce. Within this context, structure from motion (SfM) (Ullman 1979) photogrammetry and satellite synthetic aperture radar interferometry (InSAR) (Ferretti et al. 2001) represent powerful alternatives for the remote monitoring of open-pit slopes. The SfM technique allows for automatic 3D reconstruction of surfaces from highresolution imagery by applying digital photogrammetry, i.e., multi-view stereopsis algorithms that do not require background knowledge on the part of the user (James and Robson 2012). Whereas traditional softcopy photogrammetric methods require either established camera positions or the location of ground control points (GCPs) prior to scene reconstruction, the SfM approach requires neither. With SfM, the camera pose and scene geometry are reconstructed simultaneously through the automatic identification of matching features in multiple images (Westoby et al. 2012). Initial estimates of camera positions and object coordinates are refined iteratively using nonlinear least-squares minimization, multiple solutions being available in the image database (Snavely 2008). Over the last decade, SfM has been widely used to study slope instabilities related to different phenomena (e.g., de Bari et al. 2011;Lucieer et al. 2014;Jordá Bordehore et al. 2017;Sarro et al. 2018;Tomás et al. 2018), including mining-related slope instabilities (Thoeni et al. 2012;Chen et al. 2015;Haas et al. 2016;Salvini et al. 2018;Xiang et al. 2018). Satellite InSAR has been applied to monitor mass movements over wide areas with submillimeter accuracy (Crosetto et al. 2016). Although the literature contains many case studies about monitoring mining-induced ground movements using satellite InSAR, most of them focus on ground subsidence (e.g., Carnec and Delacourt 2000;Raucoules et al. 2007;Samsonov et al. 2013;Ng et al. 2017). While groundbased SAR interferometry (GB-SAR) has proven to be operational and effective for monitoring slope stability in open-pit mines (Monserrat et al. 2014), the use of satellite InSAR for this purpose is uncommon, given the decorrelation produced by strong topographical variations in active open-pit mines. Detecting and mapping open-pit slope instabilities by means of ERS and ENVISAT satellites has, however, been documented for abandoned mines in Spain (Herrera et al. 2010). The high spatial resolution and short acquisition frequency of TerraSAR-X satellite have likewise been applied to monitor slope instabilities in an active open-pit mining area in Brazil (Paradella et al. 2015). The improved acquisition frequency of Sentinel-1 satellite mission (with respect to previous C-band missions), made it possible to predict slope failure in an open-pit mine in Turkey through back analyses of InSAR and GB-SAR displacement time series (Carlà et al. 2018).
Slope stability analyses, in turn, have traditionally been performed using limit equilibrium (LE) approaches involving methods such as the ordinary method of slices (Fellenius, Bishop simplified, Spencer, and Morgenstern-Price methods) (Duncan 2003). In recent years, advanced finite element (FE) modeling approaches, e.g., using the shear strength reduction (SSR) technique (Zienkiewicz et al. 1975;Dawson et al. 1999), are gaining attention. SSR simply reduces the soil shear strength, in terms of friction angle Φ and cohesion c, until collapse occurs. The resulting safety factor (SF) is the ratio of the soil's actual shear strength to the reduced shear strength at failure. SSR presents a number of advantages over LE. Most importantly, the critical failure surface is found automatically, i.e., no assumption needs to be made in advance about its shape or location. In addition, it allows straightforward 3D analysis to be undertaken and complex pore pressure scenarios to be defined, including suction and unsaturated states. The contribution of matric suction to shear strength generally results in a significant increase in the factor of safety (Fredlund and Rahardjo 1993). For example, Rassam and Williams (1997) demonstrated that the stability of tailing dams was enhanced by nearly 30% when the matric suction contribution to shear strength was taken into account in the slope stability calculations.
This paper puts forth a methodology that combines SfM photogrammetry, satellite InSAR, and FE modeling to remotely evaluate slope instabilities in mining environments. We illustrate this approach through the back analysis of a large slope failure that occurred in January 2019 at Las Cruces, an open-pit copper mine in southwest Spain (Fig. 1).
The paper is organized as follows. The case study section describes the slope failure, the study area, the mining operation, and the lithology involved. Then, in the methodology section, we provide a description of the methods and datasets used. The photogrammetry section describes the SfM results, while the satellite InSAR section focuses on the results generated using an advanced InSAR processing chain. Further, the correlation of the InSAR results with rainfall is evaluated, and inverse velocity plots are analyzed. The modeling section offers details of 3D and 2D stability analyses of slope failure in view of different geotechnicalhydrogeological scenarios. Finally, in the wake of the remote sensing monitoring and modeling performed, and from a qualitative point of view, the most likely mechanisms triggering failure are discussed.
Case study: Las Cruces slope failure A massive landslide occurred in the early morning of 23 January 2019 at Las Cruces, an open-pit mine located 20 km northnorthwest of Seville (Fig. 1). Fortunately, the failure occurred during a shift change, and no fatalities or injuries were reported (Ecologistas en Acción Sevilla 2019). The landslide affected the north slope of the open pit, leading to the collapse of the north dump. Several pieces of mining equipment and the entrance to an investigation ramp whose construction had been completed the previous year were buried by the moving mass. Las Cruces mine is located in the southeastern part of the Iberian Pyrite Belt (IPB), which forms part of the South Portuguese Zone, the southernmost terrain of the Variscan Belt of Europe (Oliveira 1990). The region is characterized by a fairly simple stratigraphic succession that includes a~2000-m-thick basal siliciclastic sequence made up of shale and quartz arenite, which constitutes the PQ Group of Late Devonian age (Tornos et al. 2017). The transition to the overlying volcanic sedimentary complex (latest Devonian to early Carboniferous), which represents the host rock of Las Cruces deposit, is defined by a discontinuous and highly variable unit comprising mass flows, carbonate reefs, and continental sediments (Moreno et al. 1996). These volcano-sedimentary rocks constitute the Paleozoic basement, located at a depth of 150-200 m at the mine site. The sulfide deposit exploited in Las Cruces is located beneath Cenozoic sediments deposited over the ENE-WSW-trending Guadalquivir basin. This basin was developed during Neogene compression between Africa and Eurasia and was later filled by Miocene-aged marine sediments (Sanz de Galdeano and Vera 1992). Its stratigraphic sequence starts with a few meters of calcarenites hosting the so-called "Niebla-Posadas" aquifer (Miguélez et al. 2011), overlain by a 150-m thick unit of homogeneous, poorly stratified, bluish to gray, shallow marine marl, known as "blue marl" (Tsige et al. 1994), which increases in thickness towards the southeast. Quaternary fluvial sediments overlie these Tertiary deposits.
Las Cruces mine has now been operating for 9 years and is foreseen to run until late 2021 (First Quantum Minerals Ltd. 2015). Production entails six phases and four main waste dumps (north, west, south, and in-pit backfill dumps) (Fig. 2). The monitoring system at Las Cruces integrates piezometers, inclinometers, and prisms (Montero et al. 2009;Cooper et al. 2011). By the end of 2018, fourth phase was almost finished and the exposed ore from fifth and sixth phases was being mined. The open pit is roughly 1500-m long, 900-m wide, and 200-m deep. It is developed in benches 10m high with a 60°slope gradient, except for the top two benches, having a 45°gradient. The overall slope and interramp angles in the failure sector were 28°and 32°, respectively. Montero et al. (2009) established the following lithological levels, from top to bottom, at the mine site: (i) three levels of weathered marl, ranging from highly weathered to moderately weathered, with a total thickness of around 30 m; (ii) four levels of fresh marl, with consistencies ranging from very hard to soft rock and a total thickness of some 100 m; (iii) one level of sandy marl, around 5 m thick; (iv) one level of sand, with thickness ranging from 0 to 15 m; and (v) the Paleozoic bedrock. Figure 2 depicts the layout of Las Cruces, including a crosssection of the north slope displaying the lithostratigraphic sequence (adapted from Montero et al. 2009, Cooper et al. 2011, and First Quantum Minerals Ltd. 2015. Weathered marls are labeled as M1, M2, and M3 and fresh marls as L1, L2, L3, and L4. As the embankment or barrier berms were built to a height of 25 m (First Quantum Minerals Ltd. 2015), we considered a height of 30 m for the north dump. Taking into account the presence of bedding planes dipping 3°to the south within the marl layers (Cooper et al. 2011), an average dip of 3°to the south was adopted for the marls and sands.

Methodology
Structure from motion (SfM) Immediately after the failure, the available literature was collected, as well as data on orthoimages, topography, geology, hydrogeology, and geotechnical properties of the materials involved. Figure 3 shows two orthoimages of the study area before and after the failure.
A post-failure orthoimage (Fig. 3b), corresponding to January 2019, was obtained applying SfM techniques. A set of 204 photos (resolution of 3840 × 2160 pixels) was processed using the software Agisoft PhotoScan Pro (Agisoft LLC 2018). The photos were taken on 25 January 2019 (i.e., 2 days after the failure) from a private helicopter with a conventional digital camera for visual inspection purposes. As no camera positions were available for the images, five GCPs were extracted from a 2013 aerial laser scanning (ALS) point cloud and placed on the model. The 2013 ALS point cloud was directly downloaded from the Instituto Geográfico Nacional (2019). Due to the presence of the moving mass, no GCPs could be placed inside the pit. Only GCPs in areas not affected by the mining operation were selected for referencing. A post-failure DEM (Fig. 4) corresponding to January 2019, was obtained by processing these 204 photos.
In addition, a pre-failure orthoimage (Figs. 2a and 3a) and a set of ten quasi-vertical single-frame aerial images (1 m resolution) covering the study area -both of them corresponding to February 2016were downloaded from the Instituto Geográfico Nacional (2019). A pre-failure DEM (Figs. 5 and 6), corresponding to February 2016, was obtained by processing those ten single-frame aerial images likewise using Agisoft PhotoScan Pro software. Again, no camera positions were available; the GCPs used for referencing the post-failure DEM were also placed on the model, including an additional in-pit GCP placed at the bottom of the pit to ensure a reliable photogrammetric processing result in the vertical direction.
The external boundaries of the slope failure (Figs. 3 and 4) were defined by photointerpretation using the pre-and post-failure orthoimages. The source area for the initial estimation of the moving mass volume was obtained by selecting the negative elevation change values (red to yellow area in Fig. 4a and blue line fill pattern in Fig. 4b). The elevation change was calculated as the difference between the post-and pre-failure DEMs and the initial volume estimate were calculated as the mean elevation difference within the source area multiplied by its surface area. Then, the geometry of the basal shear plane (red dashed line in Fig. 4b) was deduced by comparing the two DEMs and taking into account the lithological profile of the open pit. The initial estimation of the overall volume of the instability was further refined by considering this basal shear plane profile. To this end, a smooth theoretical 3D shear surface was derived from both the basal shear plane profile and the external boundary of the failure using inverse distance weighting interpolation. The source area for the refined volume estimate was delimited by selecting the pixels where the difference between the theoretical 3D shear surface and the pre-failure DEM was negative (red line fill pattern in Fig. 4b). The refined volume estimate was calculated as the mean of this differential value within the recomputed source area multiplied by its surface area.
Advanced differential satellite synthetic aperture radar interferometry (InSAR) Ground deformation velocity maps and time series in Las Cruces were derived by processing Sentinel-1 images using the commercial SqueeSAR® (Ferretti et al. 2011) processing chain (https://site.tre- altamira.com/), which provides ground deformation velocities and time series over the study area. The first one aims to identify areas affected by ground movements, whereas the second one provides information about sudden changes or accelerations in the monitoring period, essential to determine triggering factors through back analysis. Moreover, the processing of Sentinel-1 images through SqueeSAR® provides a good picture of the temporal evolution of the displacements with great spatial coverage, thanks to the short revisit time of the satellite (6 days). Both ascending and descending data were used, covering the period from 4 January 2018 to 17 January 2019. The results were obtained by processing, respectively, 63 and 62 Sentinel-1 ascending and descending images (spatial resolution of 3 × 14 m). Thus, ground deformation velocity maps (Fig. 5) and time series (Fig. 6) could be calculated in the area of the mine for the year prior to failure. The correlation of the deformation time series with rainfall was furthermore analyzed taking into account daily precipitation records from the Agencia Estatal de Meteorología (2019).
Finally, inverse velocity calculations ( Fig. 7) were carried out using the methodology proposed in Dick et al. (2015) to derive considerations on the failure predictability through back analysis. According to the inverse velocity method, the time of failure for an instability in accelerating creep conditions corresponds to the Black line: cross-section profile trace. White box: photogrammetry processing spatial extent. White dashed box: 3D geo-mechanical modeling spatial extent. b Crosssection of the north slope and lithostratigraphic sequence (adapted from Montero et al. 2009, Cooper et al. 2011and First Quantum Minerals Ltd. 2015 point of intersection on the abscissa of the extrapolated trend in a 1/velocity vs. time plot (Fukuzono 1985;Voight 1988Voight , 1989Petley et al. 2002;Crosta and Agliardi 2003;Sornette et al. 2004;Carlà et al. 2018). In this study, the noise of the radar time series was filtered using a 5-day moving mean; inverse velocity values over ± 200 d/cm were deleted to improve the visualization. Inverse velocity was calculated using the following equation: where iv i is the inverse velocity at the date of interest t i , t i − 1 is the previous date, d i is the accumulated displacement at the date of interest t i , and d i − 1 is the accumulated displacement at the previous date t i − 1 .

Numerical modeling
Stability analyses were performed with an FE modeling approach using the SSR technique implemented in the GeHoMadrid code (Fernández-Merodo 2001;Mira 2002). Governing equations describing the hydro-mechanical coupling correspond to the static version of the Biot equation, wherein partial saturation is taken into account by assuming pore air pressure p a equal to zero, and taking the relation between effective stress σ' and total stress σ in view of a simplified version of the Bishop effective stress σ ' = σ + S r • p w m (Schrefler 1984): The following notation was used: -S: matrix defined with the kinematic strain-displacement rela- u: solid displacements p w : pore-water pressure -S r (s): degree of saturation, defined as a function of suction s through the soil water retention curve (SWRC). In that case p a = 0 and suction s = p a − p w = − p w (suction is positive for negatives values of p w ) -ρ: density of the mixture defined as ρ = (1 − n)ρ S + nS r ρ w , where n is the porosity, ρ S is the density of the solid grains, and ρ w is the water density b: external forces (gravity) The mechanical behavior of the materials was represented via an elastoplastic constitutive model in the form of an extended Mohr-Coulomb failure criterion, where the shear strength τ of unsaturated soils is defined as: and where c ′ is the effective cohesion for saturated states, σ 0 n and σ n are respectively the effective and total normal stress, ϕ ′ is the effective friction angle, and c is the apparent cohesion which includes the friction term due to suction: For the SSR algorithm, failure condition was defined using a non-convergence criterion of the Newton-Raphson nonlinear algorithm with a normalized unbalanced force limit of 10 −4 (Fernández-Merodo et al. 2014;Bru et al. 2018).
The analyzed geometry corresponds to the pre-failure DEM matching the February 2016 open-pit configuration obtained through SfM techniques (see Sect. 3.1). The precise lithology was taken into account using the layout defined by Montero et al. (2009) and plotted in Fig. 2.
The geo-mechanical properties of the thirteen considered materials were taken from Montero et al. (2009) andFirst Quantum Minerals Ltd. (2015) and are listed in Table 1. Friction angle and cohesion values correspond to drained effective strength parameters of intact materials (peak conditions). The blue marls behave as an over-consolidated clay with brittle failure behavior. Softening behavior should be taken into account and will be addressed in the discussion section. The bedding planes and ubiquitous joints with very low strength parameters (ϕ ′ = 12°a nd c ′ = 10 kPa) described by Cooper et al. (2011) were not explicitly introduced in the present global stability analysis.
The soil water retention curve (SWRC) of the blue marls provided by Boza (2014) (Fig. 8) was adopted for partial saturation. Note the large air-entry value (AEV) that provides nearly saturated conditions for values of suction up to 1 MPa.
Water pressure distribution inside the pit, within the blue marls, at the time of failure is a key aspect behind the failure mechanism. Nevertheless, it is an extremely difficult task to arrive at this quantitative determination. For one, the blue marls are fairly impermeable. The permeabilities obtained by Montero et al. (2009) range from 10 −9 to 10 −11 m/s, suggesting significant undrained behavior and a very slow consolidation process. Secondly, water pressure distribution depends on the historical hydrological boundary conditions since the beginning of the first excavation phase. Montero et al. (2009) presented a coupled hydro-mechanical analysis to predict the significant pore pressure drop in the marls due to the high undrained unload caused by the pit excavation. A calculated drop up to 1 MPa matches the piezometric data but does not reveal any pore suction detail above or below the water table. Moreover, long-term pore pressure dissipation since the excavation would depend on the artificial drainage systems installed and on the natural seasonal fluctuations at surface (evaporation and rain). Here, given the difficulties in modeling realistic historical hydrological conditions, and in order to show how even minor changes in surficial suction can influence slope stability, three different simplified water pressure distributions were considered. The first consists of a hydrostatic configuration prescribing a pore pressure on the surface equal to − 300 kPa trying to emulate a long steady-state condition with a 30 m ground water table depth (the same as the water table level measured in natural conditions prior to the pit excavation). The second one consists of a less favorable hydrostatic configuration prescribing a pore pressure on the surface equal to − 150 kPa (15m-deep water table). The third is a hypothetical worst-case scenario prescribing zero suction on the surface (fully saturated case).
For the other boundary conditions, all displacements were fixed in the bottom plane; perpendicular displacements were fixed in the lateral planes. For the applied external load, only gravity force was contemplated.
Stability was addressed in three successive analyses: -First, a simplified 3D stability analysis over the open pit was performed in order to locate the most critical slope due solely to the slope geometry of the open pit. For this purpose, we adopted a box size of 1860 m × 2200 m × 340 m centered over the open pit (white dashed box in Fig. 2). To simplify the model, (i) we considered a single, homogeneous material representative of the blue marls (averaged material properties are taken from units M and L in Table 1, ϒ average = 15 kNm-3, Φ' average = 20°, c' average = 250 kPa), and (ii) only the worst assumed pore pressure scenario (fully saturated case) was taken into account. The 3D finite element mesh, made up of 227,202 nodes and 149,325 quadratic tetrahedral elements (10 nodes per element), and the adopted hydrostatic pore pressure contours are depicted in Fig. 9.

A' A
-Secondly, detailed 2D stability analyses (plane-strain conditions) were performed over the recognized critical section (north slope), considering precise lithology (Table 1) and evaluating the three proposed simplified pore pressure distributions. The finite element mesh, with 8493 nodes and 4122 quadratic triangular elements (6 nodes per element), is shown in Fig. 10a. The three adopted hydrostatic pore pressure contours and equivalent saturations are shown in Fig. 10c and d.
-Finally, a stability back analysis was undertaken to compute mobilized friction angles Φ' mov in the basal shear plane for the three proposed simplified pore pressure distributions. In this case, the critical slip surface is not computed, but is introduced as a thin layer (deducted from the SfM analysis) in the finite element mesh (see Fig. 10b). The effective mobilized friction angle in this layer (with zero cohesion) is computed through back analysis to obtain the specified factor of safety SF = 1 (limit of stability condition).

Remote-sensed topography and geomorphic mapping
The post-failure orthoimage (Fig. 3b) was used to conduct, along with the pre-and post-failure DEMs, a detailed geomorphic analysis of the landslide (Fig. 4). We classify the slope failure as a complex movement (Varnes 1978) with a well-defined exposed shear plane clearly present in the post-failure DEM (Fig. 4a), which suggests that a primary failure initiated somewhere near the north dump barrier berm, triggering subsequent failures towards the dump. The entire moving mass consists of two distinct parts: (i) an upper part with the typical features of a rotational landslide such as head scarp, secondary scarps, and sagged crown and (ii) a mudflow-earthflow lower part, with debris on the frontal section. The external boundaries of the slope failure (red lines in Figs. 3 and 4a) were used to determine the landslide surface area (785,000 m 2 ). The pre-and post-failure DEMs obtained were then used to estimate failure volumes. The post-failure DEM yielded 0.97 m resolution and 0.86 m total error, whereas the pre-failure DEM resulted in 1.65 m resolution and 1.3 m total error. The difference between the two DEMs yielded a volume of 11.4 million m 3 , with a mean height difference of 31.4 m and a source area of 364,000 m 2 . Based on the difference between the theoretical 3D shear surface and the pre-failure DEM, estimation of the overall volume of the landslide gave 14 million m 3 , which is~1.2 times the initial estimates. The recalculated source area was 3% larger (375,000 m 2 ), and the average sliding surface depth is 6 m deeper (37.4 m).

Ground deformation maps and time series analysis
Following SfM analysis, we elaborated line of sight (LOS) velocity maps (Fig. 5). The results in ascending geometry (Fig. 5a) (Fig. 2). Once more, the south

Recent Landslides
Landslides 17 & (2020) dump exhibits the maximum negative velocity, with a value of − 23.0 cm/year; but in contrast to the ascending results, the map in descending geometry shows positive LOS velocities (movement towards the satellite) over the north slope, with a maximum velocity of 7.7 cm/year. The opposite direction of the movements measured along the two satellite LOSs over the north slope indicates that the actual displacements have a strong horizontal component. Note that the higher LOS velocities in the south dump do not necessarily imply that real deformation velocities are lower in the north slope. Since the velocities are measured along the satellite LOS, measured percentages of the real displacements depend on the main direction of movement with respect to the LOS. In the north slope, real displacements have a strong horizontal component, which is consistent with the direction of the slope, oriented southeast with a slope gradient of around 30°. For this reason, the LOS vector partially detects displacement over this slope in both ascending and descending geometry. In the south dump, the measured percentage of the real displacements is higher due to the  predominantly vertical component of the actual displacements (related most probably to consolidation processes), determined by the relatively flat geometry of the dump. In view of the resulting displacement time series, we conducted a zonal analysis for the ascending and descending geometries (Fig. 6). Five areas were selected (processing plant, north slope, south dump, north dump 1, and north dump 2 in Fig. 6a), their mean deformation velocities shown in Fig. 6c. As seen in Fig. 6b, accumulated displacements over the north slope reached − 5.7 and 4.6 cm in ascending and descending geometry, respectively, from January 2018 to January 2019. In ascending geometry, the trend shows an acceleration in mid-November 2018, and a deceleration in early January 2019. Yet in descending geometry, two acceleration events are detected in mid-January 2018, and early October 2018, in addition to two deceleration events in late March 2018, and mid-December 2018. All other areas exhibit negative accumulated displacements both in ascending and descending geometry, except for the processing plant (Fig. 6d), which shows stable trends for comparison. The south dump (Fig. 6e) is found to have an acceleration event from early March 2018 to late May 2018, a deceleration event from early September 2018 to late October 2018, and a linear trend in between. These acceleration and deceleration events would respectively correspond to the initial and final stages of the consolidation process induced by the enlargement of the dump. Finally, the north dumps ( Fig. 6f and g) show almost linear settlement trends, most likely related to the consolidation process as well; they resemble the trend of the north slope in ascending geometry (Pearson correlation coefficient between 99 and 100%), indicating a potential relationship between the deformation processes that affect the north slope and the north dumps. The events observed over the north slope in descending geometry moreover correlate with accumulated precipitation (Pearson coefficient 99%), from two rainy episodes from early January to late April 2018, and from early September to late December 2018 (Agencia Estatal de Meteorología 2019). Contrariwise, the acceleration event marked by the ascending results in the north slope -slightly delayed with respect to the last acceleration event observed in descending geometry-is poorly correlated with rainfall (Pearson coefficient 90%). To conclude, in the south and north dumps, the results show an even lower correlation with rainfall (Pearson coefficient between 64 and 83%).
Finally, an analysis of the displacement time series (both in ascending and descending geometry) using the inverse velocity approach was undertaken (Fig. 7). The inverse velocity plots indicate noisy but stable trends in the south and north dumps (Fig. 7d, e, and f), although this also applies to the plot corresponding to the north slope in ascending geometry (dark gray line in Fig. 7a). In contrast, in descending geometry, the inverse velocity plot of the north slope shows a distinctly nonlinear trend (light gray line in Fig. 7a). According to the definition of acceleration event, the two events detected in the north slope in descending geometry, as well as the acceleration observed in the north slope in ascending geometry, can be classified as temporary accelerations. Hence one must proceed with caution, in this case, when attempting to predict the time of failure using the inverse velocity method, as it assumes that the slope displaces at an accelerating rate up to the point of failure. Considering the period from 1 September to 31 December 2018 -i.e., the period corresponding to the rainfall episode prior to the failurelinear regressions of the inverse velocity plots of the north slope in ascending and descending geometry (black dashed lines in Fig. 7b and c) are of low quality (R 2 values of 0.03 and 0.28, respectively). However, linear regression of the inverse velocity plot of the north slope in ascending geometry provides, in retrospect, an excellent prediction of the time of failure (R 2 = 0.84) when considering the period from 7 October to 31 December 2018 (red dashed line in Fig. 7b). At any rate, this prediction makes sense only in the context of back analysis, as there was also a temporary acceleration in mid-January 2018 over the north slope according to the descending results, but no failure was triggered.

Stability analysis
Results of the 3D stability analysis are displayed in Fig. 11. The 3D stability analysis gave a SF equal to 1.20, indicating that all the open-pit slopes were stable for the assumed pore pressure scenario (fully saturated case). The computed critical slope was located on the northwest side of the open pit, i.e., on the north slope; plastic strain developed along a shear surface, the mobilized mass sliding along this shear surface with a typical rigid body motion. The computed failure source area (light blue polygon in Fig. 11a) is estimated to amount to 314,000 m 2 . It was verified that the other two proposed simplified pore pressure distributions indicated a similar failure mechanism.
The results of 2D stability analysis are presented in Fig. 12. For the best-case scenario, i.e., − 300 kPa pore pressure on the surface (first row in Fig. 12), the computed SF is equal to 1.28. In the intermediate scenario, − 150 kPa pore pressure on the surface (second row in Fig. 12), the computed SF is equal to 1.14. Finally, in the worst-case scenario, that is 0 kPa pore pressure on the surface (third row in Fig. 12), the computed SF is equal to 1.00. For this 2D stability analysis, three meshes were evaluated in terms of accuracy and consumed CPU time. Figure 13 illustrates the influence of the mesh size (safety factor SF and consumed CPU time vs. number of degrees of freedom d.o.f.). In order to achieve a good balance between accuracy and consumed CPU time, MESH 2 was selected for all the analyses. Finer meshes such as MESH 3 did not improve the results substantially.
In the stability back analysis, the mobilized friction angle (Φ' mov ) in the basal shear plane at failure gave 18°, 21°, and 29°, respectively for the best-case, intermediate, and worst-case pore pressure scenarios. The methodology, relying exclusively on the exploitation of remote sensed and published data, proves useful to analyze geomorphic features, failure volumes, source areas, precursory displacements, accelerations related to rainfall, triggering factors, and failure mechanisms. Regarding its limitations, the two main drawbacks of this methodology surround the uncertainty associated with available data, which can only be solved using in situ data for validation (surveying data, borehole data, inclinometers, piezometers, etc.), and the data availability itself. For this reason, the purpose of the stability analyses presented herein is essentially to discuss potential triggering and contributing factors, rather than to accurately quantify failure conditions. A landslide source area of 375,000 m 2 and a failure volume of 14 million m 3 were estimated using SfM techniques. The location and extent of the landslide source area as computed through 3D stability analysis (314,000 m 2 ) coincide fairly well with the SfM estimates. Both the SfM-and the modeling-derived results were based on a pre-event DEM from 2016 (obtained through the application of SfM techniques). Despite not representing in detail the actual 2019 pre-event geometry, it served for analysis of Las Cruces north slope, where no substantial topographical variations were observed between 2016 and 2019.
InSAR results highlight the presence of sub-decimeter rainfallcorrelated precursory displacements with a clear horizontal component on the north slope the year prior to failure. Owing to the temporary nature of the observed accelerations, predicting the time of failure by means of the inverse velocity method provided no conclusive results. Still, an excellent prediction of the time of failure (R 2 = 0.84) was obtained in retrospect, in ascending geometry, when considering the period from 7 October to 31 December 2018.
In light of the FE modeling results, it is remarkable how the simplified 3D stability analysis (Fig. 11) was able to reproduce and exactly locate the failure mechanism that acted in January 2019 (Figs. 3b and 4a). Moreover, computed displacement contours for the detailed 2D stability analysis (Fig. 12a) reveal a failure mechanism equivalent to that deduced from our SfM analysis (Fig. 4b).
A comparison of the simplified 3D and detailed 2D stability analyses also leads to noteworthy observations. For the proposed worst-case scenario (saturated case), the computed failure mechanisms are found to be similar, with respective SFs equal to 1.20 and 1.00. This finding is in line with those of previous research confirming that SF 3D > SF 2D (Duncan 2003), and can be explained by the precise lithostratigraphic geometry description adopted in the 2D model. A sand layer dipping towards the slope and near its toe, with weak resistance properties, may have had a key influence on the reported failure event. Indeed, as seen in all rows of Fig. 12b, plastic strain initiates in this material and progresses until reaching the contact between shales and tuffs. It should be noted that the bedding planes dipping towards the slope (Cooper et al. 2011) -i.e., in the worst direction from a stability standpointwere not considered in the present analysis yet would have contributed to a more unstable scenario.
As anticipated, results of the 2D stability analyses show that water pressure distribution inside the pit bore a very strong influence on the slope stability. Tensile pore stresses (suction) near the slope surface increase soil strength through apparent cohesion (see Eq. (4)) and increase slope stability. We should underline that large air-entry values (AEV) of the SWRC and almost fully saturated conditions give rise to a nearly linear increase of the friction term due to suction in Eq. (4). In the wake of their coupled hydromechanical analysis of Las Cruces pit mine, Montero et al. (2009) affirmed that "lower pore pressure distribution due to the volumetric expansion associated with the excavation of the pit allows for a more aggressive and economical slope design". In any event, quantification and prediction of suction states must be addressed very carefully in such designs. Considering the proven correlation between the accumulated displacements over the north slope as measured by InSAR techniques and the accumulated rainfall, one may envisage as a first triggering factor the loss of historical suction resulting from a pore-water pressure increase driven by rainfall. Nevertheless, given the computed safety factor SF equal to 1 for an unrealistic worst-case hydraulic scenario (no suction with 0 kPa on the surface), a complementary triggering factor may have been at work.
Because of our current knowledge of the blue marls as an overconsolidated clay with brittle failure behavior, softening behavior should be taken into account. Progressive strain-softening behavior, in full agreement with the accumulated displacements evidenced in our InSAR data, could lead to a progressive failure mechanism. According to Alonso and Gens (2006) and Gens and Alonso (2006), the mechanical strain-softening behavior of blue marls led to the outstanding progressive failure of Aznalcóllar dam foundation, just 10 km away (Fig. 1)  mobilized friction angle of the sliding basal plane of 17-18° , a value between the average peak (Φ p = 24.1°) and residual (Φ r = 11°) friction angles measured by Gens and Alonso (2006) in drained direct shear tests with constant vertical effective stresses in the range of 100-800 kPa. Aznalcóllar residual strength parameters nearly coincide with Las Cruces' bedding planes strength parameters provided by Cooper et al. (2011). As with other sliding failures in brittle clays (Stark and Eid 1997), a residual factor close to R = 0.5 seems to apply in the case of Aznalcóllar dam; as defined by Skempton (1964), R measures the degree of development of residual strength along the failure surface and is a useful parameter for characterizing progressive failure: where τ is the average shear stress at collapse (directly linked to computed mobilized friction angle) and τ p and τ r are the average peak and residual strengths at the prevailing normal effective stress along the eventual rupture surface. At Las Cruces, the performed stability back analysis gives a mobilized friction angle for the sliding basal plane equal to 18°, 21°, and 29°, respectively for the proposed best, intermediate, and worst pore pressure scenarios. These values mean approximated residual factors R equal to 0.65, 0.49, and 0.04. Apparently R is very sensitive to pore pressure and would indicate that the intermediate pore pressure scenario (− 150 kPa on the surface) is more representative of failure. Therefore, according to our stability analyses, Las Cruces failure could be explained by a delayed and progressive failure mechanisms triggered by two factors: (i) the loss of historical suction resulting from a pore-water pressure increase driven by rainfall, and (ii) the strain-softening behavior of the blue marls. These two triggers come in full agreement with the InSAR-derived accumulated displacements over the year prior to failure.
The precise and quantitative determination of pore pressure distribution inside the pit at the time of failure remains an extremely difficult open question. From a modeling point of view, it can only be addressed through real time-dependent analysis, able to characterize the undrained behavior of the blue marls and taking into account the excavation process history, the long-term pore pressure dissipation, and the natural seasonal fluctuations at the surface boundary. Such analysis, beyond the scope of the present paper, has been successfully achieved in the past in the context of cut slopes in stiff clays, explaining the delayed collapse of railways and motorways cutting slopes (Potts et al. 1997).
This work can be viewed as novel support for the application of remote sensing when monitoring and modeling slope failures in open-pit mines in those cases where data availability is limited. Although the methodology was applied here to perform post-event back analyses of mining-related landslides, it might be extrapolated to remotely identify unstable slopes and provide pre-alert warning. New remote sensing tools stand as a new paradigm to improve the design of prevention and mitigation measures in open-pit mines.