Analyzing Gradual Vegetation Changes in the Athabasca Oil Sands Region Using Landsat Data

Oil sand mining in northern Alberta/Canada in the Athabasca region is a major intrusion into the otherwise pristine natural environment. The various types of oil sands mining, transport, and processing are causing large-scale discharge of pollutants. Accordingly, this study examined the gradual changes in the physically undisturbed vegetation, that occurred from 1984 to 2021 in the Athabasca oil sands monitoring region. First, the abrupt changes were masked out with the help of auxiliary and Landsat data. Subsequently, a normalized burn ratio Landsat time-series was applied to the LandTrendr algorithm on the Google Earth Engine. In order to interpret gradual changes, measurement criteria were used to describe vegetation development, vulnerability, and variability. In addition, the spatial and temporal relationship of these to oil sand opencast mines, processing facilities, and steam assisted gravity drainage (SAGD) mines was examined. The results showed that a major part of the vegetation in the Athabasca oil sand monitoring region underwent a positive development (65.9%). However, around the opencast mines a negative vegetation development and stability within a radius of 10 km could be observed. In the surroundings of processing facilities, the development and stability of vegetation was disturbed within a radius of 2 km. Thereby the analysis of land cover classes showed that deciduous, coniferous, and mixed forest are disproportionately affected. Conversely, no negative influences on neighboring vegetation could be detected around SAGD mines. The temporal analysis showed that vegetation disturbance was most pronounced between 1990 and 2000, but recovered in recent years.


I. INTRODUCTION
I N THE last 50 years, the exploitation of the Athabasca oil sands in Northern Alberta/Canada has increased substantially. Driven by growing global demand for petroleum products and the development of advanced technologies, an extensive industry of extraction, processing, and transportation of oil-related products has grown in the boreal forest [1]. The oil sands at shallow depths are extracted by opencast mining, but in greater depths steam assisted gravity drainage (SAGD) is used [2]. Subsequently, a resource-intensive process takes place to extract the crude oil from the oil sand. This goes hand in hand with extensive land use, high water consumption, and pollution of the surroundings by exhaust gas and leaking Manuscript  tailing ponds [3]. Due to this development and natural factors, vegetation disturbances occur, caused for example by mining, oil processing and transportation, pollution, forest harvesting, urban expansion, wildfires or insect infestations [4]. These can be divided into abrupt, gradual, and seasonal changes [5]. Abrupt changes are mostly due to forest harvesting, mining or urban expansion [6]. Seasonal changes are associated with phenological period of vegetation and can mostly be attributed to climatic changes [7]. Gradual changes, on the other hand, indicate vegetation disturbance due to pollutant input or insect infestation [8]. Most studies of vegetation development tend to focus on abrupt changes [9]. Due to the strong change in spectral reflectance, this is considerably easier to detect. By analyzing gradual changes, disturbances in the vegetation can be detected at an early stage and thus facilitate ecosystem management and nature conservation by allowing preemptive remedial action [8]. Particularly in the Athabasca oil sands region, where the natural environment has been significantly altered, the physically undisturbed part of the vegetation should be monitored, as Canadian Law mandates oil companies to ensure land restoration after use [10]. In this context, there has been plenty of research on the input of pollutants and trace elements from of oil sand mining operations and their impact on the living fauna [11], [12], on rivers and lakes [13], [14], [15], as well as on the surrounding vegetation [16], [17], [18]. Furthermore, the spatial and temporal patterns of trace element deposition using sediment cores were reconstructed [19], [20]. In addition, geographic information systems were used to combine several of these studies and derive spatial structures of the pollutant input [21]. In the field of remote sensing, some studies have already tried to detect these possible gradual vegetation disturbances around the opencast mines in the Athabasca oil sand region. The first focused on land cover change via Landsat and used Advanced Very High Resolution Radiometer from the National Oceanic and Atmospheric Administration (NOAA-AVHRR) to monitor vegetation health [22]. Gillanders et al. researched the land cover development in the immediate vicinity of the opencast mines and could thus determine what kind of vegetation was affected by the mine expansion and will be in the future [23]. Latifovic and Pouliot used Landsat and MODIS data and the normalized difference vegetation index (NDVI) and normalized difference water index (NDWI) to determine gradual changes within the vegetation using linear regression between 1984 and 2012, while using auxiliary data to mask out abrupt changes [24]. In a further analysis, Pouliot and Latifovic in 2016 used temporal segmentation to classify different types This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ of vegetation disturbance and recovery in the vicinity of the opencast oil sand mines [4]. More recent remote sensing studies investigated surface displacement due to SAGD operations [25] and enhanced peatland classification with radar data [26]. Most of these studies focus on the surrounding vegetation of oil sand opencast mines and do not conduct a large-scale monitoring of the entire Athabasca oil sand monitoring (OSM) region. Therefore, little attention has been paid to the effects of other mining methods such as SAGD mines and the effects of processing facilities.
This leads to the following research questions. 1) What gradual changes can be observed in the vegetation within the Athabasca OSM region that has not been influenced by abrupt transformation processes? 2) How do these gradual changes behave in temporal and spatial terms in the vicinity of oil sand opencast mines, processing facilities, and SAGD mines? Accordingly, the objective of the study is to investigate the gradual changes and trends in vegetation between 1984 and 2021 in the Athabasca OSM region and to examine the spatial and temporal implications of different oil sand mining and processing methods [27].
The following analysis was carried out in the Google Earth Engine (GEE), which provides access to the entire Landsat archive and can improve calculation time through parallel processing [28]. In a first step, the part of the vegetation that underwent an abrupt transformation is masked out with the help of auxiliary data. The Landsat-based detection of trends in disturbance and recovery (LandTrendr) algorithm provided the basis for the investigation of the gradual changes [29]. To describe the gradual changes, the physical quantities velocity, frequency, and variance were calculated to evaluate gradual changes and describe vegetation development, vulnerability, and variability. These results were analyzed within different land cover classes and buffers created around oil sands opencast mines, processing facilities, and SAGD mines. For the temporal investigation, the mean velocity of six-year time periods within these buffers was calculated.

