35 Years of Vegetation and Lake Dynamics in the Pechora Catchment, Russian European Arctic

: High-latitude regions are a hot spot of global warming, but the scarce availability of observations often limits the investigation of climate change impacts over these regions. However, the utilization of satellite-based remote sensing data o ﬀ ers new possibilities for such investigations. In the present study, vegetation greening, vegetation moisture and lake distribution derived from medium-resolution satellite imagery were analyzed over the Pechora catchment for the last 35 years. Here, we considered the entire Pechora catchment and the Pechora Delta region, located in the northern part of European Russia, and we investigated the vegetation and lake dynamics over di ﬀ erent permafrost zones and across the two major biomes, taiga, and tundra. We also evaluated climate data records from meteorological stations and re-analysis data to ﬁnd relations between these dynamics and climatic behavior. Considering the Normalized Di ﬀ erence Vegetation Index (NDVI) and the Normalized Di ﬀ erence Moisture Index (NDMI) in the summer, we found a general greening and moistening of the vegetation. While vegetation greenness follows the evolution of summer air temperature with a delay of one year, the vegetation moisture dynamics seems to better concur with annual total precipitation rather than summer precipitation, and also with annual snow water equivalent without lag. Both NDVI and NDMI show a much higher variability across discontinuous permafrost terrain compared to other types. Moreover, the analyses yielded an overall decrease in the area of permanent lakes and a noticeable increase in the area of seasonal lakes. While the ﬁrst might be related to permafrost thawing, the latter seems to be connected to an increase of annual snow water equivalent. The general consistency between the indices of vegetation greenness and moisture based on satellite imagery and the climate data highlights the e ﬃ cacy and reliability of combining Landsat satellite data, ERA-Interim reanalysis and meteorological data to monitor temporal dynamics of the land surface in Arctic areas. (June, July, and August) and annual average temperature. The monthly precipitation and snow water equivalent were aggregated to the summer total and annual total. These monthly, summer and annual climatic data was estimated across the analyzed areas (Pechora catchment and the Pechora delta within di ﬀ erent permafrost types and major biomes) and compared with the satellite data observations on vegetation characteristics and lakes.

. Location of the study area (a), altitude (b), main biomes (c) and permafrost types (d) within Pechora catchment; Pechora delta is highlighted by a red square.
The second study area focuses on the Pechora Delta with an additional buffer zone and covers 18,133 km 2 , being a sub-region of the Pechora catchment. The delta has a length of 120 km, a total area of around 2900 km 2 and with an elevation range less than 6 m. A floodplain dominates the entire area characterized by numerous channels and lakes [40] as well as depressions occupied by polygonal peatlands and fens with peat thickness ranging from 0.5 to 5 m [41]. As permafrost develops more under convex and flat surfaces, the permafrost table is deeper in valleys that are usually dry and drained by streams, so that the area is geocryologically unstable [41]. Besides, the Pechora delta is an important breeding area and stopover site for migratory birds.

Vegetation Dynamics Analysis
The Landsat archive was utilized to select satellite images from missions 5, 7, and 8. Surface reflectance images with 30 m spatial resolution were chosen between 1985 and 2019. The peak of vegetation phenology was analyzed using only images from the summer season (July-August) [20,42,43] since the first snowfalls start from September in the Pechora catchment. 906 satellite images from Landsat 5, 172 from Landsat 7 and 126 from Landsat 8 mission were used in this study. The image collection in each year, comprising a various number of satellite images, was subjected to median temporal reduction to obtain one image per year. Every pixel in the final image holds the median value from the stack of pixels (one from each image). This process led to a satellite image covering the entire study area for 30 out of 35 analyzed years, containing the following spectral bands: blue, green, red, near-infrared, short-wave infrared 1 and short-wave infrared 2.
From Landsat spectral bands, two normalized indices were derived that illustrate the consistency of green vegetation and vegetation moisture [44].
Normalized Difference Vegetation Index (NDVI), developed by Rouse et al. (1974) [45], illustrates the greenness of vegetation, or more precisely, the magnitude of photosynthetic activity. The negative values, close to −1, correspond to water. Values close to zero generally correspond to areas without vegetation (rocks, sand). Positive low values (between 0.2 and 0.4), represent shrubs and meadows, and high values indicate the presence of forests. The formula uses the red (RED) and near-infrared (NIR) bands: Normalized Difference Moisture Index (NDMI) shows the moisture differences of the land surface, being highly correlated with the water content of the vegetation and is a good indicator of vegetation change, especially in forest applications [46]. NDMI is calculated using the near-infrared (NIR) and the short-wave infrared 1 (SWIR 1) bands, according to the following formula: The analysis was conducted using the cloud computing capabilities of the Google Earth Engine (GEE) with two aims. First, the vegetation indices were aggregated to annual averages of the summer season, resulting in charts illustrating their temporal evolution and trends between 1985 and 2019, within different geographical areas. Second, the slope coefficient of the linear regression was mapped within each pixel using the 30 image stack of vegetation indices. The maps reveal spatial general trends and dynamics of the vegetation, detecting the direction and magnitude of change over time, expressed as change of vegetation indices value per year. The linear regression method was used in many remote sensing studies, using time as the independent variable and vegetation indices as the dependent variable [47].

