Spatiotemporal Evolution of a Landslide: A Transition to Explosive Percolation

Abstract

Patterns in motion characterize failure precursors in granular materials. Currently, a broadly accepted method to forecast granular failure from data on motion is still lacking; yet such data are being generated by remote sensing and imaging technologies at unprecedented rates and unsurpassed resolution. Methods that deliver timely and accurate forecasts on failure from such data are urgently needed. Inspired by recent developments in percolation theory, we map motion data to time-evolving graphs and study their evolution through the lens of explosive percolation. We uncover a critical transition to explosive percolation at the time of imminent failure, with the emerging connected components providing an early prediction of the location of failure. We demonstrate these findings for two types of data: (a) individual grain motions in simulations of laboratory scale tests and (b) ground motions in a real landslide. Results unveil spatiotemporal dynamics that bridge bench-to-field signature precursors of granular failure, which could help in developing tools for early warning, forecasting, and mitigation of catastrophic events like landslides.

1. Introduction

Granular failure is not spontaneous. There exists a precursory failure regime in which patterns in motion emerge and evolve in space and time [1,2,3,4,5,6,7,8,9,10]. In natural field conditions, these patterns are further complicated by various multiphysics phenomena including seismic [11], rainfall [12], human activities [13,14], or a combination of these [6]. Quantitative knowledge of the spatiotemporal dynamics of these patterns is key to early detection of failure from the laboratory to the field. Unfortunately, despite the considerable efforts invested to the study of kinematics in granular systems, both of grain motions in the laboratory and ground motions in the field, a broadly accepted and useful bellwether for failure and its spatiotemporal statistical signature are still missing [3,5,7,15]. Filling this knowledge gap is crucial for developing effective diagnostic and forecasting tools for structural health monitoring (SHM) [4,8,10] and geohazard early warning systems (EWS) [12,13,16,17]. This research need has arguably become one of the most pressing issues in an era where remote sensing and imaging technologies are generating big data on motion and associated deformation continuously in real-time [10,18].

In this study, our goal is to develop a kinematic data-driven approach that can help fill this knowledge gap using the concept of critical transitions in percolation theory. Recent preliminary work [19] explored the process of percolation in time-evolving graphs embodying information on grain motions in a dense granular material during failure. That study uncovered evidence of explosive percolation in the state space of kinematics, but the study was essentially confined to the failure regime. By contrast, our interest lies in the pre-failure regime with a view towards applying percolation theory to find clues to the likely location and time of impending failure.

The concept of explosive percolation transition was first introduced a decade ago by Achlioptas, D’Souza, and Spencer (2009) [20]. In striking contrast to the continuous phase transition from classical percolation theory [21], the so-called “Achlioptas process” leads to an abrupt and discrete jump in the order parameter at the percolation transition—an “explosive” emergence of a system-spanning giant component. This dramatic change in percolation transition arises by introducing an element of competition in the Erdös-Rényi random graph construction process [21]. In the Achlioptas process, a pair of candidate edges, drawn at random, competes for selection. The winning edge is then chosen following a rule that systematically suppresses the growth of the largest cluster. In the short-term, the process promotes the growth of multiple coexisting connected clusters (or components) of comparable size. However, in the long-term the components become increasingly primed for a spontaneous merger that then precipitates an explosive, albeit delayed, critical transition to a system-spanning giant component. Although this new paradigm of percolation has attracted considerable attention in the literature, its applications to real-world systems are still being discovered [22].

Here, our interest lies not in the Achlioptas process per se, but in the outcome of this process from the perspective of a time evolving clustering point pattern in a given feature state space. In particular, recent studies have shown that collective motion develops in the precursory failure regime, where the granular body splits into subgroups, in each of which constituent grains move in near-rigid body motion [2,5,9]. When mapped to a graph where edges are assigned based on kinematic similarity, this collective motion leads to a modular structure comprising densely connected communities, in between which are bridge nodes of relatively sparse connectivity. These bridge nodes, which distinguish themselves by having high closeness centrality, can reliably predict the failure location in the early stages of the pre-failure regime. Encouraged by this, we herein adopt a connection rule that is also governed by kinematic similarity: nodes within a distance r of each other in kinematic state space are connected. Some investigations into explosive percolation have also used connection rules based on the similarity of relevant system features [23,24].