A. Study Area: The Athabasca Oil Sand Region
The study area is the Athabasca OSM region Fig. 1, which covers an area of 92 885 km 2 in northern Alberta. This is one of four regions of the Oil Sands Monitoring Program established by the governments of Canada and Alberta. The aim of the OSM is to study the environmental impacts of oil sands mining and the effects on indigenous people [27]. The Athabasca OSM region was used as study area, as this region contained different kinds of oil sands mining and large amounts of untouched vegetation. The climate of the site is characterized as subhumid with cold winters and hot summers with long wet and dry periods. The annual average temperature is 1.0 • C and total average precipitation is 418.6 mm [30]. The study area is characterized by boreal forest permeated by wetlands, lakes, and the from south to north crossing Athabasca River [31], [32]. A human footprint can be mainly found in mining activities, urban expansion in the nearby city Fort McMurray and agriculture and logging in the south.
Oil sand is a mixture of bitumen, water quartz sand and clays. Bitumen is a form of heavy crude oil, which have to be extracted from the surrounding soil materials to use it for further products [33]. Shallow oil sand can be mined by opencast mining and extracted by mixing the oil sand with hot water using natural gas [34]. However, two thirds of the deposits are too deep to mine it from the surface and thus, the SAGD technology is used [2]. Here, natural gas is mixed with water to produce steam, which is pumped into the deep-seated oil sand deposits [35]. Therefore, the viscosity of the oil sand is increased, which is then plumbed to the surface by production wells [3]. The different deposits can be seen in Fig. 1, where the dominance of SAGD technology is evident.

