Monitoring the risk of a tailings dam collapse through spectral analysis of satellite InSAR time-series data

Slope failures possess destructive power that can cause significant damage to both life and infrastructure. Monitoring slopes prone to instabilities is therefore critical in mitigating the risk posed by their failure. The purpose of slope monitoring is to detect precursory signs of stability issues, such as changes in the rate of displacement with which a slope is deforming. This information can then be used to predict the timing or probability of an imminent failure in order to provide an early warning. In this study, a more objective, statistical-learning algorithm is proposed to detect and characterise the risk of a slope failure, based on spectral analysis of serially correlated displacement time series data. The algorithm is applied to satellite-based interferometric synthetic radar (InSAR) displacement time series data to retrospectively analyse the risk of the 2019 Brumadinho tailings dam collapse in Brazil. Two potential risk milestones are identified and signs of a definitive but emergent risk (27 February 2018 to 26 August 2018) and imminent risk of collapse of the tailings dam (27 June 2018 to 24 December 2018) are detected by the algorithm. Importantly, this precursory indication of risk of failure is detected as early as at least five months prior to the dam collapse on 25 January 2019. The results of this study demonstrate that the combination of spectral methods and second order statistical properties of InSAR displacement time series data can reveal signs of a transition into an unstable deformation regime, and that this algorithm can provide sufficient early warning that could help mitigate catastrophic slope failures.