At each observed time state of the system, we characterize the evolving connectivity of the graph $G(r,t)$ as r is increased incrementally. In this context, the sudden emergence of a system-spanning component arises when r reaches a critical value that quantifies the relative motion of the subregions of the system at the given time state. We summarize this with respect to the probability profiles $p\left(r\right)$ versus r, where $p\left(r\right)$ represents the probability that a given node belongs to the giant component: $p\left(r\right)$ is close to zero (one) when $G(r,t)$ is sparsely (densely) connected.

Our strategy rests on the existence of an intermediate or precursory failure regime when the system transitions from the stable to the failure regime, which can be characterized by a kinematic field manifesting the following spatiotemporal dynamics. (a) Intracluster similarity in motion increases as grains belonging to the same cluster begin to move in near rigid-body motion, at the same time as (b) intercluster similarity in motion decreases (i.e., intercluster separation in kinematic state space grows). We test this approach using data from simulations of two classic laboratory tests on granular failure, and a real data from ground-based radar monitoring of a landslide that develops in an operational open-pit mine. Our objective is twofold: (a) to determine whether the components in the resultant network can provide an early prediction of the failure location, and (b) to establish the extent to which the transition from continuous to explosive percolation with so-called powder keg properties in the sense of [25] can be used to identify a regime change point marking the time of imminent failure.

The remainder of this paper is organized as follows. We describe the data set in Section 2 and our method in Section 3. Key results are presented in Section 4 and then discussed in Section 5. Concluding remarks are given in Section 6. A glossary of key terms and symbols used in this paper is provided in the Appendix A.

2. Data

We studied data labeled B1, B2, and Mine. These data sets are described in detail elsewhere: see the work in [26,27,28] for B1 and B2, and the work in [9] for Mine. B1 and B2 are from discrete element simulations of a well-studied assembly of initially homogeneous, densely-packed, polydisperse spherical grains, confined to planar motion. The samples, while identical initially, are submitted to two different biaxial loading conditions, resulting in two distinct failure evolution patterns. B1 (B2) is a displacement-driven compression in the vertical direction under constant confining pressure (constant volume). The failure of these initially homogeneous samples is realized through strain localization, with the ultimate persistent shear band pattern portraying: a single band extending diagonally from the bottom left corner to the top right corner of B1, and two shear bands in the shape of the letter V in B2 (Figure 1). Once fully formed, the shear band partitions the sample into near-rigid body (or underforming) components that can slide relative to each other along their shared boundary—the shear band itself. In the case of B2, the two shear bands interact and take turns in going from passive to active phase. During the passive phase, there is a stress build-up or “jamming” in the shear band; in the active phase, there is stress relief or “unjamming”, which is accompanied by the components on either side of the band undergoing relative sliding [28,29]. Shear band evolution in these samples has been studied using different techniques including those from: complex networks [27], spatial statistics [28], and dynamical systems [29]. In these past studies, particular attention was paid to quantifying energy dissipation and the correlated measure of nonaffine deformation, as well as identifying the mechanisms underlying dissipation including the birth-death evolution of structural load-bearing motifs (e.g., 3-cycles and subclasses of these, force chains, and vortices) [26,27,28,30].

The Mine data is from monitoring of a rock slope in an operational open pit mine [9]. The mine operation, location, and year of the rockslide are confidential. The rock slope stretches to approximately 200 m in length and 40 m in height (Figure 2a). The slope stability radar technology that was deployed to scan the slope and continuously measure surface movement with sub-millimeter precision is described in [13,14]. The Mine data comprise surface displacement along a line-of-sight (LOS) from the stationary ground-based monitoring station to each observed location on the surface of the rock slope at every six minutes. The monitoring period is three weeks—from 10:07 31 May to 23:55 21 June—during which time series data were gathered from 1803 pixel locations across the entire slope at high spatial and temporal resolutions. On 15 June, a rockslide occurred on the western side of the slope with an arcuate back scar and a strike length of around 120 m (red zone in Figure 2b). A global average peak velocity of 0.56 mm/min (33.61 mm/h) was recorded at 13:10 15 June: we refer to this as the time of failure (or time of landslide) $t_F$. There is a “competing slide”, a second region of instability, to the southeast region (Figure 2c). 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. In this context, this region is sometimes referred to as a false alarm in the sense that it did not eventuate into a collapse [18]. The same phenomenon is observed in laboratory tests on real sand samples: competing localized failure patterns sporadically form (e.g., transient microbands and shear bands) just before failure but these later disappear to give way to the “winning” so-called persistent pattern of shear band(s) in the failure regime [1,2,9,28,31,32].

