Tailings Dams Failures: Updated Statistical Model for Discharge Volume and Runout

This paper presents a statistical model to estimate the volume of released tailings (VF) and the maximum distance travelled by the tailings (Dmax) in the event of a tailings dam failure, based on physical parameters of the dams. The dataset of historical tailings dam failures is updated from the one used by Rico et al., (Floods from tailings dam failures, Journal of Hazardous Materials, 154 (1) (2008) 79–87) for their regression model. It includes events out of the range of the dams contained in the previous dataset. A new linear regression model for the calculation of Dmax, which considers the potential energy associated with the released volume is proposed. A reduction in the uncertainty in the estimation of Dmax when large tailings dam failures are evaluated, is demonstrated. Since site conditions vary significantly it is important to directly consider the uncertainty associated with such predictions, rather than directly using these types of regression equations. Here, we formally quantify the uncertainty distribution for the conditional estimation of VF and Dmax, given tailings dam attributes, and advocate its use to better represent the tailings dam failure data and to characterize the risk associated with a potential failure.


Introduction
The failure of tailings storage facilities can have disastrous consequences for nearby communities, the environment, and for the mining companies, who may consequently face high financial and reputational costs. Tailings are waste resulting from mining operations and are commonly deposited as slurry behind earthen or masonry dams. We refer to this form of tailing storage facility as TSF. In 2015, the breach of the Fundão TSF at Samarco mine in Minas Gerais (jointly owned by BHP Billiton Brasil and Vale S.A.) resulted in 19 fatalities, and was declared the worst environmental disaster in Brazil's history. The company entered an agreement with the Federal Government of Brazil and other public authorities to remediate and compensate for the impacts over a 15 years period. Jointly, BHP and Vale recognized a US$ 2.4 billion provision for potential obligations under the agreement [1,2]. Twenty-one company executives were charged with qualified murder, and up until September 2017 the mine had not resumed operations. The ecosystems impacts caused by a TSF failure can last for many years depending on the nature of the tailings. Samarco is in the process of restoring 5000 streams, 16,000 hectares of Permanent Conservation Areas along the Doce River basin, and 1200 hectares in the riverbanks [3]. It is estimated that the livelihoods of more than 1 million people were affected because of the failure [4]. Improvements in the design, monitoring, management, and risk analysis of TSFs are needed to prevent future failures and to estimate the consequences of a breach.
The design of tailings dams has changed significantly from the 1930s to the present [5][6][7]. Construction of the early TSFs was done by trial and error [8]. During the 1960's and 1970's geomechanical engineering started to be used to assess the behavior of the tailings and the stability of the impoundments [8]. Currently, various studies are required to approve a TSF design and increasingly the plans for remediation and closure of the impoundments have to be included in the feasibility phase. Breach assessments are now part of the requirements in the permitting process of a new TSF or an expansion in many countries. Different parameters need to be estimated while conducting these assessments [9]. These include the volume of tailings (V F ) that could potentially be released, and the distance to which the material may travel in a downstream channel, called the run-out distance (D max ). Empirical regression equations for this purpose were developed by Rico et al. [10] using historical TSF failure data, and are commonly used to characterize such failures (similar empirical relationships have been developed for dams holding water [11,12], but the lack of tailings and differences in design and construction make them inapplicable to tailings dams). However, at site conditions in the mines can vary substantially and there is considerable residual uncertainty associated with the conditional mean value estimated by these equations. In this paper, we rigorously update these regression equations using an updated data set, and characterize the uncertainty associated with the prediction. Using the uncertainty distribution for the conditional estimation of V F and D max using TSF parameters provides a better way to interpret the TSF failure data and to characterize the risk associated with a potential failure.
The calculation of V F is of particular importance for inundation analyses. Typically, TSFs are not totally emptied in case of failure (as opposed to water dams), and only a portion of the tailings are released [10]. In TSFs containing a large amount of water (supernatant pond), the breach would usually result in an initial flood wave followed by mobilized/liquefied tailings [9]. Therefore, the methods developed to estimate the released volume of water or the inundation extent from a regular dam (such as the water dam break-flood analyses methods in [13,14]), do not apply to tailings dams. Empirical equations based on past failures, dam height, and the impounded volume of tailings, are commonly used to get a first estimate of the volume of tailings that could be released and the run-out distance. In Rico et al. [10] V F is calculated using the total impounded volume (V T ) in Mm 3 as in Equation (1) and the outflow run-out distance travelled by the tailings in km (D max ) is obtained using V F and the dam's height in meters at the time of failure (H) as in Equation (2) Many investigators directly use such regression equations in a deterministic way to specify exposure. However, at site conditions vary significantly, and there is considerable uncertainty that needs to be quantified. This uncertainty increases as we consider TSF volumes that are near or beyond the range of the data included in the regression equation. Equation (1) predicts that approximately a third of the tailings in the impoundment (including water) will be the outflow volume. This approach may result in unrealistic estimates when liquefaction is a known risk as it does not take into account the tailings mass rheology (viscosity and yield stress) [9]. As Rico et al. [10] point out, some parameters contributing to the uncertainty in the predictions include sediment load, fluid behavior (depending on the type of failure), topography, the presence of obstacles stopping the flow, and the proportion of water stored in the tailings dam (linked to meteorological events or not). Therefore, it is important to account for the uncertainty in these estimates to derive a probabilistic measure of risk that also accounts for how well the regression fits in a certain range of values of the predictors.
Additional information about TSF failures since Rico et al. [10] published the above-mentioned equations is available. In this paper we update the original data used by Rico from 22 complete cases (including height, storage volume in m 3 , released volume in m 3 , and distance traveled) to 29 complete cases with data compiled in Chambers and Bowker [15]. We compare the results of the original linear Environments 2018, 5, 28 3 of 10 regressions done by Rico et al. [10] with the results using the updated dataset. A new model for the calculation of D max is proposed introducing the predictor (H f ), which is defined as: This variable was introduced to consider that the potential energy associated with the release volume, may be better related to the fractional volume released as opposed to the total volume of the TSF.
For each of the models we consider, and for the final model we recommend, we consider the uncertainty analysis of prediction. We compare the predicted intervals and observed values of V F and D max of three TSF failures across the models that were evaluated to see how well the prediction intervals fit the observed data. The indicated probability of exceedance of the observation as per each model was also assessed.