Introduction
Slope failures in the form of landslides or the collapse of engineered structures (e.g., tailings dams) pose a considerable risk to both life and infrastructure.Mitigating the risk they pose through providing an early warning relies critically upon monitoring slopes for precursory signs of instability (Intrieri et al. 2013).Conventionally, this has involved using survey monuments, inclinometers, piezometers, extensometers and ground-based radar to monitor how a rock mass is deforming, with changes in the displacement and velocity providing the most reliable indication of the stability of a slope (Carlà, Farina, Intrieri, Botsialas & Casagli 2017).Accordingly, several empirical approaches have been developed to predict the timing of a slope failure based on displacement monitoring data (Intrieri & Gigli 2016).
The majority of failure timing prediction approaches utilise the concept of accelerating (tertiary) creep theory (Saito 1969, Fukuzuno 1985), under which the behaviour of a material in the terminal stages of failure under constant stress and temperature conditions is governed by the empirical power-law (Voight 1988(Voight , 1989)): where v is the creep velocity, t is the time, and α and A are constants.It has been observed that for a wide range of slope failures α ≈ 2 (Voight 1989).Consequently, for α = 2, equation (1) has the solution: where t f is the time of failure.Equation (2) implies that approaching the time of failure, the inverse of the velocity scales as a linear function of time.Consequently, the intercept point on the time axis of an inverse velocity (1/v(t)) vs. time (t) plot corresponds to the time of failure (t f ).This approach, often referred to as the inverse velocity method, is commonly applied to surface displacement monitoring data to predict the time of slope failures, by using linear regression to extrapolate the inverse velocity trend to the point of intersection on the time axis (Carlà et al. 2018).
Although predicting the timing of slope failures can be crucial in mitigating the risk they pose, there are several factors that affect the ability to do so effectively.Firstly, the conventional groundbased monitoring techniques often provide measurements that are low density, have limited coverage and have a low or irregular temporal sampling frequency, which can make it difficult to detect the precursory tertiary creep (Carlà et al. 2019).Secondly, empirical prediction approaches like the inverse velocity method require expert user-intervention in order to derive reliable estimates of the slope failure timing.This includes filtering the velocity time-series measurements to remove instrumental and random "noise" to help determine the onset of tertiary creep (i.e., Onset Of Acceleration) and increase the fit of the linear regression line on the inverse velocity-time plot (Dick et al. 2015, Carlà, Intrieri, Di Traglia, Nolesini, Gigli & Casagli 2017).However, in the absence of a general rule regarding filtering, the selection of the optimal method is subjective and typically determined through trial-and-error (Rose & Hungr 2007).Furthermore, the inverse velocity trend often only approaches linearity in the final few weeks prior to a slope failure (Rose & Hungr 2007), therefore only enabling reliable short-term predictions that leave little time for risk mitigation measures to be implemented.
The use of satellite Interferometric Synthetic Aperture Radar (InSAR) has been increasingly recognised as an effective solution in overcoming the issue of limited surface displacement measurement coverage.The InSAR method is based upon the concept of interferometry, with precise changes in the surface elevation being derived from the interference of electromagnetic waves from two Synthetic Aperture Radar (SAR) acquisitions (Rosen et al. 2000).Accordingly, satellite InSAR has been used to measured ground motion for a wide range of applications, including volcanology (Massonnet & Feigl 1995, Pinel et al. 2014), landslides (Fruneau et al. 1996, Colesanti & Wasowski 2006), seismology (Massonnet et al. 1993, Peltzer & Rosen 1995) and various other types of geohazard monitoring (Gee et al. 2019).More recently, satellite InSAR has also been applied to slope failure prediction in relation to landslides, tailings dams and open-pit mines (Intrieri et al. 2018, Carlà et al. 2019).A particularly pertinent example concerns the 2019 catastrophic collapse of tailings Dam I in Brumadinho, Brazil.The dam collapsed causing a catastrophic mudflow that resulted in the death of more than 250 people, without any apparent risk indicators from the array of conventional methods being used to monitor its stability (Robertson et al. 2019).However, several studies have retrospectively analysed the stability of the dam preceding the collapse using satellite InSAR (Du et al. 2020, Holden et al. 2020), with some of these deriving inverse velocity-based failure timing predictions from the precursory deformation (Gama et al. 2020, Grebby et al. 2021).While successful predictions close to the actual time of failure were possible using the inverse velocity method, there remains some degree of uncertainty over the reliability of these predictions, due largely to the fitting of the linear regression line and the amount of prior warning of a potential failure that can be provided.
There have been attempts to generalise the inverse velocity method in equation ( 1) by extending the temporal dynamics of the governing mechanisms of failure to consider locally (spatially) varying dynamics.Recent examples include studies by Tordesillas et al. (2018) and Niu & Zhou (2021) and references therein.Furthermore, there is a growing body of literature exploring the potential to offer probabilistic predictions of slope failure by modifying the inverse velocity method in various ways, such as through the application of linear regression (Zhang et al. 2020), stochastic differential equations (Bevilacqua et al. 2019), quantile-regression (Das & Tordesillas 2019, Wang et al. 2020), and mixture-distribution (Zhou et al. 2020).However, with the exception of Bevilacqua et al. (2019), none of these studies explicitly consider the temporal serial correlation of the displacement or velocity fields.Additionally, almost all of these studies assume a parametric statistical distribution on the observed response (i.e., displacement or velocity).Nevertheless, rapid developments in technology for monitoring ground displacements, such as ground-based (Casagli et al. 2010) and satellite-based InSAR (Wasowski & Bovenga 2014) present the opportunity to develop much faster and sensitive precursory warnings using algorithms that are based on streaming data and require less stringent statistical assumptions about the displacement mechanism.Das & Tordesillas (2019) proposed a data-driven, non-parametric multivariate statistical methodology for near real-time monitoring of slope failures.However, Das & Tordesillas (2019) acknowledged that serial correlation inherent in such deformation signals was not modelled.Nonetheless, their framework led to a much earlier prediction of the failure than that offered by the inverse velocity method.
Accordingly, the aim of this study is to further develop the framework of Das & Tordesillas (2019) to propose a new generalised algorithm that can accommodate temporal correlation whilst retaining its non-parametric nature.One key advantage of this statistical approach is that often regular sampling may reveal a change in the slope stability regime much earlier than that observed using the mechanistic failure law.The failure dynamics obtained by the algorithm is driven entirely by properties of the data and a fundamental statistical property underpinning conventional theory of time-series analysis called second order stationarity.Consequently, the statistical nature of this approach also makes it more objective in comparison to the inverse velocity method because the subjectivity associated with manually identifying the Onset of Acceleration and optimally filtering the time-series data is circumvented.
Heuristically, a second order stationary time-series has time-invariant statistical properties of second order.It is therefore postulated that the displacement time-series of a slope, sampled during a stable epoch, must be second order stationary.However, gradual deviations from second order stationary should occur as the slope transitions into an unstable epoch.Here, it is this principle that is utilised in order develop the new algorithm, which can be broadly defined as two phases: 1. Phase 1: Detecting a change in regime -this involves first characterising (estimating) the second order statistical features of an area during a stable regime based on regularly sampled InSAR displacement time-series data.This is achieved through estimating statistical features of multiple time-series in terms of their dynamic periodograms.The stable regime (stationarity) is then defined as the entire history of the observations up until the detection of statistical variation in periodogram features; this point marks the time of regime change.
2. Phase 2: Risk classification -this involves continuously tracking the deviation of the timeseries from second order stationarity after the time of regime change detect in Phase 1.Using a combination of supervised and unsupervised learning algorithms, this deviation is described as a concave probabilistic function of time.Risk of failure warnings are then generated based on estimates of inflection points on the trajectory of this function.
In this study, the proposed algorithm is applied to satellite InSAR time-series data to retrospectively analyse the risk of the Dam I tailings dam collapse.Section 2 describes the study area and the InSAR time-series data used in the analysis.The proposed algorithm is then presented in Section 3, before the results and discussion for its application are provide in Section 4. Finally, the key conclusions, limitations and future research associated with this study are outlined in Section 5. 2 Study Area and Data Dam I (Figure 1a), at the Córrego do Feijão iron ore mine complex in Brumadinho, Minas Gerais State (southeastern Brazil), was an inactive (since 2015) upstream tailings dam that underwent a sudden and complete failure across its full 720 m width and 86 m height on 25 January 2019 as a result of a flow (static) liquefaction mechanism.The expert panel technical report by Robertson et al. (2019) attributed this to a combination of internal creep and the loss of suction induced by a period of heavy rainfall from about October 2018 to the time of the failure.The catastrophic failure of Dam I caused a mudlfow of approximately 11 million m 3 of mining waste to flow rapidly downstream, destroying properties, infrastructure and agricultural land and entering the Paraopeba River (Gama et al. 2020).Within days of the collapse, the mud covered an area of 3.13 x 10 6 m 3 (Rotta et al. 2020), ultimately claiming the more than 250 lives and affecting the whole region's ecosystem (Porsani et al. 2019).
The data used to analyse the precursory deformation and risk in this study comprises Sentinel-1 InSAR displacement time-series data for two overlapping descending orbit tracks (tracks 53 and 155), for the 17 months preceding the dam collapse.The relative line-of-sight displacement (RLOSD) time-series data covering the extent of Dam I at 20-m spatial resolution were generated using the Intermittent Small Baseline Subset (ISBAS) differential InSAR technique (Sowter et al. 2013) and made publicly available by Grebby et al. (2021).Representative displacement time-series for areas exhibiting distinctive deformation within each of the main structural elements of the dam (i.e., dam wall, tailings beach) were obtained by extracting the time-series of a subset of contiguous pixels at each locality (Grebby et al. 2021).Figure 1b) shows the locations of the six areas (labelled locations 1-6) identified as exhibiting distinctive deformation over the dam.The number of contiguous pixels for which the time-series were extracted at locations 1 through to 6 is 30, 44, 46, 25, 38 and 37, respectively.
The time-series for locations 1, 2 and 3 cover the time period (sampling window) from 19 August 2017 until 17 January 2019 (track 155), while those for locations 4, 5 and 6 covered 24 August 2017 until 22 January 2019 (track 53).In both cases, displacement measurements were derived at a 12-day interval (sampling frequency).Accordingly, the multiple individual displacement time-series and the average time-series at each of the six locations were used to construct the machine-based algorithm for characterising and precursory monitoring of the tailings dam failure.