In all the systems above, instability patterns can be observed to develop in the precursory failure regime [9,26,27,28,31]. Past studies have demonstrated that such precursory patterns can be classified into one of two types: persistent versus transient [1,2,9,28,31,32]. The patterns in the former type persists in time with only minor changes to their spatial location and shape. By contrast, the transient patterns are sporadic and disappears just before failure, leaving behind the persistent patterns through to failure. In this study, we are interested in a regime change point $t^{*}$, which marks the time when failure becomes imminent. Our hypothesis is that two conditions must hold from $t^{*}$: (a) the shape and location of failure no longer change across time as the failure pattern becomes permanently incised in the system, and (b) and relative sliding or slip among the split components of the system initiates. In what follows, we develop a method that can quantify these conditions within the framework of explosive percolation theory.

3. Methodology

3.1. Algorithm Input and Network Construction

The input to our analysis is a point pattern in an n-dimensional feature state space that evolves in time. The feature used here represents kinematics: displacements $u_{\ell }\left(t\right)$ of points $\ell =1,\dots ,N$ at time states $t=1,\dots T$. At each time state, a corresponding graph $G(r,t)$ is constructed in the network state space (NSS), where the nodes represent grains in B1 and B2 and pixels on the rock slope (Figure 3). An $N\times N$ matrix of weights is created, where the $ij$-th element describes the kinematic similarity between nodes i and j: $\Vert u_i-u_j\Vert$ for B1 and B2 where $u_i$ is a 2-dimensional displacement vector in the x and y directions; and $|u_i-u_j|$ for the Mine, where $u_i$ is a 1-dimensional displacement vector in the line-of-sight between the radar and a monitored pixel location on the slope surface. Two nodes, i and j, share a link if the $ij$-th element is less than r, where $r>0$ is a free parameter. We then study how $G(r,t)$ changes as r increases from 0 by a constant and small increment $\delta$. When r is small, only a few links exist and $G(r,t)$ is a sparsely connected graph. As r increases, an increasing number of nodes are connected, so that eventually $G(r,t)$ contains a system-spanning component. In previous studies of explosive transitions in real-world systems [23,24], the components in a network can be identified when links are added competitively. In this study we grow our network differently by using the parameter r. Consequently, our method identifies not just the components in a network but also quantifies their separation.

3.2. Algorithm Outputs

The outputs to our analysis at each $t=1,2,\dots ,T$ are (a) the order parameter $p\left(r\right)$ evaluated over a range of increasing r values, (b) the components of the network $C_k$, and (c) a critical radius $r_{c,S}$. We will show that the components $C_k$ give an early prediction of the ultimate pattern and location of failure, whereas $r_{c,S}$ can be used to detect a regime change point $t^{*}$ that is indicative of imminent failure.

Let $S\left(r\right)$ be the size of the largest connected component of $G(r,t)$. The order parameter, $p\left(r\right)=S\left(r\right)/N$, characterizes the evolution of the network by the size of the largest component. If the growth of the largest component is suppressed, multiple coexisting components with sizes similar to that of the largest component will initially form in NSS [33]. Eventually, these components will coalesce explosively across multiple discontinuous jumps in $p\left(r\right)$. Each discontinuous jump reflects the merger of a component to the largest cluster in the system, such that the number of such jumps is one less than the number of components in the system. Here, we are only interested in the major components of a network. Consequently, we limit our analysis to components with size at least $\alpha N$, where $\alpha =0.1$. We analyze the influence of the tuning parameter $\alpha$ on our results in Appendix B. The effect of a tuning parameter on the number of components that emerge in a network has also been analyzed in the context of a generalized Bohman–Frieze–Wormald (BFW) model [34] in [35].

If the network comprises multiple M components, denoted by $C_j$ for $j=1,\dots ,M$, then $d(C_j,C_k)$ measures the kinematic separation of the two components by the minimum distance between two nodes in displacement state space (DSS): one in $C_j$, the other in $C_k$. The separation of a component, $d\left(C_j\right)$, is then defined as the minimum distance of $C_j$ from another component in DSS:

