Permafrost Degradation and Vegetation Growth beyond the Polar Circle in Siberia

1. Introduction

In the Siberian Arctic, air temperature increase is about four times higher in comparison with southern areas. That may lead to significant changes in the permafrost stability, melt water runoff, and hydrothermal regime within the cryolithozone [1,2]. Warming- driving permafrost melt is an important factor that may influence trees’ growth and vegetation productivity (NPP). Potentially, permafrost degradation may improve the soils’ hydrology regime and lead to an increase in the active layer depth. However, data on this problem are scarce [3,4,5,6]. Outside the Arctic, on the southern edge of cryolithozone, permafrost melting is considered as an additional water source during seasonal water stress [7]. Air temperature increase by +1 °C, according to prognostic estimates, may lead to a decrease in 25% of the permafrost within the upper 3-meter soil layer [8]. Alongside continuous permafrost melting, melting of the ice-reach soils accompanied by thermokarst, including lake formation. However, in spite of documented permafrost melting, within some parts of the cryolithozone, including the northeast Arctic, carbon- and ice-reach soils will be stable even after 2100 yr even under a pessimistic scenario [2].

Warming at high latitudes accompanied by vegetation GPP/NPP increase, including the “greening tundra” phenomenon [9,10]. However, warming leads to an increase in wildfires’ frequency and intensity [11]. On the whole, permafrost melting, alongside elevating air temperature, is considered a primary driver of vegetation dynamics within the permafrost area [9].

Estimates of warming influence on the permafrost melting and soil moisture are based mostly on the on-ground data or satellite microwave data [3,4]. Alongside that, satellite (GRACE) gravimetry also contained information related to the permafrost melting [12]. Gravimetric data are sensitive to all variables that influence local mass anomalies, including water mass anomalies. Namely, water mass anomalies are the primary variable that influenced GRACE-derived measurements. The most important results obtained based on gravimetry are the ice sheets melting in Greenland and the Arctic [13,14]. There are also indications that gravimetry is useful in the analysis of moisture regime changes influence on the tree (Larix spp.) growth in the Arctic [10,15].

Geographically, studies of permafrost melting are located mostly within the western part of the permafrost zone in North America, whereas vast permafrost areas in Siberia and Canada are poorly studied [9]. The goals of this study: (1) permafrost melt water runoff estimation beyond the Polar Circle in Siberia; (2) water anomalies and melt water runoff influence on the larch (Larix spp.) growth and vegetation GPP/NPP in the Siberian Arctic.

Checked hypothesis: warming-driven permafrost melt promotes (a) trees’ growth increase (Larix sibirica, L. gmelinii), the dominant tree species in the Arctic, and (b) promotes vegetation productivity (GPP/NPP) increase.

2. Materials and Methods

2.1. Study Area

The study area, the Siberian Arctic, included sparse larch stands, tundra, and forest-tundra communities and is located within the continuous permafrost zone (Figure 1). The area is beyond the Arctic Circle (66.6°N–78.0°N; 60°E–180°E, 168°W–180°W). The islands and coastline were removed from the analysis. The total analyzed is ~ 3×10⁶ sq. km. All the data were transformed to the Albers equal area conic projection and the nominal 463 m pixel size corresponding to the MODIS-based NPP data.

The climate of the region is characterized by its harshness, with subarctic and sharply continental features. Long winters and brief summers are the norm, accompanied by frequent occurrences of cyclones and anticyclones. The polar night brings with it low temperatures and blizzards, while the polar day is marked by intense insolation. Additionally, cloudy weather and extended periods of drizzle are common. Precipitation is distributed extremely unevenly throughout the year, with the majority falling between late May and late September. A stable snow cover begins to form in the third week of September and melts away in the first half of June. At higher elevations, the snow depth ranges from 0.4 to 0.8 meters. In contrast, at the foothills and in flat areas, it can reach up to 8–9 m. The permafrost extends up to a depth of 1 km. Taliks (melted soils) occur under the beds of major rivers. During the summer months, an active layer up to 2 meters thick develops. The depth of seasonal freezing in soil depends on factors such as slope exposure, vegetation cover, moisture content, and soil type. The lowland tundra ecosystems are marked by a high-water content, a paucity of soil fertility, the presence of marine and accretive sediments, polygonal wetlands, as well as “baijarakhs” (permafrost-caused “hills”). In the far north, the distribution of arboreal and shrub vegetation has been observed to be particularly sparse [16,17].

The majority of tree species are presented by larch (Larix gmelinii, L. sibirica) with a small proportion of birch (Betula spp.). Within the southwest part of the area Siberian pine (Pinus sibirica Du Tour) и spruce (Picea obovata Ledeb.) occurred along the creeks and rivers mostly.