Methods
The method can be summarised in the following three stages: (i) Monitoring of the RLOSD to identify the time of regime change; (ii) Estimating state-of-the-system at the time of regime change; (iii) Sequential risk estimation and hazard classification.All the above steps build on monitoring the second order statistical properties of a time-series signal.These steps and the related concepts are described in the following subsections.

Statistical preliminaries
Conventional statistical modelling of a time-series, {Y t }, often involves modelling its second order statistical properties, its trend -slowly varying mean effect (µ t ) -the autocovariance (Cov(Y t , Y t+u )), and variance (σ 2 t ), which are defined as: Key to such modelling is the notion of second order stationarity.

Second order stationarity
Let {Y t } be a time-series that is discretely sampled at n regular time intervals {t = 1, 2, . . ., n}, such that its population mean, variance (Var(.))and auto-covariance (Cov(.))are defined as in equation ( 3).Then, Y t is second order stationary if: which are a function of time interval, (t, t + u). (4) A second order stationary time-series is a characteristic property that can be observed in the Fourier domain.

Spectral density and periodogram
If Y t is a second-order stationary time-series with the auto-covariance function, c(u), satisfying that is defined as: (5) The spectrum f (ω) is a unique signature of a second order stationary time-series marked by invariance over time.It partitions the total variance (power) of the time-series among sinusoidal event frequencies.Physicist Arthur Schuster (Schuster 1906) proposed the periodogram to estimate the spectrum of a second order stationary time-series.The discrete Fourier transform (DFT) and periodogram of observed time-series data y t , t = 1, 2, . . ., n are defined, respectively, by d(ω k ) and I(ω k ) as: where (6)