Materials and Methods
The height at time of failure (H), TSF capacity (V T ), released volume (V F ), and the distance traveled by the tailings (D max ) are inputs in Rico's equations (Equation (1) and Equation (2)). The data used in this analysis is a combination of the cases used by Rico and others compiled by Chambers and Bowker [15], including failures post 2008. Seven of the 29 incidents used by Rico et al. [10] do not have complete information or the information for volume is included in million tons, which cannot be used in the analysis without density data. These cases are in red letters in Table 1. It is important to note that the original data did not include releases as large as the ones experienced in Samarco and Mt Polley (cases 15 and 19 in Table 1), and this becomes relevant for future estimations of the potential risk of larger TSFs. The data on reported failures have variations in different sources; some of these variations are included as footnotes in Table 1. Table 1. Data in Rico et al. [10] and Chambers and Bowker (CB) [15]. Entries in red are incomplete.

No
Mine  Figure 1 shows the relationships between V F and V T , and D max with H × V F (called dam factor in [10]) using the updated dataset (plots in log scale), and D maX with H f , (Equation (3)). V F and V T show a linear relationship in the log form, while for D max , there is greater dispersion with the dam factor and H f .  Figure 1 shows the relationships between VF and VT, and Dmax with H × VF (called dam factor in [10]) using the updated dataset (plots in log scale), and Dmax with Hf, (Equation (3)). VF and VT show a linear relationship in the log form, while for Dmax, there is greater dispersion with the dam factor and Hf. Note that the Dmax vs Hf plot is much tighter than the Dmax vs the dam factor plot. CB are added points from Chambers and Bowker [15] and R are from Rico [10].
VF was estimated in the same way as in Rico et al. [10] using the predictor VT with a log-log (power) transformation and the updated data. For the estimation of Dmax three models were considered: 1. A model titled Dmax.1 which is similar in functional form to the one used by Rico et al. [10] (Dmax.1 in Table 2).  [15] and R are from Rico [10].
V F was estimated in the same way as in Rico et al. [10] using the predictor V T with a log-log (power) transformation and the updated data. For the estimation of D max three models were considered:

1.
A model titled D max .1 which is similar in functional form to the one used by Rico et al. [10] (D max .1 in Table 2).

2.
A generalized linear model (glm) with the Gaussian family using a log link function (D max .2 in Table 2). 3.
A model D max .3 which uses the new predictor H f .
The observed value of V F was used to fit the D max regressions.
The goodness of fit of each model was analyzed with residual plots, outlier identification, analysis of influential observations using Cook's distance, and computing the root mean square error(RMSE) using a 5-fold cross validation (CV). The prediction intervals and probability of occurrence of V F and D max in three historical failures was compared across models using the original and the updated datasets.