2.2. Methods

We used satellite (GRACE, MODIS) gravimetry and NPP index data, dendrochronology-derived trees’ radial growth index (GI), and field data.

2.2.1. GRACE Data

We used monthly mass grid (land) products based on GRACE and GRACE-FO remote sensing data (https://grace.jpl.nasa.gov/data/get-data/monthly-mass-grids-land) [19]. GRACE mission is supported by the NASA MEaSUREs Program. The results of GRACE measurement given in EWTA (Equivalent of Water Thickness Anomalies), extracted from the TELLUS_GRAC_L3_GFZ_RL06_LND_v04 and TELLUS_GRFO_L3_GFZ_RL06.1_LND_v04 products [20]. The data were downloaded from the EarthData geoportal (https://search.earthdata.nasa.gov) in GeoTIFF format. The monthly land mass grids contain water mass anomalies given as EWTA in meters derived from GRACE and GRACE-FO time-variable gravity observations during the specified timespan (monthly), and relative to the specified time-mean reference period (2004–2009 years).

The EWTA represents the total terrestrial water storage anomalies from soil moisture, snow, surface water (incl. rivers, lakes, reservoirs, etc.), as well as groundwater and aquifers. The spatial sampling of all grids is 1 degree (above the Arctic Circle cell is ~94×40 km²) in both latitude and longitude [19]. The time coverage of the data is 2003–2023. EWTA values are presented in meters. In this study, alongside with EWTA term, we use the terms “water anomalies” and “water runoff”, where applicable.

Gravimetric data were calculated both for a hydrological year (from September of the previous year to August of the current year inclusively) and for the vegetation growth period (JJA).

2.2.2. MODIS Data

We analyzed changes in primary productivity based on the MODIS MOD17A3HGF [21]. The annual Net Primary Productivity (NPP) data were extracted from the MOD17A3HGF from 2007 to 2023 with 463 m spatial resolution [22].

The fires were mapped using a geospatial database created by the Forest Institute of the Siberian Branch of the Russian Academy of Sciences. The database uses satellite data obtained from the NOAA/AVHRR, TERRA/AQUA/MODIS, and SNPP/NOAA-20/VIIRS platforms. To exclude anthropogenic thermal anomalies, data about settlements were used (Open Street Map, https://www.openstreetmap.org/). Using high-resolution Google Earth data, thermal anomalies associated with oil and gas fields were excluded from the analysis. To estimate pyrogenic carbon emissions in the tundra, a value of 1.13 kg C / m² described in literature concerning northeastern Siberia was used [23].

2.2.3. Climate Data Processing

We used climate data extracted from the ERA5 Land database. The ERA5 Land was created by the European Centre for Medium-Range Weather Forecasts with a spatial resolution of 0.1 degrees (~9 km). ERA5 Land data are calculated based on a reanalysis of ground measurements, satellite imagery, and land cover and relief maps [24]. Monthly data on air temperature (T2m; °C), precipitation (TP; m) and evaporation (E; m), and soil temperatures (ST; m; at a depth of 0-100 cm) were downloaded using Copernicus Climate Change Service [24]. Soil temperature was considered for a layer of 0–100 cm corresponding mainly to the rooted zone. The data were preprocessed using R-Studio (v. 23.12.01; https://posit.co/download/rstudio-desktop). To estimate the influence of climate on the EWTA and vegetation cover dynamics, we calculated maps of Spearman’s correlation coefficients (p < 0.05). These maps showed spatial heterogeneity of the climate impact on changes in EWTA and NPP during 2007–2023. VPD (vapor pressure deficit) was calculated using the DewtoVPD function realized in the plantecophys Python library, and based on the air and dew temperatures, and barometric surface pressure [25].

2.2.4. Test Sites Description and Dendrochronology Analysis

Within on-ground study sites, the dominant tree species were larch (Larix sibirica Ledeb., Larix gmelinii Rupr.) with an admixture of birch (Betula pendula Roth) and, to a lesser extent, spruce (Picea obovata Ledeb.). Shrubs represented by alder (Dushekia fruticosa Rupr.), willow (Salix spp.), dwarf birch (Betula nana L.) and Vaccinium sp. The mean forest closure was about 0.2–0.3, although in some areas it may exceed 0.5. On-ground vegetation is presented mainly by lichens, moss, and sedges.

Dendrochronological analysis was carried out at four test sites: Ary-Mas (3 sample plots, 60 model trees); Kotuy (2 sample plots, 66 model trees), Pyasino (10 sample plots, 242 model trees) and Embechime (7 sample plots, 67 model trees) (Figure 1). Wood samples were measured using the LINTAB-6 platform with an accuracy of 0.01 mm [26]. As a result, a series of radial increments for each tree (in mm) were determined. The TSAP and COFECHA programs were used to check the quality of dating [27]. We applied a detrending procedure to account for the age trend by transforming the time series of tree-ring widths to a series of unitless indices with a mean of 1.0 and relatively constant variance [28]. We averaged tree-ring chronologies of individual trees to calculate generalized indexed tree-ring chronologies. For each test site, an indexed tree-ring chronology was constructed using the ARSTAN program; detrending was performed based on the linear regression or negative exponential curve [29]

2.2.5. Assessment of the Water Mass Runoff

We used GRACE EWTA, MODIS NPP, MODIS-based fire carbon emissions, and ERA5-Land total precipitation and evaporation data to assess water mass runoff (ΔS_m). We estimated changes in total water (ΔG_m) during 2007–2023 using the following equation:

ΔG_m = G_m(2007-2009) – G_m(2021-2023),

(1)

where G_m(2007-2009) and G_m(2021-2023) are the mean EWTA during the 2007–2009 and 2021–2023, respectively. A positive value of ΔG_m means water mass runoff and negative – water accumulation. In equation (1), the beginning year 2007 was chosen since that year is a stable decreasing trend of the water anomalies (WA or EWTA) occurred (Figure 2). To smooth variability, we averaged three years of data at the beginning (2007–2009) and the end (2021–2023) of the study period.

The total gravimetry mass influenced by vegetation cover mass ($\sum V_m$) changes, mass loss caused by wildfires (carbon emissions, ∑F), and differences (P-E) between precipitation (P) and evaporation (E) changes. To account for mass changes caused by vegetation growth, we used NPP estimations extracted from the MODIS MOD17A2HGF product [21]. The sum of annual NPP (∑NPP) during 2007–2023 years was calculated. Alongside the biomass accumulation, the produced NPP was partly decayed due to biological decomposition. To estimate loss caused by decomposition, a first-order exponential decay function was applied, which estimates the percent mass remaining at a certain time elapsed since the beginning of litter decomposition (Figure 3) [30,31,32]:

$$y=e^{–k\ast t},$$

(2)

where t – time (2007–2023); k – is the decomposition rate, which for arctic tundra is about ~0.03 g g^-1 yr^-1 [31,33,34].

We used function (2) to assess the remained biomass by the following equation:

$$\sum V_m={\displaystyle \sum }_{t=2007}^{2023}\left({\mathrm{NPP}}_t\ast e^{-k\left(2023-t+1\right)}\right),$$

(3)

We estimated wildfires-caused carbon loss by the following procedure. We mapped fire events by using the geospatial fire database created by the Sukachev Institute of Forest, SB RAS. This database uses satellite fire detections obtained from the NOAA/AVHRR, TERRA/AQUA/MODIS, and SNPP/NOAA-20/VIIRS platforms. We used settlement data from Open Street Map (https://www.openstreetmap.org/) to exclude thermal anomalies of anthropogenic origin. We also excluded from further analysis thermal anomalies related to gas flares in oil and gas fields using high-resolution data from Google Earth. In the case of larch forest stands in northern Siberia, emission estimates ranged between 3 and 3.44 kgC m^-2 [35,36,37]. The final estimate of carbon emissions was based on the data of Tsvetkov [35]. To estimate pyrogenic emissions in tundra, a value of 1.13 kgC m^-2 obtained for northeastern Siberia was used [23].

We applied linear regressions to account for annual changes in differences (α) between annual precipitation (P) and evapotranspiration (P-E):

$$y=\mathsf{\alpha }\ast t+\mathsf{\beta },$$

(4)

where t is time (yrs), and β is the intercept coefficient.

Finally, the corrected water mass runoff (WR) was estimated by equation (5):

$$WR=\mathsf{\Delta }{G}_m+\sum V_m–\sum F–l\ast \mathsf{\alpha },$$

(5)

where WR is the water runoff; ∑F – carbon emissions caused by fires; $\sum V_m$ – the remained biomass; $\mathsf{\alpha }$ is delta between (P-E); l is the analyzed period length (2007–2023).

2.2.6. Statistical Analysis

We applied ISODATA unsupervised classification realized in the ERDAS Imagine software (https://hexagon.com/products/erdas-imagine) to merge pixels to clusters (zones) with similar EWTA dynamics. In total, five zones were identified.

Statistical calculations were realized using StatSoft Statistica (http://statsoft.ru), (https://www.r-project.org), and R-Studio (https://posit.co/download/rstudio-desktop/) software.

We used Pearson’s and Spearman’s correlation coefficients and linear regressions to estimate the relationship between the estimated variables. The significance of correlation coefficients and regression equations was assessed using the t-test and F-test.

The Theil-Sen estimator was applied to calculate maps of EWTA, climate, and productivity (GPP) trends. Theil-Sen estimator is a non-parametric method that fits a regression line through the median of the slopes determined by all pairs of sample points [38], and it is more robust than simple linear regressions [39]. This estimator is less sensitive to outliers, and it is more accurate than simple linear regression [40]. We used the Theil–Sen estimator realized in the Python library pymannkendall 1.4.2 (https://pypi.org/project/pymannkendall (accessed on December 01, 2022)) imported into ESRI ArcGIS Desktop 10.8.1 (https://www.esri-cis.com/ru-ru/arcgis/products/arcgis-desktop/overview). As a result, a set of raster trend maps was generated, and p-levels were estimated.

3. Results

3.1. Melting Water Runoff Raw Data (WR_r)

Water runoff dominated within the Siberian Arctic with the exception of the eastern part (the Chukotka peninsula, Figure 4).

Throughout the Arctic, the values of WR_r formed five main clusters: West and Mid Siberian (WestSib and MidSib), the Lena River (LenaRiver), East Siberian (EastSib), and Chukotka Peninsula (ChukotkaPen) (Figure 5).

Water accumulation was observed within the ChukotkaPen cluster only. On average, during 2007–2023, the mean water runoff rate (row data) was about ~3 kg/m²/y (r² = 0.67) (Figure 6).

3.2. Data used for the Water Runoff Corrections: Mapped Data of NPP, $\sum V_m$ and Fire-Caused Carbon Loss

We corrected WR raw data based on mapped values of NPP, $\sum V_m$, and fire-caused carbon loss (Figure 6 and Figure 7).

Figure 6. (a) The accumulated NPP (∑NPP) during 2007–2023. Mean ∑NPP is 2.6 kg/m² (σ = 1.1; min = 0, max = 8.5). (b) Map of the remained biomass ($\sum V_m$). Mean ($\sum V_m$ is 2.1 kg/m² (σ = 0.8; min = 0.1, max = 6.6).

Figure 6. (a) The accumulated NPP (∑NPP) during 2007–2023. Mean ∑NPP is 2.6 kg/m² (σ = 1.1; min = 0, max = 8.5). (b) Map of the remained biomass ($\sum V_m$). Mean ($\sum V_m$ is 2.1 kg/m² (σ = 0.8; min = 0.1, max = 6.6).

The mean annual fire-caused carbon loss was 8.6 (σ = 17.5) Mt C/year or 0.08 kg C/m²/year.

The percentage of the annual burned area is about 0.17±0.34 %/year (Figure 7).

3.3. Corrected Melting Water Runoff (WR)

The map of corrected melting water runoff within the total study area is presented in Figure 8a. Water runoff within given Arctic sectors is given in Figure 8b and in Table 1. Thus, the mean annual WR within the Siberian Arctic is 7.7 ± 4.4 kg/m²/yr. Maximal values occurred in the West and Mid Siberian sectors (about 9.7 kg/m²/yr), and negative ones (i.e., water accumulation) were observed in the Chukotka Peninsula sector (-2.9 ± 3.2 kg/m²/yr) (Table 1).

3.4. Melting Water Influence on the Vegetation Productivity

Within the Arctic, positive NPP trends prevailed (~15%, whereas negative trends were observed in about ~1% of the area). Positive were trends observed mostly within the western Arctic (Figure 9a). Generally, NPP within the Arctic showed a strong increasing trend (Figure 9b).

The NPP correlated with WR mostly within northern and central parts of the western Arctic sectors (WestSib and MidSib; Figure 10a). Throughout the Arctic, positive correlations were observed within 4.5% of the analyzed area, negative and insignificant ones on the 1% and 94%, respectively (p < 0.1). Although within the WestSib sector, NPP significantly correlated with water runoff (Figure 10b; r = 0.49; p < 0.05). The latter referred to the higher permafrost melting in the West Siberian plain [41]. Meanwhile, NPP correlations with WA were observed throughout the Arctic (r = 0.55, p < 0.05) (Figure A1).

Alongside melting water runoff, NPP negatively correlated with precipitation and positively with air temperature (Figures 11a,b).

3.5. Melting Water Runoff Influence on the Larch Trees Growth

Larch trees’ growth index (GI) significantly correlated with water anomalies during the growth period (JJA) (with the exception of the Pyasino site): the less water content in permafrost, the higher GI of larch trees (r = –0.64…–0.85). For the hydrological year, correlations were significant for all sites (r = –0.57…–0.76; Figure 12a). Similarly, the more water runoff (i.e., less water content in soils), the higher tree’s growth (r = –0.51…0.57).

Alongside that, larch GI positively correlated with air and soil temperature and vapor pressure deficit (VPD) whereas negatively with precipitation (Figure 12b). Note: although, for the Embechime site, correlations were insignificant due to a shorter interval of observations (until 2012 only).

4. Discussion

In the Siberian Arctic, a long-term continuous trend of water content decrease in permafrost is observed, whereas an annual melt water runoff is increasing throughout the Arctic with the exception of the Chukotka Peninsula sector (Figure 2, Figure 8). Melt water runoff (with a mean rate of 7.7 ± 4.4 kg/m^-2/yr.) entailed soil moisture decrease (and, evidently, an active layer depth increase) that leads to an increase in vegetation cover productivity and trees’ growth index (GI) (Figure 10, 1Figure 12). A similar relationship has been established for larch trees (r = 0.55). Positive correlations between NPP and melt water runoff dominated in the western part of the Arctic (over 10% of the area), which is relevant to higher permafrost melting in that area [41]. Earlier larch growth and vegetation GPP positive correlations with soil moisture decrease were reported for the Mid-Siberian Arctic [10,15]. Larch trees’ (L. gmelinii) growth stimulation by permafrost melting is found also in the mountains of northeastern China. However, that was referred to moisture deficit in that area [7,42]. Meanwhile, seasonal water stress in spring also occurred in the permafrost zone while the water supply was blocked [7,43].

Arctic soils are typically waterlogging so-called "cold wet soils". Thus, root zone temperature, along with the air temperature, limits the vegetation growth. Soil temperature is determined also by soil moisture due to the high heat capacity of water. The found trends of water runoff and, in general, a water contend decrease indicated an improvement in the hydrothermal regime in the Arctic soils. In general, a water content decrease in permafrost soils resulted in an increase of NPP (r² = –0.31; Figure A1).

The hydrothermal regime improvement, together with warming, are triggers of larch tree growth and an increase in vegetation NPP. Meanwhile, vegetation productivity increase caused negative feedback since ground cover is a thermal insulator that decreases permafrost melting [9]. Along with meltwater runoff, water "conservation" is observed in forming bogs and ponds due to thermokarst, especially in the case of ice-enriched soils [2]. Warming also stimulates solifluction, i.e., the flow of waterlogged soils on slopes [9,44]. Thus, in the Taimyr Peninsula, the frequency of landslides has increased by more than ten times in the last decade. Notably, that was accompanied by a sharp (>20 times) increase in carbon fixation [45]. In some regions, the melting of permafrost soils was stimulated by an increase in "warm rains" [9]. However, within the studied area, we found a negative precipitation influence on the larch trees and vegetation growth (Figures 11p,q,r,t, 12b).

Primarily, permafrost melting reduces waterlogging of soils, increases soil drainage and active layer depth, and enriches the root layer with oxygen and biogenic elements. Thus, in the zone of non-continuous permafrost, the thawing depth reached 2–8 m [41]. However, the estimation of the melting water runoff influence on the active layer depth needs more studies.

The main factors determining the dynamics of gravitational mass (i.e., satellite GRACE data) include, along with the mass of water (in solid and liquid state), the NPP, precipitation, and wildfires’ variations. To consider NPP decay, we used the decay coefficient (k = 0.03 g g^-1 yr^-1; equation (3)) [31,33,34]. That coefficient may differ in different phytocenoses, but relevant data is not available in the literature. The precipitation and evapotranspiration influences are taken into account by equation (5). The estimated fire-caused carbon loss was surprisingly low compared to NPP (about two orders; Table 1). However, the majority of carbon emissions occurred southward of the Arctic Circle [46].

In general, warming-driven permafrost melt stimulated trees’ growth and vegetation productivity in the Arctic which, despite elevating the burning rate, still keeps the status of the carbon sink.

5. Conclusions

Based on the satellite gravimetry, we revealed a continuous decrease in water content in the permafrost and an increase (with the exception of the Chukotka Peninsula sector) of melt water runoff (with a rate of about 7.4 kg/m^-2/yr). A warming-driven permafrost degradation promoted larch trees’ growth and an increase in vegetation NPP. In the warming climate, permafrost melt leads to improvement of the soils’ hydrothermal regime that, in turn, stimulates vegetation growth and supports the Arctic as a carbon sink despite the burning rate increase.