Sequential quantification of non-stationarity
A non-stationary time-series is characterized by a spectrum that varies over time.The empirical methods for monitoring non-stationarity are described below.

Moving window mean and variance
Exploratory investigations into deviations from second order stationarity of a time-series Y t can be performed by estimating local sample means and variances of the time-series within subsets of time windows.If partitioning the length of the time-series n into L overlapping time windows w t l , l = 1, 2, . . ., L of equal length, n L , the local sample mean and variances of Y t in window w l are defined as: 3.2.2Dynamic periodogram: evolving features of a time-series ) is defined as the local periodogram (see equation ( 6)) of the time-series Y t partitioned in each of the w t l , l : l = 1, 2..., L overlapping time windows of length n L , for all sampling locations, {s : s = 1, 2, ..., v}, across the Fourier frequencies, ω k = 2πk/n L , k = 1, 2, ..., n L .Accordingly, I s denotes the complete periodogram matrix at any arbitrary location, s, defined as: where, for simplicity, the local time window I s (ω 1 , w t l ) is denoted by I s (ω k , l).Reordering the periodogram matrix in equation (8) for the periodograms at time window w t l across all sampling locations {s : s = 1, 2, ..., v}: For a second order stationary time-series, I wt l defines the signature of multiple time-series, composed of all second order stationary time-series, I l1 would be time invariant.However, statistically evolving time-series would lead to a time varying periodogram matrix, I l1 .

Feature space: PCA periodograms
Principal components analysis (PCA) was introduced by Karl Pearson and has since become a standard unsupervised statistical learning method (Hastie et al. 2009).A common use of PCA is for reducing the dimensionality of correlated primary features in multivariable dataset.If a dataset consists of n -potentially correlated -features, the application of PCA leads to a smaller number (m < n) of uncorrelated secondary features {Z 1 , Z 2 , . . ., Z m } that comprise linear combinations of the primary (n) features.
It is well known that periodogram ordinates I(ω k ) of second order stationary time-series are asymptotically uncorrelated at Fourier frequencies, ω k = 2πk/n, k = 0, 1, 2, ..., n − 1 (Shumway et al. 2000).However, a non-stationary time-series has correlated periodograms.Here, PCA is used for two purposes: • to reduce the dimension of periodogram vectors (i.e., spectral features in each time window, w l ) retaining most of the variation in the signal; • to obtain uncorrelated secondary spectral features (see columns of equation ( 9), I wt l ) for multiple RLOSD time-series.
Mathematically, this is achieved by solving the following optimization problem using a linear combination of periodogram columns of I wt l : Subsequently, the obtained uncorrelated secondary features are used to characterize the evolution of the state-of-the-system, using multivariate statistical learning methods, cluster analysis and classification.

State-of-the-system at regime change
The first stage in the stability monitoring process is to determine a time window within which the statistical features of the ground displacement signal undergo a structural change from a period of second order stationarity -known as the time of regime change (w t0 ).A variety of statistical methods can be used to detect the transition of the system from a second order stationary to a nonstationary regime (e.g., frequentist, Bayesian).However, owing to the relatively short time-series in this retrospective analysis, an alternative approach is adopted here.Defining: as the local variance of a time-series in window, w t l .Under second order stationarity, V ar(w t l ) remains approximately the same in all windows.However, a substantial departure indicates a change in the statistical regime, which can be identified by plotting V ar(w t l ) -for the average RLOSD timeseries for each of the six locations on the dam -against time windows, w t l , l = 1, 2, . . ., L.
Inflection points on the the V ar(w t l ) plot are chosen as potential candidates for the time of regime change.From a pool of candidate inflection points, the most likely candidate for the time window of regime change (w t0 ) is determined using the principle of maximum inter-cluster variance (Das & Tordesillas 2019).To achieve this, the local periodogram in the window of regime change, I wt 0 , for all locations is estimated, as shown in equation ( 9), before using PCA to obtain a matrix of secondary spectral features, Z wt 0 , for each location s =, 1, 2, . . ., v. The feature matrix Z wt 0 encapsulates the signature (i.e., slope stability) of the state-of-the-system (i.e., the dam) at w t0 , and is derived from the average RLOSD time-series of all pixels at each location.This feature matrix is then used to partition the (six) sampling locations of the dam into a finite number (k) of clusters, C is assumed to be the baseline stable state-of-the-system immediately preceding a dynamic transition into an unstable deformation regime.In the final stage, the evolution of the system is measured using a statistic that sequentially compares all subsequent spectral features Z wt 0 against the baseline, C.

