Pollution and Climatic Influence on Trees in the Siberian Arctic Wetlands

Abstract

Siberian Arctic wetlands located within the planetary “warming hotspot” experience pronounced climate-driven vegetation cover changes. Together with warming, wetlands, which are located within the influence of Norilsk copper and nickel industry (69.35° N, 88.12° E), have been strongly influenced by industrial pollutions (sulfur dioxide mostly) since the 1940s. In addition, petroleum products release occurred in 2020 that potentially influenced vegetation vigor. We studied the combined effect of climate warming and pollution on the larch (Larix sibirica Ledeb.) and shrubs’ (Salix spp. and alder, Duschekia fruticosa) growth. Using satellite data (MODIS and Sentinel) processing, we mapped wetlands within the study area. We used on-ground survey, and applied dendrochronology, climate variables, and emissions rate analysis. We sampled woods (kerns) and, based on the tree ring analysis, generated trees and shrubs growth index (GI) chronologies. We analyzed the influence of the SO₂ emissions and eco-climate variables (air temperatures, precipitation, soil moisture, and drought index SPEI) on the larch and shrubs GI. We mapped GPP and NPP (gross and net primary productivity) and vegetation index NDVI and temporal trends of these indexes based on the MODIS-derived products. We found that chronic SO₂ influence led to larch trees GI decrease that was followed by tree mortality, which was observed until the end of 1990. Since the beginning of the 2000s, the GI of larch and shrubs has increased, which is correlated with elevated air and soil temperature and growth season prolongation, whereas excessive soil moisture negatively influenced GI. Together with that, increasing trends of vegetation indexes (GPP, NPP, and NDVI) were observed on the part of wetland within the zone of former trees’ heavy damage and mortality. The trends began mostly in 2003–2005 and were caused by emissions volume decrease and warming, together with resistant species’ (willows, graminoids, bushes, and birch) growth and invasion. We suggested that increasing productivity trends might partly be attributed to nitrogen fertilization caused by NOx emissions. Finally, we found that diesel fuel spill which happened in 2020 caused no influence on the larch, whereas some aquatic species (mosses mostly) were damaged.

1. Introduction

Northern marchlands are under a stronger influence of elevated temperatures than similar lands at lower latitudes, since the warming rate in the Arctic is more than two times higher [1]. Together, increased warming and fire activity in the permafrost area are initiating changes in the vegetation cover and species composition [2,3]. Within the Siberian Arctic, a combined warming and anthropogenic impact on the vegetation cover located within the polluted zone was caused by the Norilsk industry region. Chronical emissions (sulfur dioxide mostly) have influenced the wetland vegetation, including sparse and open forests, since the 1940s. This is a unique combination of the technogenic and climatic influence on the wetland vegetation within the completely circumpolar zone and beyond.

Norilsk city, one of the five most northward world cities, is located 300 km beyond the Polar Circle within the permafrost area of the Taymyr peninsula (Figure 1). Norilsk industry began in 1935 with the construction of nickel and copper plants. Firstly, concentrates of copper and nickel were obtained in 1939. Refined nickel, copper, cobalt, and platinoids production took place in 1942–1944. Since then, metallurgy and mining have rapidly continued to expand until the end of the 1980s [4]. Nadezhda, a huge third plant, started to work in 1979–1981. In 1985, Nadezhda production reached its planned volumes. Emissions from Nadezhda smelter are transferred over long distances due to high (250 m) smelter located on the plateau with elevation c. 200 m. The heights of the copper and nickel smokestacks are about 180 m and they are located in the lowlands. The Norilsk industry’s maximal production reached its maxima in 1989. However, in 1991, production decreased twice with minimum values in 1993. Nowadays, Norilsk industry provides, alongside nickel and copper, gold, silver, platinoids (palladium, osmium, iridium, rhodium, and ruthenium), selenium, tellurium, and sulfur [4].

An adverse strong by-effect of the Norilsk industry is air pollution by sulfur dioxide, nitrogen oxides, and heavy metals (Cu, Ni, Co). That caused dramatic decline and mortality of ground-cover vegetation (moss and lichen communities) and trees (primary larch, Larix sibirica Ledeb.) [5,6,7,8]. Thus, in 1990, larch decline and dieback were observed at the distances of 200 km and 80–100 km in the southwest direction from the Norilsk, respectively [9,10]. Alongside air pollutions, the recent release of c. 20,000 tons of diesel oil occurred in May of 2020 [11]. The spill impact on the water terrestrial biota was a subject of recent studies [12].

