Assessing Changes in Boreal Vegetation of Kola Peninsula via Large-Scale Land Cover Classification between 1985 and 2021

Abstract

The effective monitoring of boreal and tundra vegetation at different scales and environmental management at latitudes above 50 degrees North relies heavily on remote sensing. The vastness, remoteness and, in the case of Russia, the difficulty of access to boreal–tundra vegetation make it an ideal technique for vegetation monitoring in the Kola peninsula, located predominantly beyond the Arctic circle in the European part of Russia. Since the 1930s, this area has been highly urbanised and exposed to strong influence by a number of different types of human impact, such as toxic pollutions, fires, mineral excavation, grazing, logging, etc. Extensive open archives of remote sensing imagery as well as recent advances in machine learning further enable the efficient use of remote sensing methods for assessing land cover changes. Here, we present the results of mapping northern vegetation land cover and changes in it over a large territory, in time and under human impact based on remote imagery from Landsat TM, ETM+ and OLI. We study the area of about 37,000 km${}^2$ located in the central part of the Kola peninsula in the boreal, pre-tundra and tundra between 1985 and 2021 with a time interval of approximately 5 years and confirm the correlations between the human pressure and the level of vegetation changes. We put those into the perspective of year-on-year changes in the temperature and precipitation regimes and describe the recovery of the damaged original boreal vegetation (dominated by spruce) through pine and deciduous vegetation. As a by-product of this study, we develop and test an approach for the semi-automated processing and classification of Landsat images using the novel TensorFlow machine learning technique (widely spread across other disciplines) that enables high-throughput classification, even on conventional hardware.

1. Introduction

Since the beginning of industrialisation, environmental changes have become so significant and rapid that they can no longer be fully balanced and recovered by natural adaptation and connection [1,2,3]. The timespan of ecosystem degradation and recovery is now driven by the combined effects of natural processes and human activity. This is particularly the case for northern ecosystems located at latitudes above 50 degrees North as they are more vulnerable and take longer to recover due to severe climatic conditions [4,5]. Boreal ecosystems occupy about 17% of the global land surface [6] and constitute around one-third of the world forest area [7], making them an important contributor to the formation and regulation of global environmental processes through their influence on the carbon cycle, water balance and albedo changes [5].

Natural environmental factors, such as air temperature, precipitation and solar radiation, influence vegetation cover directly. These changes are generally gradual, giving vegetation time to adapt. However, boreal vegetation has contrasting degrees of stability to industrial impact, e.g., forests are more stable than pre-tundra [8]. There is generally a strong correlation between changes in boreal vegetation and the human impact on air, water and soil through industrial development, forest logging, tourism, fire, reindeer grazing, etc. [9,10,11]. Particularly in the areas of significant human impact, the extreme environmental conditions of the North lead to long vegetation recovery times [12], making such impact highly visible and measurable. The ability of boreal vegetation to react to environmental fluctuations and human impact makes it an excellent indicator for environmental changes in general and human impact on ecosystems in particular. Such changes are clearly visible in the amount of the ground phytomass, number and variety of plant species as well as tree dryness and diseases [13].

Here, we study the changes in the boreal vegetation located in a highly urbanised area that has been heavily industrialised over the last 100 years. We investigate how the boreal vegetation reflects the trends in changing climatic factors (temperature and precipitation) and how vegetation damaged through human impact recovers when the levels of influence fluctuate dramatically.

Our study focuses on the Kola peninsula located in the European Russian North, almost completely above the Arctic Circle. It covers 0.85% of the total area of Russia or 9.77% of Fennoscandia [14]. Boreal forest covers 67% of its territory. Of that, 72.5% is coniferous forests with pine and spruce [15], with birch in between, and dwarf shrubs, grass, mosses and lichen forming the field layer. Due to the narrow crowns and tree sparseness, these forests are light with a well-developed undergrowth [16]. The boreal forests of this territory are broadly divided by wetlands and lakes and do not combine into vast single-species forest areas.

High levels of industrialisation, rapid changes in economic developments at all levels and the proximity of national borders make the central part of the Kola peninsula an attractive area for studying boreal vegetation. This is why we focus our study on the area of approximately 37,000 km${}^2$ located in the central part of the Kola peninsula (Figure 1).

There had been significant economical, political and cultural changes in this area over the course of the 36 years from 1985 to 2021, driven primarily by the overall changes in the political and economical situation in Russia. This is why this study spans this date range. The long response times of boreal vegetation to external factors and long-term changes in human impact in this area make the 4–6 years horizon a long enough span to catch land cover changes in the state of the boreal forests here. Over this time horizon, we can also observe the response of the vegetation to significant ecological, economical and political changes.

Many research groups studied the industrial impact on the terrestrial and water ecosystems in this area by various methods from the late 1980s until the early 2000s [9,10,13,17,18,19,20,21]. This was made available by the opening of the Russian borders after the collapse of the USSR and financial support from various scientific organisations and local industrial enterprises in this time period and contrasts with the limited availability of studies on the state and changes in boreal vegetation in later years. This was another reason for our interest in Kola.