Classification and risk
In this final phase we derive risk thresholds based on a sequential classification methodology against the baseline feature matrix C.

Classification
The following statistical law underpinning classification error was a key observation by Das & Tordesillas (2019).Consider a classification problem where at any given point in time, t, it is desired to classify a group of observations {1, 2, . . ., n} based on features {I 1 , I 2 , . . ., I n } into a finite number of classes, C = {C 1 , C 2 , . . ., C k }.Then: which denotes the probability of classification of a randomly selected observation j = 1, 2, . . ., n into class C l at time t.For simplicity, the superscript, t, is ignored henceforth.

Risk trajectory
The classifier allocates the j th observation to class C z ; z ∈ {1, 2, . . ., k} if the class posterior probability p j , given features I j , satisfies: , and q j = 1 − p j , that is, q j = Prob(I j / ∈ C z ) z = arg max l p Il , z ∈ {1, 2, . . ., k}, the estimated class.
(13) Defining M (p j ) = p j q j , j = 1, 2, . . ., n to be the theoretical classification variation and M (p j ) as the corresponding maximum likelihood estimator, then Das & Tordesillas (2019) showed that expectation and variance of pj share a parabolic relationship: E{M (p j } = p j q j , and its variance, as shown in Figure 2. The algorithm can then estimate two slope stability warning thresholdst R and t I -based on the progression of this trajectory (Figure 2).The first milestone, t R , is the time of an emergent instability risk, which corresponds to the theoretical maximum of the trajectory at p j q j = 0.125, indicating a time of monotonic progression towards failure.This is followed by t I , the time of imminent failure.The risk warning t I is estimated as the time window when the classification variance reaches close to its theoretical maximum that is, p j q j = 0.25.On this basis, retrospective empirical estimates of t R and t I for the Brumadinho tailings dam collapse are generated.Note that estimates of population parameters (e.g.p) based on n samples p 1 , p 2 , . . ., p n are formally denoted by p. Accordingly, functions of parameters such as M (p) are formally denoted by M (p).However, for brevity we drop hat notation from the estimates in the remainder of the article.

Algorithm implementation
Algorithm 1 summarises the methodology for generating warnings of an emergent risk of failure (t R ) and an imminent risk of failure (t I ).A primary difference of this new algorithm compared with Das & Tordesillas (2019) is in the features used in monitoring.While Das & Tordesillas (2019) focused on first order properties of time-series, here the framework has is extended to include serially correlated and non-stationary time-series.The algorithm is implemented as a two-phase process using the statistical software R (R Core Team 2021).For visualization, the package ggplot2 (Wickham 2016) is used.

Estimating the time of regime change
The first phase of the algorithm is the estimation of the time of regime change, w t0 , as described in section 3.2.4.Figures 3a) and b) show the trend and moving window variance of the time-series of the mean RLOSD observed at each location of the six locations, respectively.As previously outlined, the mean RLOSD was computed as the spatial mean of all the individual time-series extracted Algorithm 1: Algorithm for slope stability risk monitoring. Result Repeat classification of locations into one of the baseline clusters for windows, w j , j = t 0+1 , t 0+2 , . . ., t 0+m ; else Declare that slope failure is imminent end for the subsets of contiguous pixels at the six locations (see Figure 1), as defined by Grebby et al. (2021).With the exception of location 1, all locations exhibit a monotonically decreasing (subsiding) displacement trend (Figure 3a).Furthermore, for locations 2-6, it is also apparent that both mean displacement and the moving window variances vary over time, which indicates non-stationary in second order statistical properties.
A key time-varying second order property is that of sample variance.An examination of Figure 3b) shows that the Brumadinho tailings dam underwent a deformation phase transition.Prior to 2 July 2018, the trajectory of spatial variance decreased in magnitude at each location, while after this date, increases in variance can be seen.The trend-adjusted moving window variances shown in Figure 3c) retain this property.This pattern is consistent across all six locations, albeit in varying magnitude.
As discussed in Section 3, the temporal invariance of the spectral density function (equation ( 5)) is the unique signature of a secondary order stationary time-series, whereas the periodogram (equation ( 6)) is the sample estimator of the spectrum.Thus, a common approach of assessing nonstationarity in a time-series is to estimate local periodograms in moving overlapping time windows and to then analyse the temporal variation of local periodograms.Periodograms of the mean adjusted RLOSD in local windows of several different lengths (n L ) were estimated.In each local window w t l , l = 1, 2, . . ., L we add one time unit (12 days) and remove the first time point of the previous window.The resulting local periodograms for each location for a local window length of 16 and overlaps of 15 time units are shown in Figure 4, where the x − axis on the plots represents the Fourier frequencies, ω k .For a time-series of length 16, it is only obtain possibly to obtain distinct estimates of spectra at 8 Fourier frequencies due to aliasing (Nason et al. 2017).Figure 4 clearly shows clear temporal variation in the local periodogram curves for each location.The evolving nature of the time-series is evident in that the periodogram ordinates appear more similar in some time windows -in terms of both intercept and slope -than in other windows.It is the 8 periodogram ordinates that are used as features in the sequential slope stability monitoring algorithm.
The next key step in Phase I of the algorithm is to then utilise the periodograms to identify a time (window) of regime change -i.e., a phase transition of the system from a stationary to non-stationary paradigm.To do this, Das & Tordesillas (2019) suggested implementing the non-stationarity metric on second order time domain properties proposed by Das & Nason (2016) to assess deviation from stationarity.Here, however, a simpler approach that is applicable to spectral methods was adopted.It is well-known that the spectrum is a orthogonal partition of the variance of a time-series over Fourier frequencies (Shumway et al. 2000).Equivalently, the moving window spectra partitions the local variance -within a particular time window -among the Fourier frequencies.When a timeseries is second order stationary, the spectrum is its unique signature and does not vary over time.This is observed as invariance in periodogram plots, over local time windows.In other words, the local variance -obtained as the sum of local periodograms at Fourier frequencies -must be uniform across all time windows.However, non-stationary time-series would have time varying (or evolving) local variance as shown in Figure 5a).
For all locations (1-6), the spatial median local variance across all contiguous pixels across moving time blocks were estimated to identify the points of inflection on the local variance trajectory (Figure 5b).Each point of inflection is considered to be a hypothetical point of regime change, w t0 .Accordingly, as shown in Figure 5b), three prospective time windows of regime change were identified, corresponding to windows 7 (2017-10-30 : 2018-04-28), 9 (2017-11-23 : 2018-05-22), and 16 (2018-02-15 : 2018-08-14).At each prospective w t0 , medoid clustering (Hu et al. 2021) of all pixels into a finite number of clustering with periodograms as features was performed, as described in section 3.2.4.The time window that led to the highest proportion of explained inter-cluster variation (with a threshold of ≈ 80%) for the fewest number of cluster partitions was chosen as w t0 (Das & Tordesillas 2019).The corresponding variations for the different numbers of cluster partitions are given in Table 1.With inter-cluster variation of 81% for 4 cluster partitions, w t0 = 16 (2018-02-15 : 2018-08-14) was selected as the optimal time of regime change.