B. Data Collection
Annual surface reflectance composites were produced for the period between 1984 and 2021 using Landsat TM/ ETM+/ OLI, which provided a continuous resolution dataset of land surface observations of 30 m to date. For this purpose, USGS Landsat surface reflectance Level 2 Collection 2 Tier 1 datasets were selected, which have already undergone atmospheric and geometric correction. For each year all available data between June and September was used (see Sections III-A and III-B). This period represents the main growing season and therefore, it is likely that the spectral contrast between stable forest patches and gradual vegetation changes is highest. The Human Footprint Inventory 2019 enhanced for OSM Region (HFI) and Historical Wildfire Perimeter Data: 1931-2020 (HWP) are the basis for masking abrupt changes (see Section III-C). The HFI consolidates 21 human footprint categories, based on more than 115 anthropogenic disturbance types, into a single integrated dataset. It was created by intersecting various data (e.g., infrastructure data), by thematic mapping based on SPOT 6 satellite imagery from 2017 to 2019 [36]. The HWP provide shapefiles about the extent of wildfires from 1931 to 2020, here only the data from 1970 onwards was used, as the areas further in the past have already gone through extensive regeneration [37]. Data extracted from the HFI were also used to analyze the spatial and temporal relationships between vegetation changes and oil mining and processing methods. In order to compare the results on different land cover classes, the forest land cover maps of Hermosilla et al. from 2019 were used [38]. From the OSM Regions dataset, the Athabasca region was selected as the study area. Derived out of the Oil Sands Project Boundaries available between 1985 and 2015, the most recent data was used, to provide an overview of the spatial distribution of different types of oil sand mining in Fig. 1. Table I summarizes the data used herein.

A. Data Preprocessing
Annual Landsat composites were generated between 1984 and 2021. To harmonize differences in spectral characteristics between ETM+, TM and OLI, the OLI image bands 2, 3, 4, 5, 6, and 7 were transformed to the spectral properties of ETM+ bands 1, 2, 3, 4, 5, and 7, respectively, using slopes and intercepts from reduced major axis regressions reported in Roy et al.   Table II [39]. Afterward, CFMASK was applied to mask cloud, cloud shadows, snow, and water [40]. All data of a year were taken to create a composite, conducted by using the median of each spectral value [24]. Then, spectral indices were calculated, which were input for the LandTrendr algorithm and for the vegetation mask described below Fig. 2(a).
To ensure that only areas that were covered with vegetation during the study period were included in the analysis, a vegetation mask was created. For this, an NDVI threshold was defined that separates vegetation from nonvegetated areas, a widely used method to measure relative vegetation abundance [41]. Therefore, a median NDVI composite was formed in 2019 between  June and September and then intersected using the land cover product of Hermosilla et al. also from 2019 [38]. The resulting histogram is shown in Fig. 3. Based on this, the threshold value for vegetation masking was set to 0.4 in order to include all vegetation classes in the analysis. Accordingly, all areas that fell below this value between 1984 and 2021 were masked out, and the resulting masks were merged Fig. 2(a).

B. Change Algorithm and Spectral Disturbance Index
The LandTrendr algorithm was selected as change algorithm to detect gradual changes within the time series. The algorithm has been developed for application to forest disturbances and has since become one of the most widely used algorithms for the detection of abrupt and gradual forest changes [42]. In addition, the implementation in GEE enables processing large amounts of data lending scalability [29].
LandTrendr is a temporal segmentation algorithm based on a yearly time series of spectral indices, to capture both gradual and abrupt changes. It simplifies the time series using a model and detects as breakpoints inflection points in the curve, where the current trend reverses or intensifies. The trends are identified using linear regression, reducing noise by smoothing out outliers without losing significant amounts of detail. At the breakpoints, a temporal segment ends and a new one begins. The algorithm outputs include a linear regression model fitted to the spectral data and the detected breakpoints [42]. Since there were no in situ data in which the fitted models could be compared with ground truth, the parameters were left at the default settings.
A Landsat time series of the normalized burn ratio (NBR) was used as an input for the LandTrendr algorithm Fig. 2(a). NBR is the normalized ratio between near-infrared and short-wave infrared (SWIR) bands. The possibilities to characterize forest structures with SWIR wavelengths are well understood [43]. Kennedy et al. showed that the NBR was most effective spectral index to detect forest disturbances in Landsat time series [42], [44], due to the high sensitivity toward moisture [45]. Despite the moisture sensitivity being an asset for understanding vegetation disturbance, changes in soil moisture could influence the estimates of vegetation development. However, areas with a high proportion of bare soil and outliers in terms of annual wetness were masked out, we expect the influence on the quantitative vegetation development estimates to be minimal.