$$d\left(C_j\right)=\underset{k=1,\dots ,M,j\ne k}{min}d(C_j,C_k).$$

(1)

The critical radius of separation is defined as the maximum separation of a component in DSS for time states, where there are at least two components: $r_{c,S}$ = ${max}_{j=1,\dots ,M}d\left(C_j\right)$ for $M\ge 2$. This maximally separated component, denoted by $C^{*}$, has the greatest kinematic separation of all other components in DSS: $d\left(C^{*}\right)=r_{c,S}$. $C^{*}$ is the last component that merges to form a system-spanning giant component and is reflected by the final discontinuous jump in $p\left(r\right)$ that occurs at $r=r_{c,S}$. To illustrate, we depict in Figure 4 the case when $M=3$. The three components, $C_b,C_g,C_r$ (labeled according to their colors blue, green, and red, respectively), form and merge at distinct r. The separation values are $d\left(C_g\right)=d\left(C_r\right)=r_1$ and $d\left(C_b\right)=r_2$. As $r_2>r_1$, $r_{c,S}=r_2$ and $C^{*}=C_b$. If $M=2$ components emerge in NSS instead, then the separation of both components are equal: $d\left(C_1\right)=d\left(C_2\right)=r_{c,S}$ for components $C_1$ and $C_2$. In this case, an additional condition that is physically justified can be imposed to determine which of the two components should be denoted as $C^{*}$.

Components in DSS correspond to groups of grains (pixels) in physical state space (PSS). The collection of pairs of grains (pixels) from distinct groups in PSS that are physically in contact with each other forms a shared boundary. This boundary is a set of continuous lines or arcs if the components in DSS partition the system into distinct clusters in PSS. Therefore, the relative motion of the two groups of grains (pixels) on opposite sides of a shared boundary is quantified by $d(C_j,C_k)$, where $C_j$ and $C_k$ are their respective components in DSS. Naturally, $C^{*}$ induces a boundary in PSS where the relative motion of the two groups of grains (pixels) on opposite sides is maximal. In the next section, we show that the regime change point to imminent failure, $t^{*}$, is governed by two conditions: (a) the spatiotemporal persistence of $C^{*}$, and (b) a sharp increase in the motion of $C^{*}$ relative to the rest of the system. Physically, the former can be interpreted as the condition when the pattern of failure becomes incised in the system, whereas the latter condition is the onset of slip or sliding of $C^{*}$ relative to the rest of the system.

4. Results

4.1. Prediction of Location of Failure for B1 and B2

We observe a precursory failure pattern in DSS in the precursory failure regime, namely, two components emerge in B1 and three in B2 (Figure 5a,c). When mapped back to PSS, these components partition the B1 and B2 samples into two and three subregions, respectively (Figure 5b,d). For both B1 and B2, the components are comparable in size, similar to the powder keg mechanism described in Friedman and Landsberg [25]. Comparing with Figure 1, it is evident that the shared boundaries of the components in the nascent stages of the precursory failure regime provide a reliable predictor of the ultimate location and pattern of failure.

In relation to previous studies, such precursory failure patterns can be classified into two types: transient versus persistent. Transient patterns form intermittently and then disappears at the onset of failure: for example, microbands [32]. By contrast, those in the latter persists in time with only minor changes to its spatial location and shape [1,2,3,5]. Of these, the persistent patterns provide a reliable early indicator of the location of the so-called ”persistent shear band” that governs the failure regime [1,2,3,5]. Note that the patterns in Figure 5 are only shown for a single time state. To establish the spatiotemporal persistence of these, we apply two measures: (a) the accumulation over a period of time states of the boundary location induced by $C^{*}$, characterized by the maximum relative motion of grains on opposite sides of this boundary, and (b) the Jaccard similarity index defined as follows,

$$J_t=\frac{|A\cap B|}{|A\cup B|}$$

(2)

where A and B represent the set of nodes of $C^{*}$ in $G(r,t)$ at two consecutive time states $t-1$ and t, respectively [36]. As $J_t$ is measured relative to the location of the maximally separated component $C^{*}$, a value of $J_t=1$ means $C^{*}$, at two consecutive time states, is exactly co-located in PSS. Note that $J_t$ depends on the existence of an $r_{c,S}$ value, as $d\left(C^{*}\right)=r_{c,S}$. Therefore, $J_t$ is only defined for time states where $M\ge 2$ components emerge in the sample.

