Multi-Year Mapping of Disturbance and Reclamation Patterns over Tronox’s Hillendale Mine, South Africa with DBEST and Google Earth Engine

This study was devised to examine the pattern of disturbance and reclamation by Tronox, which instigated a closure process for its Hillendale mine site in South Africa, where they recovered zirconium- and titanium-bearing minerals from 2001 to 2013. Restoring mined-out areas is of great importance in South Africa, with its ominous record of almost 6000 abandoned mines since the 1860s. In 2002, the government enacted the Mineral and Petroleum Resources Development Act (No. 28 of 2002) to enforce extracting companies to restore mined-out areas before pursuing closure permits. Thus, the trajectory of the Hillendale mine remains unstudied despite advances in the satellite remote sensing technology that is widely used in this field. Here, we retrieved a collection of Landsat-derived normalized difference vegetation index (NDVI) within the Google Earth Engine and applied the Detecting Breakpoints and Estimating Segments in Trend (DBEST) algorithm to examine the progress of vegetation transformation over the Hillendale mine between 2001 and 2019. Our results showed key breakpoints in NDVI, a drop from 2001, reaching the lowest point in 2009–2011, with a marked recovery pattern after 2013 when the restoration program started. We also validated our results using a random forests strategy that separated vegetated and non-vegetated areas with an accuracy exceeding 78%. Overall, our findings are expected to encourage users to replicate this affordable application, particularly in emerging countries with similar cases.


Introduction
South Africa is richly endowed with a variety of mineral resources and is among the most prominent raw material exporters in the world. In recent years, the demand for zirconium-bearing minerals, of which South Africa holds most of the world's reserves, has proliferated [1]. This has had important socio-economic implications for South Africa as the second-largest producer of this commodity (19.4%) after Australia (41.7%), together accounting for 44 Mt (61%) of the 72 Mt world total reserve base [2]. Typically, zirconium is found in association with the principal titanium minerals ilmenite and rutile, thus providing more commercial sources for zircon [3]. The extraction of these heavy sand minerals is a mature enterprise in South Africa with three currently active production mines: Richards Bay Minerals (RBM) and Tronox, both in KwaZulu-Natal (KZN) and Namakwa Sands in the Western Cape province [4]. Our focus here is mainly on the northeast coast of KwaZulu-Natal province, where Tronox KZN Sands recover these deposits by completely obliterating the existing vegetation. Within this mineral-rich coastal area, its operation consists of the currently active Fairbreeze mine and the other in the closure stage at Hillendale, which is the subject of this paper [5].
The realization of sustainable mining has become a global priority, mainly because mine closure is increasingly accepted to restore the land to a workable post-mining state, algorithm to document the effective RBM restoration in South Africa from 1984 to 2018. They also applied a random forests (RF) classifier to distinguish vegetation from nonvegetated with an overall accuracy exceeding 90%. In a similar study, Xiao et al. [37] mapped mining-induced disturbance and reclamation in Shengli Coalfield in China using LandTrendr and GEE, and their results showed a 26% reclamation rate. Furthermore, Ang et al. [38] successfully applied an RF supervised classification to detect a land cover change in mining areas in the Philippines within the GEE environment and achieved high classification accuracy. They also derived socio-economic indicators from land use and cover trends. More recently, He et al. [39] employed Landsat imagery and GEE to monitor the subsidence of water in mining landscapes on the eastern plain in China. They successfully identified subsidence and restoration with accuracies from 79% to 88% using an RF classifier. Despite these advances, several mining areas across the world, and in South Africa in particular, remain unstudied and inconclusive.
In this study, we present a case study demonstrating the restoration of the highly transformed Hillendale mine on the northeast coast of South Africa. At present, this site is in the rehabilitation phase intending to secure a mine closure permit [5], and no study has yet been conducted to examine its trajectory. Thus, our goal is to monitor the progress of this task using GEE and the DBEST algorithm based on Landsat-derived NDVI from 2001 to 2019. To the best of our knowledge, no attempt has yet been made to segment a time series of vegetation index over mining areas using DBEST, and our study is the first trial in this regard. Overall, our results are expected to encourage users to replicate this affordable application, especially in developing countries where similar cases are expected.