Lake Dynamics Analysis
The Global Surface Water dataset contains maps with the spatial and temporal distribution of surface water from 1984 to 2018 at 30 m resolution [48]. It was developed within the Copernicus program by the Joint Research Centre (JRC) and was generated by using 3,865,618 Landsat 5, 7 and 8 scenes acquired between 16 March 1984 and 31 December 2018. Each pixel was individually classified into water or land using an expert classification, and the results were compressed into a monthly archive for the entire time period [48].
From this dataset, we used transitions which provide information on the change in seasonality between the first and last year of the analyzed period and capture changes between the three classes of seasonal water, permanent water and others/no water [49]. The dataset contains raster data with mapped transitions: unchanging (stable) permanent water surfaces; new permanent water surfaces (conversion of land into permanent water); lost permanent water surfaces (conversion of permanent water into land); unchanging seasonal water surfaces; new seasonal water surfaces (conversion of land into seasonal water); lost seasonal water surfaces (conversion of seasonal water into land); conversion of permanent water into seasonal water; and the conversion of seasonal water into permanent water. These conversions refer to changes in state from the beginning and end of the time series and, therefore, do not describe the changes in the intervening years [49].
The raster dataset was resampled to a spatial resolution of 43 m, due to GEE user computation limits. The resolution was decreased stepwise from 30 m to minimize the information loss. In order to analyze only the changes related to lakes, the transitions dataset was processed by buffer analysis to remove all watercourses and changes adjacent to them. The freely accessible watercourses dataset was provided by OpenStreetMap [50]. In the absence of a lake area at the beginning of the time interval, the percentage was calculated from the total area of the region of interest. This analysis is a relevant indicator, particularly for areas with ice-rich permafrost, knowing that one important element specific to permafrost degradation is the dynamics of thermokarst lakes [19,21].

Climatic Data
The identified dynamics of vegetation greening and moisture, as well as the dynamics of lakes, were analyzed concerning climatic data from ERA-Interim atmospheric reanalysis [32] and meteorological observations at the Naryan-Mar (67.63 • N, 53.03 • E) and Cape Konstantinovski (68.55 • N, 55.42 • E) stations, located within Pechora catchment. Atmospheric reanalysis is a technique of assimilation of historical climatic data by means of numerical circulation models with the aim to remove eventual physical inconsistencies due to deficiency of meteorological observations (e.g., sparse distribution of meteorological stations or shortcomings of observational techniques). In this study ERA-Interim is used following recommendations by Lindsay et al. (2014) [51], suggesting that ERA-Interim in comparison to other reanalysis products stands out as being consistent with independent observations for applications in high northern latitudes. The spatial resolution of ERA-Interim is approximately 80 km (~0.75 • or T255 spectral grid), and various climatic parameters describing atmospheric as well as land-surface conditions are available at 3 and 6-hourly temporal resolutions. Dee et al. (2011) [32] provide a detailed description of the numerical model, data assimilation method, and input datasets used to produce the ERA-Interim reanalysis. ERA-Interim data are available at https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era-interim. In this study, we were using monthly means of 2 m air temperature, total precipitation, and snow water equivalent (SWE). These climatic parameters are spatially averaged for the discontinuous, sporadic, isolated, and permafrost free areas, as well as for essential biomes (tundra and taiga). Seasonal cycle for each region is shown on Figure 2, complemented with annual climate provided in tables for each variable. In addition to these spatially averaged regions, we have considered ERA-Interim grid points for geolocations of Naryan-Mar and Cape Konstantinovski meteorological station and compared the respective gridded data with hydrometeorological observations conducted at these two stations operated by RIHMI-WDC. Snow depth measurements are available from the RIHMI-WDC stations, while ERA-Interim contains information about SWE. To compare these two parameters, snow depths are converted into SWE by using snow density climatological values for the Russian tundra estimated by Zhong et al. [52]. Our analysis shows a good agreement of the seasonal cycle between ERA-Interim and observations ( Figure 2). In particular, the 2 m air temperature of ERA-Interim reanalysis agrees well with the measured values from the meteorological stations available in the Pechora catchment area, while precipitation and SWE show more considerable uncertainties. Small bias and high correlation between observed 2 m air temperature and corresponding ERA-Interim gridded values is a consequence of the ERA-Interim assimilation method i.e., its optimal interpolation scheme [32] that uses the observed 2 m temperature and background temperatures from the previous analysis time step to calculate current time step temperature.  [52].
The southern part of the catchment show average annual air temperature above −2 °C [36]. The amount of annual precipitation sums to 484.5 mm at Naryan-Mar and 411.3 mm at Cape Konstantinovski ( Figure 2). However, the annual precipitation varies with location (e.g., annual precipitation amounted 640.5 mm spatially averaged for the Pechora catchment area, in the taiga area 664.4 mm while in tundra 586.2 mm, according to ERA-Interim) ( Figure 2). These results are in good agreement with previously estimated by Van der Sluis et al. [37], who estimate a range of 650-750 mm/year in the taiga, and more than 1200 mm/year in the highest parts of the Ural Mountains. The Pechora catchment is a neotectonic structural lowering area characterized by peatlands favored by biogenic soils and thermokarst lakes. At the same time, bogs occur only in the sub-zone of discontinuous permafrost [36], with various shares of underground runoff in the total runoff, ranging between 10 and 60 %, the lowest underground runoff being found in the tundra forced by the presence of permafrost [37]. The landscapes of this area experienced relatively low anthropogenic-induced ecological impacts and disturbances, however, the potential for extensive environmental degradation is now greater than ever before, due to the plans for this region to become one of the primary oil and gas and other natural resources center for Russia [38,39]. Three-quarters of this area is covered by different types of permafrost (according to the permafrost zonation model of Obu et al. [33]) with approximately 4.5% of discontinuous permafrost mostly in the north-eastern part of the area, 39.5% of sporadic and 29% of isolated permafrost (Figure 1d).
The second study area focuses on the Pechora Delta with an additional buffer zone and covers 18,133 km 2 , being a sub-region of the Pechora catchment. The delta has a length of 120 km, a total On the other hand, discrepancies between ERA-Interim and observed precipitation and SWE indicate challenges related to the representation of the hydrological cycle in atmospheric reanalysis as well as limitations of observational techniques. A possible explanation for this might be that the ratio between snowfall and rainfall in reanalysis precipitation components depends on the model physics.
In contrast, precipitation measurements may be affected by wind resulting in reduced precipitation catch accuracy (e.g., raindrops and especially snow blown away from gauges typically lead to an underestimation of the actual precipitation values). There are, however, other possible explanations related to uncertainty of SWE and snow density estimated by the snow model in the process of data assimilation of atmospheric reanalysis [53], as well as random and systematic errors occurring during routine snow surveys [52]. Nevertheless, it is widely accepted that atmospheric reanalyses contain the most consistent information on historic climate. However, comparison with observation is necessary to estimate the range of possible errors. Our analysis confirms previous findings [51] that ERA-Interim data are suitable for high-northern latitudes applications.
The monthly air temperature at 2 m estimated ERA-Interim, between 1979 and 2018, was aggregated to annual summer temperature (June, July, and August) and annual average temperature. The monthly precipitation and snow water equivalent were aggregated to the summer total and annual total. These monthly, summer and annual climatic data was estimated across the analyzed areas (Pechora catchment and the Pechora delta within different permafrost types and major biomes) and compared with the satellite data observations on vegetation characteristics and lakes. Cross-correlation function (CCF) implemented in R software was used to identify time series correlations and lags between climate parameters and vegetation greenness and moisture. The correlation coefficient (R-coefficient) and statistical significance were used to reveal relationships of time series data. Within Pechora catchment and delta, the time series data tables containing average NDVI and NDMI values each year and the corresponding annual summer temperature, annual average temperature, annual summer precipitation, annual total precipitation, and annual total snow water equivalent were used to conduct CCF analysis.