In B1, the two components, $C_b$ and $C_r$ (the blue and red components in Figure 5a respectively), are separated from each other by $r_{c,S}$: $d\left(C_r\right)=d\left(C_b\right)=r_{c,S}$. Both components share a boundary along the forward diagonal in the sample. In the B1-B2 samples, our key focus lies in identifying a persistent shear band pattern in the pre-failure regime from the boundary induced by $C^{*}$ in PSS. As a result, $C^{*}$ can be chosen to be either $C_b$ or $C_r$ as they lead to the same shear band prediction in B1. In Figure 6a,b below, the boundary of $C^{*}$ in B1 is accumulated across two periods of time states: (a) $t=78$ to 91 and (b) $t=92$ to 104, determined from the $J_t$ values in B1. Figure 6c shows that $J_t$ increases close to 1 at approximately $t=78$ and persists up to failure, with a small dip to $J_t\approx 0.7$ at $t=91$. This slight decrease in $J_t$ can be observed in the difference in the accumulated shared boundary in Figure 6b relative to Figure 6a: there is a marginal clockwise rotation of the accumulated shear band location about the bottom left corner of the sample. From $t=92$ to $t=104$, the predictions become concentrated along the forward diagonal in B1, indicating that the shear band becomes incised in the sample during this time period.

In B2, $C^{*}$ is a unique component at time states where three components emerge. Similar to B1, the boundary of $C^{*}$ in B2 is accumulated over two periods of time states: (a) $t=110$ to $t=135$ and (b) $t=136$ to $t=153$ (Figure 7a,b). In Figure 7c, the $J_t$ values appear infrequently prior to $t=110$. This is succeeded by rapid fluctuations in the $J_t$ values. Of the three components present: $C_b$, $C_g$ and $C_r$ (the blue, green and red components in Figure 5c respectively), $C^{*}$ alternates between $C_b$ and $C_g$ at different time states of the pre-failure regime. This results in the formation of two distinct shear bands, with the left (right) shear band corresponding to the boundary of $C^{*}$ when $C^{*}=C_g\left(C_b\right)$. Note that both boundaries are not always present in B2. As discussed in previous studies, not all shear bands are active at the same time: they take turns going from an active state (when the shear band unjams) to a passive state (when the shear band jams) and vice versa [28]. When $J_t=0$, the component $C^{*}$ switches from one side of the B2 sample to the opposite side (i.e., from $C_g$ to $C_b$ or vice versa). The fluctuations in $J_t$ thus indicate that the order in which the three components merge frequently changes, even though the failure pattern in B2 remains in the same spatial location. Between $t=110$ and $t=135$, the similarity index fluctuates between 0 and a maximum value of only $J_t\approx 0.8$. As such, the location of each shear band shifts over time, resulting in highly varied predictions seen in Figure 7a. Additionally, microbands [32] can be observed in the bottom right corner of the sample. From $t=136$ onward, the microbands disappear leaving only the two shear bands: $J_t$ is either 0 or much closer to 1. This implies that the two boundaries of $C^{*}$ in Figure 7b now persist in the same location (with the left boundary of $C^{*}$ or left shear band unjamming more often). Equivalently, the failure pattern of two bands in the shape of the letter V has become incised in B2.

4.2. Identification of Regime Change Point for B1 and B2

We quantify the kinematic separation of the components in B1-B2 using the order parameter $p\left(r\right)$ and critical radius $r_{c,S}$. The evolution with time of the $p\left(r\right)$ vs. r profiles is shown in Figure 8a,b. The number of components that emerges at time t is evident in the number of discontinuous jumps in $p\left(r\right)$. At the early time states, when the impending failure location is first predicted, we observe that the transition to a system-spanning component occurs at a small value of r. Consequently, the critical radius $r_{c,S}$ (Figure 8c,d) is also small at these time states. Closer to failure, a long plateau in $p\left(r\right)$ develops during which increases in r merely lead to the addition of no links or intra-component links in $G(r,t)$. An inter-component link finally forms at a large value of $r=r_{c,S}$, where $C^{*}$ merges to form a system-spanning component. This increasingly delayed transition is similar to what is observed in explosive transitions [33]. The long plateau in $p\left(r\right)$ manifests as a sharp burst to a peak in $r_{c,S}$ at $t_F$. In B1, $r_{c,S}$ increases significantly at $t=92$ and remains relatively large for majority of the time states prior to the onset of failure at $t_F=104$. Similarly, in B2, an increase in $r_{c,S}$ values is observed at $t=136$ up to the onset of failure at $t_F=153$. As such, the time states $t^{*}=92$ and $t^{*}=136$ are the respective regime change points in B1 and B2. From $t^{*}$, the failure patterns become etched in the samples enabling relative sliding of $C^{*}$ with respect to its adjacent component.