Site Description
This study was carried out at the Hillendale mine operated by Tronox KZN Sands on the mineral-rich northeast coast of KwaZulu-Natal, near the towns of Empangeni and Richards Bay, South Africa ( Figure 1). The site covers 495 ha that was subjected to extensive transformation during its operation between 2001 and 2013. The procedure entailed hydraulic excavation to recover the zirconium-and titanium-bearing heavy minerals from sand dunes. This followed a gravity separation trail to produce a heavy mineral concentrate, with sand tailings reverted to recreate dunes and the bulk of the fines fraction discharged in a residue storage facility (RSF) [5].
Typically, this process completely destroys the soil and associated mineral nutrient balance for plant growth [40]. Before mining, sugarcane plantation was the dominant activity at the Hillendale site [40]. The company is dedicated to rehabilitating the site towards a self-sustaining ecosystem that is intended for sugarcane production again. This area is characterized as subtropical to tropical with hot and humid summers and warm winters, and the mean annual temperature is 21.8 • C. The area is classified as a high rainfall region by South African standards, with a mean yearly rainfall of approximately 1300 mm. The vegetation used in the rehabilitation varies across the mining area. For instance, the walls of the RSF are planted with kikuyu grass and some indigenous grass species to mitigate erosion problems and enhance the aesthetic of this landscape [41]. The Eucalyptus and Casuarina species are intended for use in the rehabilitation of the RSF through a combination of soil improvement and soil drying [41]. Moreover, while the post-mining land use is aimed at sugar cane, these tree species are currently used for stabilizing the area.