The vastness, remoteness and, in the case of Russia, the difficulty of access to boreal vegetation at latitudes above 50 degrees North make remote sensing an ideal technique for monitoring vegetation changes in these areas. It is a well-established technique for the large-scale monitoring of land cover specifically because it permits retrospective investigations of vast areas using data from extensive open archives of historical imagery. In this study, we use remote sensing to map the land cover through the semi-automated classification of Landsat TM, ETM+ and OLI satellite images at the frequency of approximately 5 years. E. Sklyar [22] first introduced this approach using Landsat imagery and the support vector machines (SVM) classification technique in 2013; however, the study area was limited to a series of smaller regions than in this study for technical reasons. The recent advances in machine learning coming with the development of the Tensorflow machine learning technique [23] enabled us to scale the investigation to a much larger area in this study.

With this study, we introduce TensorFlow [23] for the land cover analysis. Over the course of the past 10 years, there has been a substantial shift in machine learning applications across industries towards deep learning neural networks, such as TensorFlow. Neural network-based methods have undergone a strong evolution during this time and now they generally provide a higher accuracy, stability and performance for complex problems compared to more classical algorithms, such as the support vector machines (SVM) that we used earlier [22].

In this study, we develop and analyse land cover classification maps for the years between 1985 and 2021. With these maps, we aim to answer the following two questions that result from observing the changes in the climate and human impact over the period of the study:

• Whether there is a trend to the phytomass growth over the time period of the study attributable to rising temperature and precipitation levels;
• Whether there is a large-scale recovery of the damaged vegetation following the reduction in industrial activity towards the end of the 1990s.

2. Environment and Human Activity on the Kola Peninsula

The phenological phases of vegetation (that is, the seasonal cycle of plant activity) are controlled by hydro-thermal conditions. They depend mainly on the levels of air temperature, moisture, solar radiation, snow depth, the spring snow-melt length and the onset timing for cold, late summer and autumn days [24,25]. Positive correlations are observed between the warm climate and the amounts of precipitation and sunlight on the volume of phytomass as well as on the number and variety of plant species [13].

2.1. Seasonal and Year-on-Year Trends in Air Temperature and Precipitation

An analysis of the meteorological data between 2005 and 2021 shows that steady positive air temperatures, >5 °C, at 2 m above ground are observed from around May to October, peaking at about +14.80 °C on average in July (Figure 2a). July is also the month when liquid precipitations reach their maximum (Figure 2b). The snow melting around April and May (Figure 2c) contributes to the general water supply to the vegetation towards the summer months. The meteorological conditions make July the preferred time of the year for studying land cover and vegetation. The retreat of cloudiness that dominates the northern areas most of the year also makes July the most suitable month for remote sensing studies.

Based on the meteorological data, we observe statistically significant trends of increasing air temperature and precipitation in the 50-year period from 1970 to 2021. The mean air temperature (Figure 3a) rose from about –0.9 °C at the start of the 1970s to about +1.2 °C at the end of the 2010s ($\tau =0.4$ with $p=0.00005$, according to the Mann–Kendall test [26]). Interestingly, the winter temperatures increased more significantly than the summer ones. We also observe an increase in the total amount of warmth that supports vegetation growth (computed as a sum of the daily mean temperatures above a selected cutoff value). For the cutoff of +5 °C, the increase from 1970 to 2021 is around 20% based on the data by Veselov et al. [27]. The amount of the annual precipitation (Figure 3b) rose about 10%, from around 460 mm in 1970 to around 505 mm in 2021 ($\tau =0.28$ with $p=0.0057$). The increase in the average air temperature and the amount of precipitation resulted in an earlier onset of the vegetation green-up and to a longer duration of greenness. This is why we expect a higher density and productivity of forests, as well as increased volumes of phytomass towards 2021 [28,29].

Figure 2. Mean seasonal distributions of meteo data acquired at the Monchegorsk meteo station (ID 22212) between 2005 and 2021: (a) mean monthly air temperatures at 2 m above ground [30]; (b) mean monthly precipitation [31]; (c) mean monthly snow depth [32]. The date range of 2005–2021 is the one for which all three measures are consistently available for this meteo station; we use wider date ranges in other figures where only temperature or precipitation are discussed.

Figure 2. Mean seasonal distributions of meteo data acquired at the Monchegorsk meteo station (ID 22212) between 2005 and 2021: (a) mean monthly air temperatures at 2 m above ground [30]; (b) mean monthly precipitation [31]; (c) mean monthly snow depth [32]. The date range of 2005–2021 is the one for which all three measures are consistently available for this meteo station; we use wider date ranges in other figures where only temperature or precipitation are discussed.

The data in Figure 2 and Figure 3 have been collected from the meteorological station Monchegorsk (ID 22212, located at 67°58.002′, 32°52.998′) located in the forest zone outside the town of Monchegorsk. Its meteo data can be assumed to be representative of the climate conditions in the boreal forest zone of the study area.

2.2. Industrial Impact on Vegetation of the Kola peninsula