Vegetation Dynamics
In the Pechora catchment the NDVI mean value ranges temporarily from a minimum of 0.  Figure 3). However, the NDVI only started to increase after 2002, before the NDVI was almost constant or even showed a slight decrease. These results suggest a slight vegetation decline in the Pechora catchment until 2002, but then significantly developed afterwards. The results are in accordance with the average summer air temperature, which shows several colder summers before 2002, but increased frequency of warmer summer since then. Overall, the NDVI tends to follow the average summer temperature, in most cases with a delay of one year. The CCF analysis yields statistically significant values for the highest correlation coefficients of 0.39 at lag 0, and 0.35 at 1-year lag. For example, a warm summer is followed by a peak of NDVI in the next year, such as in 1989-1990, 1993-1994, 2000-2001  The spatial distribution of regression slopes of NDVI in the Pechora catchment scale is heterogeneous with an evident general greening trend. Areas showing decreases (browning) or only slight increases in NDVI are distributed mainly at higher altitudes and seem to be correlated with the regression slopes of negative NDMI ( Figure 4).
In the Pechora Delta, the NDVI does not show a notable trend over the whole period. However, the NDVI has a pronounced variability that can be separated into three distinct periods. During the first period until 2001/2002, NDVI shows an increasing trend, especially after 1987. This is followed by a decrease until 2007. After 2007, the NDVI was strongly increasing from 0.34 in 2007 to 0.51 in 2019. Like in the Pechora catchment, the vegetation dynamics in the Pechora Delta region is quite consistent with the average summer temperature development, also showing a one-year delay for a significant number of years, the highest correlation coefficient of 0.31 being reached at 1-year lag. Within the period 2001-2006, the pronounced NDVI decrease is most probably due to the lack of cloud-free satellite images, since the average summer temperature shows an increase between 2002 and 2004. Overall, the vegetation greenness is increasing in the Pechora Delta after 1987. Other studies have shown that the increase in NDVI is correlated with the rise in temperature and decrease of precipitation over the same area [54].
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 25 peak of NDVI in the next year, such as in 1989-1990, 1993-1994, 2000-2001   The spatial distribution of regression slopes of NDVI in the Pechora catchment scale is heterogeneous with an evident general greening trend. Areas showing decreases (browning) or only slight increases in NDVI are distributed mainly at higher altitudes and seem to be correlated with the regression slopes of negative NDMI ( Figure 4).  In the Pechora Delta, the NDVI does not show a notable trend over the whole period. However, the NDVI has a pronounced variability that can be separated into three distinct periods.  The same tendencies were observed in all permafrost zones of the Pechora catchment, however with some delays. In the areas with discontinuous permafrost, the NDVI started to increase noticeably in 2007, while in sporadic permafrost zones an increase was observed distinctively after 2011. In the isolated and permafrost free areas, the NDVI began to reverse its decreasing trend in 2002, coinciding with the trend in the entire Pechora catchment. The gain in NDVI between 1985 and 2019 is much more similar in permafrost areas, with 0.127 on discontinuous permafrost, 0.136 on isolated permafrost and 0.126 on sporadic permafrost, as compared to the permafrost free area, with a gain of 0.088. These results suggest a slightly different behavior of vegetation in permafrost and permafrost free areas, as the vegetation developed in the last years more substantially in permafrost areas ( Figure 5).
with some delays. In the areas with discontinuous permafrost, the NDVI started to increase noticeably in 2007, while in sporadic permafrost zones an increase was observed distinctively after 2011. In the isolated and permafrost free areas, the NDVI began to reverse its decreasing trend in 2002, coinciding with the trend in the entire Pechora catchment. The gain in NDVI between 1985 and 2019 is much more similar in permafrost areas, with 0.127 on discontinuous permafrost, 0.136 on isolated permafrost and 0.126 on sporadic permafrost, as compared to the permafrost free area, with a gain of 0.088. These results suggest a slightly different behavior of vegetation in permafrost and permafrost free areas, as the vegetation developed in the last years more substantially in permafrost areas ( Figure 5).
Throughout the entire period time, the NDVI behaved similarly over the permafrost free area and over areas with isolated and sporadic permafrost, showing more or less the same trends and evolution. However, areas with discontinuous permafrost exhibited much more variation throughout the 35 years. Except for the discontinuous permafrost zone where no significant correlation was found, the highest coefficients of correlation between NDVI time series and average summer temperature was found at 1-year lag. For example, peaks and minimums in temperature were followed by NDVI in 1986-1987, 1989-1990, 1993-1994, 1999-2000, 2000-2001  Analyzing the essential biomes of tundra and taiga, the relative differences are small in terms of NDVI evolution, as the trends for both biomes are the same as for the entire Pechora catchment. For both biomes, NDVI was slightly decreasing until 2002, followed by a strong increase with a gain in NDVI from 2002 until 2019 by 0.225 in the taiga area and 0.174 in the tundra. Over the entire period (1985-2019), the NDVI increased by 0.128 in tundra and 0.117 in taiga area. As for the cases Throughout the entire period time, the NDVI behaved similarly over the permafrost free area and over areas with isolated and sporadic permafrost, showing more or less the same trends and evolution. However, areas with discontinuous permafrost exhibited much more variation throughout the 35 years. Except for the discontinuous permafrost zone where no significant correlation was found, the highest coefficients of correlation between NDVI time series and average summer temperature was found at 1-year lag. For example, peaks and minimums in temperature were followed by NDVI in 1986-1987, 1989-1990, 1993-1994, 1999-2000, 2000-2001, 2010-2011 or 2016-2017 ( Figure 5).
Analyzing the essential biomes of tundra and taiga, the relative differences are small in terms of NDVI evolution, as the trends for both biomes are the same as for the entire Pechora catchment. For both biomes, NDVI was slightly decreasing until 2002, followed by a strong increase with a gain in NDVI from 2002 until 2019 by 0.225 in the taiga area and 0.174 in the tundra. Over the entire period (1985-2019), the NDVI increased by 0.128 in tundra and 0.117 in taiga area. As for the cases considered above, the vegetation dynamic tends to be in accordance with the average summer temperature both in tundra and taiga, showing the same delay of one year ( Figure 6). Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 25 considered above, the vegetation dynamic tends to be in accordance with the average summer temperature both in tundra and taiga, showing the same delay of one year ( Figure 6). Regarding vegetation moisture, both Pechora catchment and Pechora Delta show a clear increasing trend in the last 35 years. In the Pechora catchment the NDMI increased by 0.07 between 1985 and 2019, while in the Pechora Delta by 0.18. Unexpectedly, the NDMI trend did not show any decrease before 2011, like NDVI. Contrary to the NDVI, which was following average summer temperature rather than the annual average, the vegetation moisture dynamics seems to better concur with annual total precipitation (R-coefficient of 0.4) rather than summer precipitation (R-coefficient of 0.2), with no lag. NDMI changes are also explained by the annual snow water equivalent, showing an even higher R-coefficient of 0.45. The annual total precipitation shows an increasing trend between 1985 and 2005/2007, both in the Pechora catchment and the delta region, in accordance with the rising trend of NDMI. After the 2005/2007, the annual total precipitation slightly decreases, while NDMI shows only a weak increase (Figure 7).
Our results show that vegetation moisture illustrated by NDMI, similar to vegetation greening illustrated by NDVI, exhibited an increasing and similar trend in the last 35 years within isolated and sporadic permafrost areas, as well as in permafrost free areas. In contrast to these areas, NDMI over the discontinuous permafrost area showed, like for NDVI, much more variation and a slighter increasing trend, suggesting an influence of this type of permafrost on vegetation processes. The gain in NDMI in the last 35 years is lower in the permafrost free area (0.03), compared to the relatively similar values in the permafrost areas, with 0.063 on discontinuous permafrost, 0.075 on isolated permafrost and 0.098 on sporadic permafrost.
These results strengthen the statement above, that permafrost determines a slightly different behavior of vegetation, tending to accelerate vegetation greening and increase moisture. This assumption has meaning in the context of a warming climate, with consequences on permafrost degradation. This is a long term process that also determines an increase in soil temperature that might accelerate the plant growth as revealed by other studies [9,55]. The vegetation moisture increase before 2007 is in accordance with annual total precipitation and annual snow water equivalent data within all permafrost types. This accordance changes somewhat after 2007 when the annual total precipitation tends to decrease while NDMI still increases for all areas except for the discontinuous permafrost zone. For the latter, NDMI also tends to decrease after 2011 (Figure 8).  Our results show that vegetation moisture illustrated by NDMI, similar to vegetation greening illustrated by NDVI, exhibited an increasing and similar trend in the last 35 years within isolated and sporadic permafrost areas, as well as in permafrost free areas. In contrast to these areas, NDMI over the discontinuous permafrost area showed, like for NDVI, much more variation and a slighter increasing trend, suggesting an influence of this type of permafrost on vegetation processes. The gain in NDMI in the last 35 years is lower in the permafrost free area (0.03), compared to the relatively similar values in the permafrost areas, with 0.063 on discontinuous permafrost, 0.075 on isolated permafrost and 0.098 on sporadic permafrost.
These results strengthen the statement above, that permafrost determines a slightly different behavior of vegetation, tending to accelerate vegetation greening and increase moisture. This assumption has meaning in the context of a warming climate, with consequences on permafrost degradation. This is a long term process that also determines an increase in soil temperature that might accelerate the plant growth as revealed by other studies [9,55]. The vegetation moisture increase before 2007 is in accordance with annual total precipitation and annual snow water equivalent data within all permafrost types. This accordance changes somewhat after 2007 when the annual total precipitation tends to decrease while NDMI still increases for all areas except for the discontinuous permafrost zone. For the latter, NDMI also tends to decrease after 2011 ( Figure 8).  For the taiga and tundra biomes, vegetation moisture shows the same increasing trends within the entire Pechora catchment in the last 35 years. However, the gain in NDMI is almost double in tundra biome, with a gain of 0.108, compared to the gain of 0.059 in the taiga biome. The vegetation