C. Masking Areas Under Influence of Humans or Natural Disasters
After masking out the nonvegetated areas, only those vegetated pixels should be part of the analysis that is not affected by human-induced changes or natural disasters. For this purpose, the following areas were, thus, masked out 1) Areas under direct human influence.
2) Former wildfire areas, which show afterwards a strong positive vegetation development caused by regeneration. 3) Other areas that are suffering abrupt changes, e.g., windthrow. Pipeline right of way, seismic lines, access roads, and other oil and gas footprints going through the reclamation process may also have vegetation (likely herbs and bryoids). However, these obviously do not mean actual vegetation regeneration and should not contribute to the present analysis. Accordingly, all human influenced areas were masked out a priori using the HFI. To mask out former wildfires in the study area the HWP were used Fig. 2(a). In order to identify further areas with abrupt changes, only changes that lasted one year and exceeded a certain magnitude of change were detected by the LandTrendr algorithm. Therefore, a threshold value t had to be determined. For this purpose, the d NBR was introduced, which represents the difference in the NBR in two consecutive years. Then, 1000 sample points were randomly selected and the mean m and standard deviation std of the d NBR over all years of the study period at these points were used to calculate the threshold value t = m + 2 * std. (1) All vegetation losses within a year that were greater than this threshold were not included in the analysis [46]. Vegetation regeneration did not have to be considered in this case, as it cannot occur at such a fast rate and is therefore smoothed by the LandTrendr [42]. All resulting masks were combined and provided the basis for the following steps.

D. Analyzing Vegetation Changes
In order to evaluate gradual changes in the ecosystem, an approach of Ye et al. was adapted and modified [9]. There, the use of physical terms Velocity V , Frequency F , and Variance S was proposed for the description of vegetation change measurement criteria Fig. 2(b). In this article, the calculation and use of V in particular is expanded and changed from the original approach. In addition, a spatial and temporal investigation is carried out with regard to different types of oil sand mining and extraction.
The measured value of V describes the average annual rate of change in vegetation. Negative values indicate a loss of vegetation, while positive values signify an additional gain. In order to get an overview of the changes over time, changes were calculated in six-year periods (first and last period include seven years) by using the fitted LandTrendr data [47]. In the following, an average over the entire study period was formed from these sections where V 1 to V 7 are the average annual change rate in vegetation in one period, fit year is the fitted LandTrendr data of the displayed year, t 1 to t 7 are the time in years for the respective period and V is average over all periods. To simplify interpretation, the result was reclassified into four levels: loss ≤ −2, weak loss ≤ 0, weak gain ≤ 2 and gain > 2 The measurement criterion F describes the stability or vulnerability of the ecosystem. Low values indicate stability and high values vulnerability [48]. This is derived from the number of fitted segments generated by the LandTrendr algorithm. These can show a positive or a negative development. The results will include the total number, number of positive, and negative segments in the investigation period. A positive segment can be interpreted as a period of regeneration and a negative segment as a period of disturbance. The measurement criterion S describes the variability of vegetation development, where high values indicate high variability. This is derived from the variance of the time series. S is calculated as follows: where x i represents the source data in a pixel in year i, x i represents the mean in year i and n represents the number of years in the investigation period. In order to understand the type of vegetation regeneration within the vicinity of the mining sites, the land cover product of Hermosilla et al. [38] was used to calculate a histogram of the respective values. Only five classes were used (Broadleaf, Coniferous, Mixedwood, Wetland, Wetland-Treed), as the others only covered a negligible area within the study region.