Since the 1930s, the Kola peninsula has been the most highly industrialised and polluted region in northern European Russia. The discovery of vast deposits of ores, such as copper-nickel and apatite-nepheline ores, led to the consequent development of non-ferrous metal smelters and mines. Smelter emissions and mining constitute the most significant sources of pollution in this area [33]. These are complemented by the footprint of accompanying open quarries, spoil heaps, extensive tailing ponds and infrastructure. An extremely severe impact had been caused by the copper-nickel Severonickel smelter in Monchegorsk city and the apatite-nepheline industry in Apatity town (see Figure 1). The pollutants of the Severonikel smelter, primarily emissions of sulphur dioxide mixed with heavy metals, are highly toxic [34,35,36]. They spread easily in the atmosphere, lithosphere and hydrosphere according to prevailing winds, topography and soil properties [9,17]. Winds and ground waters transport toxic pollutants over tens of kilometres, aiding their migration between and accumulation in the components of the ecosystem [35,37,38]. Differences in the mass and mobility of heavy metals and sulphur dioxide lead to their accumulation at different distances from the source of pollution [39,40]. The areas in the proximity of the smelter suffer from high concentrations of heavy metals while areas farther away suffer from the lighter sulphur dioxide [36,41].

An excessive accumulation of toxic pollutants in plants leads to a partial or complete suppression of their physiological activity, impacting even relatively tolerant species. Trees develop skeletal branches, needle life span reduces from 6–9 years in healthy plants to 1–3 years in ones severely damaged through toxic pollutants [42], the number of plant species decreases and so does the volume of phytomass, leading to subsequent plant death [5,13,36,43]. These processes occur over years and cause a prolonged impact on the vegetation cover, leading to a mechanical destruction of the natural landscape.

At the peak of Soviet industrialisation in the 1970s–1980s, the Severonickel smelter was emitting on average 230 kt of sulphur dioxide (SO${}_2$) per year (Figure 4). With the collapse of the USSR in the early 1990s, the industrial output in the Kola peninsula dropped dramatically, including that of the smelter. In addition, both the Severonickel smelter in Monchegorsk and the Pechenganickel smelter in Nikel moved to use local ore with a lower sulphur content, resulting in a further reduction in toxic emissions (Figure 4). From 2016, only the total emission values of the Kola division smelters (across Monchegorsk and Nikel–Zapolyarniy) were available and, therefore, not shown in the figure. It is worth pointing out, however, that in 2021, those went down to 16 from 73 kt in 2020 [44].

By 2005, the sulphur dioxide emissions from the smelter dropped by more than 80%, to under 50 kt, and have remained more or less unchanged since then [46,48]. At the same time, the emissions of heavy metal from the smelter dropped by 78%, from ~4.9 kt in 1990 to ~1.1 kt in 2005 [45,46,48]. A biogeochemical analysis of samples of soil, water and vegetation in spruce dwarf-shrub-moss and pine-shrub-lichen forests [48] between 1991 and 2012 showed a similar picture of significantly decreased concentrations of sulphur, nickel and copper. According to the report by the smelter proprietors [46], the Kola division (including the Severonickel and Pechenganickel smelters) further decreased the production of copper by 27% and nickel by 16% between 2011 and 2020 while keeping the production of platinum and palladium unchanged. However, we could not find specific emissions data of the Severonikel smelter for the period between 2017 and 2021.

The forest industry contributed further damage to the vegetation of the area. Spruce and pine forests together cover approximately 70% of the total forest area of the Kola peninsula [49]. In particular, pine represented a significant commercial interest for the timber industry. Commercial logging and, specifically, clear-cutting reduced over years since the late 1990s due to the lack of profitability. Between 1993 and 1999, many earlier logging areas were seeded by spruce [50] and, since the 2000s, increasingly by pine.

With the demise of commercial logging, illegal uncontrolled clear-cutting took place and spread widely in the area. This and the development of tourism led to an additional spread of man-laid fires. Uncontrolled logging and fires resulted in the destruction of both the forest and the shrub-lichen-moss layer. The recovery of such damaged areas occurred predominantly spontaneously and naturally through some preserved shoots of young pine and spruce and the proliferation of young birch.

Given the significant reduction in toxic pollutions since the late 1990s, we expect a substantial natural recovery of the severely damaged vegetation by the plant species tolerant to the northern conditions in the area of influence of the smelters and the growth of the vegetation phytomass from the early 2000s towards 2021. We expect the recovery of the vegetation damaged by pollution to occur through pine rather than spruce due to the malnutrition of the spruce on soils damaged by pollution. Pine, with its taproot system, obtains nutrients from deep mineral soil horizons, whereas shallow spruce roots limit their nutrients to organic soil horizons, which are susceptible to damage. Spruce suffers from a deficit of nutrients in the soil, such as Ca, K, P, Mg and Mn, and the artificial addition of these elements into the soil does not lead to positive and long-term results. Such vegetation restoration experiments were attempted for areas of technogenic barrens [35,51] but did not yield notable improvements in nutrition. Generally, the long-term nutrition deficit in industrially damaged areas should limit coniferous forest recovery or lead to a recovery through pine or deciduous vegetation, typically birch. Similarly, we expect large areas of birch in the areas restoring from uncontrolled logging and man-laid fires in the 2000s and beyond.

3. Data and Methods

Given that Landsat products [52] provide high-resolution images across a series of spectral bands and are freely available for the dates between 1985 and 2021, they became a natural choice of data for this study. While we aimed to process data at a 5-year interval, cloudiness of images and gaps in Landsat product availability in some years forced us to slightly deviate from that frequency.