Satellite Data
We obtained freely available collections of Landsat scenes from the Google Earth Engine (GEE) cloud computing platform using the JavaScript code editor (https://earthengine. google.com/ Mountain View, CA, USA; accessed 25 February 2021). These 30-m Landsat 5, 7, and 8 images were filtered and cloud masked using the JavaScript code editor. Moreover, the images were already orthorectified, and Level 1 Precision Terrain (L1TP) corrected. The L1TP-corrected product fulfills both radiometric and geometric criteria set by USGS (U.S. Geological Survey), and is considered to be of high quality and usable [42]. Following Karan et al.'s [18] study, which found the normalized difference vegetation index (NDVI) to be the best indicator for detecting mine-induced vegetation changes, we computed monthly mean composites of this index from the collection of usable images from 2001 to 2019. The NDVI values were averaged for the entire mining area to get an overall picture of the mine disturbance and recovery patterns. The NDVI is computed by differencing the NIR and red band and dividing by their sum, which is expressed as follows: where values range between -1 and 1; the negative values reflect an absence of vegetation, while positive values represent vegetative areas. In our context, rehabilitated sites are associated with larger NDVI pixel values than earlier, while lower values indicate exposed soil and very low vegetation cover. However, the resultant time series was irregular because some months had missing observations due to cloud cover. Only the best available cloudless scenes were intended to be used for analysis, as cloud cover is persistent in the studied region. To resolve this problem, we emulated Halabisky et al. [43] and applied spline interpolation to fit a cubic spline to each missing NDVI value over the time series using the R statistical package "spline". We then denoised the NDVI time series using wavelets to implement our intended algorithm effectively. The denoising procedure is expected to remove the noisy data and retain as much high-quality data as possible [44]. The flowchart is illustrated in Figure 2. Before the implementation of the change algorithm, it is important to note that NDVI values were averaged for the entire Hillendale mine site, and the mean annual NDVI was used for analysis. This is because rehabilitation efforts in many mining companies are not often implemented across the entire mining area at the same time. Therefore, the averaging of the NDVI for the entire mining area is aimed at considering the overall recovery. We also analyzed patterns of disturbance and recovery in selected mining sites to demonstrate how different sections of the mines are recovering.

Change Detection Using DBEST Algorithm
In this study, we applied the Detecting Breakpoints and Estimating Segments in Trend (DBEST; [35]) program for segmenting and analyzing trend changes in the time series of vegetation indices. DBEST has two primary algorithms: the change detection algorithm and the generalization algorithm [35]. The main aim of the change detection algorithm is to detect a trend, determine the type of changes (abrupt or gradual), and estimate change timing, change magnitude, change number, and the direction of change ( Figure 3). On the other hand, the generalization algorithm simplifies the detected trend into the main features of the time series. In its operation, DBEST employs a novel segmentation algorithm that estimates the direction into linear segments using one of the three user-defined parameters, namely a generalization-threshold parameter δ, the m most considerable changes, or the threshold β for the magnitude of changes of interest for detection [35]. The cited work found DBEST to be a quick, accurate, and flexible tool that determines trend alteration estimating the timing, extent, and direction. This study used the DBEST trend generalization and change detection algorithm to detect breakpoints with the most remarkable changes in NDVI as a proxy for mining-induced disturbance and recovery. We executed the DBEST application in the R environment using the "DBEST'" package [45] (https://cran.r-project. org/web/packages/DBEST/index.html; accessed 1 March 2021).

GEE Framework for Change Monitoring Over Mining Areas
The advent of the Google Earth Engine (GEE) cloud-based platform brought unprecedented transitions in remote sensing and greatly enhanced the democratization of satellite data. The GEE leverages computational services for planetary-scale analysis and consists of petabytes of geospatial data, including a complete archive of Landsat images, climate data, and state-of-the-art algorithms with highly dense time series analysis [46]. Furthermore, the GEE eliminates many resource-based barriers to satellite data gathering, storage, and processing and hosts ready-to-use products such as NDVI [47], permitting users to devote more time to making sense of the data than on preparation. GEE further enables users to share scripts and assets [48]. These attributes have leveled the playing field for all users irrespective of their location, increasingly making GEE the preferred choice for geospatial data analysis [49]. Since its inception, several studies have widely employed GEE in various land use and land cover change applications [49], but relatively few have exploited this asset for monitoring dynamics over mining landscapes, namely only Dlamini and Xulu [36] in South Africa until now.

Validation and Accuracy Assessment
The study area was categorized into vegetation and non-vegetated classes to analyze the vegetation changes due to mining activity. For validation purposes, we used highresolution Google Earth Pro to identify reference classes for each pixel sample within the Hillendale site, where 100 samples for each category were determined. Google Earth Pro was the preferred data source to ensure the reliability of the training sample due to limited access and data procurement [50]. We randomly selected samples for each year at fixed locations that were adequately spread over the entire mining site to ensure the diversity of samples for each class, as illustrated in Figure 4. We then applied a random forests (RF; [51]) algorithm to classify vegetated and nonvegetated classes, and after that, we validated the performance of our classification within the GEE asset (https://developers.google.com/earth-engine/classifcation; accessed 28 February 2021). This process entails the degree to which a classifier correctly separates a random set of samples [52]. RF is a popular machine learning algorithm of object identification and classification within GEE [48] because of its flexibility, non-parametric nature, and ability to limit overfitting [53], and it has displayed superior performance against other classifiers across a variety of datasets [54]. For the classification, the samples were randomly split into 70% for training and 30% for validating the datasets. From the resultant matrix, we calculated a set of accuracy assessment metrics (Equations (2)-(4)). The producer accuracy (PA) is calculated as the number of vegetated pixels that were correctly classified as vegetation, B is the proportion of the mine that was incorrectly classified, user accuracy (UA) is the number of non-vegetated pixels that were incorrectly classified as vegetated, and D is the proportion of non-vegetated pixels that were classified correctly. The overall accuracy (OA) is computed by dividing the correct pixels by the total number of pixels [55]. We also computed the regions of classified areas for each year within the GEE environment.

Mann-Kendall Test
It is vital to work out the monotonic trends in the time series of any geophysical data. Here, we used the Mann-Kendall (MK) test [56][57][58] to detect the trends in the vegetation and non-vegetated time series extracted in the study area over the study period. The MK test is a non-parametric, rank-based method that is commonly used to extract monotonic trends in the time series of climate data, environmental data, and hydrological data. In this study, the MK test was used for the trend analysis of annual data series for vegetation and non-vegetated area means (in Ha). MK test statistics S can be computed using the formula: where the average value of S is E[S] = 0, and the variance σˆ2 was calculated using the following equation: where, t j is the number of data points in the jth tied group, and p is the number of the tied group in the time series. Under the assumption of a random and independent time series, the statistical function S is approximately normal distributed given that the following Z-transformation equation is used: The value of the statistic S is associated with Kendall's τ = S D where: With respect to the above-defined Z-transformation equation, in this study, there is a 5% confidence level, where the null hypothesis of no trend is rejected if |Z| > 1.96. One other important output of the Mann-Kendall statistic is the Kendall tau τ, which is a measure of a correlation that measures the strength of the relationship between any two independent variables. Furthermore, the Mann-Kendall trend method can be modified to a sequential version of the Mann-Kendall test statistic, which is called the Sequential Mann-Kendall (SQ-MK). The SQ-MK trend test method was proposed by Sneyers [59], and it can be employed to detect approximate potential trend turning points in long-term time series, as well as the dynamics of the trend along with the time frame. The SQ-MK test method produces a forward/progressive trend (u(t)) and a backward/retrograde trend (u (t)) of the time series. In a mathematical form, SQ-MK is calculated by using ranked values of y i of a given time series (x 1 , x 2 , x 3 , . . . , x n ). The magnitudes of y i , (i = 1, 2, 3, . . . ,n) are compared with y i , (j = 1, 2, 3, . . . , j−1). At each comparison, the number of cases where y i > y j are counted and then donated to n i . For effective use of this method, these time series (forward and backward trend) are plotted in the same figure. In the regions where they cross each other and diverge beyond the specific threshold (±1.96 in this study), it is regarded as a statistically significant trend. In addition, the region where they cross each other indicates the time locations where the trend turning point starts [60]. The statistic t i is is calculated using the following equation: The mean and variance of the statistic t i are given by: and: The values of the statistic u(t i ), which are standardized, are calculated using the following equation: The above equation describes the forward sequential statistic which is also called the progressive statistic. The backward/retrograde statistic values (u (t i )) are calculated at the same time, but statistic values are computed by starting from the end of the time series. For the purpose of this study, a 95% confidence level was considered. This method has been used successfully [60][61][62].

Spatiotemporal Patterns of Mine-Induced Disturbance and Recovery
It is worth noting that the extraction operations started in 2001 and ceased in 2013. Here, we began the analysis in the year 2000 as a benchmark to observe progressions of vegetation disturbance until 2013, after which rehabilitation was implemented. The results showing NDVI maps of the Hillendale mine obtained from 2000 through to 2019 are displayed in Figure 5. As expected, the Hillendale mining area was dominated by vegetation with varying NDVI signals in 2000, and this is partly attributable to various classes that were dominated by commercial sugarcane plantations. In 2001, the appearance of considerable mining activity was notable, with reduced NDVI values visible over the southern sections of the study area. From 2002 to 2013, the continuous decline of NDVI signifying progressions of mining severity is apparent. We observed a substantial improvement of NDVI from 2014 to 2019, although some bare patches were disappearing slightly through this time. According to Beukes et al. [5], the restoration process is ongoing as most activities have been accomplished. The remaining process will require some time since the growth of the Eucalyptus trees is expected to increase the vegetation coverage as they mature. This process is projected to stabilize the site to ensure that the closure requirements are met. The recovery pattern remains on a trajectory towards a vegetation-dominated site from a remote sensing perspective but is yet to reach full recovery for the NDVI. Overall, our results show that the rehabilitation of the Hillendale mine has not yet reached its final state; however, they indicate satisfactory growth and are not inconsistent with the government restoration targets, and this sets a good precedent for the currently active Tronox Fairbreeze mining operation, which has a similar background.

Overall and Site-Specific Trends of Disturbance and Recovery
We expanded the above analysis through a detailed examination of the overall NDVI time-series pattern of the Hillendale mine using the DBEST algorithm to give a complete picture of the disturbance and recovery dynamics from 2001 to 2019; the results are displayed in Figure 6. From this figure, the tendency of NDVI to form a V-shaped pattern is apparent. This pattern is expected for progressively recovering mining operations. Higher  We further inspected the site-specific pixel trajectories for three sites within the operational area, and the outcome is presented in Figure 7. Site A represents the RSF-an earth-fill embankment facility formed by the disposal of residuals after extraction and beneficiation of the minerals [63], and is characterized by very fine-textured material and very high salinity, resulting in a much-reduced drainage capacity and leading to prolonged times for drying and stabilizing the substrate [64], as illustrated in Figure 8. As expected, lower NDVI values (~0.1) are apparent throughout the course of the mining operation until a marked increase in 2013 when the rehabilitation began. The vegetation planted on the RSF is doing well and successfully fulfills key objectives: drying the dam, suppressing dust, improving conditions for sugarcane production, and enabling the rehabilitation team to draw valuable lessons [65]. While the greater section of the RSF shows a good vegetation recovery pattern, the rest of the surface is being left to vegetate naturally, which it is doing successfully [41]. Site B indicates a relatively consistent NDVI pattern throughout the study period; however, a slight increment is noticeable from 2013. Site C exhibits a clear disturbance pattern in 2003-2004, with NDVI dropping from~0.6 to 0.1 and remained low through to 2013, when a marked recovery to 0.5 transpired, forming an escarpment. A further increase of NDVI beyond 2013 is inspiring. These results show that the magnitude of disturbance is uneven within the mining site, and so is the temporal pattern of recovery rates; however, the overall picture signifies an effective reclamation.  Our results are consistent with Tronox [41] in that some sections of the mine have successful vegetation plots while some areas are less successful. The evidence highlighted by the cited work suggests that the ineffective plots may be attributable to a high water content, high sand content, or competition from other species, especially Phragmites in the RSF. The walls of the RSF are vegetated with indigenous grass to prevent surface erosion, as illustrated in Figure 8. Further, the experimental Eucalyptus and Casuarina trees established on the RSF and other mined areas show continuous success [41]. The Phragmites australis is the dominant species over the more extensive section of the RSF and is noticeably the best adapted to the conditions there [5].

Areal Extent of Mining and Reclamation between 2001 and 2019
The annual averaged non-vegetated area and vegetated area for the period from 2001 to 2019 was computed for the study area (see Figure 9). The red line represents the non-vegetated area annual mean time series, while the green line represents the vegetated area. The non-vegetated time series increases from its initial value of 63 Ha in 2001 and reaches its maximum peak of 405 Ha in 2013. From 2013 onwards, it dropped to its lowest downward point at 136 Ha in 2019. This is presumably because the Tronox mine started its mining activities in 2001, which meant bare soil began to be exposed. In 2013, the largest portion of vegetation must have been removed, which is why the area of bare soil is at its highest. The vegetation begins to drop in 2001, with its lowest point being in the year 2013 with an area of 96 Ha, presumably due to the increase in mining activities in the study area. In 2013, it began to peak, reaching a maximum of 359 Ha in 2018. For the purpose of the trend analysis of the non-vegetated and vegetated time series, MK and SQ-MK test methods were used. In terms of the MK trend test, the non-vegetated and vegetated time series were observed to have a z-score value of 0.28. This indicates an upwards but not significant trend in terms of the total time series. In order to determine the temporal dynamics of the trend, the SQ-MK test method was used. Figure 10

Accuracy Assessment
We now consider to what extent our classification procedure separated the vegetated and non-vegetated areas within the Hillendale mine site. The overall accuracies of RF classification for mine disturbance and recovery varied slightly above 0.78 across the years from 2001 and 2019 ( Figure 11). These accuracies are satisfactory and comparable with related mining studies. For example, Zhang et al. [66] achieved overall accuracies between 93% and 94% when they examined forest cover from 2006 to 2017 over mining, indicating a successful recovery in a mining region of Nanjing, Eastern China. Xiao et al. [37] also achieved high overall accuracies (0.81 and 0.84) of vegetation disturbance and reclamation from 2004 to 2019 over a surface coal mining region in the Shengli Coalfield in Inner Mongolia, China.

Discussion
The overall goal of rehabilitating mined-out areas is to restore the affected landscape conditions as closely as possible to their pre-mining reference state [67]; in our case, the Hillendale mine had a legal responsibility to restore its site to sustainable sugarcane production land use [40]. Here, we showcased the capacity of the DBEST algorithm to successfully track this development using Landsat and GEE resources. We observed an overall vegetation recovery pattern, pointing towards a successful rehabilitation, although not yet finalized. Similarly, Xiao et al. [37] used the LandTrendr algorithm within the GEE platform and successfully determined land reclamation in open-pit mining areas of a coalfield between 2003 and 2019.
Our results showed that the DBEST algorithm is highly capable of detecting breakpoints in the NDVI time series for the Hillendale mine area. The general reduction of NDVI until 2009-2011 and the recovery from 2013 reflects a continuation of the rehabilitation effort by the mine. In addition, our results reaffirmed the value of GEE-based algorithms for the rapid generation of products for tracking dynamics over mining landscapes. We found classification accuracies as high as 0.78 using the RF algorithm, and this was expected given the spectral contrast between mined-out and vegetated classes within the mining site. These results are consistent with Dlamini and Xulu [36]. The GEE-based resources have increased in use and importance for mining applications as they offer a fast and effective means for computing remotely sensed analysis. Despite this advantage, users are expected to have a basic understanding of programming.
While this approach is quick and affordable, several considerations should be taken into account. First, the use of Landsat data presents a challenge in problematic cloud regions. Second, the NDVI assessment must be placed in the proper context because it is based on NDVI values and not on-field evaluations [68]. The implication of this is that the status of reclamation is based only on vegetative cover and not on the actual vegetation species within the mining area. Therefore, there is a strong need for field verification or knowledge of the vegetation type responsible for the NDVI pattern. From an environmental perspective, the NDVI recovery results from alien plant species as natives may not be considered suitable. Our approach is disadvantaged by these methodological limitations.

Conclusions
In this study, the importance of time-series remote sensing analysis for characterizing mine-induced disturbance and restoration over the Tronox KZN Sands Hillendale mining site between 2001 and 2019 has been illustrated. Our results reveal a decreasing NDVI trend when the operation started in 2001 and deepens to the lowest point over the period 2009 to 2011, with a subsequent recovery coinciding with the cessation of operations and rehabilitation efforts noticeable from 2013 onwards. However, different sites are progressing at varying levels, and this is mainly due to the growth of diverse vegetation types within the site, which include sugarcane, reeds, and the Eucalyptus tree species. Furthermore, our random forests classification was capable of separating vegetated and non-vegetated areas within the affected area with accuracies exceeding 78%.
From this analysis, we conclude that the DBEST algorithm is a valuable tool for detecting patterns of disturbance and restoration over surface mining environments, especially Hillendale. Our study showed uneven recovery patterns for selected sites within the mining site; however, the overall picture indicates an effective restoration. Moreover, our affordable approach can better extract mine-induced disturbance and restoration patterns with a high degree of fidelity and be replicable in many areas undergoing intense vegetation transformations, particularly in data-poor countries with limited resources. We suggest, however, that a complete translation of the DBEST algorithm code into a GEE platform could enhance its application as it would be more easily accessible for the broader community.

Data Availability Statement:
The datasets used and/or analyzed in this study are available from the corresponding author on reasonable request.