Spatiotemporal slope stability analytics for failure estimation (SSSAFE): linking radar data to the fundamental dynamics of granular failure

Impending catastrophic failure of granular earth slopes manifests distinct kinematic patterns in space and time. While risk assessments of slope failure hazards have routinely relied on the monitoring of ground motion, such precursory failure patterns remain poorly understood. A key challenge is the multiplicity of spatiotemporal scales and dynamical regimes. In particular, there exist a precursory failure regime where two mesoscale mechanisms coevolve, namely, the preferred transmission paths for force and damage. Despite extensive studies, a formulation which can address their coevolution not just in laboratory tests but also in large, uncontrolled field environments has proved elusive. Here we address this problem by developing a slope stability analytics framework which uses network flow theory and mesoscience to model this coevolution and predict emergent kinematic clusters solely from surface ground motion data. We test this framework on four data sets: one at the laboratory scale using individual grain displacement data; three at the field scale using line-of-sight displacement of a slope surface, from ground-based radar in two mines and from space-borne radar for the 2017 Xinmo landslide. The dynamics of the kinematic clusters deliver an early prediction of the geometry, location and time of failure.


Symbols and their meaning
The magnitude of the relative displacement of two grains i and j (or two pixels) linked in  Total flow leaving the source node q G Directed network G * Undirected link-capacitated network H(X(t)) Entropy of the clustering assignment X(t) I(X(t); X(t − 1)) Mutual information between X(t) and X(t − 1) i, j Grain (pixel) in N which corresponds to a physical contact e k Sink node l Link in T l m Link of minimum weight on the path joining nodes u and v in T LOS Line-of-sight ℓ Observation point (grain or pixel) NMI Normalized mutual information N Physical contact network n Number of nodes in N , or G * O (n − 1) × 3 matrix whose entries are all zeros O [i, : ] i th row of O p Pixel with the highest rate of movement in at the time of failure PFR Precursory failure regime, mesoregime q Source node R min Radius of the smaller of the two contacting particles S Silhouette score for the whole body s(i) Silhouette score for a node i SSR Slope stability radar T Final time stage T Gomory-Hu tree t Time stage t F , t F B , t F 1 , t F 2 , t F X Time of failure for a general case, Biax, Mines 1, 2 and Xinmo, respectively t * , t * B , t * 1 , t * 2 , t * X Regime change point for a general case, Biax, Mines 1, 2 and Xinmo, respectively V Set of nodes v, u Node in G * W, W ′ Two disjoint components of N corresponds to the removal of edges in B(F ) or removal of a link l in T w(l) Non-negative weight on link l in T w(l m ) Capacity of the minimum cut separating u and v in N X (t) Clustering assignment at time stage t x(e) Flow of link e Natural and engineered slopes, composed of granular materials like rocks, concrete and soil, can maintain their structural integrity even as damage spreads. But there is a tipping point, beyond which damage can propagate to cause catastrophic failure with little to no apparent warning signs [1][2][3] . The landslides in Xinmo (China, 2017) and the dam collapse in Brumadinho (Brazil, 2019) are recent reminders of the devastating impact of slope failure on human lives and livelihoods, infrastructure, and the environment [1][2][3][4][5] . Here the term failure is used from an operative point of view and is the moment when the slope totally or partially collapses, displaying a paroxysmal acceleration and a disintegration of the mobilized material. A critical frontline defense against these hazards is large-scale monitoring and analysis of slope movement using remote sensing technologies 1,3,6-12 . Some of these measurements have now reached spatial and temporal resolutions (e.g., Slope Stability Radar 7 ) which enable direct connections to be made to the fundamental deformation and failure of granular materials 1,3,6,11 . Nevertheless, there are significant challenges to overcome before the full potential of these data assets can be harnessed for geotechnical risk assessment and hazard management 11,12 . One of, if not, the biggest challenge lies in the analysis and interpretation of monitoring data with respect to the underlying micromechanics and dynamics of deformation in the precursory failure regime (PFR) [6][7][8][9][10] . Here we address this challenge by formulating a holistic framework for spatiotemporal slope stability analytics for failure estimation (SSSAFE). SSSAFE is physics-based and bears explicit connections to the micromechanics and dynamics of ductile to brittle failure in granular solids (e.g., [13][14][15][16] and references therein). A hallmark of SSSAFE is its detailed characterization of the spatiotemporal coevolution of the preferred pathways for force and damage in PFR using kinematic data. As highlighted in various reviews 3,6,[9][10][11][12] , scant attention has been paid to the spatiotemporal dynamics of landslide deformation, with existing approaches in landslide forecasting and early warning systems (EWS) falling into one of two categories: (a) spatial analysis of an unstable slope to estimate the location and geometry of a landslide 17 , or (b) temporal analysis of ground deformation of single measurement points exhibiting tertiary creep, to deliver a short-term forecast of the time of failure 1,10,18 . The former partially relies on expert judgment (e.g. the choice of the failure criterion and the method of analysis 6  www.nature.com/scientificreports/ and on in situ data (depth of the lithologies and of the water table, resistance parameters of the rock or soil) that always bear a certain level of uncertainty and representativeness bias 19 .
In temporal analysis, the inverse velocity (INV) theory originally proposed by Fukuzono 20 is the most widely applied method for prediction of the time of collapse in the terminal stages of PFR. This approach has no spatial aspect and depends on assumptions which motivate areas for improvement in forecasting 1,7 . Being based on the inverse value of a derivative parameter, this method is heavily affected by noise, especially when the velocity is not particularly high. While this can be addressed by smoothing the data using a moving average 10,18 , this comes at the cost of diminished sensitivity to important changes in the acceleration trends due to short-terms events, including: surface boundary conditions (e.g., civil engineering and mining works 6 ), variations in the trigger factors of slope instability (e.g., rainfall 2 , seismic 21 , mining blasts 7,8 ), and the inherently complex mechanical interactions between different parts of the slope. Concurrent sites of instability may also interact and induce stress redistributions that lead a landslide to "self-stabilize" [22][23][24] . Efforts 1,7 to improve the INV approach give prima facie evidence to suggest that more accurate forecasts can be achieved when the spatial characteristics of slope displacements are incorporated in the temporal analysis of monitoring data.
Accordingly, recent work focused on the spatiotemporal evolution of landslide kinematics in PFR in two case studies using: (a) ground-based radar data of a rockfall in an open pit mine (Mine 1) where two sites of instability emerged, leading one to self-stabilize before the larger one collapsed; and (b) satellite-based Sentinel 1 radar data (Xinmo) of the catastrophic collapse in Xinmo, which led to 83 fatalities 4,24-28 . Guided by lessons learned from the physics and dynamics of granular failure, these delivered a reliable early prediction of the location and geometry of the failure region 4,24,25,27,28 , as well as regime change points in PFR 4,25,26,28 . In this study, we build on these efforts to develop a holistic data-driven framework which eliminates the uncertainties associated with a postulated stress-strain model for the slope, yet holds explicit connections to the first principles of fracture and failure mechanics of heterogeneous and disordered granular solids (e.g., [13][14][15][16] and references therein). To do this, we adopt a transdisciplinary approach which integrates network flow theory of granular failure [13][14][15][16] and mesoscience [29][30][31] . Given the novelty of this formulation from several fronts, the next section gives a brief review of the relevant developments which, woven together, form the basis of SSSAFE.
Precursory dynamics of granular failure across system levels and scales (a) Preferred paths for transmission of damage versus force. In complex systems, not all paths for transmission are created equal. Some are preferred over others. Experimental studies into the transmission of force and energy in natural and synthetic granular media (e.g., sand, photoelastic disk assemblies) and associated discrete element simulations have shown that the mesoregime of PFR is governed by the coupled evolution of two dominant mechanisms [13][14][15]32,33 . The first comprises the preferred paths for force transmission (mechanism A): a set of system-spanning paths that can transmit the highest force flow along direct and shortest possible routes through the system. Distinct force chains (Fig. 1) can be readily observed to form along these percolating paths, in alignment with the major principal stress axis [34][35][36] . The second are the preferred paths for damage (mechanism B), where cracks and/or shearbands emerge. Note that the term damage is broadly used here to mean the separation of two grains in contact, bonded or unbonded.  Arguably the best manifestation of the coupled evolution between force and damage can be observed in deforming photoelastic disk assemblies 33,36 . Here one can readily observe forces continually rerouted to alternative pathways as damage spreads (Fig. 1). This scenario is similar to traffic flows where vehicles are diverted to alternative routes when a road is closed off for repairs or other incidents. Following this analogy to road networks, grain contact networks similarly give rise to emergent flow bottlenecks. Prior network flow studies have shown that these sites, which are highly prone to congestion, ultimately become the preferred paths for damage in the failure regime 13,14 . Counter to intuition, however, the bottlenecks do not generally coincide with the location of damage sites in the nascent stages of PFR, which is when ideally predictions should be made to allow enough time to enact mitigative measures. Instead, a process which can be described as a compromise-in-competition between the preferred paths for force and damage develops, which effectively shields the bottleneck from damage. Specifically, force congestion in the bottleneck is relieved by complex stress redistributions that redirect forces to other parts of the sample, where damage can be accommodated with minimal reduction to the system's resistance to failure (or global force transmission capacity). This may explain why current failure detection methods, which rely on damage sites in the early stages of PFR for spatial clues on where catastrophic failure ultimately forms, sometimes suffer high false positive rates in laboratory 37 and field levels 11 .
(c) The principles of mesoscience. To account for the compromise-in-competition among preferred transmission paths and simultaneously 'jump scale'-from laboratory to field-we integrate the network flow approach 13,14 with the principles of mesoscience (Table 1, Fig. 2). Pioneered by Li and co-workers [29][30][31] in the area of chemical and process engineering, mesoscience has enabled the upscaling of models of gas/solid-particle flow systems from laboratory to industrial scale. Mesoscience is predicated on the concept of a compromise-in-competition between at least two dominant mechanisms in a so-called mesoregime of a complex system. In the simplest case of two competing mechanisms (A and B), the mesoregime mediates two limiting regimes; each is governed by one dominant mechanism, A (B) in the A-dominated (B-dominated) regime, which is formulated as an extremum. Li et al. 29 argues that, while the classical single objective optimization formalism applies to each limiting regime, Table 1. Mechanisms underlying the strength and failure of granular materials compete in the precursory failure regime (PFR).

Regime (emergent structures)
A-Preferred paths for force (force chains) B-Preferred paths for damage (cracks, shear bands) . This chart is analogous to the mesoscience perspective depicted for gas-or solid-particle flow systems 31 . www.nature.com/scientificreports/ the compromise-in-competition in the A − B mesoregime necessitates a multiobjective optimization approach. Results from prelude studies 13,14 , employing a dual objective network flow analysis, corroborate this view. Moreover, opposing trends manifest as the system evolves from one limiting regime to the other ( A → A − B → B and vice versa), consistent with the mesoscience principles (Fig. 2, Table 1). In laboratory tests where detailed analysis of underlying mechanisms are possible, the B-dominated failure regime is characterized by bursts to a peak in all the indicators of stored energy release and dissipation, including: kinetic energy, dissipation rate, population of buckling force chains and their supporting 3-cycles, average values of local nonaffine motion, grain velocity and rotation 16,[38][39][40] . By contrast, at the opposite extreme, in the A-dominated stable regime, all of these quantities are negligibly small. In PFR, these opposing tendencies compromise and give rise to spatiotemporal dynamical patterns 4,[25][26][27] .
(d) Clustering patterns in the kinematics characterize the mesoregime PFR. The compromise-in-competition between force transfer and damage paths in PFR gives rise to collective motion or kinematic partitions: mesoscale clusters where constituent members move collectively in near rigid-body motion 4,[25][26][27] . Interestingly, Li and coworkers also observed particle clusters in the mesoregime of gas/solid-particle flow systems, and conjectured that these emerge from particles tending to minimize their potential energy, while the gas tries to choose a path of least resistance through the particle layers [29][30][31] . Analogously, in the systems studied here, damage favors the path of least resistance to failure -the path that forms the common boundaries of kinematic clusters 13,14 .
(e) Dynamics of kinematic clusters provide early prediction of failure across scales. A complex network analysis of individual grain motions in sand in laboratory tests 41 and of surface ground motion in a slope (e.g., Mine 1) 24 has shown that the impending failure region develops in between subregions of transient but high kinematic similarity early in PFR. Moreover, the spatiotemporal dynamics of these clusters can deliver a reliable change point t * from which such partitions become incised in the granular body, giving rise to their near relative rigid body motion: for example, when the active 'slip region' of a slope begins to detach and accelerate downslope from a relatively stationary region below; or when parts of a rock mass on either side of a developing crack undergo relative slip. That is, persistent partitions in kinematics space forewarn of impending partitions in physical space 4,[25][26][27] . In a parallel effort 28 , the computational challenges of embedding knowledge of kinematic clustering in a stochastic statistical learning model from high-dimensional, non-stationary spatiotemporal time series data were overcome, with displacement and velocity trends and the failure region of Mine 1 successfully predicted more than five days in advance.
(f) Establishing a connection to first principles fracture and failure mechanics for granular solids. Relative motions at the grain-grain level were used to study the coevolution of force and damage propagation in a network flow analysis -with explicit connections to the most popular fracture criteria, starting with Griffith's theory for crack propagation 13,14 . The emerging flow bottlenecks for force and energy, proven to be the paths of least resistance to failure, were found to deliver an accurate and early prediction of the location of shear bands and macrocracks that ultimately develop in the failure regime. The question that now arises is: Can a combined mesoscience and network flow approach detect the bottlenecks and kinematic clusters from radar-measured surface ground motion data and, if so, how can their spatiotemporal evolution be used to deliver an early prediction of a likely place and time of failure?
Here we answer this question and demonstrate our approach through SSSAFE. Based solely on kinematic data for input, SSSAFE first applies the network flow model to identify and characterize the emerging kinematic clusters in PFR, and then uses their dynamics to deliver an early prediction of where and when failure is likely to develop. Different from recent past work [24][25][26][27] which adopt an essentially pattern-mining approach, SSSAFE rigorously predicts the path of least resistance to failure in a manner consistent with the fundamental failure micromechanics and dynamics across different system levels and scales. Four systems are analyzed: a standard laboratory test (Biax); and three rock slopes, man-made slopes Mine 1 and Mine 2 and a natural slope Xinmo. The input kinematic data to SSSAFE comprise individual grain displacements in Biax, and radar line-of-sight displacement data gathered from ground-based radar (Mine 1 and Mine 2) and space-borne radar (Xinmo).

Data
The input data to our analysis consist of the following system properties at each time state t = 1, 2, . . . , T : (a) coordinates of observation points ℓ = 1, 2, . . . , L ; (b) displacement vector recorded at each point d 1 We have four data sets (Fig. 3). The first is Biax, a well-studied simulation of granular failure in a standard laboratory test in which an assembly of polydisperse spherical grains is subjected to planar biaxial compression 16,38-40 . Here each point ℓ is a moving grain and the vector d ℓ is two-dimensional. The sample begins to dilate at around t = 50 . Collective buckling of force chains initiate at around t = 98 , giving way to a brief period of strainsoftening and the development of a single shear band along the forward diagonal of the sample. This shear band becomes fully formed at t = 104 , referred to as the time of failure t F B (Fig. 3a). From this point on, the sample exists as two clusters, in each of which constituent grains move collectively as one: two 'solids' in relative rigidbody motion along their common boundary, viz. the shear band. Details of this simulation and mechanisms underlying its bulk behavior in the lead up to and during failure are provided elsewhere 16,[38][39][40] .
Three large field scale data sets, corresponding mainly to the pre-failure regime of three rock slides, are examined 22 . In these case studies, for ease of illustration of SSSAFE's practical implementation in EWS, we refer to t = 1, 2, . . . , T as the monitoring period, but keep in mind that this may constitute only a subset of the total time period of the actual monitoring campaign that was undertaken for these slopes. Mine 1 and Mine 2 are from monitoring data of a rock slope in two different open pit mines using ground-based SSR-XT -3D real aperture radar 7,8 (Fig. 3b,c). The mine operation, location and year of the rock slides are confidential. However, we have all the information needed for this analysis. Each observation point ℓ is a grid cell or pixel, ranging in size from www.nature.com/scientificreports/ 3.5m x 3.5m to 7m x 7m, in a fixed grid. The vector d ℓ is 1D, which corresponds to the displacement along a line-of-sight (LOS) between the radar and the point ℓ on the slope surface. Mine 1 is a slope constructed from dumped loose rock, derived from blasted and/or excavated overburden (Fig. 3b). The monitored domain stretches to around 200 m in length and 40 m in height. Movements of the rock face were monitored over a period of three weeks: 10:07 May 31 to 23:55 June 21. Displacement at each observed location on the surface of the rock slope was recorded at every six minutes, with millimetric accuracy. This led to time series data from 1803 pixel locations at high spatial and temporal resolutions for the entire slope. A rock slide occurred on the western side of the slope on June 15, with an arcuate back scar and a strike length of around 120 m. Mine 1 reached peak pixel velocity of around 640 m/yr. Considering a precautionary correction for the radar line-of-sight, this falls in the moderate velocity category 42 and corresponds to an evacuation response 22 . The time of collapse t F 1 occurred at around at 13:10 June 15, close to when the global average peak velocity of 33.61 mm/hr was reached. There is a competing slide: a second region of instability, to the east (encircled area, Fig. 3b). This region intermittently developed large movements, but the instability was somehow arrested and movement slowed down the day before the collapse of the west wall 24,25 . In this context, this region is sometimes referred to as a false alarm in the sense that it did not eventuate into a collapse 11 . While in many cases "tertiary creep" ends with a total or partial failure, it is also possible, like in Mine 1, that the whole landslide or a part of it finds a new equilibrium 10,22 . There are many possible reasons for this, such as a reduction of the destabilizing forces through stress redistributions or the geometric configuration of the sliding surface, which slows down and ultimately arrests the whole or part of the landslide body.
Mine 2 is a rock slope of an open cut mine dominated by intact igneous rock that is heavily structured or faulted by many naturally occurring discontinuities (Fig. 3c). A slope stability radar scanned the section of the rock face for displacement for approximately 6 days from 15:39 August 19, until 07:05 August 25, each scan taking approximately 6 minutes, again with millimetric precision. Measurements at 5394 pixel locations were taken every 6 minutes giving high spatial and temporal resolution for the entire domain, measuring 1280 m wide and around 224 m high. A rock slide occurred on the Southern wall at 03:00 August 25; we refer to this as the time of failure www.nature.com/scientificreports/ t F 2 for the rest of this paper. The area that failed, measuring approximately 135m wide and 145m high, moved over 1 million tonnes of debris. Mine 2 reached peak pixel velocity of 2.8 m/day, which is classified as rapid 42 . The Xinmo landslide (Fig. 3d) can also be classified as a rock slide-debris avalanche 22. It is composed of metamorphic sandstone intercalated with slate, which detached on June 24, 2017 and hit Xinmo village (Maoxian, China, 32 • 03 ′ 58 ′′ N, 103 • 39 ′ 46 ′′ E) causing 83 deaths and destroying 64 houses. The analyzed data set is focused only on the original source area that was located near the crest of the mountain ridge north of Xinmo village, at an altitude of 3431 m a.s.l.. As this source moved along the slope, it entrained new rock material and reached an estimated volume of 13 million m 3 and a terminal velocity of 250 km/h 43 . The site was not actively monitored at the time, but displacement data obtained from Sentinel-1 constellation, that takes periodical interferometric acquisitions of the area, have been retrospectively analyzed to determine if a forewarning would have been possible 44 . The data used consisted of 45 SAR images in C-band (6.5 cm wavelength), at 5 m × 14 m spatial resolution, acquired along the descending orbit (incidence angle of 40.78) and spanning from 9 October 2014 to 19 June 2017 (that is five days before the failure). The pixels are of size 5m × 14 m. Data, covering an area of 460 km 2 , were elaborated with the SqueeSAR algorithm 45 and comprised more than 130,000 measurement points. Xinmo reached peak pixel velocity of around 27 mm/yr in the tertiary creep phase, which is very slow 42 but later reached a terminal velocity of 250 km/h during failure 43 .

Method
The core components of our proposed spatiotemporal slope stability analytics for failure estimation (SSSAFE) framework are summarized in Table 2 and Fig. 4. The key idea is to model the transmission of force in each studied system in a way that accounts for the coupled evolution of the preferred pathways for force and damage, and to use this model to predict the emerging kinematic clusters in the mesoregime PFR. To achieve a consistent formulation across different system levels, we model force transmission as a flow through a network. At the core Table 2. Combined mesoscience and network flow formulation behind SSSAFE. A compromise-incompetition between mechanism A (preferred paths for force) and mechanism B (preferred paths for damage) governs the mesoscale at the laboratory level and field level.  www.nature.com/scientificreports/ of this formulation is a set of optimization problems on a network in accordance with network flow theory and mesoscience principles. We emphasize that our implementation of this model is confined only to finding the preferred paths for damage which represent the common boundary of the kinematic clusters. Detection and characterization of the preferred paths for force are outside the scope of this investigation. Such paths have been characterized for different laboratory samples, including concrete (e.g., 14,15,33 ).
Core components of SSSAFE. The core components of SSSAFE are implemented in three consecutive steps, with Steps 2 and 3 delivering respectively the predictions on the likely location and time of failure. Recent work 13 shows explicit connections between the formulation below and the most popular fracture criteria of fracture mechanics, beginning with Griffith's theory for crack propagation. The method proposed here, consistent with Griffith's theory, was found to provide the most accurate and early prediction of failure in PFR for a range of ductile to quasi-brittle laboratory samples.
Step 1: Construct the flow network F . Forces are transmitted along physical connections. Hence, the construction of the flow network F begins with an undirected network N that represents the physical connectivity of the system: the grain contact network in Biax; the proximity network in the three slopes, Mine 1, Mine 2 and Xinmo, where pixels within a distance d of each other are connected. Each node of N represents a grain (pixel), while each link in N represents a grain-grain contact (pixel-pixel connection). The links in N vary with loading history in Biax, but is fixed across the monitoring period for the three slopes.
Next, N is transformed to a directed network G = (V , A) , where V, A are the set of nodes and set of arcs respectively. That is, each link connecting nodes i ∈ V and j ∈ V in N is represented by a pair of symmetric arcs e ∈ A : one from i to j and another vice versa. Given this symmetry, we will use the symbol e to also denote a link. Every link in G is then assigned a non-negative capacity c(e) which corresponds to the maximum flow value that the given link can support. Since the model concerns force transmission, c is thus the force that must be overcome to break the connection: the strength or resistance to failure of the grain-grain contact or pixelpixel connection. Given that what is measured reliably in both laboratory and field levels is motion -not forces or stresses -we express this capacity c in terms of the motions of the connected elements. Hence, the contact capacity function c is given by where | − − → u ij | is the magnitude of the relative displacement of two grains (or two pixels) linked in N . Note that since we are only interested in the flow bottleneck 14,46 , what is important in this analysis are the relative values of the link capacities and not their absolute values. Indeed, for this purpose, the model for the link capacity need not be in units of force, as previously shown (e.g., 14,15,33 ). Consequently, in Equation (1), we set the capacity to be such that the higher the relative motion of grains (pixels) linked in N , the less stable is the connection and in turn the lower is its corresponding capacity c (see Fig. 5). Finally, a direction for the flow is dictated by a pair of artificial nodes q and k called the source and the sink of G. The quadruple F = (G, c, q, k) is called a flow network. In the Biax sample, the natural choice for the source q and sink k are the top and bottom walls so that the direction of flow is in alignment with the direction of the applied vertical compression (and major principal stress axis) of the sample. For Mines 1 and 2 and Xinmo, however, there is no obvious choice for the source-sink pair to direct www.nature.com/scientificreports/ the flow since the loading conditions are unknown for these slopes. To address this problem systematically, we construct the Gomory-Hu tree (GHT) 47 for the network N , as explained in Step 2.
Step 2: Find the kinematic clusters from the bottleneck of F . The bottleneck of F , B(F ) , is given by the cut of F with the least capacity. Any cut of F , Ŵ , is a set of links in N which, if disconnected, represents a literal cut of F into two disjoint components {W, W ′ } of V such that no flow can be transmitted from source q ∈ W to sink k ∈ W ′ . Thus, any cut Ŵ contains all arcs emanating from a node in W and terminating on a node in W ′ .
Physically, a cut Ŵ may be thought of as a virtual crack of the studied granular body or domain whose connectivity is described by N . Physical disconnection of the contacts associated with the links in Ŵ would thus result in a literal system-spanning crack which splits the body into two disjoint pieces. The capacity of Ŵ is defined as c(Ŵ) = e∈Ŵ c(e). Following Equation (1), this represents the total force flow that must be overcome to disconnect every link in Ŵ.
Here we are interested in finding that cut with the least capacity -the so-called minimum cut, also known as the bottleneck B(F ) . Thus, the capacity of the bottleneck B(F ) represents the global failure resistance, F * : the minimum amount of force flow needed to overcome the resistance of the connected links B(F ) to break apart and split the granular body into two disjoint pieces. Note that this analysis does not preclude a body from splitting apart into more than two pieces: in such cases, one can repeat the same analysis described here for each piece to obtain further subpartitions. In the cases studied here, this is unnecessary as the studied systems split apart essentially into two components with the bottleneck being their shared boundary. In Biax, the bottleneck B(F ) predicts the location of the shear band that forms in the failure regime. In the case of Mines 1 and 2 and Xinmo, B(F ) predicts the boundary of the landslide. As time to failure draws near, we expect motion in the components to become increasingly coherent and near-rigid-body resulting in kinematic clustering. The active cluster in PFR, denoted by , distinguishes itself by manifesting an increasing downward motion (viz. increasing trend in cumulative displacement and velocity) due to gravity, while the stable cluster remains relatively stationary.
To find the bottleneck of Biax at each time, we solve the Maximum flow -Minimum cut (MFMC) problem on F , following earlier work 14,33 . This is a two stage calculation. Stage 1 solves the Maximum flow problem to find the global flow capacity, F * , the maximum flow that can be transmitted through N given its topology and link capacities. More formally, given a flow network The flow value that solves the above is the maximum flow F * . Once F * is established, we move to Stage 2 to solve the Minimum Cut Problem.
The Minimum Cut Problem of F = (G, c, q, k) is the cut Ŵ min such that The above is typically solved using the Ford-Fulkerson algorithm 46 . This exploits the well known max-flow mincut theorem which states that the maximum flow possible F * is the capacity of the minimum cut or bottleneck 48 . Using this theorem and Equations (1) -(5), we can now directly relate the conditions on where and when catastrophic failure occurs to the bottleneck. That is, catastrophic failure occurs when the force flow exceeds the resistance to breakage of all the links in the bottleneck. Where the system physically breaks apart is given by the bottleneck itself, which is the limiting shared boundary of the kinematic clusters.
As earlier mentioned in Step 1, the uncontrolled and unknown loading conditions of Mine 1, Mine 2 and Xinmo pose a challenge for the assignment of an appropriate source and sink nodes to direct the flow. We address this issue by constructing the Gomory-Hu tree (GHT) 47 for the network N . In what follows, we outline this procedure briefly for completeness. Full details are described elsewhere 49 . Let G * = (N, c) be an undirected, linkcapacitated network. For every pair of nodes u, v in G * , the GHT stores information on the minimum u-v cut of G * that separates u and v. For a network with n nodes, there are a total of n(n − 1)/2 possible source-sink pairs, each with a corresponding minimum cut. However, the construction of the GHT shows that the minimum cuts for some pairs of nodes are identical. In fact, the GHT contains information on exactly n − 1 distinct minimum x(e). www.nature.com/scientificreports/ cuts 47 corresponding to a set of n − 1 explicit source-sink pairs. One could infer all the remaining implicit source-sink pairs from this set using the GHT, as illustrated in Fig. 6. Formally, the GHT is defined as follows. In Fig. 6, we illustrate an example contact network N with n = 9 pixels, its corresponding Gomory-Hu tree T and a table summarizing the outcome of removing a link in T . There are 36 possible source-sink pairs. T contains 8 explicit source-sink pairs (column 1, Fig. 6c). Removing link l, connecting nodes u and v in T , gives two distinct components W, W ′ : these correspond to the kinematic clusters of N when the edges in the minimum cut separating the source-sink pair u and v are removed. All other source-sink pairs and their corresponding minimum cuts can be inferred from T. Observe this global minimum cut is biased towards highly imbalanced cuts where one component is significantly smaller than the other in terms of the number of member nodes. Such highly imbalanced partitions may correspond to the smaller component having only one or at most a few pixel locations out of thousands or more. As such, this may not provide a complete summary of emerging partitions that lead to catastrophic failure. Larger partitions that span the system, where the part that dislodges from the rest of the slope constitutes a sizeable portion of the slope, are of interest. Accordingly, we introduce a cut ratio ρ which is defined as the number of nodes in the smallest to largest component upon removal of a link in T . Hence, in this example, if we are interested in the smallest component containing at least 3 pixels, we find the minimum cut such that 0.3 ≤ ρ ≤ 1 . This yields the cut that corresponds to the removal of link (2,5) in T ; the explicit source-sink pair is (u = 2, v = 5) as before. By inspecti n g T , w e c a n s e e t h a t s o u r c e -s i n k p a i r s correspond to For Mines 1 and 2 and Xinmo, the number of nodes in G * and possible source-sink pairs are, respectively: 612, 186, 966; 5394, 14,544,921; 610, 185,745. It is thus computationally expensive to enumerate every source-sink pair. However, we can use T to capture a failure event such that the minimum cut identifies a failure area that is no smaller than a prescribed fraction of the studied domain size. That is, we can remove the link in T with the minimum weight such that 0.3 < ρ ≤ 1 to obtain two corresponding cluster components W, W ′ . It is easy enough to check the partitions, and ensure that parts of landslide boundaries that are close to the boundary of the monitored domain (recall the top left boundary of the rockfall in Mine 1 (Fig. 3b) are also captured. This is done simply by checking the cases where 0.03 ≤ ρ ≤ 0.3 , which identify partitions that lead to the smaller cluster being as small as 3% of the number of nodes in the larger cluster (Fig. 7). In summary, the key output from Step 2 is the bottleneck B(F ) , the path of least resistance to failure which separates the active cluster that will likely collapse from the rest of the slope, and the slope's failure resistance F * .
Step 3: Characterize the cluster dynamics. As depicted in Fig. 4, at each time state up until the current time t, we find the flow bottleneck B(F ) , its associated clusters and the failure resistance. We can use this historical information to characterize the dynamics of the cluster motions as the monitoring advances in time. Here we are interested in one of the defining aspects of granular failure, namely, collective motion. As time advances towards the failure regime, we quantify the extent to which: (a) intracluster motions become increasingly coherent and similar-at the same time as intercluster motions become more and more different (separated in kinematic state space); and (b) the predicted clusters no longer change in member elements, suggesting that the pattern of impending failure has become physically incised in the system. To do this, we compute the Silhouette score S 50 to quantify the quality of the clustering pattern obtained from the network flow analysis, coupled with an information-theoretic measure of Normalized Mutual Information (NMI) 51 to quantify the temporal persistence of the clustering pattern.
The Silhouette score S ∈ [−1, 1] gives an overall measure of the quality of clustering 50 . It is the global average of s(i) which measures how similar is a given node i to the other nodes j in its own cluster (cohesion) compared to the nodes in the other clusters (separation): here a(i) is the average distance in the displacement state-space from i to all other nodes in the same cluster, and b(i) is the average of the distances from i to all points in the other cluster. As shown in Figure 8, a good clustering pattern (high S) is one where the nodes in the same cluster exhibit very similar features (nodes are tightly packed in feature state space hence small a(i)); while nodes from different clusters have very different features (red nodes are well separated from blue nodes in feature state space hence large b(i)). As a general guide, values below 0.2 suggest essentially no clustering pattern was found, while the closer S is to one, the more compact are the individual clusters while being more separated from each other. Given the studied feature is motion, an here I(X(t); X(t − 1)) is the mutual information between X(t) and X(t − 1) and H(.) is the entropy of the corresponding clustering assignments. NMI ∈ [0, 1] : 0 means there is no mutual information, as opposed to 1 where there is perfect correlation or similarity, between the clusters at t and t − 1 . Intuitively, NMI measures the information that the clustering assignments X(t) and X(t − 1) share: the higher the NMI, the more useful information on the clustering pattern is encoded in X(t − 1) that can help us predict the clustering at the next time state X(t).
In summary, based on the results from Steps 2 and 3, we can identify a regime change point t * from which the failure resistance F * drops close to its minimum of zero, as S rises and/or levels above 0.2, while NMI stays close to 1. For all t ≥ t * , a prediction on the landslide region is given by , the active or fastest moving cluster. Thus, in general, the mere fact that a landslide experiences deformation does not translate into a prediction of the likely location of an impending failure from SSSAFE. A clear regime change point t * must be identified. That said, in the event that monitoring commenced after t * , such that the active and fast moving cluster no longer changes during the studied time states, then we can expect S to remain high above 0.2, while NMI stays close to 1. In other words, even though no further regime change point may be detected over the studied time states, the temporal persistence of high values for both S and NMI should still serve as a good indicator that is a high risk area that is prone to failure.
In addition, for landslides exhibiting tertiary creep deformation, the time of failure t F can be predicted by performing a linear regression analysis with a continuous and overlapping rolling time window for the inverse mean velocity of for t ≥ t * . Thus with respect to detecting the time of failure t F , the main advance achieved here is that SSSAFE not only obviates the need to subjectively select a pixel to implement the Fukuzono INV analysis 20 but also ensures the INV analysis takes into account the spatiotemporal and coupled evolution of force and damage pathways in PFR. While a comprehensive SSSAFE-INV analysis of different landslides manifesting tertiary creep deformation is outside the scope of this work, we nevertheless illustrate in the next section an implementation of this strategy for detecting t F from a few selected time windows for each slope studied here.

Results and discussion
In all of the systems studied, SSSAFE uncovers three dynamical regimes over the course of the monitoring campaign, consistent with a compromise-in-competition between force and damage (Figs. 9-12). In Biax, the global mean velocity steadily rises in PFR, before a sudden burst to a peak in the failure regime (Fig. 9a). Simultaneously, the opposite trend can be observed in the time evolution of the system's resistance to failure F * , which decreases progressively as damage spreads in PFR, eventually dropping to its minimum value close to zero at stage 80 (Fig. 9b). Extensive published studies of this sample has shown that columnar force chains at stage 80 have lost considerable lateral support in the region of impending shear band, due to dilatancy 39,40,52 . While force redistributions around force chains continually occur during this period (recall Fig. 1), ultimately, the degradation in the region precipitates collective force chain buckling at the peak stress ( t = 98 , Fig. 9a inset), culminating in a fully developed shear band at t F B = 104 , when kinematic clusters move in almost relative rigid-body motion (recall Fig. 5a). These events were previously observed in various types of sand and photoelastic disk assemblies 15,33,41 . A consistent dynamics emerges in the time evolution of NMI and S in Fig. 9c. From stage t * B = 80 , S rises from around 0.5 before levelling off at the start of the failure regime at t F B = 104 ; NMI stays close to 1 from t * B = 80 . These trends imply a recurring bottleneck B(F ) as evident in 80 ≤ t ≤ 110 of Fig. 9d, such that grains on either side progressively move collectively as one in opposite directions (Fig. 5a).
At the field scale using radar data, SSSAFE delivers qualitatively similar trends for Mines 1 and 2 and Xinmo. The presence of large fluctuations in Mines 1 and 2 (Figs. 10-11a-c) is not surprising given these mines were operational with blasting, pumping, transport and drilling works taking place at various times over the course of the monitoring period. Like in Biax, the failure resistance of Mine 1 drops close to zero as early as around t * 1a = 69 (Figs. 10b), with corresponding rises in NMI and S towards 1 (Figs. 10c), even though the rock slope (H(X(t))H(X(t − 1))) ; www.nature.com/scientificreports/ appears intact with near-zero global mean velocity (Fig. 10a). Note that we present this failure resistance in log scale to ease identification of a subsequent second drop which constitutes a second regime change point at t * 1b = 3322 = 12 : 14 June 14. The rapid drop of the failure resistance at t * 1a = 69 suggests that internal cracks and shear bands have started to propagate internally, as delineated by the preferred paths of damage (black points in Fig. 10d). Damage spread, possibly initiating even before the time period covered by the data analyzed here, is to the extent that the capacity for force transfer between adjacent connected material points (pixels) along the recurring bottleneck B(F ) to the west, has significantly reduced as early as t * 1a = 69 , even though there are still many remaining connections in the rock slope that keep it manifestly intact: Supplementary Movie Mine1 clearly shows that the west wall persists as an active area from the beginning of the monitoring campaign.
As time advances towards failure over the period 69 ≤ t < 3322 (Fig. 10b-c), the interaction between the two regions of instability leads to an initial decline in S while NMI stays close to 1 due to the persistence of the west wall cluster, the site that eventually collapses. But the day before the collapse, S sharply rises from t * 1b . This rise in S suggests that the clustering pattern has now become incised in the slope to the extent that the clusters are now essentially undergoing relative motion along their common boundary, as accelerates [24][25][26][27] . The change point t * 1b improves on earlier work using a pattern mining approach which detects the time of imminent failure to be one to two hours later: t = 3350 = 14 : 53 June 14 26 and t = 3333 = 13 : 16 June 14 25 . Results from the SSSAFE-INV analysis using the mean velocity of for t ≥ t * 1b (Fig. 10e) corroborates the change point t * 1b ; here is the actual time of failure t F 1 minus the predicted time of failure from the SSSAFE-INV analysis. To support this, and in keeping with current pixel-based INV analysis 10,18 , we also add corresponding results from analysis of the velocity of the fastest moving pixel p in (Fig. 10e-inset). Recall that the novelty of SSSAFE-INV lies in the establishment and use of in an INV analysis for t ≥ t * , as noted in the earlier section on the core components of SSSAFE. The three linear regression fits (from three time windows) in Fig. 10e serve only to illustrate this procedure. We envisage that the implementation of SSSAFE-INV in practical EWS would involve the use of rolling and overlapping time windows, since this process allows for continuous updates in the prediction of t F 1 , concomitant with incoming new data on slope displacement.
In Mine 1, multiple sites of instability interact mechanically in PFR. Our method can reliably identify and differentiate these regions (Fig. 10). The west wall where catastrophic failure occurs can be distinguished early in PFR by the temporal persistence of the predicted landslide boundary (black points) in this area, in contrast to the eastern corner where this boundary only occasionally appears (Supplementary Movie Mine1). This intermittent dynamics in the latter is due to redundant force pathways which the system exploits to relieve the build up  www.nature.com/scientificreports/ of stress in B(F ) along the west wall by diverting the forces and damage there to alternative paths, including to the competing slide to the east. In laboratory samples undergoing quasi-brittle failure 13 , the interaction between competing cracks manifest in the form of stress redistributions along the preferred force pathways between the bottleneck where the macrocrack ultimately forms and the competing crack which later undergoes structural arrest (self-stabilize). The same is observed in Mine 1: note the concentration of black points in the area between the actual failure region to the west and the competing slide to the east in 0.3 < ρ ≤ 1 , 69 ≤ t < 3322 of Fig. 10d. This compromise-in-competition continues until all such paths are exhausted, t = t * 1b = 3322 , from which time B(F ) remains fixed and becomes primed for uncontrolled crack propagation, along the landslide boundary ( 3322 ≤ t ≤ 3600 , Fig. 10d). Mine 1 provides a good example of why early prediction of failure rests crucially on methods that can account for the spatiotemporal compromise-in-competition between force and damage pathways. Essentially, the ultimate effect of stress redistributions is to delay failure, since any damage to B(F ) leads to a reduction in F * . But there is an undesired concomitant which is the considerable uncertainty they pose for early prediction of failure, given damage is rerouted and concentrated elsewhere -away from the region of impending failure in PFR 13 .
In Mine 2, the rock slope is dominated by intact igneous rock that embodies many natural joints or faults. There is the A-dominated stable regime over the first day of the monitoring period 1 ≤ t < 221 , where the global mean velocity fluctuates initially around 0, as F * portrays a decreasing trend (Fig. 11a,b). Trends in both NMI and Silhouette coefficient S suggest that no substantial clustering structure in the kinematics developed on the first day: NMI fluctuates between 0 and 1 while S remains below 0.5 (Fig. 11c). The system embodies redundant pathways to divert stresses away from the area of impending failure ( 1 ≤ t < 221 , Fig. 11d). In the A-B dominated mesoregime of PFR, S progressively increases to 1, implying the emergence of collective motion ( 221 ≤ t ≤ 1370 , Fig. 11d). As failure draws near, intracluster motions become coherent and near rigid-body, while intercluster motions become separated (Fig. 11c inset), as the cluster corresponding to the location of impending failure accelerates (Fig. 11e). We see these trends are precisely mirrored by the NMI of the clusters (Fig. 11c). Note that the landslide boundary, shown at t = 221 in Fig. 2 actually appears as early as t = 104 and persists up until t = 178 , which explains the high NMI scores. However, the kinematic clusters undergo a short period of change during 179 ≤ t < 221 which may reflect any number of perturbations on the mine site, including blasting. Around the same time interval, large fluctuations can also be observed in S. Close to and during the B-dominated stable regime, S flattens out close to 1, indicative of a strong clustered motion. Altogether, the evidence from F * , S and NMI marks a regime change point at t * 2 = 221 = 13 : 39 August 20, which is just over 4 days prior to the collapse on t F 2 = 1315 = 03 : 00 August 25. Results from the SSSAFE-INV analysis also supports the progressive evolution to collapse at t F 2 (Fig. 11e). For the Xinmo landslide, two key trends are evident from SSSAFE: (a) the regime change point on t * X = 26 = August 23, 2016 which is 10 months in advance of the actual time of collapse from when the active area became fixed in the location that later became the rock avalanche source 11 (Fig. 12a-e); and (b) the development of cracks in a smaller competing failure zone above and to the east of the actual rock avalanche source (black diamonds in 0.03 ≤ ρ ≤ 0.3 , 1 ≤ t < 26 of Fig. 12e, Supplementary Movie Xinmo) as early as 2015. Thus results here corroborate prior findings 1,4,11,43,53 that the rock avalanche was generated from the upper part of the slope (Fig. 12a inset), near the mountain crest, and only afterwards did it entrain old landslide deposits along its way, which showed no sign of movement and played a role only in the propagation and not in the triggering phase. In the months immediately following t * X , the velocities recorded were slow 42 and predisposed to self-stabilization 22 . Interestingly, the shift in from the eastern to the western flank of the source area match the reconstruction proposed by Hu et al. 53 , who attribute the triggering of the rock avalanche to an initial rockfall on June 24, 2017 at 05:39:07 (local time) and impacting the bedrock in the source area, where a crack network was present but still locked by persistent rock bridges. While the long creep behavior of the source area demonstrates that the failure was the result of a long-term process, instead of the sudden outcome of an external impulse, it is credible that the trigger rockfall hypothesized by Hu et al. 53 could have imposed the final stresses needed to overcome the failure resistance (capacity c in Equation 1) of the remaining connections in the recurring bottleneck B(F (t)), ∀t ≥ t * X . That B(F ) persisted in the same location from August 23, 2016 strongly suggests a progressive degradation in rock strength all along this path, with the antecedent prolonged rainfall 43 likely aiding this condition and rendering B(F ) increasingly poised for uncontrolled crack propagation in the lead up to the failure event on June 24, 2017 (Fig. 12b). For the SSSAFE-INV analysis, we select only two time windows to illustrate this procedure since there are fewer time states in Xinmo compared to Mines 1 and 2. Around 40 days prior to the collapse, began to manifest a linear trend in the temporal evolution of its inverse mean velocity, which delivered a prediction for t F X that is a day later than the actual collapse (� = −1). SSSAFE offers a lead time of a day to weeks. This constitutes sufficient forewarning to undertake evacuation and other response actions 54 . SSSAFE takes only a few tens of seconds per time state to generate predictions on the likely region of failure: 30 seconds for Mine 1 and Xinmo, and 50 seconds for Mine 2, on a standard laptop computer with 8 cores 1.30 GHz CPU. Thus a prediction can be returned before the next measurement even for the most advanced radar technology (e.g., 1-5 minutes). At this rate, a reasonable number of time states (e.g., 30 consecutive time states would take at most 30 minutes) to establish robustly the dynamics of the region of interest for the purposes of identifying t * and t F . Finally, we note the following limitations of SSSAFE. In its current form, SSSAFE is designed to support slope stability monitoring efforts by identifying: (a) the likely landslide location , (b) a regime change point t * from which can be confidently established, and (c) for landslides manifesting tertiary creep in t ≥ t * , the likely time of failure t F . As SSSAFE builds on knowledge of the precursory dynamics of catastrophic failure in granular systems, the problem of identifying if and when such failure can be expected for landslide types where the active site moves slowly for hundreds of years and rarely fail catastrophically 55  We demonstrate how SSSAFE can be applied to identify emergent kinematic clusters in the early stages of the precursory failure regime for four case studies of catastrophic failure: one at the laboratory scale using individual grain displacement data; and three slopes at the field scale, using line-of-sight displacement of a slope surface, from ground-based and space-borne radars. The spatiotemporal dynamics of the kinematic clusters reliably predicts where and when catastrophic failure occurs. The clusters share a common boundary along the path of least failure resistance. Here we found this path to precisely locate the impending, shear band in the laboratory sample and the landslide boundary in the natural and man-made slopes. The regime change point is marked by intracluster (intercluster) motions becoming very similar or rigid-body (separated) which, in turn, induces a spatial pattern of physical partitions that become invariant in time through to failure. Our findings illuminate a way forward to rationalize and refine decision-making from broad-area coverage monitoring data for improved geotechnical risk assessment and hazard mitigation. To that end, ongoing efforts are focused on the extension of SSSAFE to a probabilistic platform 28 that incorporates uncertainty systematically for various slopes and relevant scenario projections.

Data availability
The data that support the findings of this study are available from GroundProbe Pty Ltd but restrictions apply since these data were used under license for the current study.Appendix

A. The biaxial compression test data
The laboratory scale data set Biax is a well studied data set from a simulation of polydisperse spherical grains confined to planar biaxial compression test (Fig. 3) 39 . It is governed by a classical DEM model 34 , modified to incorporate a moment transfer to account for rolling resistance and thus capture the effects of non-idealized particle shapes. A combination of Hooke's law, Coulomb's friction, and hysteresis damping is used to model the interactions between contacting particles. An initially isotropic assembly of 5098 spherical particles is prepared and subject to biaxial compression with motion constrained to the plane. The ensemble is subjected to constant confining pressure, with a coefficient of rolling friction of µ r = 0.02 . The initial packing fraction is 0.858. A rolling resistance and a sliding resistance act at the contacts, both of which are governed by a spring up to a limiting Coulomb value of µ|f n | and µ r R min |f n | , respectively, where f n is the normal contact force and R min is the radius of the smaller of the two contacting particles. A summary of all the interaction parameters governing the contacts and other quantities relevant to this study is presented elsewhere 39 .