Released Volume of Tailings
The presence of the new data points in the updated dataset has an effect in the uncertainty of V F .1 as seen in the R 2 , standard error and 5-fold cross validation results in Table 3. In Samarco, 58% of the tailings contained at the Fundão dam were released (P.15 in Table 1), whereas in the incident of the Gypsum tailings dam in Texas only 1.2% of the contained tailings were released (P.17 in Table 1). The Gypsum tailings dam incident (P.17) was not included in the updated regression (Equation (4)) because it was a minor release, different from the characteristics of the rest of the dataset (identified as an outlier with high influence, Table 3), and had a strong effect in the normality of the residuals. The final regression equation with the updated dataset excluding P.17 is, R 2 = 0.887; standard error: 0.315. Using Rico's notation, this may be expressed as V F = 0.332 × V T 0.95 . However, since the model is fit in terms of log(V F ) it is important to consider the uncertainty of prediction in terms of log(V F ), and then transform the prediction intervals to real space to determine the proper uncertainty intervals for V F . Tests for the residuals from the fit provided by Equation (4) indicate that a Gaussian distribution Environments 2018, 5, 28 6 of 10 cannot be rejected for the residuals at the 5% level (Shapiro-Wilk test p-value = 0.1161). The prediction intervals are then computed at the desired significance level for log(V F ) and then transformed to real space for V F . Table 4 shows the prediction of V F and 90th prediction interval for Samarco, Mt. Polley and the Gypsum TSF using Equation (4). These cases were selected based on the influence they have in the regression of V F with the updated data ( Table 3). The prediction intervals are wider and the predicted mean is larger for Samarco and Mt. Polley when the original dataset is used. In the case of the Gypsum tailings dam, the observed value is not within the prediction interval in neither regression because the volume released was so small and the probability of exceeding it was very high (the same was observed even when that data point was included in the regression). As mentioned in the introduction, the volume of released tailings will vary greatly depending on the type of failure and the composition of the tailings but for a first estimate of potential damage, the scenarios within the prediction interval are useful to assess a range of potential consequences. The upper prediction interval for V F should be constrained by the total volume of tailings available since in all the cases presented in Table 4, the 95th percentile is more than what is physically possible (more than the total impounded are released). In this case, finding the probability associated with totally emptying the dam would be a better approach for risk estimation. We developed an online pp that has the capability of computing the probability of exceeding a value of V F specified by the user (available at https://columbiawater.shinyapps.io/ShinyappRicoRedo/). In this manner, the uncertainty around V F can be considered when estimating D max . The app also provides Q5, Q50 and Q75 of V F . From Table 1, P.10, P.13, P.28, and P.29 were nearly or totally emptied.

Run-Out Distance (D max )
The analysis of D max .1 (the original model by Rico et al. [10]) using the updated and original datasets, shows that the uncertainty increases when the new data points are introduced; R 2 is reduced and the cross validated error increases ( Table 5). The Samarco failure is an influential observation in all the D max regressions except for D max .3 ( Table 5). The distance traveled by the tailings reached 637 km, although more than 90% of the tailings stayed within 120 Km of the dam, the rest were transported in the Doce river all the way to the Atlantic Ocean [15]. The Bonsal TSF (P.7) is also identified as an influential observation (Table 5), this incident had a large D max (close to 1 km) compared to the released volume of tailings (the released volume of tailings was approximately 0.01 m 3 ). The failure mode reported at P.7 was overtopping so it is likely that the tailings dam had a large proportion of water at the time of failure. P.12 also appears as an influential observation; in that case the distance traveled by  Based on the results from Table 5 (R 2 , and 5-fold CV) and the analysis of the residual plots, the best model found was D max .3, which uses the new predictor H f and has the form: Residual standard error = 0.658. Using Rico's notation, this may be expressed as D max = 3.04 × H f 0.545 Table 6 includes the results of the prediction and prediction intervals using the original model (D max .1) with the original and the updated datasets, and includes the results of D max .3 for three of the influential observations. The uncertainty distributions are obtained as before by considering the residuals associated with Equation (5), testing for Gaussian structure (Shapiro-Wilk test p-value = 0.388), and then computing the prediction intervals at the desired significance levels, and transforming them to real space. Figure 2 has examples of the prediction intervals obtained from D max .1 O and D max .3 U compared to the observations. From the results in Table 6 and Figure 2 is evident that the additional data used to fit the regressions of D max dramatically reduces the uncertainty bounds in the prediction of large events such as Samarco and Mt. Polley, although in smaller D max events such as Bonsal, Los Frailes or Omai, the uncertainty is similar or it increases. The probability of exceeding the observed run-out distance shown in Table 6 was calculated transforming the observation to the log space, and evaluating its location in the distribution associated to the prediction interval (t distribution). This is exemplified in Figure 3. Table 6. Predicted values (in km) using all the data points for training the models (using the observed V F ). The importance of considering the uncertainty distribution around the regression of Dmax, rather than using the conditional mean directly is illustrated by the examples in Table 6. For the Samarco incident, the mean value of the predicted Dmax is 174 km using Dmax.3, while the predicted 5th (95th) percentile is 10 (2933) km. The Dmax reported from the actual failure was 637 km downstream, which based on the uncertainty distribution associated with the regression equation, has a probability of exceedance of approximately 22% using model Dmax. 3. In this case the tailings were deposited directly in the Doce River [4], transporting the tailings all the way to the Atlantic Ocean, whereas for other TSF failures an immediate river receptor may not be there, limiting the travel distance. Consequently, if just the conditional mean of the regression equations such as those developed by Rico et al. [10] is used, then one would be rather poorly informed as to the range of potential consequences of a failure. For a probabilistic risk evaluation then, for Samarco, the concern would have been the greater than 600 km impact with a 22% chance rather than the modest 174 km indicated by the regression. It is important to highlight that the observed VF was used to fit all the Dmax regressions but in reality this value might not be known prior to a failure. Therefore, the uncertainty of the estimation of VF from Equation (4) has to also be considered in the predictions of Dmax, which will increase the uncertainty. This is also an issue with Rico's approach reported in Equations (1) and (2), and it is acknowledged in their paper. In the online app we developed, this is addressed and Dmax can be calculated taking into account the uncertainty around VF in a two-step model.