E. Spatial and Temporal Distribution of Gradual Changes
In order to analyze the spatial and temporal distribution of the described measurement criteria and to understand the impacts on the vegetation, the surroundings of oil sands opencast mines, processing facilities and SAGD mines were evaluated Fig. 2(c). The data basis is formed by the HFI; this consists of shapefiles where the attribute FEATURE_TY indicate the anthropogenic disturbance type. For the opencast mines, the attributes MINES-OILSANDS, TAILING-POND, and RIS-MINING-OILSANDS were selected, what corresponds to areas with bare and/or low vegetated ground and low human density for the purpose of oil sands mining and tailings ponds. In place of oil processing facilities, the attribute OIL-GAS-PLANT [36] was selected, under which different industrial sites used for oil production were summarized. For SAGD mines, the attribute WELL-BIT was selected, what represents bitumen well pads. The shapefiles were checked for correctness by visual interpretation using satellite images from Google Earth. To quantify the spatial distribution buffers of different sizes (0.25 km, 0.5 km, 1 km, 2 km, 3 km, 4 km, 5 km, 10 km) were used and afterwards the mean and the standard deviation (stdDev) of the measurement criteria of all underlying pixels were calculated. Further, to examine the temporal distribution of the gradual change, the mean of the measurement criteria V of all pixels inside the buffer within the six-year periods were calculated.

F. Distribution of Abrupt Changes
By masking out areas under human influence, previous wildfire areas, nonvegetated areas and areas undergoing abrupt changes, an area of 51 766 km 2 was excluded from the analysis. Fig. 4 shows the distribution of the different masks used. The mask of forest fires Fig. 4(b) demonstrate the greatest extent, particularly evident in the peripheral areas in the north, east, and west. In addition, there is the mask of human impact Fig. 4(a), which is dominant in the north-east. The nonvegetated and abrupt changes mask Fig. 4(c) and (d) mostly coincide with the previously mentioned ones. The nonmasked area of 41 119 km 2 , which is subsequently analyzed, is mainly located in the center of the study area.

IV. RESULTS
Hereafter, the results of the analysis of the measurement criteria Velocity (V ), Frequency (F ), and Variance (S) are described. Afterward, the characterisation of the spatial and temporal distribution of gradual changes is presented.

A. Analyzing Annual Change Rate
The velocity (V ) was used to describe the average annual vegetation change over the past 37 years. The spatial and temporal distribution of V is shown in Figs. 5 and 6. The spatial distribution of V is presented in Fig. 5, which shows that the negative trends are particularly concentrated in the north-eastern part around the opencast mines. In contrast, positive trends dominate in the western and southern parts. The proportion of positive and negative vegetation changes over the entire study period can be seen in Fig. 6(a). While positive changes clearly dominate (65.9%), a substantially smaller proportion (34.1%) showed a negative trend. A clear trend can be observed in the temporal distribution of V in Fig. 6 parts b to h. From 1984 to 2000, negative trends increased, especially in the center and north-east of the study area. From 2000 onward, this trend was reversed, especially in the center of the study area. It is noteworthy that negative trends continue to persist in the north-east.

B. Analyzing Stability and Vulnerability
The Frequency (F ) is an indicator for the stability and vulnerability of the vegetation Fig. 8. Here, low values represent the stability and high values vulnerability of the ecosystem, which is derived form the number of fitted segments by the LandTrendr algorithm. An overview of the spatial distribution of the number of all segments is given in Fig. 8(a). Fig. 8(d) provides the corresponding histogram. A large part of the area is characterized by low F values. Moreover, the spatial distribution does not show any visible trends, except for minor local variations. The number of segments broken down by loss and gain is shown in Fig. 8(b) and (c). It can be observed that a considerable amount of the pixels contain only one loss segment, whereas the amount of pixels containing two or more gain segments is substantially higher. Fig. 8(d) shows that most of the pixels contain one segment, then the number of pixels decreases with each further segment, with the exception of three segments, where an increase can be seen. Fig. 8(d) shows also the number of pixels per land cover class, whereby it becomes evident that the Coniferous, Wetland-Treed, and Broadleaf classes has a disproportionate number of pixels with three segments.

C. Analyzing Variability
The variance (S), where high values refer to a larger variability in the vegetation development, is depicted in Fig. 9. The spatial distribution of S in Fig. 9(a) shows that the majority of the study area is dominated by low S, interspersed with isolated areas of high S. High S is remarkably frequent in the center, south and north-east around the opencast mines. These observations can be confirmed by the variance within land cover classes in Fig. 9(b): A small part shows values close to zero S, while the majority of pixels exhibit low S, followed by monotonically decreasing number of pixels with increasing S. There are no noticeable differences between the classes, except for the wetland class, which shows a much later increase and therefore a higher overall variance value.