We hypothesize that current climate change might mitigate adverse pollution influence on the trees and shrubs growth.

The purpose of this study is an analysis of trees growth response to the combined influence of air pollution (SO₂ emissions) and climate change on the trees growing in arctic wetlands. In addition, we studied the influence of the diesel oil spill on the tree growth.

We seek the answers to the following questions:

-

What was the impact of SO₂ emissions on the larch (Larix sibirica Ledeb.) trees and shrubs’ (Salix spp. and alder, Duschekia fruticosa) growth?

-

How has climate warming influenced larch trees and shrubs’ growth?

-

What are the temporal trends of the gross (GPP) and net (NPP) primary production and vegetation index NDVI (Normalized Difference Vegetation Index) within the zone of pollution influence?

2. Materials and Methods

2.1. Study Area

The study domain is located within the northern wetlands, surrounded by offshoots of the Putorana plateau (Figure 1). This is the zone of continuous permafrost with multiple lakes. Sparse and open forests formed by larch, spruce (Picea obovata Ledeb.), and birch (Betula spp.). Bushes are presented mostly by willows (Salix spp.), alder (Duschekia fruticosa Pouzar), and Betula spp. Stony tundra is occupied by lichens, moss bushes, swamp ledum, and alder communities. The climate is sharply continental, with below zero (−9 ÷ −10 °C) annual air temperature (recorded minimum was −57 °C).

2.2. Field Studies

Over 2021–2022, we made a field survey within the area with heavily damaged forests (“impact zone”) and background (control) territories. Twenty-five temporary test plots (TP; approximately 0.5 ha each one) were established around the Pyasino Lake (13 and 12 TP within impact and control zones, respectively; Figure 2). Tree and stand characteristics (stem diameter, tree height, stand density, regeneration number, ground cover, and soil type) were determined. For dendrochronological analysis, we randomly selected trees within each TP. We extracted individual tree cores (small, pencil-sized pieces of wood oriented from bark to stem center perpendicular to the tree’s stem) with an increment borer from about 20 trees to represent each TP. Cores were taken at DBH height (1.3 m above ground-line) or lower (to root collar) from the stem. In the analysis we used 228 larch cores (195 survived and 33 dead trees), 19 alder cores, 59 cores of willow species from the impact zone, and 169 larch cores from the control zone.

2.3. Dendrochronological Analysis