Lake Dynamics
In the Pechora catchment, 161,119 permanent lakes were identified for the analyzed period (1984-2018). Their unchanging water area of 592,289 ha represents 1.47% of the total area of the Pechora catchment. The permanent lakes have a density of 0.4 lakes/km 2 , an average area of 3.7 ha, a minimum area of 0.19 ha, and the largest lake covers 18,736.3 ha. The minimum area corresponds to the resolution of the resampled Global Surface Water dataset, representing one pixel.
The area of permanent lakes experienced a net loss of 2685.1 ha during 1984-2018, with a gain of 4277.8 ha through the conversion of land into permanent lakes, but a loss of 6962.9 ha through the transformation of permanent to seasonal lakes ( Figure 10).
The stable water area of seasonal lakes was 200,101 ha (0.49%), while the number of seasonal lakes was 255,703. Therefore, an average seasonal lake area was 0.79 ha. The smallest seasonal lake covers 0.19 ha and the largest, 10,011.9 ha, with a density of 0.63 lakes/km 2 . The seasonal lake area showed the highest changes with a net gain of 158,986.5 ha (0.397%). Increases in this area were associated with the loss of land (0.51%) and permanent lakes (0.054%), while decreases were induced by land expansion (0.13%) and permanent lake expansion (0.037%) ( Figure 10). The net gain of seasonal lakes could be related to annual snow water equivalent, which increased by 150 mm between 1984 and 2018.
In the Pechora Delta, permanent lakes cover an area of 100,089.4 ha (5.52% of the Pechora Delta study area) between 1984 and 2018. Counting 24,037 permanent lakes, the average density reaches 1.33 lakes/km 2 , more than three times higher than lake density in the Pechora catchment. The average area of permanent lakes is 4.2 ha, with a minimum of 0.19 ha. The largest permanent lake in Pechora Delta is also the largest in the entire Pechora catchment, with an area of 18,736.3 ha. The same change tendencies as for the entire catchment emerged in the Pechora Delta, with a net loss of permanent lake area of 2393.9 ha (0.13%) (Figure 11).