For performing the classification, we used the state-of-the-art Tensorflow deep learning method [23], which we integrated into the in-house research software that provided image pre-processing, classification, result post-processing, etc. (see Supplementary Materials for links).

3.1. Landsat Product Entities

We collected and processed passive multi-spectral, multi-temporal Landsat TM, ETM+ and OLI products, acquired with the pixel resolution of 30 m. We used Landsat Level-2 Science Products from Landsat Collections 2 that contain atmospherically corrected bottom-of-atmosphere (BoA) reflectance data [52]. The comparative analysis of Landsat Level-2 Science Products and Level-1 Data Products showed improvements of surface reflectance products over all spectral bands for Landsat-8 OLI Level-2 and Landsat-7 ETM+ Level-2 [53], which makes Level-2 Products a good choice for this study.

There are only a couple of periods in this area that demonstrate low cloudiness [32]. First, it is the period from January to April, when the sun elevation is low, ground illumination is poor and vegetation is in the dormant phenological phase. Then, it is the period of July and partially August, when the sun elevation is at its peak, providing the best possible illumination, and when vegetation is at the peak of phenological maturity. For high-quality Level-2 products, the solar zenith angle should be below 76° [53], which is guaranteed in the summer months.

We have collected images for July or early August.

Table 1. Criteria for selecting Landsat images. Here, VNIR stands for visible and near-infrared portion of the electromagnetic spectrum and corresponding Landsat bands.

Parameter	Value	Value for the Training Set
Path	187 and 188	186–190
Row	012 and 013
Sensor bands	VNIR
Cloudiness	<30%	<10%
Sun elevation	High	High
Time of year	July to beginning August
Years	1985–2021	1985–2010

lists essential criteria which we used to select the most suitable images for this study.

Table 2. Downloaded Landsat product entities used for collecting training and testing data for the classifier.

Spacecraft ID	Sensor ID	Path/Row	Date	Entity ID
Landsat 5	TM	188/012	28 July 1986	LT05_L2SP_188012_19860728_20200918_02_T1_SR
Landsat 5	TM	190/012	13 July 1993	LT05_L2SP_190012_19930713_20200914_02_T1_SR
Landsat 7	ETM+	188/012	26 July 2000	LE07_L2SP_188012_20000726_20200918_02_T1_SR
Landsat 7	ETM+	186/012	28 July 2000	LE07_L2SP_186012_20000728_20200918_02_T1_SR
Landsat 5	TM	187/012	9 July 2005	LT05_L2SP_187012_20050709_20200902_02_T1_SR
Landsat 5	TM	190/011	25 July 2009	LT05_L2SP_190011_20090725_20200827_02_T1_SR

and

Table 3. Downloaded Landsat product entities used for land cover classification.

Spacecraft ID	Sensor ID	Path/Row	Date	Entity ID
Landsat 5	TM	188/012	9 July 1985	LT05_L2SP_188012_19850709_20200918_02_T1_SR
Landsat 5	TM	188/013	9 July 1985	LT05_L2SP_188013_19850709_20200918_02_T1_SR
Landsat 5	TM	188/012	23 July 1990	LT05_L2SP_188012_19900723_20200915_02_T1_SR
Landsat 5	TM	188/013	23 July 1990	LT05_L2SP_188013_19900723_20200916_02_T1_SR
Landsat 5	TM	188/012	8 August 1996	LT05_L2SP_188012_19960808_20200911_02_T1_SR
Landsat 5	TM	188/013	8 August 1996	LT05_L2SP_188013_19960808_20200911_02_T1_SR
Landsat 7	ETM+	187/012	17 July 1999	LE07_L2SP_187012_19990717_20200918_02_T1_SR
Landsat 7	ETM+	187/013	17 July 1999	LE07_L2SP_187013_19990717_20200918_02_T1_SR
Landsat 5	TM	187/012	9 July 2005	LT05_L2SP_187012_20050709_20200902_02_T1_SR
Landsat 5	TM	187/013	9 July 2005	LT05_L2SP_187013_20050709_20200902_02_T1_SR
Landsat 5	TM	187/012	10 July 2011	LT05_L2SP_187012_20110710_20200822_02_T1_SR
Landsat 5	TM	187/013	10 July 2011	LT05_L2SP_187013_20110710_20200822_02_T1_SR
Landsat 8	OLI	187/012	10 July 2017	LC08_L2SP_187012_20170710_20200903_02_T1_SR
Landsat 8	OLI	187/013	10 July 2017	LC08_L2SP_187013_20170710_20200903_02_T1_SR
Landsat 8	OLI	187/012	5 July 2021	LC08_L2SP_187012_20210705_20210713_02_T1_SR
Landsat 8	OLI	187/013	5 July 2021	LC08_L2SP_187013_20210705_20210713_02_T1_SR

list Landsat product entities selected for this study, which satisfy these criteria. Those have been downloaded free of charge from Landsat Missions at the USGS and EROS (see Data Availability Statement for links), unpacked into per-entity directories and used as is by the algorithm described below.