Tree cores were mounted on wooden backing. Cores and disks were finely sanded and treated with contrast powder to enhance visualization of ring boundaries. The measurements were carried out on a LINTAB-6 platform with an accuracy of 0.01 mm. The cross-dating quality was assessed by COFECHA software (https://www.ldeo.columbia.edu/tree-ring-laboratory/resources/software; [13] (accessed on 1 December 2022)). From the measured series, we constructed a growth index (GI) chronology for each TP from raw tree-ring width in mm to unitless GI values by a detrending method with program ARSTAN ([14]; https://www.ldeo.columbia.edu/tree-ring-laboratory/resources/software (accessed on 1 December 2022)). We removed low frequency trends by fitting tree-ring series with a negative exponential curve or line with a negative slope. For statistical analysis, the GI chronology was normalized by conversion to Z-scores with a mean of zero and standard deviation of 1.0.

2.4. Eco-Climate Variables

Tree GI was analyzed against the pollution (SO₂) emissions, monthly air and soil temperatures, precipitation, SPEI drought index (standardized precipitation-evapotranspiration index), soil moisture (volumetric soil water layer depth until to 7 cm), and EWTA (equivalent water thickness anomalies, or terrestrial water storage). Air temperature, precipitation, and soil moisture were obtained from ERA-5 database (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land-monthly-means?tab=overview, accessed on 15 October 2022; spatial resolution 0.1° × 0.1° [15]).

Drought index SPEI is the difference between precipitation and potential evapotranspiration [16]. Decreased SPEI indicated increased drought severity. We calculated SPEI based on ERA5-Land data by R-script (https://github.com/sbegueria/SPEI (accessed on 1 December 2022)). SO₂ emissions were obtained from scientific literature and official reports [5,17,18].

The EWTA is the product of gravimetric measurements made by the GRACE satellite with on-ground resolution 1° × 1° [19]. This is characteristic of the moisture regime and represents total terrestrial water storage. EWTA data were obtained from NASA-JPL database (https://podaac-opendap.jpl.nasa.gov/opendap/hyrax/allData/tellus/L3, accessed on 15 October 2022).

GPP data (C, kg/ha) were obtained from the Terra/MODIS product MOD17A2H version 6.0 (https://lpdaac.usgs.gov/products/mod17a2hv006/, accessed on 15 October 2022). Data temporal and spatial resolutions were 8 days and 500 m, correspondingly [20]. Mosaics of GPP (as well as NPP and NDVI) composites were compiled based on median values from 8-day composites over summer periods. NPP annual values were obtained from the Terra/MODIS product MOD17A3HGF version 6.0 (https://lpdaac.usgs.gov/products/mod17a3hgfv006/, accessed on 15 October 2022). The database contained annual NPP values with a spatial resolution of 500 m [21]. NDVI were calculated based on data from the Terra/MODIS products MOD09Q1 и MOD09A1 version 6.1 (https://lpdaac.usgs.gov/products/mod09q1v006/ [22], accessed on 15 October 2022 and https://lpdaac.usgs.gov/products/mod09a1v006 [23], accessed on 15 October 2022; spatial resolution 250 m).

Wind rise was generated based on wind direction hourly data that were obtained from the nearest of the Norilsk stations (station 23077 and 23078; coordinates 69.33° N/87.95° E and 69.33° N 88.25° E, correspondingly).

2.5. Method of Groundcover Mapping

The created land-cover classification was based on the analysis of MODIS summer scenes obtained in 2022. Seven MODIS spectral bands were used. Bands 3–7 were resampled from 453 m spatial resolution to 231 m using the non-linear downscaling algorithm [24,25]. We used the Support Vector Machine classification method [26]. Training samples were generated based on the ground data, high-resolution scenes and aero-photography, and topographic maps. The accuracy of the classification map was estimated using the error matrix and kappa statistics [27]. The overall accuracy of the classification is 92% (kappa = 0.88). Further land-cover classes along the river valleys were identified based on the GMTED digital elevation model (https://www.usgs.gov/core-science-systems/eros/coastal-changes-and-impacts/gmted2010, accessed on 15 October 2022) and the dataset of water bodies Global Water Surface Extent [28]. The spatial resolution of the created map is ~250 m, corresponding to medium-scale maps of M 1:500,000 [29].

2.6. Statistical Analysis

Four sets of statistical analysis were applied to the GI and environmental variables. Pearson correlation analysis identified significant relationships between the GI and individual climate variables. Spearman’s correlation coefficient was used to identify the relationship between GI and sulfur dioxide values. The piecewise (also known as segmented or broken-stick) regression analysis was used to detect breakpoints in the time-series of larch GI, GPP, and NPP [30]. We used StatSoft Statistica (http://statsoft.ru (accessed on 1 December 2022)) and IBM SPSS Statistics Base V27 (https://www.ibm.com/analytics/spss-statistics-software (accessed on 1 December 2022)) software for statistical analysis.

We used the Theil–Sen regression to determine tendencies (trends) in dynamics of NDVI, GPP, NPP, and climate. The Theil–Sen estimator is a non-parametric method for fitting a line to sample points in the plane by choosing the median of the slopes of all lines through pairs of points and it is based on the Kendall tau correlation coefficient [31,32]. This estimator is less sensitive to outliers, and it is more accurate than the simple linear regression [33]. We used the Theil–Sen estimator realized in the Python library pymannkendall 1.4.2 (https://pypi.org/project/pymannkendall (accessed on 1 December 2022)). The source code of the pymannkendall was adopted to Python 2.7 and used in the ESRI ArcGIS software to apply it to analyze multiband raster datasets represented by time-series of NDVI, GPP, NPP, and climate variables for 2000–2021. As a result, raster datasets of regression slopes and p-values were determined. Additionally, for NDVI, GPP, and NPP, a floating starting year (2001–2017) of the most significant trends were calculated with the last year, fixed, and equal to 2021.

3. Results

3.1. Tree and Test Plot Characteristics

Tree species are represented by Larix sibirica, Picea obovata, and Betula sp. Larch the dominating tree species in the region. Larch tree characteristics varied, including mean tree height (6–10 m), DBH (diameter at breast height, i.e., 1.3 m) (13–16 cm), and age (15–435 years at DBH height; median age 95 years). The shrub layer was composed of alder, willow, and birch (Duschekia fruticose, Salix spp. and Betula spp.). The low shrubs and herbaceous ground cover layer contained Ledum palustre L., Vaccinium vitis-idaea L., Vaccinium uliginosum L., Carex spp., and Equisetum arvense L.

3.2. Climate Variables

Within the study area, northwest winds prevailed during the summer and southeast ones prevailed in winter (Figure 3a). Summer and winter air temperatures have been increasing since 1990 (R² = 0.13–0.34; Figure 3b). Soil summer moisture has been decreasing since 2000 (R² = 0.63; Figure 3c). A strong decreasing trend experienced minimal values soil moisture in July, which indicated an acute drought frequency increase (R² = 0.81; Figure 3c). Terrestrial moisture content (GRACE-derived June values) has been decreasing since 2000 (R² = 0.20; Figure 3d). Meanwhile, precipitation values have stagnated.

Maps of eco-climate variables within the study area show their spatial distribution within the study area (Figure 4). The annual temperature varied from c. −3 °C (in wetlands) to −15 °C (in mountains). Positive temperature trends were observed on the major part of the lowlands (Figure 4a,b). Annual precipitation varied from c. 600 mm/year (wetlands) to 1000 mm/year (mountains) (Figure 4c). Maximal values of atmospheric drought (indicated by drought index SPEI) were observed in the area southeast of Norilsk (Figure 4d). The highest soil temperatures are located around the Pyasino Lake and southeast of Norilsk (Figure 4e). Minimal soil moisture content was observed southward of Pyasino Lake (Figure 4f).

3.3. Pollution Sources and Emissions Transfer

The locations of pollution sources (Nadezhda, copper, and nickel plants) are presented in Figure 5. Pollution transfer is controlled by wind direction and speed and terrain topography. Although the main directions of pollution transfer (according to wind rose, Figure 3a) are southeast and northwest (Figure 5b), pollution might be transferred to the other directions (Figure 5a,c) or stay mostly within smelters location in the case of anticyclone-caused still weather (Figure 5d). Under certain metrological conditions, part of the sulfur dioxide emissions extends to the middle and northern parts of the territories adjacent to Lake Pyasino. However, no significant effect of SO₂ emissions on the radial increment of larch trees growing in this area was found.

3.4. Trees’ Growth Index Responses to Pollution Impact and Climate Variables

The dynamics of larch growth index (GI) in the impact and control zones until 1950 was highly synchronous (correlation coefficient r = 0.9; Figure 6). In the second half of the 20th century, that synchronicity decreased (r = 0.79) due to changes mainly in the low-frequency component. Since 2000, a GI increase has been recorded for trees in the impact zone mostly, and during the last decade (compared to the same pre-warming period of 1990–1999) larch GI increased by 1.6 times; for comparison, within the control zone, GI increased 1.2 times (Figure 6b).

The dynamics of radial increment (a decrease in GI until the second half of the 1990s with a subsequent increase at the beginning of the 21st century) and the dynamics of SO₂ emissions are in antiphase (an increase and subsequent decrease in emissions; Figure 7). Although the local minima of the radial increment during this period correspond to similar minima recorded for other wood species and are likely to be of a climatic origin, the increase in the radial increment in the impact zone (low frequency component; Figure 6b,c) coincides in time with a decrease in SO₂ emissions. According to dendrochronological analysis, the onset of tree mortality dates back to the late 1950s. The maximum mortality of larch trees was observed in 1980–1990 against the background of an increase in sulfur dioxide emissions (Figure 8).

Comparative analysis of the GI and sulfur dioxide emissions dynamics for the period 1985–2020 (i.e., for period of reliable SO₂ emissions data) indicates significant dependence of GI on emissions (r = −0.56 for larch and r = −0.41 for alder, p < 0.05; Figure 9). No significant correlations between GI and sulfur dioxide emissions were found for the larch trees in the control zone and for the willow in the impact zone.

Correlation analysis revealed significant stimulating effects of the air temperature increase in July and inhibition effects of precipitation (July–August) and soil moisture on trees’ GI (Figure 10). There are no significant differences in the reaction of trees’ GI on climatic changes between the impact and control zones. It should be noted that the observed prolongation of the growth period, due to increased temperature in June and September, stimulates radial growth increase. Meanwhile, there is significant negative correlation between GI and the spring warmth (temperature increases in May; Figure 10a).

The multiple regression analysis indicates the volume of SO₂ emissions and air temperature as the main factors influencing the radial increment of larch trees in the impact zone. At the same time, the air temperature in June–July has a stimulating effect on radial increment, and the May temperature has an inhibitory effect (Equation (1)).

GI = 0.65 × T_JJ − 0.43 × SO₂ − 0.35 × T_MAY − 0.04,

(1)

where GI—growth index, T_JJ—June–July temperature, T_MAY—May temperature, SO₂—values of SO₂ emissions. Explained dispersion: R² = 0.7. Period: 1985–2020.

In the control, the radial increment does not depend on sulfur dioxide emissions and is determined by air temperature and soil moisture (“cold soils effect”; Equation (2)):

GI = 0.04 + 0.47 × T_JJ − 0.38 × P_JUNE − 0.31 × SM_MAY,

(2)

where GI—growth index, T_JJ—June–July temperature, P_JUNE—June precipitation, SM_MAY—soil moisture. Explained dispersion: R² = 0.55. Period: 1985–2020.

Shrub species (willow and alder) also have a positive response to climate warming (Figure 11). In addition to significant positive correlations with July temperature, the growing season is prolongation indicated by positive correlations with September and June temperatures. As in the case of larch, negative correlations with May temperatures are observed. Similar to larch, an inhibitory effect of excessive soil moisture on radial increment is also observed for bushes (Figure 11d).

3.5. GPP, NPP, and NDVI Dynamics

GPP, NPP, and NDVI maps indicated the minimal values of GPP, NPP, and NDVI in the areas of the main air pollution impact (Figure 12a–c). Meanwhile, in the areas where increasing trends of vegetation productivity were observed (although in the most affected zone (westward from Norilsk), trends were mostly zero (Figure 12d–f). The start of the positive trends dated mostly to 2003–2005 (Figure 12d–f).

4. Discussion

Norilsk industry emissions caused the biggest large-scale decline and mortality of forests beyond the Arctic circle. Tree growth in the studied northern wetlands was influenced by both air emissions (SO₂ mostly) and the warming climate. The beginning of industrial emissions referred to the end of 1930, whereas significant “emissions-growth index” correlations has been observed since 1950. The latter coincided with the first reported larch forests’ decline and mortality [9]. The growth decrease was continuous until the end of 1990, when the tree mortality rate reached its maxima. In the following period (since 2000), larch growth experienced an increase, which was attributed to the combined effect of pollution load decrease and air and soil warming. SO₂ emissions (the main pollutant) decreased from c. 2.5 M t/year in 1980 to c. 1.6 M t/year at present due to technology improvements, together with nickel plant closure (2016). Some contribution to the growth positive trends may be related to the CO₂ fertilization effect. However, that effect is low in cold habitats [34]. In the “background zone”, the larch growth increase correlated with air temperature only, although certain wind directions emissions partly covered that area, which indicated that emissions load was within the larch compensatory range.

Since the end of 1990, elevated air temperatures in the beginning (June) and in the end (August–September) stimulated larch growth, whereas July temperatures’ influence on the growth decreased. For comparison, during the previous period (1950–2000), larch GI was mostly dependent on the July temperatures (R² = 0.54), whereas June and August–September influences were insignificant. A similar effect of low growth response to July temperature was described for the Siberian pine growth within mountain forest-tundra ecotone [35]. Thus, with warming larch growth limitation by July, temperatures shifted to temperatures in the beginning and the end of growth period. Meanwhile, early (May) warms led to a GI decrease. That effect was caused by “warming provocation” of the phloem and cambium activity, while soil water was not still available with consequent needles and twigs desiccation, as well as the damage by late frosts. In extreme cases, larches may flash needles prior to snowmelt.

Opposite to temperature, soil moisture increases negatively influenced larch growth due to the “cold soils effect”, i.e., low soil temperatures slow down root activity with consequent tree growth decreases (Figure 10d and Figure 11d). This effect is most pronounced in waterlogged habitats. The other negative factor is a limited root oxygen supply within poorly drained soils. However, growth limitation by water availability has also been described in permafrost zones [36,37,38]. Larch may experience a moisture deficit at the beginning of the growth period when elevated air temperatures are sufficient for photosynthesis, while soils are still mostly frozen [37].

For bushes (willow and alder), we observed similar larch growth and temperature shifts. Notably, that GI of willow more strongly responded to the June and September temperatures in comparison with larch and alder. That effect may be caused by non-leaf photosynthesis, i.e., photosynthesis occurred in the bark tissues. For example, for aspen, which also belongs to the Salicaceae family, non-leaf (bark) photosynthesis produced 15–20% of the total photosynthetic production [39]. In a harsh environment, willow species may grow even in the leafless form. Meanwhile, the significance of non-leaf photosynthesis for shrubs and trees surviving in the harsh habitat needs further studies. Notably, the growth index of Salix species did not significantly correlate with SO₂ emissions, even within the zone of strong emissions impact that indicates a high resistance of these species to SO₂ influence. Thus, willow species may be used for reclamation of the polluted territories.

Alongside tree growth increase, increasing trends of vegetation indexes (GPP, NPP, and NDVI) were observed in part of the wetland. Interestingly, the high increase was located mainly within the area of larch trees mortality that occurred during 1960–1990. Vegetation productivity increase was located mostly within areas with lower soil water content, whereas in wetlands with high moisture content (the majority of the background territories), GPP and NPP trends were mostly insignificant. The beginning of those trends is dated mostly by 2003–2005 and is associated with emissions decrease, warming, and invasion and growth of resistant species (birch, willow and other bushes, and graminoids). As noted by [40], the main direction of the phytocenosis successions is the substitution of the forests and moss–lichen associations to the grass and bushes communities. An additional factor of polluted area greening is warming-driven larch growth increase. Although [8] reported zero larch regeneration within the area of total larch mortality, we observed a minor amount of larch seedling in the area of former tree mortality. Certainly, regeneration establishment and, more generally, species composition is a dynamic process that needs further study.

Alongside the above-discussed factors, pollution-origin nitrogen fertilization may also promote positive trends in vegetation productivity. Annual Norilsk NOx emissions varied within 9–25 thousand tons/year. Studies of pollution-caused nitrogen fertilization indicated an increase of tree growth (e.g., [41,42]).

Finally, we found that dispel fuel spill’s influence (which occurred in 2020) on the tree species was insignificant. Meanwhile, all moss species and some vascular ones are highly intolerant to petroleum pollution, whereas aquatic and riverside grasses are petroleum-resistant [12].

5. Conclusions

Sulfur dioxide emissions, air temperature, soil moisture, and temperature are the main constraints of the larch growth in the studied northern wetlands. SO₂ influence led to larch growth decrease from 1950 until the end of 1990 when tree mortality reached its maxima. Since 2000, larch growth has been increasing due to combined effect of air pollution decrease and climate warming. With warming, larch and bushes growth limitation by July temperatures switch to the temperatures in the beginning of the growth period. Meanwhile, early spring warmth led to GI decrease due to “metabolism provocation” while soil water was not available. In some years, even “larch trees greening” phenomenon occurred while the background was still snow-covered. Consequently, this caused twigs and needles desiccation, together with damage by the late frosts.

During recent decades, increasing trends of GPP, NPP, and NDVI were also observed on parts of the wetland. Notably, high productivity trends were located partly in the zone of former trees’ mortality. Positive trends were associated with emissions volume decrease, climate warming, invasion of resistant species (willows, graminoids, bushes, and birch), and increased larch tree growth. In addition, soil fertilization due to NOx emissions may potentially stimulate vegetation productivity increase. Finally, we found that the fuel spill that happened in 2020 had no influence on the larch, whereas some aquatic species (mosses mostly) were damaged.