Risk of failure warnings
Using local periodograms for a window length of 16 as features, all time-series were subsequently classified into the class labels of w t0 at all subsequent points of time w t0+1 , w t0+2 , . ... Mis-classification (pq) is treated as an indication of the evolution of the state-of-the-system.At each time window, Consequently, the algorithm identifies a risk of imminent failure of the dam at least a month prior to the actual event.This is comparable with the 40 days of advanced warning based on the inverse velocity method by Grebby et al. (2021).If monitored using this approach in near real-time, this would have provided adequate warning for detailed site investigations or risk mitigation measures to have been implemented.Moreover, the warning of an emergent risk is even earlier, at least five Table 1: Explained intra-and inter-cluster variations for different numbers of cluster partitions, 2, 3, 4, 5, at various hypothetical times of regime change, wt 0 = 7 (2017-10-30 : 2018-04-28), wt 0 = 9 (2017-11-23 : 2018-05-22), and wt 0 = 16 (2018-02-15 : 2018-08-14).months prior to the collapse, coinciding with an initial phase of accelerating deformation highlighted by Gama et al. (2020).This is considerably earlier than the earliest reliable indication of an emerging failure 51 days prior based on the inverse velocity method (Grebby et al. 2021).Although, with hindsight, this deformation could correspond to a period of secondary creep rather than tertiary, at the time, such a warning would still have been useful in initiating crucial follow-up assessments of the stability of the dam.The discussion thus far has focused on the results associated with the estimation of local periodograms over time windows of length 16.In order to assess the effect of the length of the moving window on the regime change and risk estimation, several other window lengths were used in the computations.These included periodogram estimates for the same time-series over moving windows of length 24 and 32, with varying overlap lengths.
Comparing Figures 7 and 8 with Figure 4, it can be observed that all three window estimates show temporal variation (and clustering) of local periodogram curves.The primary differences are two-fold: • longer windows allow for a broader spectrum of variation.That is, when choosing local periodograms for time-series of length 32 it is possible to estimate (low frequency) periodic components at twice the length of that allowed by time window of length 16 (1 time unit = 12 days).
• Accordingly, the local window (numbers) of clustering vary.
Overall, the analysis reveals that all time window lengths estimated the maximum variance to coincide with the Fourier frequency that corresponds to the period of July-August 2018.This confirms that the window length over which the local periodograms are estimated have negligible effect on the prospective time windows of regime change and, hence, the subsequent timing of the risk of failure warnings.However, this may not necessarily always be true, and so a more objective window selection approach would be beneficial in generalising the algorithm for application to a broader range of cases.
The key to objectively optimising window selection is a methodology that can estimate periodograms at all frequencies that have significant variance components, for given a slope failure process.In other other words, such a method must allow estimation of power, inherent in the RLOSD time-series, at the smallest typical frequency or longest time period typical of that type of slope failure, conditional on the sampling frequency.Examples of window width estimation (as also a window kernel) include those from within the statistical sciences and signal processing (equivalent to bandwidth), such as density estimation (Green & Silverman 1993, Bowman 1984), trend estimation (Fried 2004), and periodogram estimation (Ombao et al. 2001).As demonstrated by many of these examples, a rigorous and automatic approach can be developed using the resampling methodology known as cross-validation.One potential solution, based on mechanistic apriori knowledge of prevalent Fourier spectra of the RLOSD time-series for a given type of slope failure, would be to use a cross-validation-based decision criteria to select the optimal window length by minimizing a loss function which will ensure that the spectra at the lowest relevant frequency can be estimated.