Collecting data for training a land cover classifier is a laborious and time-consuming process. Nevertheless, it is common practice to collect training data from the same image that is being analysed. This approach is wasteful and not directly applicable for classifying images for other areas or dates. In this study, we wanted to develop and test an approach that can scale, that can be applied to a broad and remote territory and a long date range. For this, we separated not just the locations but the dates and images used for collecting training data (Table 2) from dates and images used to compute results (Table 3). Generally, this permits reusing the classifier with other images across the European North for a wide range of dates.

3.2. Data Processing

Having acquired the Landsat images and having attributed specific areas to distinct classes used to train the classifier based on field descriptions and general knowledge of the study area, the problem of assessing land cover changes via the classification of Landsat images includes the following steps:

1

Defining the set of input variables for the classifier and output land cover classes;

2

Extracting the training and testing sets of input variables from atmospherically corrected Landsat images selected for training/testing;

3

Training the TensorFlow-based classifier and assessing its accuracy using the testing set;

4

Performing the full-image classification of Landsat images;

5

Filtering and subsampling of resulting classification images to the area of interest;

6

Mapping land cover classification and computing statistics;

7

Computing change maps between pairs of years as well as contingency tables for class changes between years.

The implementation of this algorithm relied heavily on existing commonly used libraries for geo-image processing and classification, namely libgdal (linkable shared library of GDAL [54]), libtensorflow (linkable shared library of TensorFlow [23]) and libgeotiff (linkable shared library of GeoTIFF [55]), all installed via Homebrew [56] on Apple OS X.

The algorithm itself is a collection of code (hereafter research software) that glues together the above libraries for the specific problem and study area (e.g., it makes use of the specific projection of Landsat images for this area). It was implemented as a free and open-source application (links provided in the Supplementary Materials section) developed from scratch as a part of this study, using the Go language and relying only on free and open-source dependencies. It was developed, tested and executed on Apple OS X, using conventional Apple M1 (arm64) hardware. With minor tweaks to the mechanism of sourcing the dependent libraries, the code can be compiled and run on any Linux or Unix operating system. Detailed instructions on compiling it and rerunning the analysis are provided in the source code itself.

Given that the problem in this study is not of an exploratory nature but has a well-defined sequence of operations, a highly efficient statically typed language is better suited than much less efficient dynamically typed languages, such as Python or R, typically used for machine learning. C++ is the industry standard for performance; however, it is difficult to master and comprehend and is highly error prone. This is why we have chosen Go which is almost as efficient as C++, equally compiled and statically typed, yet extremely eloquent and easy to master. It is even simpler than Python in terms of syntax and provides type checking at compile time, leading to fewer errors.

3.2.1. Class Definitions, Input Variables and Training Data

Landsat Level-2 is highly comparable across satellites: characteristic histograms of bands 1–5, and 7 in Landsat-5/7 and the corresponding bands 2–7 in Landsat-8 demonstrated a high degree of comparability as shown in Figure 5. These histograms were collected across a large, fixed representative area (with respect to land cover class variability) across Landsat images spanning a number of years to compensate for any biases in selecting a specific image, area or date. These diagrams are not intended as an accurate representation of any particular image, area or date but provide the means to validate that, on average, Landsat-5, Landsat-7 and Landsat-8 provide comparable response for comparable land cover in this region. This permitted us to use all typical bands and indices: bands 1–5 and 7 in Landsat-5/7 (bands 2–7 in Landsat-8); indices NDVI [57], NDWI [58], NBR [59] and NBR2 [60].

Landsat Level-2 products provide reflectance values roughly in the range [0, 1]. The indices, that are normally expected to centre around 0.0, were offset by +0.5 to bring them roughly into the range of [0, 1]. This was performed to equalise the contributions of input variables when distances between classes are computed by the algorithm and the cost function is applied.

Class definitions were based on field descriptions over multiple years and the expert knowledge of the area. Training and test data were collected via mapping Earth coordinates of specific class annotation to Landsat band values in specific years and specific preselected Landsat images. The initial set of 43 classes (defined through a variability of original field descriptions) was reduced in the final iteration to 11 classes describing different degrees of vegetation damage along with a number of non-vegetation classes.

Training and test data for these 11 classes were collected from four Landsat-5 and Landsat-7 images for years 1986, 1993, 2000, 2005 and 2009, listed in Table 2, mapping site field descriptions into pixel values. The field descriptions were collected in years 1993, 1994, 1999–2009 and 2011. Where possible, the sites were located in large homogeneous areas of forest, wetland, tundra, lakes, etc., but also in homogenous areas of human impact through fires, technogenic barrens, industrial waters, etc. A few classes could not be collected in homogeneous areas and those show much smaller numbers of pixel values. We used the closest in date satellite image to the date of the field site description to collect the corresponding pixel values. The span of years and different satellite systems were used to reduce systemic biases in the training data. Some classes were represented much more frequently in the field data than others. To avoid overfitting the model due to non-uniform distribution of class sizes, each class was reduced to a maximum of 40,000 px values using a uniform random sampler. Of these, the smaller of the following two was used to put aside classifier testing data: 20% of available pixel values, or the maximum of 3000. This resulted in the following datasets for training the model and assessing its accuracy given in

Table 4. Classes used in the Tensorflow classification of Landsat images along with their class value in the resulting TIFF file, class description and legend colour. The last two columns provide the number of pixel values used for training and testing the classifier.