Lake Dynamics
In the Pechora catchment, 161,119 permanent lakes were identified for the analyzed period The stable area of 39,418 seasonal lakes is 39807.8 ha (2.2%), leading to an average lake area of approximately 1 ha. The smallest seasonal lake in the Pechora Delta is 0.19 ha and the largest covering 670.9 ha, with an average density of 2.17 lakes/km 2 . Throughout the analyzed period, seasonal lakes increased by 16,457.9 ha (0.91%) (Figure 11), which is also most probably related to the increase of annual snow water equivalent by 110 mm. The stable water area of seasonal lakes was 200,101 ha (0.49%), while the number of seasonal lakes was 255,703. Therefore, an average seasonal lake area was 0.79 ha. The smallest seasonal lake covers 0.19 ha and the largest, 10,011.9 ha, with a density of 0.63 lakes/km 2 . The seasonal lake area showed the highest changes with a net gain of 158,986.5 ha (0.397%). Increases in this area were associated with the loss of land (0.51%) and permanent lakes (0.054%), while decreases were induced by land expansion (0.13%) and permanent lake expansion (0.037%) (Figure 10). The net gain of seasonal lakes could be related to annual snow water equivalent, which increased by 150 mm between 1984 and 2018.
In the Pechora Delta, permanent lakes cover an area of 100,089.4 ha (5.52% of the Pechora Delta study area) between 1984 and 2018. Counting 24,037 permanent lakes, the average density reaches 1.33 lakes/km 2 , more than three times higher than lake density in the Pechora catchment. The average area of permanent lakes is 4.2 ha, with a minimum of 0.19 ha. The largest permanent lake in Pechora Delta is also the largest in the entire Pechora catchment, with an area of 18,736.3 ha. The same change tendencies as for the entire catchment emerged in the Pechora Delta, with a net loss of permanent lake area of 2393.9 ha (0.13%) (Figure 11). The stable area of 39,418 seasonal lakes is 39807.8 ha (2.2%), leading to an average lake area of approximately 1 ha. The smallest seasonal lake in the Pechora Delta is 0.19 ha and the largest covering 670.9 ha, with an average density of 2.17 lakes/km 2 . Throughout the analyzed period, seasonal lakes increased by 16,457.9 ha (0.91%) (Figure 11), which is also most probably related to the increase of annual snow water equivalent by 110 mm. The highest percentage of permanent lakes from the total area of permafrost type (4.1%) and the highest permanent lake density (1.1 lakes/km 2 ) were found for discontinuous permafrost, followed by sporadic and isolated permafrost, with the lowest values in permafrost free regions (0.03% and 0.01 lakes/km 2 respectively). Isolated permafrost and permafrost free regions feature fewer lakes with higher lake average size (4.5 ha), while discontinuous and sporadic permafrost areas are covered by much more permanent lakes with a smaller size (3.6-3.7 ha) ( Table 1). The latter is the dominant waterbody type in the Arctic region [56]. Area (ha) Figure 11. Changes regarding permanent and seasonal lakes in Pechora delta.
The stable area of 39,418 seasonal lakes is 39807.8 ha (2.2%), leading to an average lake area of approximately 1 ha. The smallest seasonal lake in the Pechora Delta is 0.19 ha and the largest covering 670.9 ha, with an average density of 2.17 lakes/km 2 . Throughout the analyzed period, seasonal lakes increased by 16,457.9 ha (0.91%) (Figure 11), which is also most probably related to the increase of annual snow water equivalent by 110 mm.
The highest percentage of permanent lakes from the total area of permafrost type (4.1%) and the highest permanent lake density (1.1 lakes/km 2 ) were found for discontinuous permafrost, followed by sporadic and isolated permafrost, with the lowest values in permafrost free regions (0.03% and 0.01 lakes/km 2 respectively). Isolated permafrost and permafrost free regions feature fewer lakes with higher lake average size (4.5 ha), while discontinuous and sporadic permafrost areas are covered by much more permanent lakes with a smaller size (3.6-3.7 ha) ( Table 1). The latter is the dominant waterbody type in the Arctic region [56]. Stable seasonal lakes show different behavior depending on the permafrost type. The highest percentages occur for sporadic (1%) and discontinuous permafrost (0.8%) with reversed density values (1.8 lakes/km 2 in discontinuous permafrost, and 1.2 lakes/km 2 in sporadic permafrost). The isolated permafrost seems to favor the development of larger seasonal lakes (mean area of 0.9 ha), but it is covered by a lower percentage with a much lower lake density ( Table 2). The results reveal that discontinuous and sporadic permafrost type seems to favor the development of a much denser lake network, covering a much higher percentage than isolated and permafrost free regions, regardless of lake type. Permanent lake area decreased by 0.01% (1520 ha) and 0.12% (2144.5 ha) across sporadic and discontinuous permafrost, respectively. The isolated and permafrost free region exhibited positive net changes regarding permanent lakes, but with much lower amounts, of 0.01% (863.3 ha) and 0.001% (110 ha), respectively ( Figure 12). Stable seasonal lakes show different behavior depending on the permafrost type. The highest percentages occur for sporadic (1%) and discontinuous permafrost (0.8%) with reversed density values (1.8 lakes/km 2 in discontinuous permafrost, and 1.2 lakes/km 2 in sporadic permafrost). The isolated permafrost seems to favor the development of larger seasonal lakes (mean area of 0.9 ha), but it is covered by a lower percentage with a much lower lake density (Table 2). The results reveal that discontinuous and sporadic permafrost type seems to favor the development of a much denser lake network, covering a much higher percentage than isolated and permafrost free regions, regardless of lake type. Permanent lake area decreased by 0.01% (1520 ha) and 0.12% (2144.5 ha) across sporadic and discontinuous permafrost, respectively. The isolated and permafrost free region exhibited positive net changes regarding permanent lakes, but with much lower amounts, of 0.01% (863.3 ha) and 0.001% (110 ha), respectively ( Figure 12). Seasonal lakes exhibited more pronounced changes, with higher gains and losses over all types (Figure 12). The net gain of seasonal lakes is also in accordance with the annual snow water equivalent, which increased by 200 mm across discontinuous permafrost and 130 mm across sporadic permafrost between 1984 and 2018.
The analysis of lake changes by permafrost types reveals that sporadic and discontinuous permafrost impose a much higher dynamics of lake area across both lake types.
The analysis of the two major biomes shows that in the tundra region a much denser stable lake network exists with higher cover percentage (0.99 lakes/km 2 and 3.8% for permanent lakes, 1.38 lakes/km 2 and 0.8% for seasonal lakes). The tundra region also has larger permanent lakes (average size of 3.9 ha compared to 2.9 ha in taiga), and smaller seasonal lakes (0.55 ha compared to 1.37 ha in taiga) (Tables 3 and 4). Lakes in the tundra region exhibited much higher gains and losses, both for permanent and seasonal lakes. In the tundra region, the net loss of permanent lakes is 0.05% (6153 ha), while in taiga biome a net gain of 0.01% (3462.2 ha) occurred.
As expected, seasonal lakes show much higher net changes, with larger expansions and losses in the tundra region compared to the taiga biome. The net changes of seasonal lakes are positive both in tundra and taiga biomes, with 0.26% (34,159.3 ha) and 0.46% (124,812.6 ha) respectively. Like previous results, the seasonal lake expansion is in accordance with annual snow water equivalent, which increased by 110 mm across the tundra region, and 170 mm across taiga. It is worth to mention that while in taiga biome, both permanent and seasonal lakes showed positive changes, in tundra biome the permanent lakes reduced their area ( Figure 13). Seasonal lakes exhibited more pronounced changes, with higher gains and losses over all types ( Figure 12). The net gain of seasonal lakes is also in accordance with the annual snow water equivalent, which increased by 200 mm across discontinuous permafrost and 130 mm across sporadic permafrost between 1984 and 2018.
The analysis of lake changes by permafrost types reveals that sporadic and discontinuous permafrost impose a much higher dynamics of lake area across both lake types.
The analysis of the two major biomes shows that in the tundra region a much denser stable lake network exists with higher cover percentage (0.99 lakes/km 2 and 3.8% for permanent lakes, 1.38 lakes/km 2 and 0.8% for seasonal lakes). The tundra region also has larger permanent lakes (average size of 3.9 ha compared to 2.9 ha in taiga), and smaller seasonal lakes (0.55 ha compared to 1.37 ha in taiga) (Tables 3 and 4).  Lakes in the tundra region exhibited much higher gains and losses, both for permanent and seasonal lakes. In the tundra region, the net loss of permanent lakes is 0.05% (6153 ha), while in taiga biome a net gain of 0.01% (3462.2 ha) occurred.
As expected, seasonal lakes show much higher net changes, with larger expansions and losses in the tundra region compared to the taiga biome. The net changes of seasonal lakes are positive both in tundra and taiga biomes, with 0.26% (34,159.3 ha) and 0.46% (124,812.6 ha) respectively. Like previous results, the seasonal lake expansion is in accordance with annual snow water equivalent, which increased by 110 mm across the tundra region, and 170 mm across taiga. It is worth to mention that while in taiga biome, both permanent and seasonal lakes showed positive changes, in tundra biome the permanent lakes reduced their area ( Figure 13).