As $r_{c,S}$ corresponds to the maximum kinematic separation of components in DSS, it quantifies the extent to which the clusters of grains in PSS move in opposing directions as near-rigid-bodies along the predicted shear band location. Therefore, it is not surprising that this measure of relative sliding between components correlate with energy dissipation in B1 and B2, as quantified in [26,28]. In B1, the two sharp bursts in $r_{c,S}$ and $t=94$ or $t=103$ represent two major slips between the two components along the predicted shear band location in B1. In B2, multiple peaks in $r_{c,S}$ between $t^{*}$ and $t_F$ suggest regular slips between the components: these occur along each of the two shear bands at different time states. We also note that previous studies of B1 and B2 [28] have shown that the shear band in B1 is less dense with lower average grain coordination number at failure than those in B2. This results in a higher degree of relative rigid-body-like motion between the two parts of the sample which explains the higher peak $r_{c,S}$ values in B1 compared to B2.

4.3. Prediction of Location of Failure for the Mine

Unlike the controlled and idealized loading conditions of the B1–B2 samples, the operational open pit mine was subjected to various perturbations such as mine blasting and weather changes. However, details of these external stimuli were not available. In this context, this study explores the extent to which meaningful information on impending failure can be uncovered from monitoring displacement data alone using the concept of explosive percolation. At this juncture, it is instructive to consider an additional kinematic feature state space, a 1-dimensional velocity state space (VSS), that can be easily derived from the given displacement data. A motivation for this is that velocity is considered widely in forecasting the time of failure in landslide monitoring [18]. Note that in the B1–B2 samples, the constant time intervals means that the velocities and displacements of the grains are proportional. Therefore, DSS is equivalent to VSS. This is not the case for the Mine since measurements at distinct pixel locations on the rock slope are taken along a line-of-sight: displacements are calculated with respect to a reference position while velocities are calculated with respect to the previous time state.

As our aim is to predict the location of the rockslide (Figure 2), we are interested in $M=2$ components with the expectation that the failure region is that which exhibits anomalously high cumulative motion and rate of motion close to the time of failure. Thus, we assign $C^{*}$ to be that component with the higher mean cumulative displacement (velocity) in DSS (VSS). The predicted failure region from $C^{*}$ in DSS is shown in red at representative time states in Figure 9 (left column). In the early stages of the precursory regime, the prediction of the failure zone from DSS captures much of the west wall that ultimately collapses, along with the false alarm region to the southeast corner. As time advances, the red zone along the west wall persists, whereas that to the southeast corner appears intermittently. These two unstable zones interact which complicates the dynamics of the precursory failure regime [9,31]. As noted in [9], instabilities from the west shift to the middle, closer to the southeastern region when this is active. In fact, for most of these time states, these unstable zones combine to form one of the emerging components (e.g., at $t=1500$ and $t=3200$ in Figure 9 (left column)), consistent with Figure 2c. On the other hand, the location of $C^{*}$ in VSS is uninformative and corresponds to a group of pixels that is scattered in PSS (red zone in Figure 9 (right column)) at all time states in the early pre-failure regime. The group of pixels do not form a continuous separating boundary that is indicative of impending failure. The feature state space VSS can give misleading results in the early stages of the precursory regime, as transient false alarm regions can distort rates of motion by increasing at intermittent time periods only to slow down and stabilize later close to failure. This can be further compounded by using the spatially averaged velocity of pixels to identify $C^{*}$ in the early stages of the pre-failure regime.