Time of regime change windows
As previously shown in Figure 5b), three inflection points were identified and then subjected to cluster analysis (Table 1) to estimate the optimal time window of regime change ( ŵt0 ).As a result, w t0 = 16 (2018-02-15 : 2018-08-14) was selected.However, to assess the effect of the time window of regime change on the risk of failure warnings, estimates of the risk thresholds, t R and t I , were computed for the two alternative prospective time windows of w t0 = 7 (2017-10-30 : 2018-04-28) and 9 (2017-11-23 : 2018-05-22).Figures 9 and 10 show the time-series of the median and interquartile range (IQR) for the risk statistic (pq) and risk thresholds t R and t I for w t0 = 7 and w t0 = 9, respectively.As for w t0 = 16, it can be seen that both alternative candidates for the time of regime change also predict an imminent risk of failure (t I ) well in advance of the actual collapse of the dam.For w t0 = 7, an imminent risk of failure is detected in the time window 11 March 2018-7 September 2018, while for w t0 = 9 this risk is detected during the period 16 April 2018-13 October 2018.These risk warnings precede that for the optimal window for the time of regime change (w t0 = 16) by 2-3 months and the actual collapse by several months.However, as also observed by Das & Tordesillas (2019), these earlier time windows of regime change (w t0 = 7 and w t0 = 9) produce maximum values for the risk statistic (pq) that are less than that for both w t0 = 16 and the theoretical maximum of 0.25.Consequently, the risk warnings generated for the earlier alternative time windows are less accurate -that is, more prone to false positive -than those for the optimal time of regime change window of 16.
Although in this case the cluster analysis was able to identify the optimal time of regime change window from three prospective candidates, a more robust approach could be more effective in addressing the issue of identifying a unique solution.The problem of detecting a time of regime change is similar in essence to that of detecting a structural change point in a time-series.With regards to change point detection, there are numerous potential approaches from both the frequentist and Bayesian perspectives that require future investigation (Killick et al. 2012, Cho & Fryzlewicz 2012, 2015, Ghassempour et al. 2014).Alternately, detecting ŵt0 can also be considered analogous to detecting anomalies in streaming data, as in Hill et al. (2009) and Dereszynski & Dietterich (2011).