Discussion
The Landsat images proved to be suitable for assessing the vegetation response to the increasing air temperature since they provide a fair temporal resolution and spatial coverage [20]. Our Landsat NDVI analysis over the Pechora catchment indicates that this sensitive Arctic region experienced a gradual moderate greening trend, mainly since 2002. NDVI values of the entire catchment split by permafrost categories and main vegetation type displayed positive increasing trends in all areas, including the Pechora Delta. Greening trend values of NDVI/year are similar to others captured elsewhere in the northern regions [20,27,43]. Thus, the NDVI has a mean decadal increase of 0.031 in the Pechora catchment over the last 35 years, whereas in the Lena Delta an increase of 0.035 was calculated between 1999 and 2014 [20]. A similar gradual moderate greenness trend (0.003 unit NDVI/year) over 1984-2012 was derived for Canada and Alaska in the US [27]. For the Arctic Mackenzie Delta, a mean decadal NDVI increase of 0.031 was calculated for the 1985-2018 period [43]. The results mentioned above are in agreement with recent findings by Epstein et al. (2018) [57], revealing an overall greening trend throughout most of the Arctic tundra. According to this report, in Eastern Siberia, North Slope of Alaska and the southern zones of the Canadian tundra the NDVI increase is stronger.
Moreover, the spatial distribution of NDVI trend slope reveals the prevalence of regions with extensive and moderate greening of the Pechora catchment, with only a few exceptions ( Figure 4). Thus, the north-eastern extremity and the Ural Mountains are associated with the lowest increasing rate of NDVI (0.022 unit NDVI/decade). In contrast to these areas underlain by discontinuous permafrost, the central and southern lowlands of the catchment, dominated by the taiga biome are characterized by a strong greening trend (0.032 unit NDVI/decade). A pronounced greening tendency was also observed in the floodplains in the Pechora Delta (0.030 unit NDVI/decade). At the catchment scale, the distribution of the NDVI trend slope generally decreases as the altitudes rise. Thus, above 500 m the increasing rate of NDVI is relatively small (0.023 unit NDVI/decade), whereas between 200 and 500 m a maximum of 0.033 unit NDVI/decade was calculated.
Arctic greening is generally associated with the lengthening of the growing season [26], changes in land cover type [58], or increases in vegetation productivity [59] as a consequence of warming atmosphere temperatures [60]. In this study, we show that the average summer air temperature is a good proxy for vegetation productivity in the Pechora catchment. Based on plot-based evidence for vegetation transformation or tree-ring analysis, other studies [25,58,61] also revealed increased productivity over the Arctic tundra vegetation due to recent summer warming. Analyzing the relationships between summer maximum NDVI and climatic factors (temperature and precipitation) across four biomes in northern West Siberia, Miles et al., (2019) [62] highlighted strong relationships between the peak annual value of NDVI and June-July surface air temperature. Based on similar findings, other studies [63,64] assumed a recent increase in canopy abundance, height, and radial growth due to increasing summer warmth. Along coastlines, the sustained warming pattern is linked to the decrease of sea ice concentration and, according to Bhatt et al. (2010) [65], the sea-ice retreat is a driving factor for positive NDVI trends. The Pechora Sea located in the south-eastern corner of the Barents Sea is seasonally ice-covered, experiencing the largest fraction of ice loss in March throughout the Northern Hemisphere [66].
Recent studies showed that the most prominent increase of NDVI might be found in the southern tundra and northern taiga [67] where the abundance of relatively tall vegetation leads to greater productivity and a decreased albedo compared with low vegetation of the typical arctic tundra, composed by mosses, lichens, herbs, and small shrubs [68]. In the Pechora catchment, 'tall shrub tundra' dominated by dwarf birch (Betula nana L.), lower shrubs and forest-tundra prevail [69] and this might explain the positive NDVI trend. According to the Circumpolar Arctic Vegetation Map [34], the greatest part of the Pechora taiga belongs to bioclimate subzone E and is dominated by shrubs greater than 40 cm tall. This subzone has undergone significant positive trends in the peak annual value of NDVI, both in Eurasia and North America in the last three decades [67]. The north-eastern part of the catchment is characterized by the occurrence of erect dwarf-shrub tundra, which belongs to subzone D. Here, small shrubs (less than 40 cm tall) occupy most of the wet areas, whereas lichen-rich dwarf shrub tundra is common on the sandy soils of the Pechora Sea shore. In this area underlain by discontinuous permafrost, the mean decadal NDVI increase is only 0.022, considerably lower than on sporadic and isolated permafrost area or permafrost-free terrain. Based on our findings it seems that dwarf-shrub tundra in the north-eastern part of the Pechora catchment responds in a different way to increasing summer warmth and has lower productivity. Similar inconsistencies were reported in different Arctic sites, where vegetation productivity was not tightly linked with summer air temperature variations [70].
In the discontinuous permafrost regions, significant differences in productivity and thus, high spatio-temporal variations of NDVI are associated with acidic versus nonacidic soils [63]. The dynamics of NDVI can also be attributed to the composition of vegetation species assemblages since in the discontinuous permafrost zone small shrubs mainly occur in a matrix with lichen-covered bedrock [34]. Despite the overall greening illustrated by the positive NDVI trend throughout most of the Pechora catchment, signs of browning have also been observed, particularly at higher altitude. Several localized hot-spots of browning were also noticed near the river channels, similar to those mentioned in the Lena Delta region [20]. These negative NDVI trends correspond to local-scale processes and might be linked with an increase in bare ground cover as a result of permafrost thaw [71], the occurrence of retrogressive thaw slumps [72], or disturbances of vegetation due to selective foraging and trampling by large herbivores [73]. In the Pechora catchment, the NDMI trend correlates well with the increase of NDVI. For NDMI a mean decadal increase of 0.014 is calculated for the entire Pechora catchment, with the highest increase in the Pechora Delta (0.022 unit NDMI/decade) and the smallest among the areas underlain by discontinuous permafrost (0.004 unit NDMI/decade). The NDMI analysis also revealed strong wetting in the Lena Delta, where a decadal increase of 0.0271 was determined [20]. In the Pechora catchment, as the altitudes increase, the NDMI slope trend decreases considerably. In the Ural Mountains, above 500 m, a low increase was observed (0.004 unit NDMI/decade), whereas between 0 and 100 m the increase was nearly four times larger (0.015 unit NDMI/decade). Similar moderate wetting trends characterize the landscapes with sporadic and isolated permafrost as well as permafrost free regions (between 0.013 and 0.016 unit NDMI/decade). Insignificant differences in the NDMI trend were observed when comparing tundra (0.015 unit NDMI/decade) and taiga (0.014 unit NDMI/decade) biomes. Related to lake dynamics in the Pechora catchment and Pechora Delta, an overall slight decreasing trend has been observed in permanent lakes for the last 35 years. The results indicate a net loss of nearly 27 km 2 in the surface area of water in permanent lakes, which is insignificant compared to other arctic regions. Similar findings of disappearing or shrinking lakes were reported for lakes located in different parts of the Arctic [74][75][76]. In Canada, a net reduction of 6700 km 2 in lakes area was reported for a territory 25 times larger [75], whereas in Central Siberia large lakes (greater than 40 ha) declined by 11% [74]. For comparison in the Pechora catchment, the most significant decreasing (0.12%) was observed in areas underlain by discontinuous permafrost. In sporadic and isolated permafrost regions and in permafrost free areas the majority of lakes remained relatively stable. Intensive gross lake area loss (−12 to −29%) was reported in discontinuous permafrost regions in Western Siberia [77]. Pokrovsky et al. (2013) [78], also reported rapid decreasing size and water levels (approximately 50 cm) of the thermokarst lakes in the discontinuous permafrost zone within Western Siberia, after the anomalously hot summer in 2012. Analyzing four transects in Alaska, western Canada, Western Siberia, and East Siberia, Nitze et al. [16], concluded that lake area loss is dominant in discontinuous permafrost regions. In Western Siberia, a net lake loss of nearly 8% due to drainage and drying was observed [16]. However, in continuous permafrost regions, a general lake expansion pattern is dominant and was described in central Yakutia [19], Western Siberia [77], East Siberia [16]. The reduction in the surface area of water in permanent lakes is coincident with increased evaporation, due to a warming atmosphere [79]. Vanishing lakes in high-latitude regions are also caused by the permafrost thawing, which increases the subsurface drainage [80]. Also, permafrost thawing may release more nutrients into the warmer soils, stimulating the aquatic vegetation and further contribute to ponds shrinkage [81]. The expansion of seasonal lake area is most likely related to the overall increase of annual snow water equivalent, followed by snow melting peaks, increased water flow and budget in streams that connect the lakes and provide their recharge as mentioned in other studies with test sites in European Russia [21].
Regarding data quality and the confidence of the results, the usage of various Landsat sensors is known to be suitable for such applications [20], however with minor issues. As Nitze and Grosse, 2016 [20], we also used Landsat data at the highest preprocessing level, specifically surface reflectance Level-1 data, subjected to radiometric and atmospheric corrections. The radiometric corrections result in a high normalization of the signal among various sensors with 5% or less uncertainty between TM and ETM+ [82]. Therefore, the very low radiometric differences among sensors could hardly influence the trend analysis. Besides, Nill et al. [43] indicate that the trend analysis of NDVI has a high level of confidence as recording low temporal variations. The issue of the stripe pattern and the corresponding noise, recorded by other studies, has been overcome in our study by using only Landsat 7 satellite images acquired before the failing of the Scan Line Corrector (SLC).