To quantify the persistence of $C^{*}$ and distinguish between the persistent and transient regions, we introduce two similarity measures in our analysis: $J_t^{\left(1\right)}$ measures the similarity of the component $C^{*}$ in DSS across consecutive time states, whereas $J_t^{\left(2\right)}$ measures the geographical similarity of $C^{*}$ obtained separately from each kinematic feature state space, DSS and VSS, at the same time state. Note that in the terminal phase of the precursory regime, we can expect the failure location to portray not only higher cumulative displacements across time, but also significantly higher velocities compared to the stable zone. Moreover, as discussed earlier in Section 2, previous studies have shown that the transient patterns fade away close to failure leaving behind the persistent pattern that prevails through to failure [9,26,27,28,31]. Thus, altogether, these conditions imply that $C^{*}$ from DSS should become invariant in time while simultaneously converging to $C^{*}$ from VSS—resulting in a rise in both $J_t^{\left(1\right)}$ and $J_t^{\left(2\right)}$ to their maximum of 1.

In Figure 10, $J_t^{\left(1\right)}$ initially fluctuates between 0 and 1 reflecting the intermittent formation of the instability to the southeast: recall $C^{*}$ in Figure 9 (left column) most notably from $t=2800$ to 3000. Eventually, at $t=3333$, we observe a final increase in $J_t^{\left(1\right)}$ with only the pattern along the western wall remaining in the prediction. This transition at $t=3333$ is unambiguously confirmed in the sharp increase of $J_t^{\left(2\right)}$ at this time. The regime change point at $t=3333$ can be visualised in Figure 11. From $t=3333$, the transient pattern to the southeast disappears from $C^{*}$ in DSS (just as microbands to the bottom right corner in B2 have faded away (Figure 7a,b). The region along the west wall consolidates and slowly moves leftward before converging to the actual failure location at $t_F=3568$. The failure location is also present in $C^{*}$ from VSS, and coincides with that from DSS. In summary, our results here suggest that although we can predict the failure location early in the pre-failure regime using DSS, the more reliable indicator of imminent failure is the colocation in PSS of the critical components $C^{*}$ from DSS and VSS.

4.4. Identification of Regime Change Point for the Mine

The separation of the predicted failure location $C^{*}$ from the rest of the rock slope is summarised by the $p\left(r\right)$ vs. r profiles and the evolution of $r_{c,S}$ over time in both DSS and VSS. In both $p\left(r\right)$ vs. r graphs (Figure 12a,b), the value of r preceding the formation of a system-spanning component increases in the lead up to failure, with a longer plateau in $p\left(r\right)$, similar to our observations in B1 and B2. Compared to the evolution of $p\left(r\right)$ with r in DSS, the rightward shift in the graphs in VSS is more distinct. This is consistent with our earlier result in Figure 10, where the predicted failure location $C^{*}$ from VSS provided a better signal of imminent failure. A sharp increase in the values of the two $r_{c,S}$ values is also observed at $t=3333$ (Figure 12c,d). This time state coincides with the increase in geographical similarity in the predicted failure location using $J_t^{\left(1\right)}$ and $J_t^{\left(2\right)}$, and reaffirms that $t^{*}=3333$ is the regime change point of the Mine. The increase in $r_{c,S}$ at $t^{*}$ can be interpreted in a similar fashion to B1 and B2, where the rigid-body-like motion between the failure and stable regions increases along the arcuate boundary of the landslide, the border between the stable zone and the slip region that collapses.

Altogether, our results suggest that the time when failure becomes imminent is $t^{*}=3333$. This corresponds to a date–time of 14 June, 13:16 p.m., approximately 24 h before the actual time of failure of 15 June, 13:10 p.m.

5. Discussion

We uncover a novel relationship between the evolution of the order of a percolation transition and the emergence of precursory signatures of localized failure in granular systems. This relationship proves robust for both laboratory and field scale conditions. We started with kinematic data across many times states from the stable regime through to failure: individual grain motion from simulations of laboratory tests on the one hand, and ground motion from real landslide monitoring on the other. At each time state, each data set forms a point pattern in the displacement state space. In the precursory failure regime, time-evolving clusters of densely packed points in displacement state space emerge, separated by low-density regions where points are relatively far apart from each other. This clustering pattern reflects the development of system partitions in physical state space, in which constituent grains move in unison, or, in the case of the rock slope, similar surface motions manifest at distinct pixel locations. We characterize the development of this collective motion in the corresponding network state space through the lens of explosive percolation. In particular, we tracked the time evolution of the network components as they approach the critical transition to explosive percolation and measured their corresponding critical radius $r_{c,S}$, the maximum separation among all the components in the network. In physical state space, these components give an early prediction of the ultimate pattern of failure for these systems, as previously studied using other techniques [1,2,9,28,31,32]. However, a time is eventually reached in the later stages of the precursory failure regime when the shared boundaries of these components persist in the same location across time. In systems with multiple shared boundaries (due to more than two components), different boundaries may manifest at different times but their respective locations do not change. During this period, a regime change point $t^{*}$ is reached when $r_{c,S}$, which quantifies the extent of relative sliding, undergoes a sudden and sharp increase to a peak. Thus, $t^{*}$ corresponds to when the partitions become incised in the physical state space, thereby precipitating the onset of slip between components along their shared boundaries. Moreover, the order parameter profile $p\left(r\right)$ versus r at $t^{*}$ distinguishes itself from previous time states by a prolonged delay in percolation transition, as evident in a relatively long plateau preceding the jump.

With respect to findings from prior studies of the fully confined laboratory samples, $t^{*}$ coincides with the onset of collective buckling of columnar force chains, which rapidly culminate in the full development of the shear band at $t_F$ [26,29]. During $t^{*}\le t\le t_F$, the “broken pieces” of the granular body undergo relative sliding along their shared boundaries—the shear band (see also [3]). This explains the burst to a peak in the kinematic separation $r_{c,S}$. As shown in past studies of these systems, this time interval coincides with similarly sharp peaks in various measures of energy dissipation including: fluctuating kinetic energy, dissipation rate, average grain rotations, average local nonaffine deformation, and numbers of collapsing load-bearing structures (force chains, 3-cycles), and attendant vortices [26,28,29,30,32]. To test the robustness of the patterns, we quantified the spatiotemporal persistence of the shared boundary between the predicted components at consecutive time states in the lead up to the time of failure $t_F$.

In the field scale data on cumulative displacement for the rock slope, essentially the same trends were observed when compared to those from laboratory simulations. The interaction of the two unstable regions, the west wall and the southeast region, is reflected in the spatial variability of the predicted components for much of the precursory pre-failure regime. The instability along the west wall laterally shifts back and forth, but is mostly pulled towards the middle of the slope when the competing slide to the southeast is active. The complex dynamics of this interaction was recently characterized in [31]. To confirm this and the patterns uncovered, we carried out additional tests by repeating our analysis on the velocity data set, which is the information used to forecast the time of failure in landslide monitoring [15]. At the time of imminent failure $t^{*}$: the critical components of both state spaces converge to, and remain in, the same location in the physical state space, while their corresponding critical radii $r_{c,S}$ simultaneously undergo an abrupt increase to a peak. Moreover, the time of imminent failure is consistent with that recently uncovered in [31] using a combined statistical and machine learning approach. This highlights the need for future research into the optimal combination of kinematic properties for input data in forecasting the location and time of failure.

6. Conclusions

The concept of explosive percolation is used to develop a data-driven approach for characterization and prediction of granular failure from kinematic data. We test our approach using two data sets: one comprises individual grain motions from simulations of classical laboratory tests in which the granular material is driven to failure, the other from ground motions from radar monitoring of a developing landslide. Results demonstrate that our approach can provide an early prediction of the location of failure and help identify the time of imminent failure for both laboratory and field conditions. Accordingly, this work bridges a gap between fundamental knowledge gained from micromechanical studies of granular failure and the application of this knowledge in forecasting impending landslides. In the case of micromechanics, this study has also opened a path for future research into the influence of key material properties (such as grain shape and size distribution) on the overall strain localization pattern that demarcates the near-rigid body moving parts of the granular mass at failure. In the case of landslide forecasting, work is underway to explore the optimal combination of kinematic properties to analyze for a range of different landslides and to extend this approach to a probabilistic framework, which includes information on known triggers of landslides such as rainfall and seismic events.

Finally, our approach is general and can be employed to summarize the salient elements of a spatiotemporal dynamics of a system that develops a modular network structure or clustering pattern in a relevant feature state space. At variance with standard clustering algorithms (e.g., k-means used in [31]), the method proposed here does not need a pre-defined number of clusters for each point pattern. Also, the analysis can be extended to characterize the spatiotemporal dynamics in high dimensional feature state space and/or the spatial convergence of clustering structures derived from multiple lower dimensional feature state spaces.