Class Value	Forest Zone	Non-Forest Zone	Legend Colour	Training Pixels	Testing Pixels
0	Missing data	ditto	Transparent
1	Cloud	ditto		32,000	3000
2	Clean water: lake, river or shadow	ditto		32,000	3000
3	Water with sediments: mostly industrial	ditto		7775	1943
4	Non-vegetated: technogenic barren, residential or industrial area	Stone tundra		32,000	3000
5	Burnt area: mostly new	ditto		3632	907
6	Sparse tree/shrub: mostly deciduous	ditto		32,000	3000
7	Wetland	—		4235	1058
8	Coniferous: pine or 40–60% damaged spruce	—		26,574	3000
9	Coniferous: spruce	—		32,000	3000
10	Deciduous: birch, willow	—		32,000	3000
11	—	Tundra shrub, dwarf shrub, lichen		32,000	3000

.

See Supplementary Materials for information on downloading data (including training and testing datasets) and on accessing the source code for software (directory data of the research software).

3.2.2. TensorFlow Classification

In this study, we use the TensorFlow machine learning algorithm for Landsat image classification based on deep learning neural networks [23] and the Keras library [61] that provides a uniform interface for defining input for deep learning neural networks.

We use three Keras layers [61], yielding a categorical classification of 11 classes described above. The model is trained with the Python TensorFlow library and then used in the Go code for the full-image classification. Full-image classification has been performed for all Landsat images listed in Table 3.

The classifier confusion matrix is shown in Figure 6. The accuracy of the model was assessed over the testing dataset and, after fine-tuning the classifier parameters, constituted 63% in the cross-validation over 11 classes (up from 60% with the default parameters). Most classes are accurately reproduced in the training dataset (with an accuracy of >60%). The highest levels of confusion are observed between wetland and deciduous vegetation and between spruce and pine coniferous vegetation. Both misclassifications are expected due to similar spectral properties of wetland and sparse deciduous cover, and coniferous cover, respectively.

The part of the algorithm that defines and trains the TensorFlow model still relies on Python and the tensorflow library in Python, installed via conda [62], which is also available via Homebrew. It is executed only once rather than for every image and, therefore, does not constitute a performance bottleneck.

See Supplementary Materials for information on downloading data (including the serialised classifier and full classification images) and on accessing the source code for software (directory classification implements training the model and performing full image classification).

3.2.3. Trimming and Filtering of Classification Images

To focus the analysis only on the area with the highest quality of data across all years (low cloudiness, no image artefacts, etc.), all images have been trimmed to a quadrilateral defined by the following four coordinates in the sinusoidal projection EPSG:32636: top-left (7668849, 479358), top-right (7620132, 592293), bottom-right (7341672, 464008) and bottom-left (7392573, 358017).

To further reduce noise in classification results, a weighted circular filter was applied to trimmed classification images, similar to the approach originally introduced by Huang et al. [63]. For every pixel, the $5\times 5$ filter block (Figure 7a) was centred at the pixel and the total contributions of all the classes surrounding the pixel were computed using the weights favouring the actual pixel value and the nearest surroundings. The class with the highest weight was selected to represent the resulting pixel value. Typical noise reduction can be seen in the example in Figure 7c, obtained over data in Figure 7b taken from one of the classification images.

See Supplementary Materials for information on downloading data (including trimmed and filtered results) and on accessing the source code for software (directories trim and filter implement image trimming and filtering, correspondingly).

3.2.4. Computing Change Maps

To assess changes in land cover since 1985, in years 2017 and 2021, change maps are computed by comparing pixel values between the classification maps for 2017 and 2021 and the classification map for 1985. Wherever the class value has not changed between the two years, the corresponding class of the pixel was preserved. Wherever the class value changed, a special value marking that there was a difference was assigned. The resulting maps show the scale of all changes on the background of stable land cover classes, in particular, water, stone tundra, etc.

4. Results and Discussions

To assess the qualitative and quantitative changes in the boreal vegetation of the central part of the Kola peninsula over 35 years, we computed eight land cover classification maps for the years 1985, 1990, 1996, 1999, 2005, 2011, 2017 and 2021 (Figure 8).

These maps were obtained by processing the total of 16 multi-temporal Landsat images, applying the algorithm described above and combining the classification results of two images (rows 012 and 013) for each year to cover the study area in full.

A direct assessment of the quality of the classification results over all the years and such a large territory is prohibitively expensive. We have to rely on reference maps and our local knowledge. We have used a variety of environmental maps (specifically vegetation maps, landscape maps, old-growth forests maps, the map of metal content in bioindicator plants and others) compiled in different years to validate our classification maps. We also used our field data from the area, which was not used for training the classifier, as well as our expert knowledge about the area, its development and natural evolution for the general quality assessment. While we cannot guarantee the accuracy of each individual pixel, the general land cover distribution correlates well with the reference maps and our local knowledge.

4.1. Land Cover Change Detection

We start with the general land cover change detection, looking at which areas changed between 1985 and 2017/2021 and how they changed.

Figure 9 shows two land cover change maps computed between the years 1985 and 2017 (Figure 9a) and the years 1985 and 2021 (Figure 9b). The computation employed the maps in Figure 8a for 1985, Figure 8g for 2017 and Figure 8h for 2021.