D. Spatial Distribution of Gradual Changes
The spatial distribution of the measurement criteria V , F, and S is represented by the mean and standard deviation within buffers around the oil sand opencast mines, processing facilities and SAGD mines is shown in Fig. 10. Within V Fig. 10(a), a positive dependency of proximity can be observed for opencast mines and processing facilities; this is clearly more pronounced for opencast mines and the negative influence persists up to a distance of 10 km, while for processing facilities a negative dependency can only be seen up to a distance of 2 km. In the case of SAGD mines, a negative dependency is visible up to a distance of 1 km. There seems to exist no dependency between  proximity and the development of F Fig. 10(c), but it is striking that F is considerably lower in the vicinity of SAGD mines than at processing facilities and opencast mines. S shows a negative dependence of proximity for opencast mines Fig. 10(e), processing facilities and SAGD mines. The negative dependency is most pronounced for opencast mines. The stdDev shows a decline with increasing distance in all cases Fig. 10(b), (d), and (f). Here, opencast mines show the highest stdDev followed by processing facilities and SAGD mines.

E. Temporal Distribution of Gradual Changes
The temporal distribution of gradual changes in the surroundings of oil sand opencast mines, processing facilities, and SAGD mines is shown in Fig. 11. In the vicinity of the opencast mines Fig. 11(a), the development of enveloping vegetation is consistently negative between 1990 to 2010, although a slight positive dependency can be observed, which is particularly strong between 1990 and 2000 at short distances. The periods from 1984 to 1990 and from 2010 to 2021 show a positive vegetation development, but no specific dependence. In the vicinity of processing facilities Fig. 11