Conclusions
Conventional approaches to predicting slope failures are based on the inverse velocity method.However, predictions made using this method can be unreliable for risk management owing to its highly subjective nature and the fact that the deterministic relationships of this exponential law often materialises close to the occurrence of a failure, therefore offering limited prior warning of a potential catastrophe.In this study, an alternative risk identification approach is proposed, which incorporates theories from spectral analysis of time-series to extend a data-driven sequential monitoring algorithm developed by Das & Tordesillas (2019).This statistical-based algorithm is more objective and can be applied to any displacement time-series, sampled at regular time intervals, accounting for serial correlation.
The efficacy of this new algorithm has been demonstrated through application to satellite InSAR displacement monitoring data associated with the 2019 Dam I Brumadinho tailings dam collapse.The algorithm analyses second order properties to show that the tailings dam transitioned into an unstable epoch around July 2018, several months before the actual collapse.This was evident in the moving window variance plots and further corroborated using the periodogram-based local variance plots; the latter of which accounts for serial correlation within each displacement time-series.Whereas most slope stability monitoring analysis only considers the trend or first order moments of the deformation signal, the results of this study shows that investigating second order statistical properties is critical.
The optimum time of regime change window, w t0 = 16 (2018-02-15 : 2018-08-14), was identified based on cluster analysis and then subsequently used to retrospectively estimate risk of failure warnings for Dam I.The algorithm detected an emerging risk of failure during the period 27 February 2018-26 August 2018 and an imminent (maximum) risk of collapse of the tailings dam during 27 June 2018-24 December 2018.The latter warning comes at least several weeks prior to the actual dam failure on 25 January 2019, while the first sign of an emerging risk is evident at least five months prior to the failure.The amount of advanced warning for the risk of imminent failure provided by the algorithm is comparable to previous analysis based on the inverse velocity method, although considerably longer in terms of reliably detecting signs of an emerging failure.This is ultimately due to the statistical nature of the algorithm making it more objective and robust than empirical approaches, which also means that it can be readily applied to other slope failures.Although the risk warnings generated here are based on retrospective analysis, application of the algorithm on near real-time monitoring data could have provided sufficient early-warning to enable detailed site investigations or more urgent risk mitigation measures to have been implemented in order to advert a humanitarian and environmental catastrophe.Overall, the results of this study further attests the role that satellite InSAR can have in integrated dam monitoring and early warning systems.
Further work is required to make the algorithm more widely applicable to analysing the risk associated with any potential slope failure scenario.This would primarily include the implementation of more objective or automated approaches for selecting the optimal time window length and time window of regime change.Additionally, the algorithm may also be further augmented by considering both spatial and temporal auto-correlation, in order to develop a spatio-temporal risk statistic.This would enable not just a temporal warning of a risk of failure, but also highlight the specific location of the instability.

Figure 1 :
Figure 1: The Brumadinho (Brazil) study area.a) Satellite image of Dam I pre-and post-collapse, b) Relative LOS displacement maps for tracks 155 and 53 prior to the collapse, with locations 1-6 indicated.

Figure 2 :
Figure 2: Relationship between the average classification uncertainty and its variance.The three vertical lines correspond to the three risk thresholds in time; wt 0 : time of regime change, tR: time of emergent risk, and tI : time of imminent risk.

Figure 3 :
Figure 3: Spatial trend and moving window variances at six sampling locations.a) Time-series of spatial sample mean of relative LOS displacement (RLOSD).b) Moving window spatial sample variance (MWV) of mean RLOSD.c) Moving window spatial sample variance (MWV) of the detrended mean RLOSD.

Figure 4 :
Figure 4: Local periodogram (window length of 16) of detrended mean RLOSD in successive time windows for each location.

Figure 5 :
Figure 5: Local variance as sum of periodogram ordinates.a) Pixel-wise local variance of mean RLOSD at successive time windows.b) Median local variance of mean RLOSD at successive time windows.

4. 3
Effect of time windows on risk of failure warnings 4.3.1 Local periodograms

Figure 7 :
Figure 7: Local periodogram (window length of 24) of detrended mean RLOSD in successive time windows.

Figure 8 :
Figure 8: Local periodogram (window length of 32) of detrended mean RLOSD in successive time windows.
: Phase I: Feature estimating state of the system if w t < w t0 then Estimate local periodogram, I(ω k , w t ); Estimate a non-stationarity metric S(w t ) based on I(ω k , w t ); Classification and risk At j = t 0+1 classify all locations into C class labels; p l (t) = class posterior probability, q l