We have computed maps for two different target years to obtain a more general picture of the vegetation transformation. Looking at the changes up to these two years next to each other, we can make a better overall picture of the land cover changes over the 30-year period without biasing our interpretation too much to one outcome that could be affected by the climatic conditions of the specific year.

The year 2021 was visibly warmer and generally wetter than average, and we see a strong change in the land cover compared to 1985. This is because various species exhibit a different onset, duration and end of their vegetation phenological phases, which also vary yearly depending on climatic conditions [64,65]. At the same time, substantial alterations took place in the industry over five years from 2017 to 2021 and further reduced the industrial output that could result in a visibly higher fraction of changed land cover due to the vegetation recovery in 2021 compared to 2017.

The extensive spatial changes in the vegetation state since 1985 are observed for both 2017 and 2021 in Figure 9. These changes affect mostly the forest zone, while the vegetation in the latitudinal tundra zone has changed only slightly. First, the changes due to climatic factors are less visible in the tundra vegetation cover. Then, no notable industrial activity in the tundra zone affects the lichen cover. The forest vegetation is susceptible to both climatic factors and industrial activity.

In the southern part of the study area, we know of extensive forest logging areas in the 2010s. These areas are shown as changed in 2017. Here, spruce was partially planted and partially left for natural recovery after forest logging activities, and it recovered substantially towards the 2020s, which is also partially observed in the change maps (where the changed state of 2017 is replaced by the original spruce class, matching 1985 in some areas). The distribution of coniferous forests (both pine and spruce) in our maps in Figure 8, and the partial preservation of spruce throughout the years (observed in Figure 9), correlates well with the distribution of old-growth forests of the Murmansk region mapped in 1999 [66].

In the central and northern parts of the study area, we observe the disappearance of old-growth forests between 2017 and 2021. They show up partially unchanged in Figure 9a for 2017 but are no longer seen in Figure 9b for 2021. This corresponds to the recent forest logging in those areas.

To better understand how these areas have changed, that is, which classes converted to which on average, we then computed the confusion matrices for the land cover class transformations between 1985 and 2017 (Figure 10a) and 2021 (Figure 10b).

4.2. Recovery of Damaged Vegetation

Vegetation reflectance depends on the chemical characteristics of the plants (e.g., plant age and ability for water intake) and their structural properties (e.g., leaf or needle thickness and crown density). For example, smaller crowns of coniferous trees lead to lower canopy reflectance [67]. These parameters also determine the sensitivity and tolerance of different plant species to environmental changes [68,69,70]. For example, lichen is more sensitive to air quality than dwarf shrub (Empetrum-Vaccinium) vegetation, while the grass is relatively stable towards air pollutions [71]. This way, the loss of more sensitive plant species and changes in the canopy density due to natural and human factors make it challenging to separate the damaged and undamaged vegetation based purely on spectral reflectance properties. One of the outcomes is the domination of the class ‘Coniferous: pine or 40–60% damaged spruce‘ that we observe in Figure 9. This class captures a wide range of partially damaged vegetation or vegetation restoring through pine. We see how the spruce forest, burnt areas, deciduous with birch, etc., are converting to damaged/pine over time, highlighting that the areas were damaged and are recovering.

Given the concentration of the human impact in the boreal zone (the white quadrilateral in Figure 8), this zone had seen the highest impact of the industry and toxic pollutions. Consequently, the disturbed areas here (originally spruce forests) had been recovering through deciduous plant species and pine. In both 2017 and 2021, we consistently observe the recovery of the burnt areas through deciduous vegetation and pine. The burnt areas here are of human origin, either in places of dry forest, damaged by industrial emissions, or in areas of the recreation close or near the lakes or rivers.

Since the 1980s, significant reductions in industrial air emissions have led to a decrease in the concentrations of sulfur and heavy metals in coniferous needles but not in soil [72]. As expected, this has led to the recovery of damaged spruce forests through pine and deciduous vegetation. Further similar effects are expected in 2021 and beyond due to the Severonickel smelter closing the copper processing unit for modernisation in the spring of 2020. It is expected that sulphur dioxide emissions will drop down to 7–15 times lower than the level of sulphur dioxide in the 1980s.

The quantitative confirmation of the damaged spruce recovery over time through deciduous vegetation and pine can be clearly seen in Figure 11b. The Mann–Kendall test [26] shows a relatively steep downward trend for spruce ($\tau =-0.5$) of moderate confidence ($p=0.11$) compensated by moderate upward trends for deciduous vegetation ($\tau =+0.29$) and pine ($p=0.38$), each of a low confidence (p-values of 0.38 and 0.54, respectively). The low confidence is due to a small number of observations and visible outliers in the middle of the series.

The average yearly air temperatures increased by more than 2 °C between 1985 and 2021 (Figure 3). This is a general worldwide trend across the Arctics and sub-Arctic, which led to the growth of the primary productivity and expansion of shrubs and trees [24,73]. The summer of 2021 was the warmest through all the studied years, and its spring was one of the warmest in the set (Figure 12b). At the same time, the winter of 2021 was one of the coldest, second only to 1985, with generally low precipitation the year round (Figure 12b). Due to these factors, the deciduous vegetation was flourishing in 2021, much more substantially than the coniferous one, and we clearly see both in the classification maps in Figure 8h and in the outlier value of the trend in Figure 11a. An additional potential contribution to the volume of the deciduous vegetation in 2021 could be that the copper department at the Severonikel smelter closed for modernisation in the spring of 2020 and the reduction in emissions became noticeable a year later.