V. DISCUSSION
In this study, gradual changes in physically undisturbed vegetation were identified using LandTrendr algorithm and measurement criteria V , F, and S. To achieve this, abrupt changes were masked out using different approaches. The most prevailing abrupt changes were forest fires and mining activities. Similar results for abrupt changes were obtained by Pouliot and Latifovic in 2016, where the gradual changes showed that insect infestations played an important role [4]. Gradual changes can also increase with stress from forest disturbance due to clearing and expansion of industrial activities [4]. As these could not be masked out due to the lack of in situ data, insect infestation may therefore be part of the gradual changes detected in this study. Another unnoticed aspect that influences the gradual changes is climate. The climate normalization of a gradual trend analysis of the vegetation by Latifovic and Pouliot in 2014 showed that the climate or climate change only slightly influences gradual changes [24].
The investigation of V showed that the vegetation in undisturbed parts of the Athabasca OSM was largely subject to a positive development, although it must be acknowledged that some of this detected regeneration may have been contributed by invasive species. However, as satellite imagery from Landsat was used here, and it was not possible to get corresponding high spectral or spatial resolution data, the specific contribution of invasive species could not be estimated [49], [50]. The detected negative development was spatially concentrated in the area around opencast mines in the northeast, which was additionally confirmed by the spatial investigation carried out in Section IV-D. These demonstrated a negative influence on vegetation up to a distance of 10 km and beyond. It also showed that processing facilities have a negative impact on the surrounding vegetation at distances of up to 2 km, whereby the analysis of the land cover classes showed a disproportionate impact on deciduous, coniferous, and mixed forest. Conversely, no negative influence could be detected in the vegetation surrounding SAGD mines using the methods proposed herein. One of the reasons for this could be that land disturbances from SAGD mining are typically observed off-site and vegetation impacts are likely displaced upstream (see [51] for further details). As only the vegetation changes in the immediate vicinity of mining sites was investigated herein, the upstream and off-site vegetation disturbances could not be quantified. These observations were further supported by the measurement criteria of F and S, which indicate variability and stability of vegetation development. No dependence between proximity and the value of F was found, but it is worth noting that both the surroundings of opencast mines and processing facilities showed substantially higher values and thus weaker stability in the vegetation than SAGD mines. On the other hand, S decreased with increasing distance in all studied areas; this could be attributed to a smoothing effect with a larger amount of pixels within larger buffers. SAGD mines and processing facilities showed a very similar development, while opencast mines showed a notably higher value of S. In summary a negative vegetation development prevailed in the vicinity of opencast mines and that this was accompanied by high variability and unstable vegetation. Similar results in vicinity to opencast mines were also observed by Latifovic et al. [22] and by Latifovic and Pouliot [24]. Indeed, it would be valuable to validate the results, even if only for a small part of the study area, but unfortunately, no suitable datasets could be found. In future, however, curating field-based and high-resolution (e.g., drone/aerial photography) vegetation regeneration observations could help validate such studies, contributing to the development of satellite-based systematic vegetation monitoring systems. In relation to the negative vegetation development around the opencast mines and processing facilities, numerous studies have investigated the impact of atmospheric deposition of acid [16], [17], [52], polycyclic aromatic hydrocarbons (PAHs) [53], [54] or trace metals [19], [20] emitted by oil sand mining and processing in the surrounding vegetation. These are mostly released into the environment by dust particles and increased concentrations are related to the vicinity to opencast mines and lead to an impact on vegetation growth and species distribution [13], [55]. Another factor is contamination of the surrounding water system. For example, large quantities of mostly toxic trace metals emitted during the oil sand processing are released into the rivers during the snow melt [14], [34].
It has also been identified that the tailing ponds, which are supposed to absorb incoming waste water, are leaking it into the groundwater [56]. Furthermore, the impact of oil sands mining on the climatic conditions in the vicinity of open pit oil sands mining was also investigated, revealing that the average temperature increased by 1.2 • C compared to the regional average [57]. This combination of trace metals, PAHs, contaminated water and changes in the regional climate could have led to the negative impacts found around opencast mines and processing facilities. Ground movement (heave/subsidence) due to SAGD operations could also impact the spatial and temporal distribution of vegetation changes. However, as the heave and subsidence is limited to a few centimetres per year and stop after a few years [58], [59], it has been assumed that this has limited influence on the vegetation. Indeed, this also presents an interesting opportunity for future research, although it is out of scope for this particular study.
The temporal development of V in Sections IV-A and IV-E showed that between 1990 and 2000, the majority of the vegetation in the study area underwent a negative development. The periods from 1984 to 1990 and from 2010 to 2021, on the other hand, clearly showed a positive vegetation development. In the period from 1984 to 1990, this can be explained by more distant mining activity, which has steadily expanded since then. However, the regeneration of vegetation development in the period from 2010 to 2021, could be attributed to improvements to mining and upgrading technologies, what leads to a reduced trace metal input in recent years [19], [20].

VI. CONCLUSION
This study investigated gradual changes in physically undisturbed vegetation within the Athabasca OSM region between 1984 to 2021. For this purpose, the LandTrendr algorithm was applied to a Landsat time series and three different measurement criteria were calculated from its results (Velocity V , Frequency F, and Variance S). The results were then examined in the vicinity of oil sand opencast mines, processing facilities, and SAGD mines. A major part of the vegetation in the Athabasca OSM region showed a positive development (65.9%). However, around the opencast mines in the northwest, a negative vegetation development and stability was observed, with strong disturbances noticeable up to a radius of 10 km. In the surroundings of processing facilities, the development and stability of the vegetation was disturbed within a radius of 2 km, while around SAGD mines, no effect on vegetation was detected. Thereby the analysis of the land cover classes showed that deciduous, coniferous, and mixed forests are disproportionately affected. The negative vegetation development around oil sand opencast mining and processing facilities was attributed to the input of various pollutants and trace metals emitted. Whereas the impacts of SAGD mines on vegetation are typically off-site and upstream and therefore could not be quantified in this study. Furthermore, it was evident that vegetation disturbance was most pronounced between 1990 and 2000, which eventually recovered.
This study showed the viability of assessing gradual vegetation changes using GEE with the LandTrendr algorithm, as it allows to quickly and efficiently work on large areas and long time series. Future research should continue to analyze vegetation disturbances resulting from oil sand mining operations, especially SAGD mines.