Conclusions
Based on remote sensing images, we analyzed trends of summer vegetation greenness (NDVI) and moisture (NDMI) in the Pechora catchment between 1985 and 2019. We found a dominant greening and a consistent moistening of vegetation at the Pechora catchment scale and in the Pechora Delta, however also small areas with a browning trend have been identified. The evolution of NDVI is similar in tundra and taiga biomes, but its increase is slightly higher in the tundra area. Considering NDVI over different permafrost types, a different pattern is only found on the discontinuous permafrost terrain, where NDVI has a higher variability compared to permafrost free, isolated and sporadic permafrost areas that showed a notable increase of vegetation greenness during the analyzed period. One important finding of this study is that in the Pechora catchment, the remote sensing identified dynamics and trends of vegetation greenness are consistent with the evolution of the average summer temperature with a one year lag, while vegetation moisture is following annual total precipitation and annual snow water equivalent without lag. The respective climate parameters were estimated by the ERA-Interim reanalysis and observed at meteorological stations. ERA-Interim reanalysis data was found to be reliable over the Pechora catchment since it showed a very good agreement with the measured meteorological data.
At the Pechora catchment scale, an overall decrease in the area of permanent lakes was observed, but an expansion of seasonal lakes, with distinct behavior depending on the underlying terrain. The permanent and seasonal lakes showed a higher dynamic in areas of permafrost as compared to permafrost free terrain. Regarding the distribution of lakes in major biomes, there is a more evident change in the tundra as compared to taiga areas. The seasonal lake expansion is in accordance with the overall increase of annual snow water equivalent, estimated by ERA-Interim reanalysis across all study areas and different zones.
The general consistency between the indices of vegetation greenness and moisture based on satellite imagery and the climate data highlights the efficacy and reliability of combining Landsat satellite data, ERA-Interim reanalysis and meteorological data to monitor temporal dynamics of the land surface in Arctic areas.