The field observations that we collected in this area between 2002 and 2019 also confirm the steady recovery of the damaged vegetation, even in the cold boreal climate. This is helped by the upward trends in the annual year air temperature shifting more than 2 °C between 1985 and 2021 (Figure 3a) and the relatively warm winters. The weather reports confirm the general trend to warmer winters (Figure 12a).

We observe a slight increase in the overall phytomass over the years, proxied by the fraction of classes Deciduous, Coniferous: pine and Coniferous: spruce. One can see this in a relatively steady increase in their total volume in Figure 11a. However, the increase is not statistically significant to qualify as a trend. Similarly to our earlier research [75], we observe vegetation recovery not via the original dominant spruce species but via the pine and deciduous vegetation. The recovery of the damaged vegetation is still a widely discussed and controversial topic [76,77].

4.3. Metodology Considerations

To improve the quality of the classification results under the conditions when different types of vegetation could yield similar reflectance patterns, we used areas with homogeneous land cover of the size no less than $90\times 90$ m${}^2$ (or $3\times 3$ pixels in Landsat images with the 30 m resolution). Further improvements are achieved through indices computed through a non-linear combination of spectral bands, such as NDVI and NBR, used in this study. The digital elevation models and thematic GIS layers can further improve the separation of the land cover classes but are difficult to incorporate alongside the spectral reflectance data.

Cloudiness is another significant limiting factor for selecting satellite imagery in the northern regions. For accurate land cover change detection, we need images with low cloudiness in the study area (<10%); otherwise, the error in the pixel distribution by classes is significant and it is difficult to detect the year-on-year change trends clearly. Beyond the study area, the image may exhibit the cloudiness coverage of up to 30%. At the same time, the number of cloudy days is very high in the North throughout the year; generally, every second day of the month is cloud, which leads to a shortage of appropriate images in global image archives. To collect sufficient datasets for year-on-year land cover monitoring, we had to expand the time interval by a couple of weeks around the end of July, the most probable period of the vegetation maturity stage, to the whole month of July and beyond.

For a reliable year-on-year land cover change detection, we focus on the maturity stage when the vegetation is at its growing-season peak (where four phenological stages are normally distinguished: dormancy, green-up, maturity and senescence [78]). An analysis of the phenological data shows that the maturity stage in boreal forests is reached between the middle of July and the beginning of August, which is also the warmest time in this area. This is a very short period, which complicates the selection of the satellite imagery for the aims of vegetation change detection.

Due to the short vegetation period and the high cloudiness throughout the year, it is impossible now to relate the vegetation state with the inter-seasonal fluctuations of climatic factors. It is recommended to use a continuous set of remote imagery for such studies, ideally for a defined date at the end of July every year.

5. Conclusions

Based on the Landsat multi-temporal data over a 35-year period, we observe the increase in deciduous vegetation (birch, willow and aspen) of different live forms (trees and shrubs) and the increase in pine in the central parts of the Kola peninsula in the 2010s replacing the original spruce damaged through the toxic pollution of smelters, mining, fires and other industrial activities in the years before. Deciduous plant species have spread in places of mechanical and chemical damage.

On the local level, we observe the greening of a technogenic barren and significantly damaged vegetation by more tolerant plant species. On the regional level, we do not see strong dependencies between the seasonal fluctuations of the climatic factors and vegetation cover. There is a lack of continuous data for such an analysis.

The use of Landsat images in large-scale studies spanning multiple years, while generally straightforward, poses some limitations due to the changes made in the sensor band ranges. Some bands differ insignificantly, but we had to exclude band 5 (Landsat-5/7)/band 6 (Landsat-8) due to significant differences in the band values that would not map to the same classes.

On the technical side, TensorFlow has proven to be an extremely versatile, easy-to-use and scalable technique. It yielded a classification accuracy above 60% with its default set of parameters and could be fine-tuned up to 63%. The classifier could be trained on conventional hardware in the order of minutes to classify a complete Landsat image of approximately 8000 × 8000 px. The combination of Go, TensorFlow and GDAL has also proven to be extremely efficient and scalable. In this study, we managed to process dozens of Landsat images in a matter of minutes on conventional hardware. In our earlier studies, we used the combination of R, SVM and GDAL, which required two orders of magnitude more time to proces a similar amount of information.

The use of the Go as a glue language to bring the end-to-end analysis together has proven to be efficient due to its simplicity, compile time validation and speed. Other languages that scientists are proficient in, e.g., Python, can equally be used instead. However, interpreted languages, such as Python and R (in contrast to compiled languages, such as Go or C++), generally present a higher barrier for users to reach comparable levels of execution speed and correctness. It was straightforward to bring in the GDAL, TensorFlow and (for experimentation) SVM libraries to our Go code, which might be less straightforward with other languages, such as the R Language for Statistical Computing. In our earlier studies, we used R to bring together all the analysis, which required much more resources and was much slower overall.