Conclusions
The empirical equations developed by Rico et al. [10] to estimate the volume of tailings released in a tailings dam failure and the run-out distance of the tailings were reviewed. An updated dataset provided information on dam failures that happened after the Rico et al. [10] paper was published and includes cases of dams with larger storage capacity and height than the points in the original dataset. The introduction of the new data points in the regression reduces the uncertainty of the prediction of large failure incidents such as the one occurred in Samarco.
An improved model to estimate the run-out distance is proposed. The model uses the predictor Hf that considers the potential energy associated with the released volume as opposed to the whole tailings impoundment volume. The model proposed has a better linear fit than the original model when using the updated dataset. The updated model to calculate VF is presented in Equation (4), and the new model to estimate Dmax in Equation (5). We recommend using the app we provide which contains the equations (available at https://columbiawater.shinyapps.io/ShinyappRicoRedo/). Since we recommend using the uncertainty distribution for each "prediction" it is easiest for the user to use our web app. As data on other failures becomes available, it can be brought into the app and the model can then be automatically updated. The importance of considering the uncertainty distribution around the regression of D max , rather than using the conditional mean directly is illustrated by the examples in Table 6. For the Samarco incident, the mean value of the predicted D max is 174 km using D max .3, while the predicted 5th (95th) percentile is 10 (2933) km. The D max reported from the actual failure was 637 km downstream, which based on the uncertainty distribution associated with the regression equation, has a probability of exceedance of approximately 22% using model D max .3. In this case the tailings were deposited directly in the Doce River [4], transporting the tailings all the way to the Atlantic Ocean, whereas for other TSF failures an immediate river receptor may not be there, limiting the travel distance. Consequently, if just the conditional mean of the regression equations such as those developed by Rico et al. [10] is used, then one would be rather poorly informed as to the range of potential consequences of a failure. For a probabilistic risk evaluation then, for Samarco, the concern would have been the greater than 600 km impact with a 22% chance rather than the modest 174 km indicated by the regression. It is important to highlight that the observed V F was used to fit all the D maX regressions but in reality this value might not be known prior to a failure. Therefore, the uncertainty of the estimation of V F from Equation (4) has to also be considered in the predictions of D max , which will increase the uncertainty. This is also an issue with Rico's approach reported in Equations (1) and (2), and it is acknowledged in their paper. In the online app we developed, this is addressed and D max can be calculated taking into account the uncertainty around V F in a two-step model.

Conclusions
The empirical equations developed by Rico et al. [10] to estimate the volume of tailings released in a tailings dam failure and the run-out distance of the tailings were reviewed. An updated dataset provided information on dam failures that happened after the Rico et al. [10] paper was published and includes cases of dams with larger storage capacity and height than the points in the original dataset. The introduction of the new data points in the regression reduces the uncertainty of the prediction of large failure incidents such as the one occurred in Samarco.
An improved model to estimate the run-out distance is proposed. The model uses the predictor H f that considers the potential energy associated with the released volume as opposed to the whole tailings impoundment volume. The model proposed has a better linear fit than the original model when using the updated dataset. The updated model to calculate V F is presented in Equation (4), and the new model to estimate D max in Equation (5). We recommend using the app we provide which contains the equations (available at https://columbiawater.shinyapps.io/ShinyappRicoRedo/). Since we recommend using the uncertainty distribution for each "prediction" it is easiest for the user to use our web app. As data on other failures becomes available, it can be brought into the app and the model can then be automatically updated.
This paper emphasizes that these are empirical regression equations with significant uncertainty about the mean. Some investigators directly use such regression equations in a deterministic way to specify exposure. However, at site conditions vary significantly (rheology, water content, failure type, etc.), and even with the log-log regressions presented here, there is considerable uncertainty that needs to be quantified. It is important to account for the uncertainty in these estimates to derive a probabilistic measure of risk that also accounts for how well the regression fits in a certain range of values of the predictors.